Work fluctuations in the active Ornstein–Uhlenbeck particle model

An active particle is a physical entity able to transform energy from the environment or an internal reservoir into directed motion [1–4]. These particles are the fundamental constituents of active matter, a special class of out-of-equilibrium systems which have taken centre stage of statistical mechanics in the last few years. Nature offers already a plethora of examples of active matter systems, as colonies of microorganisms [1, 5, 6], living cells [1, 2, 7, 8], swarms, schools and flocks [1, 2, 9, 10], but active particles can also be produced artificially [3, 11–14]. From a theoretical standpoint, the interest towards active systems was fuelled by the display of intriguing collective properties, both on average, such as motility-induced phase separation [15–19] and at the level of fluctuations [20–24]. In the present paper, we focus on the problem of a single active particle free or interacting with an external potential. We take the perspective of thermodynamics of trajectories [25] and study the statistics of path-dependent observables: studies of this kind might reveal interesting properties even for the simplest example of close-to-equilibrium dynamics [26], while offering the possibility for a systematic treatment of arbitrarily far-from-equilibrium systems such as active systems [27].

The dynamics of an (overdamped) active particle under the action of an arbitrary potential U( r (t), t) and immersed in a thermal bath with friction coefficient γ and temperature T is prescribed by a stochastic equation of motion of the following form:

$$\begin{equation}\dot{\boldsymbol{r}}(t)=\boldsymbol{v}(t)-{\gamma }^{-1}\nabla U(\boldsymbol{r}(t),t)+\sqrt{2{D}_{\mathrm{T}}}\boldsymbol{\xi }(t),\end{equation} \tag{ 1 }$$

where r (t) is the position of the particle and v (t) its self-propulsion velocity. The constant D_T = γ⁻¹ k_B T denotes the translational diffusion coefficient, while ξ (t) is a zero-mean and delta-correlated white noise. The statistics of v (t) depends on the specifics of the self-propulsion mechanism and several possibilities have been considered in the literature: run-and-tumble particles, where the modulus of v (t) is constant while the direction changes at random Poissonian times, or active Brownian particles, where the direction changes as a Brownian motion, to name but a few examples [24, 28–30]. In this paper we adopt the active Ornstein–Uhlenbeck particle (AOUP) model, where v (t) itself is an Ornstein–Uhlenbeck process [31, 32]. The AOUP model retains the fundamental property of active particles, i.e. persistence of motion due to self-propulsion, and has enjoyed recent analytical insights regarding collective properties [32, 33] and energetics [27, 34]. Besides, the AOUP model can also represent passive tracers immersed in a bath of active particles, such as a bacterial bath [35, 36].

The observable we choose to characterise the dynamics of the system of interest is the energy injected by the self-propulsion force, or active work,

$$\begin{equation}{W}_{\tau }={\int }_{0}^{\tau }\mathrm{d}t\enspace \gamma \boldsymbol{v}(t)\cdot \dot{\boldsymbol{r}}(t).\end{equation} \tag{ 2 }$$

For a single free active particle, the active work can be identified with the heat flowing into the thermal bath [37] and is proportional to the entropy production [20, 38, 39], provided the self-propulsion velocity is assumed to be even under time reversal [23, 27, 34]. The active work is a natural observable of interest for the stochastic thermodynamics of active systems [27, 40–42], since it measures how efficiently active driving is converted into motion. Moreover, the fluctuations of the active work have been shown to be deeply connected to structural and dynamical properties of active systems [20, 22, 43] and studying their large deviations provides a pathway to control the collective behaviour of systems of active particles [21, 23, 44–46].

We study the asymptotic fluctuations of W_τ via the rate function of the active work,

$$\begin{equation}I(w)=\underset{\tau \to \infty }{\mathrm{lim}}\enspace \frac{1}{\tau }\enspace \mathrm{log}\left(\mathrm{P}\mathrm{r}\mathrm{o}\mathrm{b}.\left\{{W}_{\tau }/\tau =w\right\}\right),\end{equation} \tag{ 3 }$$

for an AOUP which is free or confined by a potential. On the one hand, our work is motivated by some recent results [20–23] showing that, for a variety of active particle models, interaction-related effects such as the formation of clusters [47–49] and the consequent drag against the direction of the active force [20, 22] induce singularities in the active work distribution. On the other hand, we are interested in characterising the free AOUP case and compare it to the corresponding passive problem, where the rate function of the energy injected by the thermal noise displays singular behaviour [26]. The main results of our analysis are summarised below.

• We compute I(w) exactly for a free AOUP in arbitrary dimension and find it to display two asymptotically linear tails but no non-analyticities, at variance with the power injected by uncorrelated thermal noise in absence of activity [26]. This result is based on the evaluation of sub-exponential contributions which are responsible for the singularity in the case of [26].
• We estimate numerically I(w) for an AOUP moving in a harmonic and two different anharmonic potential wells, and discuss the phenomenology of rare trajectories coming from lower- or higher-than-average fluctuations. In the harmonic case, we find, within numerical accuracy, the presence of linear tails, and a similar phenomenology to that of the free AOUP. In the anharmonic potentials, we do not find singular behaviour of the rate function in spite of accumulation effects at the well boundaries [50].
• We use analytical results for the free-particle case and numerical results both for the free and the confined problem to discuss the validity of a fluctuation relation [51] for the active work. We find it to hold for free (exactly) and harmonically confined (numerically) AOUP with respect to the bath temperature, independently of the strength of the active force. For an anharmonically confined particle at high activity, we find instead that a fluctuation relation is satisfied with respect to a different temperature, whose value agrees with the effective temperature obtained from the fluctuation–dissipation relation [52–54].

The remainder of the paper is organised as follows. In section 2, we consider a free AOUP and show the calculation of the active work rate function. Our derivation, shown in subsections 2.1 and 2.2, is based on a path-integral calculation of the scaled cumulant generating function (SCGF), which yields the rate function under Legendre–Fenchel transform (LFT). Such an approach yields also preasymptotic corrections to the cumulant generating function, which are generally required in order to guarantee that the rate function coincides with the LFT of the SCGF. In subsection 2.3 we provide two independent numerical estimates of the rate function that can be directly compared to our theoretical calculations: one obtained by a direct measure of I(w) from numerical integration of the equations of motion, the other obtained by estimating first the generating function of W_τ 's cumulants with a cloning algorithm [55–57] and then performing a Legendre–Fenchel transformation (see appendix A for details). The problem of a confined AOUP is discussed in section 3. Three different confining potentials are considered: a harmonic potential (subsection 3.1), a 'stiff' potential growing as r¹⁰, with r the distance of the particle from the origin, and a circular-shaped rigid barrier of fixed radius whose borders are modelled trough a WCA potential (subsection 3.2). Section 4 is devoted to the study of fluctuation relations for the active work in three settings: the free-particle case, the harmonically confined AOUP and the AOUP confined by the 'stiff' anharmonic potential. Finally, in section 5, we report the conclusions and final comments of our study.

Let us begin by recalling the details of the overdamped AOUP model. The model consists of two equations, for the d-dimensional position r (t) and self-propulsion velocity v (t) of the AOUP. Regarding the parameters appearing in the equations, we follow the convention of [50] and write

$$\begin{equation}\dot{\boldsymbol{r}}(t)=\boldsymbol{v}(t)-{\gamma }^{-1}\nabla U(\boldsymbol{r}(t),t)+\sqrt{2{D}_{\mathrm{T}}}\boldsymbol{\xi }(t),\end{equation} \tag{ 4a }$$

$$\begin{equation}\dot{\boldsymbol{v}}(t)=-{\gamma }_{\mathrm{R}}\boldsymbol{v}(t)+\sqrt{2{D}_{\mathrm{R}}^{\prime }}\boldsymbol{\eta }(t),\end{equation} \tag{ 4b }$$

where η (t) and ξ (t) are two zero-mean, unit-variance independent white noises, i.e.

$$\begin{equation}\left\langle {\xi }_{i}(t){\xi }_{j}({t}^{\prime })\right\rangle =\left\langle {\eta }_{i}(t){\eta }_{j}({t}^{\prime })\right\rangle ={\delta }_{i,j}\delta (t-{t}^{\prime })\quad \forall \enspace i,j=1,\dots ,d.\end{equation} \tag{ 5 }$$

The damping coefficient γ_R controls the exponential decay of correlations in v (t). Intuitively, γ_R can be written as (d − 1)D_R, where d is the spatial dimension and D_R the rotational diffusion coefficient of the direction of self-propulsion. As a result, the correlations of the self-propulsion satisfy

$$\begin{equation}\left\langle \boldsymbol{v}(t)\cdot \boldsymbol{v}(t)\right\rangle \stackrel{\vert t-{t}^{\prime }\vert \to \infty }{\to }\frac{d{D}_{\mathrm{R}}^{\prime }}{(d-1){D}_{\mathrm{R}}}\enspace {\text{e}}^{-(d-1){D}_{\mathrm{R}}\vert t-{t}^{\prime }\vert }.\end{equation} \tag{ 6 }$$

The coefficient $d{D}_{\mathrm{R}}^{\prime }/(d-1){D}_{\mathrm{R}}$ appearing in the right-hand side of equation (6) is the square of the typical modulus of the self-propulsion velocity, which can be thought of as the ratio between a typical self-propulsion force F_a and the mobility γ. In other words, the parameter ${D}_{\mathrm{R}}^{\prime }$ is fixed by $d{D}_{\mathrm{R}}^{\prime }=(d-1){D}_{\mathrm{R}}{({F}_{a}/\gamma )}^{2}$. The translational diffusion coefficient D_T is given by D_T = γ⁻¹ k_B T and the rotational diffusion coefficient is also proportional to the thermodynamic temperature k_B T. For a disk-shaped particle in d = 2, for instance, D_R = 3D_T/σ², with σ the particle diameter. In general, D_T and D_R could be considered as independent parameters, incorporating both thermal and active fluctuations: our choice implies a purely thermal origin, so that the only source of departure from equilibrium is the self-propulsion velocity. To sum up, the free parameters of the model are the temperature T, the typical magnitude of the active force F_a and the particle mobility γ. Additional parameters, required in order to specify the potential, will be introduced in section 3. In the remainder of this section, we set the potential U( r (t), t) to 0 and compute the large deviations of the active work (equation (2)) for the free AOUP.

Let us then turn to the main focus of this manuscript, the asymptotics of the probability-density-function Π(w) of the active work,

$$\begin{equation}{\Pi}(w)=\left\langle \delta ({W}_{\tau }-\tau w)\right\rangle \asymp {\text{e}}^{-\tau I(w)},\end{equation} \tag{ 7 }$$

where W_τ denotes the active work as a random variable and w the specific realisations, scaled by the observation time τ. The average here is performed over realisations of the stochastic noises ξ (t) and η (t); the symbol ≍ denotes equality of the large-τ limit on the logarithmic scale [58, 59]. By introducing, as it is customary, the Laplace representation of the delta function, $\delta (x)=(1/2\pi i){\int }_{-i\infty }^{i\infty }\enspace \mathrm{d}\lambda \enspace {\text{e}}^{\lambda x}$, we can write

$$\begin{equation}{\Pi}(w)=\frac{1}{2\pi i}{\int }_{-i\infty }^{i\infty }\enspace \mathrm{d}\lambda \enspace {\text{e}}^{-\tau \lambda w}\left\langle {\text{e}}^{\lambda {W}_{\tau }}\right\rangle \equiv \frac{1}{2\pi i}{\int }_{-i\infty }^{i\infty }\enspace \mathrm{d}\lambda \enspace {\text{e}}^{-\tau \lambda w}\hat{{\Pi}}(\lambda ),\end{equation} \tag{ 8 }$$

where we have introduced the generating function of W_τ 's cumulants (CGF), $\hat{{\Pi}}(\lambda )$. Because of the exponential factor in the integrand, once $\hat{{\Pi}}(\lambda )$ is obtained, the asymptotic of Π(w) can be estimated with a saddle-point expansion. We now turn to the computation of the cumulant generating function $\hat{{\Pi}}(\lambda )$. The saddle-point expansion is performed in subsection 2.2. Finally, in subsection 2.3, we compare our analytical prediction with the result of direct simulations of the model and a refined estimate obtained with a biased sampling of the model trajectories.

2.1. The cumulant generating function

In this section we compute the generating function of W_τ 's cumulants,

$$\begin{equation}\hat{{\Pi}}(\lambda )=\left\langle {\text{e}}^{\lambda {W}_{\tau }}\right\rangle ,\end{equation} \tag{ 9 }$$

for the free AOUP, following the approach of [26]. The average appearing in equation (9) is performed with respect to the realisations, or paths, of the noises ξ (t) and η (t) affecting position and self-propulsion speed in equation (4). Since, by equation (4a),

$$\begin{equation}{W}_{\tau }=\gamma {\int }_{0}^{\tau }\mathrm{d}t\enspace \boldsymbol{v}(t)\cdot \dot{\mathbf{r}}(t)=\gamma {\int }_{0}^{\tau }\mathrm{d}t\enspace \boldsymbol{v}(t)\cdot \boldsymbol{v}(t)+\gamma \sqrt{2{D}_{\mathrm{T}}}{\int }_{0}^{\tau }\mathrm{d}t\enspace \boldsymbol{v}(t)\cdot \boldsymbol{\xi }(t),\end{equation} \tag{ 10 }$$

the path probabilities required to perform the average are those of ξ (t) (denoted $\mathcal{P}[\boldsymbol{\xi }(t)]$) and v (t) (denoted $\mathcal{P}[\boldsymbol{v}(t)\vert {\boldsymbol{v}}_{0}]{p}_{0}({\boldsymbol{v}}_{0})$, conditioned to the initial condition v ₀ with initial probability p₀( v ₀)). These path probabilities refer to time intervals of length τ, but we omit τ from the symbols so as to ease the notation. The relevant path probabilities are given by [60],

$$\begin{equation}\begin{aligned} \mathcal{P}[\boldsymbol{\xi }(t)]& =\mathrm{exp}\left\{-\frac{1}{2}{\int }_{0}^{\tau }\mathrm{d}t\enspace \boldsymbol{\xi }(t)\cdot \boldsymbol{\xi }(t)\right\}, \\ \mathcal{P}[\boldsymbol{v}(t)\vert {\boldsymbol{v}}_{0}]& =\mathrm{exp}\left\{-\frac{1}{4{D}_{\mathrm{R}}^{\prime }}{\int }_{0}^{\tau }\mathrm{d}t\enspace {(\dot{\boldsymbol{v}}(t)+{\gamma }_{\mathrm{R}}\boldsymbol{v}(t))}^{2}+\frac{d}{2}{\gamma }_{\mathrm{R}}\tau \right\}, \end{aligned}\end{equation} \tag{ 11 }$$

where the time-integrals are to be intended according to the Stratonovich discretisation scheme [61]. The corresponding measures, $\mathcal{D}[\boldsymbol{\xi }(t)]$ and $\mathcal{D}[\boldsymbol{v}(t)\vert {\boldsymbol{v}}_{0}]$ include also the proper normalisation factor. $\mathcal{P}[\boldsymbol{v}(t)\vert {\boldsymbol{v}}_{0}]$, in particular, is the path probability of the self-propulsion process with fixed endpoints v ₀ and v (τ) ≡ v _τ . Denoting with p( v _τ , τ| v ₀, 0) the transition probability of the Ornstein–Uhlenbeck process [61], the normalisation factor is thus fixed by the following condition,

$$\begin{equation}\int \mathcal{D}[\boldsymbol{v}(t)\vert {\boldsymbol{v}}_{0}]\mathcal{P}[\boldsymbol{v}(t)\vert {\boldsymbol{v}}_{0}]=p({\boldsymbol{v}}_{\tau },\tau \vert {\boldsymbol{v}}_{0},0).\end{equation} \tag{ 12 }$$

We can now unwind the average symbol of equation (9), so as to have

$$\begin{equation}\begin{aligned} \hat{{\Pi}}(\lambda )& =\int \mathrm{d}{\boldsymbol{v}}_{0}\enspace \mathrm{d}{\boldsymbol{v}}_{\tau }\enspace \mathcal{D}[\boldsymbol{\xi }(t)]\mathcal{D}[\boldsymbol{v}(t)\vert {\boldsymbol{v}}_{0}]\mathcal{P}[\boldsymbol{v}(t)\vert {\boldsymbol{v}}_{0}]\mathcal{P}[\boldsymbol{\xi }(t)]\enspace {p}_{0}({\boldsymbol{v}}_{0}) \\ & \quad \times \mathrm{exp}\left\{\lambda \gamma \left({\int }_{0}^{\tau }\mathrm{d}t\enspace \boldsymbol{v}(t)\cdot \boldsymbol{v}(t)+\sqrt{2{D}_{\mathrm{T}}}{\int }_{0}^{\tau }\mathrm{d}t\enspace \boldsymbol{v}(t)\cdot \boldsymbol{\xi }(t)\right)\right\}, \end{aligned}\end{equation} \tag{ 13 }$$

which satisfies the normalisation condition $\hat{{\Pi}}(0)=1$.

Because of the Stratonovich convention, the rules of standard calculus can be used to simplify the right-hand side of equation (13):

$$\begin{equation}\begin{aligned} \hat{{\Pi}}(\lambda )& ={\text{e}}^{\frac{d{\gamma }_{\mathrm{R}}\tau }{2}}\int \mathrm{d}{\boldsymbol{v}}_{0}{p}_{0}({\boldsymbol{v}}_{0})\enspace \mathrm{d}{\boldsymbol{v}}_{\tau }\enspace \mathcal{D}[\boldsymbol{v}(t)\vert {\boldsymbol{v}}_{0}]\enspace \mathcal{D}[\boldsymbol{\xi }(t)]\mathrm{exp}\left\{-\frac{{\gamma }_{\mathrm{R}}}{4{D}_{\mathrm{R}}^{\prime }}({\boldsymbol{v}}_{\tau }^{2}-{\boldsymbol{v}}_{0}^{2})\right\} \\ & \quad \times \mathrm{exp}\left\{-\frac{1}{4{D}_{\mathrm{R}}^{\prime }}{\int }_{0}^{\tau }\mathrm{d}t\left({\dot{\boldsymbol{v}}}^{2}(t)+{\alpha }^{2}{\boldsymbol{v}}^{2}(t)\right)-\frac{1}{2}{\int }_{0}^{\tau }{\left(\boldsymbol{\xi }(t)-\sqrt{2{D}_{\mathrm{T}}}\lambda \gamma \boldsymbol{v}(t)\right)}^{2}\right\}, \end{aligned}\end{equation} \tag{ 14 }$$

where we have set ${\alpha }^{2}={\gamma }_{\mathrm{R}}^{2}-4{D}_{\mathrm{R}}^{\prime }\lambda \gamma (1+\lambda \gamma {D}_{\mathrm{T}})$. After the linear change of variables ${\boldsymbol{\xi }}^{\prime }(t)=\boldsymbol{\xi }(t)-\sqrt{2{D}_{\mathrm{T}}}\gamma \lambda \boldsymbol{v}(t)$ the integral over the thermal noise can be performed—it equals 1 due to normalisation. The remaining terms can all be written as products over the spatial components of v (t), so that the integral, due to the independence of the components of an Ornstein–Uhlenbeck process, factorises. By isotropy, each of the factors yields the same contribution. We will, in addition, make use of the following identity for a one-dimensional Ornstein–Uhlenbeck process v(t) and $\alpha \in \mathbb{R}$ [62],

$$\begin{equation}\begin{aligned} & \int \mathcal{D}[v(t)\vert {v}_{0}]\mathrm{exp}\left\{-\frac{1}{4{D}_{\mathrm{R}}^{\prime }}{\int }_{0}^{\tau }\left[\dot{v}{(t)}^{2}+{\alpha }^{2}v{(t)}^{2}\right]\right\} \\ & \qquad ={\left(\pi \frac{4{D}_{\mathrm{R}}^{\prime }}{\alpha }\enspace \mathrm{sinh}(\alpha \tau )\right)}^{-\frac{1}{2}}\enspace \mathrm{exp}\left\{-\frac{\alpha }{4{D}_{\mathrm{R}}^{\prime }}\frac{({v}_{\tau }^{2}+{v}_{0}^{2})\mathrm{cosh}(\alpha \tau )-2{v}_{0}\enspace {v}_{\tau }}{\mathrm{sinh}(\alpha \tau )}\right\}, \end{aligned}\end{equation} \tag{ 15 }$$

and set p₀(v₀) to the stationary probability of the self-propulsion process,

$$\begin{equation}{p}_{0}({v}_{0})=\sqrt{\frac{{\gamma }_{\mathrm{R}}}{2\pi {D}_{\mathrm{R}}^{\prime }}}{\text{e}}^{-\frac{{\gamma }_{\mathrm{R}}}{2{D}_{\mathrm{R}}^{\prime }}{v}_{0}^{2}}.\end{equation} \tag{ 16 }$$

Therefore, after performing the integrals over v ₀ and v _τ , we obtain the following expression for the cumulant generating function,

$$\begin{equation}\begin{aligned} \hat{{\Pi}}(\lambda )& ={\text{e}}^{\frac{d{\gamma }_{\mathrm{R}}\tau }{2}}{\left(\frac{{\gamma }_{\mathrm{R}}^{2}+{\alpha }^{2}}{\alpha {\gamma }_{\mathrm{R}}}\enspace \mathrm{sinh}(\tau \alpha )+2\enspace \mathrm{cosh}(\tau \alpha )\right)}^{-\frac{d}{2}} \\ & =\frac{{\text{e}}^{\tau \frac{d}{2}\left({\gamma }_{\mathrm{R}}-\alpha \right)}}{{\left(\frac{1+{\text{e}}^{-2\tau \alpha }}{2}\right)}^{d/2}{\left(2+\frac{{\gamma }_{\mathrm{R}}^{2}+{\alpha }^{2}}{{\gamma }_{\mathrm{R}}\alpha }\enspace \text{tanh}(\tau \alpha )\right)}^{d/2}}. \end{aligned}\end{equation} \tag{ 17 }$$

Let us pause briefly and comment on the last result: the cumulant generating function consists of a factor which grows exponentially with τ and a sub-exponential prefactor. The exponential factor gives rise to the SCGF,

$$\begin{equation}g(\lambda )=\underset{\tau \to \infty }{\mathrm{lim}}\enspace \frac{1}{\tau }\enspace \mathrm{log}\hat{{\Pi}}(\lambda ),\end{equation} \tag{ 18 }$$

which reads, for this specific problem,

$$\begin{equation}g(\lambda )=\frac{d}{2}\left({\gamma }_{\mathrm{R}}-\sqrt{{\gamma }_{\mathrm{R}}^{2}-4{D}_{\mathrm{R}}^{\prime }\gamma \lambda (1+\gamma {D}_{\mathrm{T}}\lambda )}\right).\end{equation} \tag{ 19 }$$

This expression agrees with the earlier result of [39] in d = 2 and generalises it to any spatial dimension. Importantly, we also notice that if the parameters γ_R and ${D}_{\mathrm{R}}^{\prime }$ were free to take any value, instead of being constrained by the relations expounded at the beginning of the section, then the spatial dimension would only appear as an overall multiplicative factor in g(λ). The implication is that, for the free AOUP, increasing the spatial dimension only increases the rate at which the probability Π(w) of active work fluctuations concentrates around the mean, while it does not affect the overall shape of the rate function.

In general, the leading exponential behaviour of $\hat{{\Pi}}(\lambda )$ provides enough information for the calculation of the rate function via LFT [58, 59]. However, a complex subexponential prefactor such as the one appearing on the right-hand side of equation (17) might result in additional non-analyticities of $\hat{{\Pi}}(\lambda )$ which still affect the leading exponential behaviour [26]. This property is clear when the inverse Laplace transform in equation (8) is computed with a saddle-point method: denoting with F(λ) the subexponential prefactor, such that $\hat{{\Pi}}(\lambda )=F(\lambda ){\text{e}}^{\tau g(\lambda )}$, equation (8) can be cast in the following form,

$$\begin{equation}{\Pi}(w)=\frac{1}{2\pi i}{\int }_{-i\infty }^{i\infty }\mathrm{d}\lambda \enspace F(\lambda )\enspace {\text{e}}^{-\tau \left[\lambda w-g(\lambda )\right]}.\end{equation} \tag{ 20 }$$

Not only F(λ) concurs to the subexponential corrections to the asymptotics of Π(w), but also its non-analyticities might pose severe limitations to the deformation of the integration contour. Therefore, it is worth including the subexponential prefactor in the saddle-point estimation of the integral above, which is carried out in the next section.

2.2. Saddle-point estimation of the rate function

In order to simplify the discussion, we introduce the following rescaled variables,

$$\begin{equation}\begin{aligned} {\tilde{\gamma }}_{\mathrm{R}}& =\tau {\gamma }_{\mathrm{R}},\qquad \tilde{\lambda }=\frac{4{D}_{\mathrm{R}}^{\prime }\gamma }{{\gamma }_{\mathrm{R}}^{2}}\lambda , \\ \tilde{\alpha }& =\frac{\alpha }{{\gamma }_{\mathrm{R}}}=\sqrt{1-\tilde{\lambda }-A{\tilde{\lambda }}^{2}},\quad \text{with}\;A=\frac{{\gamma }_{\mathrm{R}}^{2}{D}_{\mathrm{T}}}{4{D}_{\mathrm{R}}^{\prime }}. \end{aligned}\end{equation} \tag{ 21 }$$

In terms of the new variables,

$$\begin{equation}g(\lambda )=\frac{d}{2}{\gamma }_{\mathrm{R}}(1-\tilde{\alpha }),\end{equation} \tag{ 22a }$$

$$\begin{equation}F(\lambda )={\left(1+{\text{e}}^{-2\tilde{\alpha }{\tilde{\gamma }}_{\mathrm{R}}}\right)}^{-d/2}{\left(1+\frac{1}{2}\left(\tilde{\alpha }+\frac{1}{\tilde{\alpha }}\right)\mathrm{tanh}({\tilde{\gamma }}_{\mathrm{R}}\tilde{\alpha })\right)}^{-d/2}.\end{equation} \tag{ 22b }$$

Although the original variables were different, in terms of the rescaled variables we get exactly the same result as [26], hence the same considerations apply here as well. In particular, the square root in the definition of $\tilde{\alpha }$ introduces two branch points at

$$\begin{equation}\tilde{\lambda }={\tilde{\lambda }}_{1/2}=-\frac{1\pm \sqrt{1+4A}}{2A}.\end{equation} \tag{ 23 }$$

Assuming $\tilde{\alpha }\enspace > \enspace 0$ for positive real arguments of the square root (i.e. ${\tilde{\lambda }}_{1}\enspace < \enspace \tilde{\lambda }\enspace < \enspace {\tilde{\lambda }}_{2}$) is equivalent to considering two branch cuts on the real axis: one originating in ${\tilde{\lambda }}_{1}$ and running towards negative values, the other originating in ${\tilde{\lambda }}_{2}$ and running towards positive values. Other possible singularities are the poles of F(λ) (equation (22b)). It is straightforward to check that the denominator of F(λ) has purely imaginary roots $\tilde{\alpha }=iy$, with y satisfying

$$\begin{equation}\left(y-\frac{1}{y}\right)\mathrm{tan}\left({\tilde{\gamma }}_{\mathrm{R}}y\right)=2.\end{equation} \tag{ 24 }$$

In terms of $\tilde{\lambda }$, these points are located at

$$\begin{equation}{\tilde{\lambda }}_{1/2}(y)=-\frac{1\pm \sqrt{1+4A\left(1+{y}^{2}\right)}}{2A},\end{equation} \tag{ 25 }$$

thus are covered by the two branch cuts (see figure 1). One concludes that, at variance with the passive case examined in [26], the sub-exponential pre-factor of $\hat{{\Pi}}(\lambda )$ does not induce additional singularities in the complex plane. Let us notice that the overall picture can be extended to all dimensions d, although for odd dimensions the roots of the denominator of equation (22b) generate branch points rather than poles.

Zoom In Zoom Out Reset image size

Knowing the analyticity of $\hat{{\Pi}}(\lambda )$ we can now proceed with the calculation of Π(w), whose integral expression we repeat here for clarity,

$$\begin{equation}{\Pi}(w)=\frac{1}{2\pi i}{\int }_{-i\infty }^{i\infty }\mathrm{d}\lambda \enspace F(\lambda )\enspace {\text{e}}^{-\tau \left[\lambda w-g(\lambda )\right]}.\end{equation} \tag{ 26 }$$

The procedure is the following: first, one must identify the stationary points of the function at the exponent, which are all saddle points in the complex plane. Then, by Cauchy's theorem, the integration contour can be deformed so that it passes through the highest of the saddle points. In particular, the deformed contour is chosen so that it crosses the saddle point along the direction of steepest descent, so that only the portion of the deformed contour close to the saddle point contributes to the integral as τ → ∞ [63].

In the present problem, the exponent of the integrand is, from equation (20), f(λ) = λw − g(λ). It is convenient to consider again the rescaled variables, i.e.

$$\begin{equation}f(\tilde{\lambda })={\gamma }_{\mathrm{R}}\left(\tilde{w}\tilde{\lambda }-1+\tilde{\alpha }(\tilde{\lambda })\right),\end{equation} \tag{ 27 }$$

with $\tilde{w}={\gamma }_{\mathrm{R}}w/(4{D}_{\mathrm{R}}^{\prime }\gamma )$. The condition ${f}^{\prime }(\tilde{\lambda })=0$ is satisfied by

$$\begin{equation}{\tilde{\lambda }}_{\pm }^{(s)}=-\frac{1\pm \sqrt{1+4A\frac{4{\tilde{w}}^{2}-1}{4{(A+\tilde{w})}^{2}}}}{2A}\end{equation} \tag{ 28 }$$

and it is straightforward to check that the highest saddle point is ${\tilde{\lambda }}_{-}^{(s)}$ for $\tilde{w}\enspace \geqslant \enspace 0$, ${\tilde{\lambda }}_{+}^{(s)}$ for $\tilde{w}\enspace < \enspace 0$. It is also worth noticing that the saddle point lies within the interval $\left[{\tilde{\lambda }}_{1},{\tilde{\lambda }}_{2}\right]$ for every value of $\tilde{w}$, so that the integration contour can always be deformed so as to pass by the saddle point without encountering any non-analyticity of the integrand.

Although analytical expressions for the steepest descents curves are not available, it can be shown that they cross the saddle point in the direction of the imaginary axis. Therefore, after performing a quadratic expansion of the exponent $f(\tilde{\lambda })$ around the saddle point, the following asymptotic expression is obtained for Π(w),

$$\begin{equation}{\Pi}(w)\asymp \frac{F({\tilde{\lambda }}^{(s)})}{2\pi }{\left(\frac{2\pi }{\tau \mathcal{C}}\right)}^{1/2}{\text{e}}^{-\tau {\gamma }_{\mathrm{R}}\left[\tilde{w}{\tilde{\lambda }}^{(s)}-1+\tilde{\alpha }({\tilde{\lambda }}^{(s)})\right]}.\end{equation} \tag{ 29 }$$

Here $\mathcal{C}$ denotes the modulus of the second derivative of f(λ) along the steepest descent path—i.e. the imaginary direction in the complex plane—and ${\tilde{\lambda }}^{(s)}(\tilde{w})$ coincides with ${\tilde{\lambda }}_{+}^{(s)}$ for negative $\tilde{w}$'s and with ${\tilde{\lambda }}_{-}^{(s)}$ for positive $\tilde{w}$'s. From the exponential factor on the right-hand side of equation (29) we can finally extract a definition, for the rate function of the active work's fluctuations,

$$\begin{equation}I(w)=\frac{d{\gamma }_{\mathrm{R}}}{2}\tilde{I}\left(\frac{2}{d}\frac{{\gamma }_{\mathrm{R}}}{4{D}_{\mathrm{R}}^{\prime }\gamma }w\right),\qquad \tilde{I}(\tilde{w})=\sqrt{\frac{1+4A}{A+{\tilde{w}}^{2}}}\left(\frac{A+{\tilde{w}}^{2}}{2A}\right)-1-\frac{\tilde{w}}{2A},\quad A=\frac{{\gamma }_{R}^{2}{D}_{\mathrm{T}}}{4{D}_{\mathrm{R}}^{\prime }},\end{equation} \tag{ 30 }$$

which is shown in figure 2 (left) for a specific choice of the system parameters. At variance with the rate function of the power injected by the thermal noise on a Brownian particle, which becomes linear abruptly at a certain critical value of w [26], the rate function of equation (3) is only asymptotically linear. Regarding the strong asymmetry of fluctuations about the mean, however, the phenomenology is similar to the passive case: positive fluctuations of the Ornstein–Uhlenbeck noise generate trajectories with a large velocity, resulting in an even higher energy uptake at later times. Due to this positive feedback, higher-than-average fluctuations of the active work are significantly more likely than lower-than-average ones, consistently with the slope of the right branch of equation (30) being smaller in modulus than the slope of the left branch.

Zoom In Zoom Out Reset image size

2.3. Numerical estimates via direct sampling and cloning algorithm

The rate function of W_τ —or any dynamical observable—can also be estimated numerically in two different ways. The first one involves the direct sampling of a large number of values for the active work, obtained by solving the dynamics of the problem for several initial conditions and realisation of the stochastic forces, then applying the definition equation (3) to the empirical distribution function. An empirical estimate is in fact all that can be obtained for the general confined AOUP problem, considered in the next section. The second method we employ is the population Monte Carlo, or simply cloning algorithm [55–57], where one samples the SCGF and then obtains the rate function using the LFT. This method is especially useful to access values of the active work in the tails of the distribution, which typically require a number of samples exponentially large in τ. The details of both numerical integration methods are reported in appendix A. Here we only report our choice of parameters: we set k_B T = 0.05, γ = 10 and σ = 1, such that D_T = 0.005 and D_R = 0.015. Moreover we set the self-propulsion force to F_a = 20, such that the Péclet number, computed as Pe = σF_a /k_B T = 4 × 10² [50], is large.

Figure 2 shows the empirical rate functions on the right, obtained from the histogram of the active work using equation (3), for different values of τ as reported in the key. As soon as the observation time τ approaches 2 × 10³, the rate function converges to the exact result. Figure 3 (left) shows instead a numerical estimate of the free-particle SCGF for the same set of parameters considered in figure 2 and increasing number of copies N_c. Notice that the closer λ to zero, the faster the convergence of the numerical estimate to the theoretical curve. In fact, far away from the 'typical' value λ ≃ 0, which corresponds to an unbiased (or almost unbiased) evolution, the estimate of g(λ) depends crucially on one of the clones achieving some rare fluctuation. Besides, it is interesting to notice that choosing a value of λ which lies outside of the domain of g(λ) results in a break down of the numerical scheme. Specifically, while the numerical estimate of g(λ) for λ close to the left limit of the domain does not overlap with the analytical curve, it is still sensible. Conversely, when applied out of the domain the algorithm returns either diverging or highly oscillating values. This feature is worth remembering when studying problems without a known analytical solution, such as the cases reported in section 3. The rate function can be finally obtained with a Legendre–Fenchel transform (LFT), provided the SCGF is smooth and steep [58, 59]. Note that the accuracy of the LFT depends on the sampling density of the SCGF itself. For instance, figure 2 (right) shows a comparison between the LFT of the SCGF and the rate function obtained via direct measurement: in order to achieve high accuracy around the minimum (shown in the inset of the figure) we have performed a denser sampling of the cloning SCGF around λ = 0. The cost of a denser sampling is mitigated by the fact that, close to λ = 0, convergence of the SCGF is achieved with a comparatively low number of clones (cf figure 3). The same sampling criterion will be used in the presence of confining potentials in order to compare the LFT with the rate function measured from direct numerical integration of the dynamics.

Zoom In Zoom Out Reset image size

Figure 3 (right) shows the typical trajectories of the clones for a few different values of the bias λ. These can be thought of as a representation of the trajectories of a model which is biased so as to have a λ-dependent average active work w(λ), with w(λ) given by the saddle-point condition w(λ) = g'(λ). At λ ≃ 0 (trajectory starting on the blue dot in the figure), the particle follows the typical AOUP dynamics—w(0) coincides with the unbiased average active work. For λ large and negative, the trajectory (starting on the green and red dots in the figure) produces an active work which is much smaller than the free-AOUP average. Specifically, for the one starting from the green dot, with λ ≃ −10, the active work is close to zero, whereas the one starting from the red dot produces a negative active work. The trajectories themselves hint at the mechanism by which non-typical active works are produced: a vanishing active work, for instance, can be produced by having a sequence of atypically small kicks from the active Ornstein–Uhlenbeck noise, or by having the thermal delta-correlated noise to act in opposition with the active noise: both cases result in the particle moving much less than on average, as it is the case for the corresponding trajectory displayed in figure 3. A negative active work, instead, requires the particle's velocity to be consistently opposite to the active noise, which can be realised by having the thermal noise not only opposite to but also larger in magnitude than the active noise. The result is again a trajectory which displays little displacement from its initial position. It is interesting to notice that the latter type of trajectory is only allowed in the presence of thermal noise: letting the thermal diffusion coefficient D_T to vanish would cause the whole branch of g(λ) with negative derivative to disappear, signalling that only positive values of the active work are allowed.

We now apply the numerical techniques discussed in subsection 2.3 to AOUPs subject to different confining potentials. For each case, we present a numerical estimate of the SCGF obtained with a cloning algorithm and compute the rate function as a LFT. We then compare such rate functions with those estimated directly from the logarithm of the empirical distribution of the active work. In addition, by looking at biased trajectories, we explore the mechanisms leading to active work fluctuations in the various potentials considered.

It is worth recalling the phenomenology of the fluctuations of the injected power for a confined Brownian particle, as it provides a reference frame for a better understanding of the confined AOUP problem. The confined Brownian problem was first examined in [26], where the author concludes that an external potential has only pre-asymptotic effects on the fluctuations of the injected power, therefore the rate function remains equal to that of the free-particle problem. In simple terms, the asymptotic fluctuations of the injected power are not affected by the presence of a confining potential. This picture is consistent with the intuitive explanation of the general shape of the rate function, which is caused by the positive feedback between fluctuations of the noise and energy uptake of the particle mentioned at the end of subsection 2.2. Much more recent than its passive counterpart, the problem of a confined active particle is currently under the scrutiny of the active matter community [50, 64–67]. It is well understood, for instance, that the two timescales of the problem—one related to the confining potential and the other to the persistent noise—compete in the determination of the steady-state distribution of the particle's position. Thus, depending on the parameters of the problem, in steady state the active particle will be either pushing against the slope of the confining well or fluctuating in the middle of the well in the low persistence limit.

Within the specific context of the AOUP model, it is interesting to notice the singularity of the harmonic case with respect to generic confining potentials. In this case, as in the free problem, the dynamics satisfies detailed balance [31, 68], resulting in an effective equilibrium regime. As a result, the harmonically confined AOUP does not display the 'pushing' phase in steady state for any value of the parameters. Such peculiar aspect of AOUPs is reflected in the fluctuations of the active work. In fact, our study reveals and highlights the differences in the phenomenology of active work fluctuations between the harmonic case (subsection 3.1) and cases where the particle is confined by nonlinear confining potentials (subsection 3.2), at variance with the passive problem where the confining potential has little-to-no effect [26].

3.1. Harmonic potential

The AOUP of this section interacts with a confining harmonic potential U( r ) = k r²/2, so that the equations of motion in the overdamped limit are

$$\begin{equation}\dot{\boldsymbol{r}}(t)=\boldsymbol{v}(t)-{\gamma }^{-1}k\boldsymbol{r}(t)+\sqrt{2{D}_{\mathrm{T}}}\boldsymbol{\xi }(t),\end{equation} \tag{ 31a }$$

$$\begin{equation}\dot{\boldsymbol{v}}(t)=-{\gamma }_{\mathrm{R}}\boldsymbol{v}(t)+\sqrt{2{D}_{\mathrm{R}}^{\prime }}\boldsymbol{\eta }(t).\end{equation} \tag{ 31b }$$

As mentioned in the previous paragraph, the harmonic potential is somewhat special for the AOUP, as it results in linear equations of motion. Therefore, the position and active force processes r (t) and v (t) are Gaussian, like the noises ξ (t) and η (t). Furthermore, the problem can be shown to have an effective formulation as an equilibrium problem [31] and the steady-state distribution of the AOUP position is a Gaussian distribution centred at the bottom of the potential well. From the perspective of active work fluctuations, there does not seem to be a fundamental difference with respect to the free case of section 2, apart from an obvious reduction of the average active work. In fact, after solving equation (31a) for r (t), it is straightforward to show that (assuming vanishing position and self-propulsion speed at t = 0 for simplicity)

$$\begin{equation}\left\langle \boldsymbol{v}(t)\cdot \dot{\boldsymbol{r}}(t)\right\rangle =\left\langle \boldsymbol{v}(t)\cdot \boldsymbol{v}(t)\right\rangle -{\gamma }^{-1}k{\int }_{0}^{t}\mathrm{d}{t}^{\prime }\enspace \left\langle \boldsymbol{v}(t)\cdot \boldsymbol{v}({t}^{\prime })\right\rangle \stackrel{t\to \infty }{\to }\frac{d{D}_{\mathrm{R}}^{\prime }}{{\gamma }_{\mathrm{R}}}\frac{{\gamma }_{\mathrm{R}}}{{\gamma }_{R}+{\gamma }^{-1}k}.\end{equation} \tag{ 32 }$$

The average active work (per unit time) is obtained by multiplying the expression above by γ: as k > 0, the harmonic average is always smaller than the free average $\gamma d{D}_{\mathrm{R}}^{\prime }/{\gamma }_{\mathrm{R}}$.

The cloning estimate of the SCGF g^harm(λ) is shown in figure 4 (left), together with some representative trajectories (right). For the sake of comparison, the analytical SCGF of the free case equation (19) is also shown in the left panel (black solid line). The rate function resulting from the LFT of g^harm(λ) is shown in figure 5 and compared to the direct numerical estimate. The trajectory starting from the green dot in figure 4 (right) is an example of typical trajectory (λ = 0). The figure also reports a trajectory with larger-than-average active work (λ > 0, starting from the blue dot in the inset) and one with smaller-than-average active work (λ < 0, red dot). On a qualitative level, the trajectory with a positive bias only looks more persistent than the average trajectory, whereas the one with a negative bias remains closer to the initial condition r ₀ = 0, as it was the case with negatively biased trajectories of the free problem. However, the cloning estimate of the harmonic SCGF reveals two important differences with respect to the free problem, where the SCGF g^free(λ) is defined in an interval $\left[{\lambda }_{1}^{\text{free}},{\lambda }_{2}^{\text{free}}\right]$ and it is steep, i.e. the derivative diverges at the boundaries of the domain. On the one hand, for large positive λ's the algorithm returns diverging results, supporting the hypothesis that also the harmonic SCGF diverges at a certain $\lambda ={\lambda }_{2}^{\text{harm}}$. In addition, the derivative of the harmonic SCGF in 0 is much smaller than in the free case and it does not vary significantly for small positive λ's. Therefore, our best estimate does not rule out the possibility of a finite derivative of g^harm(λ) at the boundary of the domain. We will comment further on this point later, when discussing rate functions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

On the other hand, for λ large and negative, our numerics do not reveal any sign of divergence, although a higher number of clones might be required for estimating the actual values of g^harm(λ). In other words, what we observe here is analogous to what observed in the free problem for λ slightly bigger than ${\lambda }_{1}^{\text{free}}$, where, because of the large magnitude of λ, the largest number of clones we can afford does not grant the coincidence of numerical estimate and analytical prediction. We conclude that, at variance with the free AOUP case, the SCGF of the active work might not diverge for large negative λ, thus having a domain of the form $(-\infty ,{\lambda }_{2}^{\text{harm}}]$.

We have also obtained an independent estimate of the rate function itself, by sampling the active work in direct numerical solutions of equation (31). With respect to higher-than-average fluctuations, the difference between a steep g(λ), the derivative of which diverges at the boundary of the domain, and a non-steep one is the following: the rate function associated with a steep SCGF is only asymptotically linear, with the slope of the linear branch coinciding with the domain boundary λ₂, whereas a non-steep SCGF results in a rate function I(w) which is exactly linear after a threshold w*. The threshold coincides with the limiting slope of the SCGF, i.e. w* = g'(λ₂). The rate function for the harmonically confined AOUP is shown in figure 5 and it does, within numerical uncertainty, approach a linear branch right after the minimum. Let us nevertheless stress that our estimate cannot exclude a very steep rise of g^harm(λ)'s derivative very close to ${\lambda }_{2}^{\text{harm}}$, which would result again in a steep SCGF. Also the left branch of the numerical rate functions appears linear and there are no reasons to expect a different phenomenology of lower-than-average fluctuations with respect to the free problem. However, the finiteness of the cloning SCGF indicates that there might be differences very far in the tails, although such regions cannot be accessed by the numerical techniques at hand. The proposed scenario seems to be robust also with respect to variation of k, taking into account that as k is increased the average active work decreases (equation (32)) and fluctuations become generally rarer, as shown in figure 5 (right).

3.2. Anharmonic potentials

In this section we let the AOUP interact with two anharmonic, radially-symmetric confining potentials. The dynamics obeys equation (39): in the first case we examine the potential U_stiff( r ) = k_stiff r ¹⁰/10, which we refer to as the 'stiff' potential. In such confining potential the AOUP displays a different steady-state behaviour depending on the Péclet number Pe. Three representative cases are shown in figure 6, where all the parameters but k_B T are fixed, and we consider the presence and also the absence of the thermal noise ξ (t). For high temperatures and low Péclet, on the left and centre, the steady-state distribution of the AOUP position is peaked at the origin, whereas for low temperatures and high Péclet, on the right, it displays an annular peak with a finite radius. We study the large deviations of the active work in the regime of parameters resulting in an annular steady-state distribution for the position in presence of the thermal noise, as in figure 6 (right), and report here the corresponding results. In the other regime, where the steady-state distribution of the position is peaked at the origin, one would only observe progressively symmetric rate functions, as the thermal contribution to active work fluctuations begins to dominate over the active contribution.

Zoom In Zoom Out Reset image size

Also in this case we provide an estimate of the SCGF g^stiff(λ) obtained via cloning and an estimate of the rate function obtained by sampling the empirical distribution of the active work over several independent numerical solutions, following the method described in subsection 2.3. The results are shown in figure 7, with the SCGF on the left and rate function on the right. As in the free (subsection 2.3) and harmonic (subsection 3.1) case, the cloning algorithm returns unphysical results for λ larger than a certain threshold ${\lambda }_{2}^{\text{stiff}}$, indicating that the domain of g^stiff(λ) might have an upper extremum. In contrast with the harmonic case, however, our numerical estimates are compatible with a steep SCGF, i.e. with g^stiff(λ) reaching an infinite derivative at ${\lambda }_{2}^{\text{stiff}}$. Direct numerical estimates of the rate function indeed do not show any linear branch, supporting the hypothesis of a steep SCGF. For λ large and negative the cloning algorithm behaves as in the harmonic problem: there are no divergences but a higher number of clones would be required for a quantitative estimate of the SCGF.

Zoom In Zoom Out Reset image size

We have also studied the large deviations of the active work for an AOUP confined in a circular well, which we have modelled with a Weeks–Chandler–Andersen (WCA) potential on the difference between the distance of the AOUP from the origin r and the radius R of the circular well, i.e.

$$\begin{equation}{U}_{\text{circle}}(\mathbf{r})=\begin{cases}0\quad & \text{for}\enspace r< R-{2}^{\frac{1}{6}}\sigma \\ {U}_{\text{LJ}}(R-r)-{U}_{\text{LJ}}({2}^{1/6}\sigma )\quad & \text{for}\enspace r\geqslant R-{2}^{1/6}\sigma \end{cases}\end{equation} \tag{ 33 }$$

with U_LJ(x) the Lennard-Jones potential,

$$\begin{equation}{U}_{\text{LJ}}(x)=4{\epsilon}\left[{\left(\frac{\sigma }{x}\right)}^{12}-{\left(\frac{\sigma }{x}\right)}^{6}\right].\end{equation} \tag{ 34 }$$

We set both the spatial (σ) and energy () scales of the Lennard-Jones potential to 1. In practical terms, the AOUP moving in the potential U_circle( r ) is free until its distance from the origin reaches R − 2^1/6, then it is pushed back by the ascending branch of the Lennard-Jones potential. Also for this potential the steady-state distribution of the AOUP position is either annular or peaked at the origin depending on the Péclet number and we focus on the regime where it is annular. Our estimates of SCGF and rate function are shown in figure 8 for R = 5σ. The features of SCGF and rate function are similar to those observed with the stiff potential, therefore we will not comment any further.

Zoom In Zoom Out Reset image size

Let us close the section by examining the biased trajectories for the anharmonic potentials considered in this section, reported in figure 9 (left panel for the stiff potential, right panel for the circular well). With respect to the harmonic confinement, here the interaction with the potential is crucial for the realisation of both typical and atypical trajectories: this is particularly evident when looking at the biased trajectory of the circular-well problem (right panel of the figure). The typical trajectories (starting from the green dots in the figure) mirror the annular structure of the steady-state distribution of the position. From a dynamical perspective, the AOUP performs a persistent motion along the angular direction, switching between the clockwise and anticlockwise directions. These trajectories still produce a positive active work, although smaller than in the free case. The trajectories with negative λ (starting from the red dots in the figure) are those that produce approximately vanishing active work and they do so by having an unusually persistent active noise, so that the particle sits at a distance from the centre such that the active force and potential restoring force are balanced while performing thermal fluctuations. By contrast, a larger-than-average active work can be produced when the active noise is less persistent than usual. Having, for instance, a sudden change in the direction of the active noise when the particle is pushing against the potential slope would cause the particle to revert its direction and slide down that same slope, thus gaining velocity. It is interesting to notice that the mechanisms leading to higher- or lower-than-average active work in the presence of an anharmonic radially symmetric confining potentials are opposite with respect to those at work in the harmonic and free problems: here highly correlated active noise is required to produce a small active work and anti-correlated active noise results in large active work, whereas, in the free problem, unusually correlated active noise results in higher-than-average active work and anti-correlations cause a reduction of the active work.

Zoom In Zoom Out Reset image size

Fluctuation relations have been formulated as general theorems fixing the symmetry of the distribution of observables Ω(t) which depend on the system trajectory. These theorems were first proved in the context of non equilibrium diffusion processes, where equilibrium systems are perturbed by external fields, and for jump processes that satisfy time reversal invariance and markovianity [69–73]. Later on, fluctuation relations received a more systematic formulation [51], and have also been considered in the framework of active matter, using the entropy production as significative thermodynamic observable [32–34, 40, 74].

A (stationary) detailed fluctuation theorem, to which we will refer in the following as fluctuation relation (FR), is said to be satisfied when the time-independent stationary distribution Π(ω) of the functional observable ω = Ω/τ satisfies the identity

$$\begin{equation}\frac{1}{\tau }\enspace \mathrm{log}\left(\frac{{\Pi}(\omega )}{{\Pi}(-\omega )}\right)=c\omega ,\end{equation} \tag{ 35 }$$

with c a constant depending in general on the system and on the observable considered. When the stationary distribution of ω(t) satisfies also a large deviation principle, i.e. Π(ω) ≡ e^−τI(ω) with I(ω) rate function, the FR takes the following asymptotic form as τ → ∞,

$$\begin{equation}I(\omega )-I(-\omega )=c\omega .\end{equation} \tag{ 36 }$$

As such, equation (36) lacks of the subexponential contributions appearing in general also in the stationary configuration distribution and negligible in the large time limit. The validity of an asymptotical FR can actually be traced down to the following property of the SCGF [69],

$$\begin{equation}g(\lambda )=g(c-\lambda ).\end{equation} \tag{ 37 }$$

In this section we check the validity of the FR for the active work (per unit time) distributions in the free-particle case and in the presence of confining potentials. Note that the direct check of the FR through the rate function, using equation (36), requires the underlying numerical distribution to be significantly sampled for negative values. This condition is not satisfied for all the cases considered so far in this work, thus we resort to simulations with a reduced magnitude of the activity (parameters values are indicated in the caption of figures 10 and 11). Alternatively, we provide also an indirect check of the FR by using the SCGF data shown figures 3 and 4, obtained with the previous choice of parameters trough the cloning algorithm (see appendix A), in equation (37).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the free-particle case, a simple calculation shows that the SCGF equation (19) satisfies the symmetry property equation (37) with c = −T⁻¹, in fact

$$\begin{equation}g\left(-\frac{1}{\mathrm{T}}-\lambda \right)=\frac{d}{2}\left({\gamma }_{\mathrm{R}}-\sqrt{{\gamma }_{\mathrm{R}}^{2}-4{D}_{\mathrm{R}}^{\prime }\gamma \left(1+\gamma {D}_{\mathrm{T}}\lambda \right)\lambda }\right)=g(\lambda ).\end{equation} \tag{ 38 }$$

Indeed, figure 10(c) shows that the SCGF g(λ) taken from figure 3 and g(−λ − T⁻¹) overlap well, apart from the last points on the left-hand and on the right-hand side of the figure which are affected by the convergence problems detailed in subsection 2.3. Therefore, we can assert that an FR for the free-particle active work is satisfied with slope c = −T⁻¹. A direct proof of the validity of this theorem can also be obtained using the rate function equation (30): after some algebra, we find the result I(w) − I(−w) = −w/T, valid for every choice of parameters. For the sake of completeness, figures 10(a) and (b) provide numerical support to this result. In particular, figure 10(a) shows the rate functions of the active work evaluated through direct sampling setting a lower activity, as indicated in the caption, while figure 10(b) shows the FR at various times computed using rate functions from (a).

The equivalence between the constant c and the negative inverse bath temperature holds for a large class of models and dynamical observables [51]. Related to these observables there is the entropy production (where c = 1), which, let us remark, is proportional to the active work defined in equation (2) provided the self-propulsion velocity v (t) is considered even under time reversal transformation [23, 27]. In general, one can define a fluctuation-relation temperature T_FR such that $c=-{T}_{\text{FR}}^{-1}$, and compare it with other quantities such as the kinetic temperature T_kin, defined from the equipartition theorem, and the effective temperature T_eff(t) defined from the fluctuation–dissipation theorem [52–54]. Hereafter we consider the value of T_eff(t) measured in the long time limit, where it reaches a constant value, so we drop the explicit time dependence (see appendix B for more details). As reported in the appendix, we can compute these quantities analytically in the free AOUP case, and find that their values grow quadratically with the Péclet number with different functional forms, making them distinguishable with respect to each other and to the bath temperature. For the free AOUP, we have that T_FR = T.

Concerning the confined AOUP, we resort to numerics to measure the T_FR for some representative cases and compare it with the values of T, T_kin and T_eff. For the harmonic potential, we checked directly the FR using a lower activity than that reported in figures 4 and 5, in order to be able to sample also negative values of the active work (figure 10(d)) and evaluate equation (35). From the results reported in figure 10(e) it is evident that, similarly to the free-particle case, the slope of the curves is compatible with T_FR = T. We remark that in this low-activity case the bath temperature T = 0.05 is essentially indistinguishable from the kinetic and effective temperatures T_kin ≃ 0.0500 and T_eff ≃ 0.0505, whose exact asymptotic expressions are reported in appendix B. In the high-activity case we can estimate T_FR by using the SCGF of figure 4. Figure 10(f) shows that the SCGF satisfies the same symmetry property of the free-particle case g(λ) = g(−λ − T⁻¹) near the function minimum, so that T_FR = T = 0.05. In the high-activity case the three temperatures are now perfectly distinguishable (T = 0.05, T_kin ≃ 0.1804 and T_eff ≃ 87.1199), and thus we conclude that only the bath temperature is compatible with T_FR.

For the stiff potential U_stiff(r) = k_stiff r¹⁰/10 we report in figure 11 the results of the analysis of the FR for two significative choices of parameters, giving rise respectively to a center-peaked stationary position distribution (low-activity) and to an annular one (high-activity), as described in subsection 3.2. In figures 11(a) and (b), the choice is such that the position distribution of the confined AOUP is peaked at the centre, as in figure 6 (centre). In this case, as it can be seen from figure 11(b), an FR is again satisfied with a slope compatible with T_FR = T = 0.05, which is indistinguishable from T_kin ≃ 0.0499(±0.01%), but lower than T_eff ≃ 0.0717(±0.01%) (both temperatures are measured numerically in this case). For higher activity (figures 11(c) and (d)), the position distribution of the confined AOUP is annular with finite radius, as in figure 6 right panel. In this case, the FR reported in (d) is again satisfied, but with a T_FR ≃ 0.4887(±4.9%) much larger than the bath temperature T = 0.05. Here, a comparison with T_kin ≃ 0.0444(±0.001%) and T_eff ≃ 0.4775(±0.01%) shows that T_FR = T_eff within numerical accuracy. This difference between low and high activity can be justified by the fact that in the former case the particle dynamics takes place near the potential minimum, while in the latter case it does not.

In this paper we have studied the large deviations of the time-averaged power injected by the self-propulsion force—or active work—for an AOUP interacting with a confining potential. We have examined, in particular, four cases: the free particle case, without any confining potential; the harmonic case, with a quadratic potential; two anharmonic potentials, namely a 'stiff' potential growing as r ¹⁰ and a circular well modelled with a WCA potential.

For the free-particle case, we have obtained the rate function of the active work in general dimension d exactly from saddle-point calculations of the inverse Laplace transform of the CGF. The corresponding SCGF was recently computed for d = 2 in [39] by solving a tilted eigenvalue problem. The rate function of the free problem can also be obtained from the SCGF via LFT. This might lead to erroneous results when the cumulant generating function displays diverging sub-exponential contributions [26]. Our analysis takes such contributions into account, but we show that, at variance with the corresponding passive problem [26], the interplay between thermal and active noise causes the non-analyticities stemming from sub-exponential contributions to be covered by those of the SCGF, showing that estimating the rate function as a LFT of the SCGF yields the correct result.

For the confining potentials, we estimated the RF and the SCGF numerically. For a quadratic confining potential, we find the development of linear tails and a similar phenomenology of the free-particle case. These tails suggest the presence of a typical condensation transition due to a sticking of the saddle-point mechanism [75–79], which would lead to a rate function growing linearly with the active work after a certain threshold—further analytical and numerical investigations will be needed to demonstrate this point.

For a general nonlinear confining potential, the AOUP will push against the slope of the potential in steady state when the Péclet number is sufficiently high. We characterized the SCGF of the active work for two representative nonlinear potentials and found that they do not show any singular behaviour despite accumulation effects at the boundaries. We found that lower-than-average fluctuations of the active work correspond to particles pushing against the potential and are realised with highly correlated kicks always pointing against the direction of the potential slope. Higher-than-average fluctuations correspond instead to particles which move in circles for most of the observation time. It would be of interest, in this respect, to extend the analysis to ensembles of interacting AOUPs. Such systems are known to display similar macroscopic average properties to those of system of interacting active Brownian particles [23], thus it is a natural question whether the analogy extends to features related to rare fluctuations [20, 21].

Finally, we have studied the fluctuation relation for active work fluctuations in three different situations: the free-particle case, the harmonically confined AOUP and the AOUP confined by the stiff anharmonic potential. In particular, we defined a typical temperature T_FR stemming from the fluctuation relation slope. For the free-particle case, we proved analytically that the fluctuation relation is satisfied with T_FR = T, and supported the results with numerical simulations. For the harmonically confined AOUP, we found numerically that the fluctuation relation is satisfied with T_FR = T, and that this result is independent of the activity. For the anharmonically confined AOUP, we found that low-activity AOUP satisfies the fluctuation relation with T_FR = T, while high-activity AOUP have a much higher value of T_FR than the bath temperature, which is equal to an effective temperature T_eff estimated independently using the fluctuation–dissipation theorem. This increase in the effective temperature corresponds to giving rise to an annular position distribution instead of a center-peaked one when increasing the activity.

These results could be connected in some way to recent works showing that in the AOUP model the entropy production rate is different from zero only when the particle is under the action of a more-than-quadratical potential [32–34]. In order to better understand this possible connection, the role of the effective temperature and the general behaviour of active systems, a more in-depth analysis of fluctuation theorems for active systems could be of interest. We then leave as a future possible work further extensions of our results to analytical and numerical checks of fluctuation relations.

G G acknowledges MIUR for funding (PRIN 2017WZFTZP).

The first method we applied consists in performing numerical integration in two dimensions of the equations

$$\begin{equation}m\ddot{\boldsymbol{r}}(t)=-\gamma \dot{\boldsymbol{r}}(t)+\gamma \boldsymbol{v}(t)-\nabla U(\boldsymbol{r}(t),t)+\gamma \sqrt{2{D}_{\mathrm{T}}}\boldsymbol{\xi }(t),\end{equation} \tag{ 39 }$$

which includes the inertial term $m\ddot{\boldsymbol{r}}(t)$, with m the mass of the particle and in which different potentials expressions are considered. The equations are integrated using the Vanden-Eijnden–Ciccotti algorithm [80]. The parameter choice is different, as indicated in the figures caption in each section, but in general we choose the mass m = 1 and the friction coefficient γ in such a way that the ratio m/γ ≪ 1 and the system can always be considered as effectively overdamped for timescales larger than the inertial time t_I = m/γ = 0.1. We use as integrating timestep dt = 0.01, and evolve the system up to a maximum observation time of order τ ∼ 10³, or 10⁵ dt, computing the active work directly from the definition equation (2). The empirical distributions of the active work are obtained by integrating the dynamics for N_c = 10⁵ independent realisations of the stochastic forces, with initial conditions sampled from equilibrated systems evolved for τ_eq = 10² ≫ t_I .

The second method we use to estimate the rate function is the cloning algorithm, an algorithm based on importance sampling ideas and described briefly as follows. A large number N_c of copies, or clones, of the system is evolved simultaneously using the same integration scheme described in the previous paragraph. Copies are then cloned or pruned depending on the value of the active work so as to generate a biased ensemble with a lower- or higher-than-usual average active work. Specifically, we divide the entire simulation time interval τ in M subintervals such that τ = MΔt, with Δt = 1 (corresponding to 100 timesteps). We compute the active work of each copy after every time interval Δt and we register the values ${\Delta}{W}_{t}^{a}={W}_{t+{\Delta}t}^{a}-{W}_{t}^{a}$, with a the copy index, together with the weighted sum over all copies ${G}_{t}(\lambda )={\sum }_{b}{\text{e}}^{\lambda {\Delta}{W}_{t}^{b}}$. λ is the biasing parameter, coinciding with the argument of the cumulant generating function. Thus an integer score n^a is assigned to each copy according to

$$\begin{equation}{n}_{a}=\lfloor \frac{{\text{e}}^{\lambda {\Delta}{W}_{t}^{a}}}{{G}_{t}(\lambda )}{N}_{\text{c}}+\zeta \rfloor ,\end{equation} \tag{ 40 }$$

where ζ is a random number uniformly distributed in [0, 1] and ⌊.⌋ denotes the lower integer part. For λ positive, trajectories with a higher ${\Delta}{W}_{t}^{a}$ will have a higher score. Each copy a with positive score is then cloned so that it appears n_a times, whereas copies with n_a = 0 are pruned. After this step, copies are deleted or cloned at random in order to keep the number of copies constant. After repeating this evolve-and-clone procedure for the M steps, the SCGF can be estimated as

$$\begin{equation}{g}_{\tau }(\lambda )=\frac{1}{\tau }\enspace \mathrm{log}\hat{{\Pi}}(\lambda )\simeq \frac{1}{(M-1){\Delta}t}\sum\limits _{t=0}^{M-1}\mathrm{log}\left(\frac{{G}_{t}(\lambda )}{{N}_{\text{c}}}\right).\end{equation} \tag{ 41 }$$

In this appendix we recall how the effective and kinetic temperatures can be defined on the basis of the fluctuation–dissipation relation and of the equipartition theorem (see e.g. [52–54]). Concerning the effective temperature, the starting point is just the equilibrium fluctuation–dissipation relation

$$\begin{equation}2T\chi ({t}^{\prime },t)={{\Delta}}^{2}({t}^{\prime },t),\end{equation} \tag{ 42 }$$

where

$$\begin{equation*}{{\Delta}}^{2}({t}^{\prime },t)=\langle {[\mathbf{r}(t)-\mathbf{r}({t}^{\prime })]}^{2}\rangle \end{equation*}$$

is the total mean square displacement and

$$\begin{equation*}\chi ({t}^{\prime },t)={\int }_{{t}^{\prime }}^{t}\mathrm{d}{t}^{{\prime\prime}}\sum\limits _{\alpha =1}^{d}{R}_{\alpha \alpha }({t}^{{\prime\prime}},t)\end{equation*}$$

is the integrated linear response, with

$$\begin{equation*}{R}_{\alpha \beta }({t}^{\prime },t)={\left.\frac{\delta {\langle {r}_{\alpha }(t)\rangle }_{h}^{\lambda }}{\delta {h}_{\beta }^{\lambda }({t}^{\prime })}\right\vert }_{{h}_{\beta }^{\lambda }=0}\end{equation*}$$

linear response of the system, α, β dimension indices, d dimension of the system and ${h}_{\beta }^{\lambda }({t}^{\prime })$ an external perturbation depending on the parameter λ. The idea is to exploit equation (42) to define in the active out-of-equilibrium system the time dependent effective temperature

$$\begin{equation}{T}_{\text{eff}}({t}^{\prime },t)=\frac{{{\Delta}}^{2}({t}^{\prime },t)}{2\chi ({t}^{\prime },t)}.\end{equation} \tag{ 43 }$$

The kinetic temperature is instead defined exploiting the velocity fluctuations and the equipartition theorem. From the latter, we can in fact write for each degree of freedom i of a system that

$$\begin{equation*}\frac{1}{2}m\langle {\dot{\mathbf{r}}}_{i}^{2}(t)\rangle =\frac{1}{2}{k}_{\text{B}}{T}_{\text{kin}}(t),\end{equation*}$$

where k_B is the Boltzmann constant and ⟨...⟩ denotes an ensemble mean. From this expression one simply obtains the kinetic temperature definition as

$$\begin{equation}{T}_{\text{kin}}(t)=\frac{m\langle {\dot{\mathbf{r}}}_{i}^{2}(t)\rangle }{{k}_{\text{B}}}.\end{equation} \tag{ 44 }$$

In the free-particle case, we can use the expressions for the mean square displacement and the mean value of the squared velocity reported in [81], and the computed integrated response function

$$\begin{equation*}\chi ({t}^{\prime },t)=\frac{t-{t}^{\prime }}{\gamma },\end{equation*}$$

with γ the bath friction coefficient. We find that in the limit t → ∞

$$\begin{equation*}{T}_{\text{kin}}(t)\to T+{\left(\frac{{F}_{a}}{m}\right)}^{2}\frac{{D}_{\mathrm{R}}^{\prime }}{{\gamma }_{\mathrm{R}}}\frac{1}{\gamma }\frac{1}{\left({\gamma }_{\mathrm{R}}+\frac{\gamma }{m}\right)}\end{equation*}$$

and

$$\begin{equation*}{T}_{\text{eff}}(t)\to T+{\left(\frac{{F}_{a}}{m}\right)}^{2}\frac{{D}_{\mathrm{R}}^{\prime }}{{\gamma }_{\mathrm{R}}}\frac{1}{{\gamma }^{2}{\gamma }_{\mathrm{R}}},\end{equation*}$$

with t' = 0. In these formulas both T_kin(t) and T_eff(t) reach a constant value for large enough times and have both a quadratic dependence on the Pe but with a different functional form, such that they can be distinguished with respect to each other and from the bath temperature T when Pe is sufficiently high.

We can provide the expression for T_eff(t) and T_kin(t) also in the harmonically confined case (U(r) = k r²/2). Using the expressions for the mean square displacement and the mean value of the squared velocity reported in [81] for a harmonically confined AOUP, and the integrated linear response

$$\begin{equation*}\chi ({t}^{\prime },t)=\frac{1-{\text{e}}^{-\frac{\gamma }{k}({t}^{\prime }-t)}}{k},\end{equation*}$$

with γ the bath friction coefficient and k the elastic constant of the harmonic potential, we find that in the limit t → ∞

$$\begin{equation*}{T}_{\text{kin}}(t)\to T+{\left(\frac{{F}_{a}}{m}\right)}^{2}\frac{{D}_{\mathrm{R}}^{\prime }}{{\gamma }_{\mathrm{R}}}\frac{1}{\gamma }\frac{1}{\left({\gamma }_{\mathrm{R}}+\frac{\gamma }{m}+\frac{k}{m{\gamma }_{\mathrm{R}}}\right)}\end{equation*}$$

and

$$\begin{equation*}{T}_{\text{eff}}(t)\to T+{\left(\frac{{F}_{a}}{m}\right)}^{2}\frac{{D}_{\mathrm{R}}^{\prime }}{{\gamma }_{\mathrm{R}}}\frac{1}{\gamma }\frac{1+\frac{\gamma }{m{\gamma }_{\mathrm{R}}}}{\left({\gamma }_{\mathrm{R}}+\frac{\gamma }{m}+\frac{k}{m{\gamma }_{\mathrm{R}}}\right)},\end{equation*}$$

with t' = 0. Similarly to the free-particle case, both temperatures reach a constant value for large times.

In the anharmonically confined case (U_stiff(r) = k_stiff r¹⁰/10), we estimated T_eff(t) and T_kin(t) numerically. The former was estimated directly using the definition equation (44). The latter was estimated from the definition equation (43), measuring independently Δ²(t) and χ(t). Δ²(t) was measured as the x-component mean square displacement, while for χ(t) we applied a constant force h in the x-direction. We chose h = 0.2 in order for the force to be small enough to remain in the linear regime and high enough to overcome the large fluctuation effects. In figure 12 (left), we report χ(t) and Δ²(t) in a high-activity case with F_a = 10.0, γ = 100, T = 0.05. Note that these functions reach a constant value after a few permanence times ${\gamma }_{\mathrm{R}}^{-1}$. In figure 12 (right), we report instead T_eff(t). Notice that for times smaller than the persistence time, $t\ll {\gamma }_{\mathrm{R}}^{-1}$, T_eff(t) = T, with T the bath temperature. For $t\gg {\gamma }_{\mathrm{R}}^{-1}$, instead, we reach a fitted constant value of T_eff(t) ≃ 0.4775(±0.01%). T_kin is instead constant over time and has a fitted value of T_kin ≃ 0.0444(±0.001%). We also estimated the two temperatures for a low-activity case, with F_a = 0.5, γ = 10, T = 0.05. In this case the numerical measurements yield T_kin ≃ 0.0499(±0.01%) and for large times T_eff ≃ 0.0717(±0.01%).

Zoom In Zoom Out Reset image size