Perturbative Calculation of Quasi-Potential in Non-equilibrium Diffusions: A Mean-Field Example

Abstract

In stochastic systems with weak noise, the logarithm of the stationary distribution becomes proportional to a large deviation rate function called the quasi-potential. The quasi-potential, and its characterization through a variational problem, lies at the core of the Freidlin–Wentzell large deviations theory (Freidlin and Wentzell, Random perturbations of dynamical systems, 2012). In many interacting particle systems, the particle density is described by fluctuating hydrodynamics governed by Macroscopic Fluctuation Theory (Bertini et al., arXiv:1404.6466, 2014), which formally fits within Freidlin–Wentzell’s framework with a weak noise proportional to $1/\sqrt{N}$, where N is the number of particles. The quasi-potential then appears as a natural generalization of the equilibrium free energy to non-equilibrium particle systems. A key physical and practical issue is to actually compute quasi-potentials from their variational characterization for non-equilibrium systems for which detailed balance does not hold. We discuss how to perform such a computation perturbatively in an external parameter $\lambda$, starting from a known quasi-potential for $\lambda =0$. In a general setup, explicit iterative formulae for all terms of the power-series expansion of the quasi-potential are given for the first time. The key point is a proof of solvability conditions that assure the existence of the perturbation expansion to all orders. We apply the perturbative approach to diffusive particles interacting through a mean-field potential. For such systems, the variational characterization of the quasi-potential was proven by Dawson and Gartner (Stochastics 20:247–308, 1987; Stochastic differential systems, vol 96, pp 1–10, 1987). Our perturbative analysis provides new explicit results about the quasi-potential and about fluctuations of one-particle observables in a simple example of mean field diffusions: the Shinomoto–Kuramoto model of coupled rotators (Prog Theoret Phys 75:1105–1110, [74]). This is one of few systems for which non-equilibrium free energies can be computed and analyzed in an effective way, at least perturbatively.

Appendices

Appendix 1: Smooth Stationary Solutions of the Shinomoto–Kuramoto Model at $T=0$

In Sect. 4.3.4, we have discussed simple delta-like stationary solutions of the McKean–Vlasov dynamics for the Shinomoto–Kuramoto model at $T=0$. They described the particles accumulated at the fixed point of single particle drift. These are however not the only stationary solutions for vanishing temperature and also smooth solutions are present, which, however, are linearly unstable in the ferromagnetic model considered in this paper ($J>0$). We briefly discuss such smooth stationary solutions at $T=0$ and their stability in this appendix.

For $T=0$, Eq. (4.36) implies that

$$\begin{aligned} \rho _{inv}(\theta +\theta _0)=\frac{\sqrt{1-\left( \frac{y}{F}\right) ^2}}{2\pi \left( 1-\frac{y}{F}\sin \theta \right) } \end{aligned}$$

(6.1)

provided that $0\le y<F$ and the self-consistency relation (4.43) holds taking now the form

$$\begin{aligned} \left( \frac{y}{J}\right) ^2+\bigg (\frac{_{1-\sqrt{1-\left( \frac{y}{F}\right) ^2}}}{^{\frac{y}{F}}}\bigg )^{2}=\left( \frac{h}{J}\right) ^2. \end{aligned}$$

(6.2)

The left hand side of the last equation grows monotonically from 0 at $y=0$ to $(F/J)^2+1$ at $y=F$. Hence there is a unique solution for $0\le y<F$ if and only if $F^2+J^2> h^2$. We infer that there is a unique stationary smooth solution at $T=0$ when the latter constraint holds and none otherwise.

To examine the stability of such solutions, we examine the eigenfunctions $\delta \rho _{\mu }(\theta +\theta _0)$ with vanishing integral of the linearized Fokker–Planck operator $R_{\rho _{inv}}$ of (4.44) at $T=0$. We may write $\delta \rho _{\mu }(\theta +\theta _0) =\partial _{\theta }\delta f_{\mu }(\theta )$, where $\delta f_{\mu }$ is periodic on $[0,2\pi ]$. The eigenequation takes the form

$$\begin{aligned}&\big (R_{\rho _{inv}}\partial _{\theta }\delta f_{\mu }\big )(\theta )\nonumber \\&\quad =-\partial _{\theta }\Big ((F-y\sin \theta )\, \partial _{\theta }\delta f_{\mu }(\theta ) - Y_1\,\rho _{inv}(\theta +\theta _0)\sin \theta + Y_2\,\rho _{inv}(\theta +\theta _0) \cos \theta \Big )\nonumber \\&\quad =\mu \,\partial _{\theta }\delta f_{\mu },\,\ \end{aligned}$$

(6.3)

where

$$\begin{aligned} Y_1=J\int _0^{2\pi }\sin \theta \,\,\delta f_{\mu }(\theta )\,d\theta ,\qquad Y_2=-J\int _0^{2\pi }\cos \theta \,\,\delta f_{\mu }(\theta )\,d\theta . \end{aligned}$$

(6.4)

We may now define a new operator S such that

$$\begin{aligned} \partial _{\theta }\big (S\,\delta f\big ) =R_{\rho _{inv}}\partial _{\theta }\delta f_{\lambda } . \end{aligned}$$

(6.5)

Explicitly,

$$\begin{aligned} \big (S\,\delta f\big )(\theta )=-(F-y\sin \theta )\,\partial _{\theta }\delta f(\theta ) - Y_1\rho _{inv}(\theta +\theta _0)\,\sin \theta + Y_2\rho _{inv}(\theta +\theta _0) \,\cos \theta .\nonumber \\ \end{aligned}$$

(6.6)

Clearly, the diagonalization of $R_{\rho _{inv}}$ on the subspace of functions $\delta g$ with vanishing integral is equivalent to the diagonalization of S.

A direct calculation shows that the 3-dimensional subspace $V_3$ spanned by functions

$$\begin{aligned} f_1(\theta )=\rho _{inv}(\theta +\theta _0)\,\sin \theta ,\qquad f_2(\theta )=\rho _{inv}(\theta +\theta _0)\,\cos \theta ,\qquad f_3(\theta )=\rho _{inv}(\theta +\theta _0)\nonumber \\ \end{aligned}$$

(6.7)

is invariant for S. In particular, the representation of S on such subspace is given by the matrix

$$\begin{aligned} \left( \begin{array}{ccc} J\,\frac{1-\sqrt{1-\big (\frac{y}{F}\big )^2}}{\big (\frac{y}{F}\big )^2} &{} F &{} J\,\frac{1-\sqrt{1-\big (\frac{y}{F}\big )^2}}{\frac{y}{F}} \\ -F &{} J\Bigg (1-\frac{1-\sqrt{1-\big (\frac{y}{F}\big )^2}}{\big (\frac{y}{F}\big )^2}\Bigg ) &{} -y \\ 0 &{} -y &{} 0 \end{array} \right) . \end{aligned}$$

(6.8)

A little algebra shows that the above matrix has one zero eigenvalue which corresponds to the eigenvector $(-y,0,F)$, i.e. to the constant function

$$\begin{aligned} \delta f(\theta )\,=\,\rho _{inv}(\theta +\theta _0)\,(F-y\sin {\theta })\,=\, \frac{\sqrt{F^2-y^2}}{2\pi } . \end{aligned}$$

(6.9)

In this case $\delta g=\partial _\theta \delta f=0$. Consequently, the zero eigenvalue does not occur for the operator $R_{\rho _{inv}}$. The other two eigenvalues of the matrix (6.7) are given by

$$\begin{aligned} \lambda _\pm \,=\,\frac{J}{2}\,\pm \,\sqrt{\left( \frac{J}{2}\right) ^2\,-\,J^2 \frac{_{\sqrt{1-\big (\frac{y}{F}\big )^2}}}{^{\Big (1+\sqrt{1-\big (\frac{y}{F}\big )^2}\Big )^2}}\,-\,F^2+y^2} . \end{aligned}$$

(6.10)

of which at least one has positive real part if $J>0$. It follows that for $J>0$, the stationary solution (6.1) is linearly unstable.

Appendix 2: Algorithm for the Taylor Expansion of the Quasi-Potential Close to a Stationary State

We discuss briefly in this Appendix how the Taylor expansion of the quasi-potential of the Shinomoto–Kuramoto model may be numerically calculated. This will be possible, in principle, at all the orders of the Taylor expansion.

Let us first concentrate on the quadratic term $\mathcal {F}^{(0)}$. In order to explicitly calculate kernel $\varphi ^0(\theta ,\vartheta )$, we proceed similarly as in Sect. 4.5.3 and consider the eigenfunctions of operators $R_{\rho _{inv}}^\dagger$ and $R_{\rho _{inv}}$ acting in space $H_0$:

$$\begin{aligned} R_{\rho _{inv}}^\dagger u_{k}=\alpha _k\,u_{k}\qquad \text {and}\qquad R_{\rho _{inv}} v_{k}=\overline{\alpha }_k\,v_{k} . \end{aligned}$$

(7.1)

Assuming $\alpha _k$ to be different, we impose the orthogonality relations (4.99). Eigenfunctions $u_k$ form a basis of $H_0$ but they are not orthogonal. Similarly for $v_k$. The kernel of the quadratic term of the quasi-potential may be represented as

$$\begin{aligned} \varphi ^0(\theta ,\vartheta )=\sum _{k,l}\Phi _{kl}\,v_k(\theta )\, \overline{v_l(\vartheta )}\, \end{aligned}$$

(7.2)

where $\Phi _{kl}$ is defined by

$$\begin{aligned} \Phi _{kl}= \int _{0}^{2\pi }\overline{u_k(\theta )}\,d\theta \int _{0}^{2\pi } \varphi ^0(\theta ,\vartheta )\,u_l(\vartheta )\,d\vartheta =\langle u_k|\Phi u_l \rangle . \end{aligned}$$

(7.3)

The problem is thus reduced to the calculation of $\Phi _{kl}$. Now, from Eq. (4.122), we can easily compute the matrix elements of $\Phi ^{-1}$ in the basis $u_k$. Indeed, from Eq. (4.122),

$$\begin{aligned} (\Phi ^{-1})_{kl}\equiv \langle u_k|\Phi ^{-1}u_l\rangle = - \frac{2k_BT}{\alpha _k^*+\alpha _l}\int d\theta \, \overline{(\partial _\theta u_k)(\theta )}\,\rho _{inv}(\theta )\, (\partial _\theta u_l)(\theta )\,. \end{aligned}$$

(7.4)

Note, however, that the matrix $(\Phi ^{-1} )_{kl}$ is not the inverse of $\Phi _{kl}$ because the basis formed by $u_k$ is not orthonormal. This problem can be handled with simple linear algebra. We introduce the matrix of scalar products

$$\begin{aligned} P_{kl}=\langle u_k|u_l\rangle \end{aligned}$$

(7.5)

whose inverse will be denoted by $(P^{-1})_{kl}$. Then, we construct another matrix

$$\begin{aligned} (B^{-1})_{kl}\equiv \sum _i(P^{-1})_{ki}\,(\Phi ^{-1})_{il}\, \end{aligned}$$

(7.6)

with inverse $B_{kl}$. Finally,

$$\begin{aligned} \Phi _{kl}=\sum _i P_{ki}\,B_{il} . \end{aligned}$$

(7.7)

Once $\Phi _{kl}$ is known, it is straightforward to write the kernel $\varphi ^0$ in the real space using Eq. (4.126).

We conclude by observing that, once $\varphi ^0$ is known, one could also numerically evaluate the higher order kernels $\varphi ^n$ of the Taylor expansion by solving Eq. (4.126) using the fact that the basis $\Phi v_k$ is composed of eigenstates of operators $K^{(0)\dagger }_r=-\Phi \,R_{\rho _{inv}}\Phi ^{-1}$, see Eqs. (4.123) and (7.1).