Physically consistent numerical solver for time-dependent Fokker-Planck equations

We present a simple thermodynamically consistent method for solving time-dependent Fokker-Planck equations (FPE) for overdamped stochastic processes, also known as Smoluchowski equations. It yields both transition and steady-state behavior and allows for computations of moment-generating and large-deviation functions of observables defined along stochastic trajectories, such as the fluctuating particle current, heat, and work. The key strategy is to approximate the FPE by a master equation with transition rates in configuration space that obey a local detailed balance condition for arbitrary discretization. Its time-dependent solution is obtained by a direct computation of the time-ordered exponential, representing the propagator of the FPE, by summing over all possible paths in the discretized space. The method thus not only preserves positivity and normalization of the solutions but also yields a physically reasonable total entropy production, regardless of the discretization. To demonstrate the validity of the method and to exemplify its potential for applications, we compare it against Brownian-dynamics simulations of a heat engine based on an active Brownian particle trapped in a time-dependent quartic potential.

Figure 1

Sketch of the configuration space discretization used for the numerical solution of the two-dimensional overdamped Fokker-Planck equation (1). The black points mark the states inside the domain $[x_{-},x_{+}]\times [y_{-},y_{+}]$, and the colored points form the boundary (see Sec. 3b). The black full arrows depict the allowed transitions with the “bulk” transition rates (13) or (18) (horizontal transition) and (14) or (19) (vertical transitions). The red boundary is reflecting and thus the particles cannot cross the red states (hence no red arrows). The blue boundary is absorbing and thus particles can leave the system from these sites (depicted by one-way dashed blue arrows). The states in the corners of the domain require two boundary conditions. In the figure, we impose reflecting boundary condition in the $y$ direction (depicted by red circumferences of the points) and periodic boundary conditions in the $x$ direction (depicted by green interiors of the points). The periodic boundary allows the particles to leave the system in the $x$ direction. The leaving particles reenter the system at the opposite side of the domain, as depicted by the green dot-dashed arrows.

Figure 2

Sketch of the phase-space discretization used for the numerical solution of the two-dimensional overdamped Fokker-Planck equation (1) in case of the driven active particle (Sec. 7). The meaning of the arrows and point colors is the same as in Fig. 1.

Figure 3

Active particle confined to a single dimension and driven by the quartic potential (66) of time-dependent strength $k\left(t\right)$.

Figure 4

The parameters of the microscopic heat engine consisting of the periodically driven active particle depicted in Fig. 3, during one period of the cyclic driving protocol. Boxes comprise maximum and minimum values of the corresponding variables during the cycle.

Figure 5

Comparison of observables for the periodically driven active particle, depicted as a function of time during one limit cycle, as computed from BD simulations (symbols) and the MNM (lines): (a) averages $〈x^2〉, 〈xcos\theta 〉$, and $〈x^4〉$ and (b) marginal probability density $\xi (x,t)$ for the particle position $x$ at five time instants $t$ during the cycle. The PDF $\xi (x,0)$ at the initial time 0 (dashed blue line) and $\xi (x,1)$ at the final time 1 (full orange line) in the first panel of (b) coincide, because the system operates in the limit cycle as described by Eq. (71).

Figure 6

Application of the MNM to moment-generating functions (MGFs): (a) the MGFs ${\chi }_w$ and ${\chi }_q$ for work and heat, (b) the first 11 raw moments, and (c) the first 10 cumulants [c)] of the net extracted work $-w$ and net supplied heat $q$ per cycle, as calculated from the MGFs depicted in panel (a) ($*$) and from corresponding BD simulations of $200\times {10}^6$ trajectories (box plots).

Figure 7

Box plots of relative differences between moments (77) [panels (a)] and cumulants (78) [panels (b)] for work (left) and heat (right) as computed using MNM and BD, respectively [see Figs. 6 and 6].

Figure 8

Large-deviation limit of heat and work distributions. BD simulations of the cumulative distributions ${\rho }_q$ and ${\rho }_w$ of the net heat $q$ supplied ($◯$) and the net work $-w$ extracted ($\square$) over 100 cycles show that the work distribution has converged to the common limiting form obtained from the MNM (superimposing lines, with the vertical line indicating the average), while the heat distribution has not (left panel). The Legendre-Fenchel transformed logarithms of the MGF of heat and work (right panel) elucidate the unequal convergence toward the common large-deviation function (LDF) $J_q\left(q\right)\sim J_w(-w)$ with the number of cycles $N=30$, 50, 100. While the work distribution ($\square$) has already converged for $N=30$, the heat distribution ($◯$) keeps evolving (top to bottom).