Stochastic thermodynamics of Brownian motion in a flowing fluid

For simplicity, we examine two-dimensional Brownian dynamics in a fluid with time-independent flow. Generalization to three-dimensional time-dependent flow is straightforward. The velocity field of the fluid is

$${\bm{v}}({\bm{x}})=v_{x}(\bm{x})\,\hat{\bm{e}}_{x}+v_{y}(\bm{x})\,\hat{\bm{e}}_{y},$$

(1)

where $\hat{\bm{e}}_{x},\hat{\bm{e}}_{y}$ are respectively the unit vectors along $x$ and $y$ directions, and $\bm{x}=x\,\hat{\bm{e}}_{x}+y\,\hat{\bm{e}}_{y}$. Assuming that the Brownian particle is further confined by an external potential $V(\bm{x})$, its motion can be described by the following over-damped Ito-Langevin equations:

$$\begin{split}-\gamma(dx-v_{x}dt)-\partial_{x}Vdt+\sqrt{2\gamma T}\,dW_{x}&=0,\\ -\gamma(dy-v_{y}dt)-\partial_{y}Vdt+\sqrt{2\gamma T}\,dW_{y}&=0,\end{split}$$

(2)

where $\gamma$ is the friction constant, $T$ is the temperature, and $dW_{x},dW_{y}$ are the standard Wiener noises, which have the following basic properties:

	$\displaystyle\langle dW_{i}\rangle$	$\displaystyle=$	$\displaystyle 0,$		(3a)
	$\displaystyle\langle dW_{i}dW_{j}\rangle$	$\displaystyle=$	$\displaystyle dt\,\delta^{ij}.$		(3b)

Note that the first term in each of Eqs. (2) is the friction force multiplied by $dt$.

Equations (2) can be rewritten into:

$\displaystyle dx^{i}+\frac{T}{\gamma}(\partial_{i}U^{0}-\varphi_{i}^{0})dt=\sqrt{\frac{2T}{\gamma}}dW_{i},$

(4)

where $U^{0},{\bm{\varphi}}^{0}$ are given respectively by

	$\displaystyle U^{0}(\bm{x})$	$\displaystyle=$	$\displaystyle\beta\,V(\bm{x})+C^{0},$		(5a)
	$\displaystyle{\bm{\varphi}}^{0}(\bm{x})$	$\displaystyle=$	$\displaystyle\beta\gamma\,{\bm{v}}(\bm{x})=\varphi^{0}_{i}(\bm{x})\hat{\bm{e}}_{i},\quad$		(5b)

where $C^{0}$ is an irrelevant normalization constant. We do not need to distinguish superscripts from subscripts because we will only use Cartesian coordinate systems. Equation (4) is a special case of the following covariant nonlinear Ito-Langevin equation with non-conservative forces [21] (with all repeated indices summed over):

$$dx^{i}+L^{ij}(\partial_{j}U^{0}-\varphi_{j}^{0})dt-\partial_{j}L^{ij}dt=b^{i\alpha}dW_{\alpha}(t),$$

(6)

where all variables are even under time reversal, and the $2\!\times\!2$ matrices $L^{ij}$ and $b^{i\alpha}$ are given respectively by

	$\displaystyle{\bm{b}}$	$\displaystyle=$	$\displaystyle\sqrt{\frac{2T}{\gamma}}\begin{pmatrix}1&0\\ 0&1\end{pmatrix},$		(7)
	$\displaystyle{\bm{L}}$	$\displaystyle=$	$\displaystyle{\bm{L}}^{T}=\frac{T}{\gamma}\begin{pmatrix}1&0\\ 0&1\end{pmatrix}=\frac{1}{2}{\bm{b}}{\bm{b}}^{T}.$		(8)

Since both $T$ and $\gamma$ are constants, $\partial_{j}L^{ij}$ in Eq. (6) vanishes identically. We shall call $U^{0}$ and $\varphi_{i}^{0}$ respectively the generalized potential and the non-conservative force ¹¹1Strictly speaking, the non-conservative force is defined as $T\bm{\varphi}$ in Ref. [21]..

The Langevin equation (4) is equivalent to the following covariant Fokker-Planck equation (FPE):

$\displaystyle\partial_{t}p-\partial_{i}\frac{T}{\gamma}(\partial_{i}+\partial_{i}U^{0}-\varphi_{i}^{0})p=0,$

(9)

which can also be written in the form of:

$\displaystyle\partial_{t}p+\partial_{k}j_{k}=0,$

(10)

where $j_{k}$ is the probability current:

$\displaystyle j_{i}=-\frac{T}{\gamma}(\partial_{i}+(\partial_{i}U^{0})-\varphi_{i}^{0})p.$

(11)

It is easy to see that the following transformation:

	$\displaystyle U^{0}$	$\displaystyle\rightarrow$	$\displaystyle U=U^{0}+\psi,$		(12a)
	$\displaystyle\varphi_{i}^{0}$	$\displaystyle\rightarrow$	$\displaystyle\varphi_{i}=\varphi_{i}^{0}+\partial_{i}\psi,$		(12b)

leaves the combination $\partial_{i}U^{0}-\varphi_{i}^{0}$ invariant, and hence also leaves the Langevin equation (6) and the Fokker-Planck equation (9) as well as the probability current (11) invariant. Inspecting Eqs. (5a) and (5b), we see that the transform (12) may be understood as a simultaneous change of the external potential and the flow field that preserves the Brownian motion. We shall call it a gauge transformation. A particular decomposition of the combination $\partial_{i}U-\varphi_{i}$ into $\partial_{i}U$ and $\varphi_{i}$ shall then be called a gauge.

The most convenient gauge is the Gibbs gauge [21], where $U$ is related to the NESS via

$\displaystyle p^{\rm ss}(\bm{x})=e^{-U(\bm{x})}.$

(13)

Substituting this back into Eq. (11), we find the NESS probability current is then given by

$\displaystyle j_{i}^{\rm ss}(\bm{x})=\frac{T}{\gamma}e^{-U(\bm{x})}\varphi_{i}(\bm{x}).$

(14)

which is proportional to $\varphi_{i}$. The fact that $\varphi_{i}$ is non-vanishing characterizes the non-equilibrium nature of the NESS. Substituting Eq. (14) into the steady state FPE $\nabla\cdot{\bm{j}}^{\rm ss}=0$, we obtain the Gibbs gauge condition:

$$\displaystyle\partial_{i}\varphi_{i}-(\partial_{i}U)\varphi_{i}=0,$$

(15)

which, using Eq. (12), can be further rewritten into:

$\displaystyle\partial_{i}(\varphi_{i}^{0}+\partial_{i}\psi)-(\varphi_{i}^{0}+\partial_{i}\psi)\partial_{i}(U^{0}+\psi)=0.$

(16)

In Sec. II.3, we solve this nonlinear differential equation for the case of simple shear flow and harmonic confining potential, and use Eqs. (12) to determine $U,\varphi_{i}$ for the Gibbs gauge.

In the Gibbs gauge, the Langevin equation and the FPE are given by

	$\displaystyle dx^{i}+\frac{T}{\gamma}(\partial_{i}U-\varphi_{i})\,dt$	$\displaystyle=$	$\displaystyle\sqrt{\frac{2T}{\gamma}}dW_{i},\quad\quad$		(17a)
	$\displaystyle\partial_{t}p-\frac{T}{\gamma}\partial_{i}(\partial_{i}+\partial_{i}U-\varphi_{i})\,p$	$\displaystyle=$	$\displaystyle 0.$		(17b)

Equation (17a) will be called the Gibbs representation of the Langevin dynamics, whereas Eqs. (2) and (4) will be called the natural representation of the Langevin dynamics.

We now define the adjoint Langevin dynamics, which is related to the original dynamics (17) via the following transform in the Gibbs gauge:

$\displaystyle U^{\rm Ad}=U,\quad\varphi_{i}^{\rm Ad}=-\varphi_{i}.$

(18)

Using Eqs. (13) and (14), we see that the adjoint process has the same NESS pdf and opposite NESS probability current as the original process:

	$\displaystyle p^{\rm Ad,ss}(\bm{x})$	$\displaystyle=$	$\displaystyle e^{-U(\bm{x})}=p^{\rm ss}(\bm{x}).$		(19)
	$\displaystyle j_{i}^{\rm Ad,ss}(\bm{x})$	$\displaystyle=$	$\displaystyle-\frac{T}{\gamma}e^{-U(\bm{x})}\varphi_{i}(\bm{x})=-j_{i}^{\rm ss}(\bm{x}).$		(20)

Just as the original dynamics, the adjoint dynamics can also be realized by many different combinations of flow field and confining potential, each characterized by a pair $\{U^{0,\rm Ad},\varphi_{i}^{0,\rm Ad}\}$ that is related to $\{U^{\rm Ad},\bm{\varphi}^{\rm Ad}\}$ via a gauge transformation:

	$\displaystyle U^{\rm 0,Ad}$	$\displaystyle\rightarrow$	$\displaystyle U^{\rm Ad}=U^{\rm 0,Ad}+\psi^{\rm Ad},$		(21a)
	$\displaystyle\varphi_{i}^{\rm 0,Ad}$	$\displaystyle\rightarrow$	$\displaystyle\varphi_{i}^{\rm Ad}=\varphi_{i}^{\rm 0,Ad}+\partial_{i}\psi^{\rm Ad},$		(21b)

which is the counterpart of Eqs. (12). The gauge function $\psi^{\rm Ad}$ is arbitrary. The most convenient choice is:

$\displaystyle\psi^{\rm Ad}$

$\displaystyle=$

$\displaystyle-\psi.$

(22)

Substituting this back into Eqs. (21), we may express $\{U^{0,\rm Ad},\varphi_{i}^{0,\rm Ad}\}$ in terms of $\{U^{\rm Ad},\bm{\varphi}^{\rm Ad},\psi\}$. Combining these results with Eqs. (18), we obtain:

	$\displaystyle U^{0,\rm Ad}$	$\displaystyle=$	$\displaystyle U^{0}+2\psi;$		(23a)
	$\displaystyle\varphi_{i}^{0,\rm Ad}$	$\displaystyle=$	$\displaystyle-\varphi_{i}^{0}.$		(23b)

Substituting these into Eqs. (5), we find the confining potential and the flow field for the adjoint process:

	$\displaystyle V^{\rm Ad}({\bm{x}})$	$\displaystyle=$	$\displaystyle T\,(U^{0}+2\,\psi-C^{0})$		(24a)
		$\displaystyle=$	$\displaystyle V({\bm{x}})+2\,T\,\psi({\bm{x}}),$
	$\displaystyle\bm{v}^{\rm Ad}({\bm{x}})$	$\displaystyle=$	$\displaystyle-\frac{T}{\gamma}\bm{\varphi}=-\bm{v}({\bm{x}}).$		(24b)

Hence the flow field of the adjoint dynamics is the opposite of that of the original process.

The Gibbs representation of the adjoint Langevin dynamics can be obtained by using Eq. (18) in Eq. (17a), whereas the natural representation of the adjoint Langevin dynamics can be obtained by using Eqs. (23) in Eq. (4), or, equivalently, by using Eqs. (24) in Eqs. (2).

For numerical simulations (to be detailed in Sec. V), we shall only consider a harmonic confining potential and a simple shear flow:

	$\displaystyle V(\bm{x})$	$\displaystyle=$	$\displaystyle\frac{K}{2}(\bm{x}-\bm{x}_{0})^{2},$		(25a)
	$\displaystyle{\bm{v}}(\bm{x})$	$\displaystyle=$	$\displaystyle{y\,\zeta}\,\hat{\bm{e}}_{x}.$		(25b)

Using Eqs. (5) we find

	$\displaystyle U^{0}$	$\displaystyle=$	$\displaystyle\frac{\beta K}{2}(\bm{x}-\bm{x}_{0})^{2}+C^{0},$		(26a)
	$\displaystyle{\bm{\varphi}}^{0}(\bm{x})$	$\displaystyle=$	$\displaystyle{\beta\gamma}\zeta\,y\,\hat{\bm{e}}_{x},$		(26b)

Note that $\zeta$ has the dimension of inverse time, and $\gamma\,\zeta$ is even under time reversal. The dimension of $\bm{\varphi}$ is inverse of length, and hence also even under time reversal. The shear flow and the confining potential are illustrated in Fig. 1, together with a contour line of the NESS pdf.

We define a dimensionless parameter $\epsilon\equiv{\gamma\,\zeta}/{K}$, which characterizes the relative importance of the flow field compared with the confining potential. For colloidal particles in shearing fluid under typical experimental conditions, this parameter is expected to be much less than the unity, hence we expect that the flow field only leads to a small perturbation of the equilibrium distribution of the Brownian particle. Therefore we may solve Eq. (16) by expanding $\psi$ in terms of $\epsilon$, and subsequently use Eqs. (12) to find $U$ and $\bm{\varphi}$. Since the calculation is rather straightforward, we skip all details and directly present the second order results:

	$\displaystyle\psi$	$\displaystyle=$	$\displaystyle\frac{\beta K\epsilon}{2}(-xy-xy_{0}+x_{0}y)$		(27a)
		$\displaystyle+$	$\displaystyle\frac{\beta K\epsilon^{2}}{8}(-x^{2}+y^{2}+2x_{0}x+2y_{0}y),$
	$\displaystyle U$	$\displaystyle=$	$\displaystyle C+U^{0}+\psi,$		(27b)
	$\displaystyle\bm{\varphi}$	$\displaystyle=$	$\displaystyle\frac{\beta K\epsilon}{2}\left[(y-y_{0})\hat{\bm{e}}_{x}+(-x+x_{0})\,\hat{\bm{e}}_{y}\right]$		(27c)
		$\displaystyle+$	$\displaystyle\frac{\beta K\epsilon^{2}}{4}\left[(-x+x_{0})\hat{\bm{e}}_{x}+(y+y_{0})\,\hat{\bm{e}}_{y}\right],\quad$

where the constant $C$ is such that Eq. (13) is normalized. We shall not need the concrete expression for $C$. The reader may verify directly that Eqs. (27a) does satisfy the Gibbs gauge condition (15) up to $o(\epsilon^{2})$.

If $\epsilon$ is not small, Eqs. (27) are not good approximations. Nonetheless, it is easy to see from Eq. (16) that, due to the quadratic nature of $U^{0}$ and the linear nature of $\varphi^{0}_{i}$, $\psi$ is quadratic in ${\bm{x}}$. Consequently, $U$ is also quadratic in ${\bm{x}}$, and $\bm{\varphi}$ is linear in ${\bm{x}}$. Contour lines of $U$ are therefore all ellipses, one of them being illustrated in Fig. 1. We can therefore set

$\displaystyle U=Ax^{2}+Bxy+Cy^{2}+Dx+Ey+F,$

(28)

and numerically find all coefficients $A,B,C,D,E$. $F$ is determined by the condition of normalization. The numerical method is explained in App. A.1.

We define the fluctuating internal energy as

$\displaystyle H(\bm{x};\lambda)\equiv TU({\bm{x}};\lambda)=T(U^{0}+\psi),$

(29)

where $U({\bm{x}};\lambda)$ is the generalized potential in Gibbs gauge, to be found by solving Eq. (15). The NESS distribution Eq. (13) can then be rewritten as

$\displaystyle p^{\rm ss}(\bm{x};\lambda)=e^{-\beta H(\bm{x};\lambda)}=e^{-U(\bm{x};\lambda)}.$

(30)

To realize such a NESS, one only needs to hold the shear flow and the control parameter $\lambda$ fixed for a sufficiently long period of time. The equilibrium free energy $F(\lambda)$ is defined as

$\displaystyle F(\lambda)\equiv-T\,\log\int_{\bm{x}}e^{-\beta H}=-T\,\log 1=0,$

(31)

where in the second step we used the fact that $p^{\rm ss}(\bm{x})$ is properly normalized. Hence our definition (29) of energy guarantees that the equilibrium free energy vanishes identically. This leads to certain simplification of the fluctuation theorems, as we shall see below.

We define differential heat at the trajectory level as

	$\displaystyle d\bar{}\hskip 1.00006pt{\mathscr{Q}}$	$\displaystyle\equiv$	$\displaystyle d_{\bm{x}}H-T\,\bm{\varphi}\circ d{\bm{x}}$		(32)
		$\displaystyle=$	$\displaystyle T\,(d_{\bm{x}}U-\bm{\varphi}\circ d{\bm{x}}),$

where $d_{\bm{x}}H$ means differential of $H$ due to variation of $\bm{x}$, and $\circ$ is the stochastic product in Stratonovich’s sense:

$\displaystyle\bm{\varphi}\circ d{\bm{x}}\equiv\varphi_{i}({\bm{x}}+d{\bm{x}}/2)\,dx^{i}.$

(33)

The differential work at the trajectory level is defined as

	$\displaystyle d\bar{}\hskip 1.00006pt{\mathscr{W}}$	$\displaystyle\equiv$	$\displaystyle d_{\lambda}H+T\,\bm{\varphi}\circ d{\bm{x}}.$		(34)
		$\displaystyle=$	$\displaystyle T\,(d_{\lambda}U+\bm{\varphi}\circ d{\bm{x}}),$

where $d_{\lambda}H$ is the differential of $H$ due to variation of $\lambda$.

The above defined heat and work can be decomposed into a housekeeping part and an excess part:

	$\displaystyle d\bar{}\hskip 1.00006pt{\mathscr{Q}}$	$\displaystyle=$	$\displaystyle d\bar{}\hskip 1.00006pt{\mathscr{Q}}^{\rm hk}+d\bar{}\hskip 1.00006pt{\mathscr{Q}}^{\rm ex},$		(35)
	$\displaystyle d\bar{}\hskip 1.00006pt{\mathscr{W}}$	$\displaystyle=$	$\displaystyle d\bar{}\hskip 1.00006pt{\mathscr{W}}^{\rm hk}+d\bar{}\hskip 1.00006pt{\mathscr{W}}^{\rm ex},$		(36)

where $d\bar{}\hskip 1.00006pt{\mathscr{Q}}^{\rm hk},d\bar{}\hskip 1.00006pt{\mathscr{W}}^{\rm hk}$ are respectively housekeeping heat and housekeeping work, whereas $d\bar{}\hskip 1.00006pt{\mathscr{Q}}^{\rm ex},d\bar{}\hskip 1.00006pt{\mathscr{W}}^{\rm ex}$ are respectively excess heat and excess work, defined as

	$\displaystyle d\bar{}\hskip 1.00006pt{\mathscr{Q}}^{\rm hk}$	$\displaystyle\equiv$	$\displaystyle-T\,\bm{\varphi}\circ d{\bm{x}},$		(37a)
	$\displaystyle d\bar{}\hskip 1.00006pt{\mathscr{Q}}^{\rm ex}$	$\displaystyle\equiv$	$\displaystyle T\,d_{\bm{x}}U,$		(37b)
	$\displaystyle d\bar{}\hskip 1.00006pt{\mathscr{W}}^{\rm hk}$	$\displaystyle\equiv$	$\displaystyle T\,\bm{\varphi}\circ d{\bm{x}}=-d\bar{}\hskip 1.00006pt{\mathscr{Q}}^{\rm ex},$		(37c)
	$\displaystyle d\bar{}\hskip 1.00006pt{\mathscr{W}}^{\rm ex}$	$\displaystyle\equiv$	$\displaystyle T\,d_{\lambda}U.$		(37d)

The first law of thermodynamics at the trajectory level is given by either of the following two forms:

$\displaystyle\begin{split}dH&=d\bar{}\hskip 1.00006pt{\mathscr{Q}}+d\bar{}\hskip 1.00006pt{\mathscr{W}}=d\bar{}\hskip 1.00006pt{\mathscr{Q}}^{\rm ex}+d\bar{}\hskip 1.00006pt{\mathscr{W}}^{\rm ex}.\end{split}$

(38)

Note that the housekeeping heat and housekeeping work exactly cancel each other.

Heat and work at the ensemble level can be obtained by averaging the corresponding quantities at the trajectory level over both noises and the pdf of ${\bm{x}}$. To obtain a well-defined continuum limit, these differential quantities at the ensemble level must be computed up to the first order in $dt$. It is important to remember that the Wiener noises are square root of $dt$, and hence according to Eq. (17a), $d{\bm{x}}$ contains parts scaling with $\sqrt{dt}$. Consequently, we need to expand these differential quantities up to the second order of $d{\bm{x}}$, to keep all terms linear in $dt$.

As an example, let us compute the excess heat at the ensemble level:

$\displaystyle d\bar{}\hskip 1.00006ptQ^{\rm ex}$

$\displaystyle=$

$\displaystyle\langle\!\langle d\bar{}\hskip 1.00006pt{\mathscr{Q}}^{\rm ex}\rangle\!\rangle,$

(39)

where $\langle\!\langle\,\cdot\,\rangle\!\rangle$ means double average over Wiener noises and over probability distribution of $\bm{x}$. First we use Eq. (37b) to expand $d\bar{}\hskip 1.00006pt{\mathscr{Q}}^{\rm ex}$ up to the second order in $d{\bm{x}}$:

$\displaystyle d\bar{}\hskip 1.00006pt{\mathscr{Q}}^{\rm ex}=T\partial_{i}Udx^{i}+\frac{T}{2}\partial_{i}\partial_{j}Udx^{i}dx^{j},$

(40)

where all products are in Ito’s sense. We now use the Langevin equation (17a) to express $d{\bm{x}}$ in terms of $dt$ and $d\bm{W}$. All terms in the form of $dt^{2}$ and $dtdW_{i}$ can be neglected, since they are higher order than $dt$. Then we average over the Wiener noises $d\bm{W}$. Finally we multiply the result by the pdf $p({\bm{x}};t)$ and integrate over ${\bm{x}}$, and obtain the differential excess heat at the ensemble level:

$\displaystyle d\bar{}\hskip 1.00006ptQ^{\rm ex}=-\frac{T^{2}\,dt}{\gamma}\left\langle\partial_{i}U(\partial_{i}U-\varphi_{i})-\partial_{i}\partial_{i}U\right\rangle,$

(41)

where $\langle\,\cdot\,\rangle$ means average over the pdf $p({\bm{x}},t)$ of ${\bm{x}}$:

$\displaystyle\langle\,\cdot\,\rangle=\int_{\bm{x}}\cdot\,p({\bm{x}},t).$

(42)

Similarly, the differential housekeeping heat at the ensemble level is given by

$\displaystyle d\bar{}\hskip 1.00006ptQ^{\rm hk}=-\frac{T^{2}\,dt}{\gamma}\left\langle\varphi_{i}(\partial_{i}U-\varphi_{i})-\partial_{i}\varphi_{i}\right\rangle.$

(43)

Further using the Gibbs gauge condition (15) we may rewrite the above result as

$\displaystyle d\bar{}\hskip 1.00006ptQ^{\rm hk}=-\frac{T^{2}\,dt}{\gamma}\int_{\bm{x}}p(\varphi_{i})^{2}\leq 0,$

(44)

which is non-positive definite.

The differential housekeeping and excess work at the trajectory level can be similarly computed:

	$\displaystyle d\bar{}\hskip 1.00006ptW^{\rm hk}$	$\displaystyle=$	$\displaystyle-d\bar{}\hskip 1.00006ptQ^{\rm hk}=\frac{T^{2}\,dt}{\gamma}\int_{\bm{x}}p(\varphi_{i})^{2}\geq 0,$		(45)
	$\displaystyle d\bar{}\hskip 1.00006ptW^{\rm hk}$	$\displaystyle=$	$\displaystyle T\int_{\bm{x}}p\,d_{\lambda}U$		(46)

The total heat and work at the ensemble level are then the sum of the corresponding housekeeping parts and excess parts:

	$\displaystyle d\bar{}\hskip 1.00006pt{Q}$	$\displaystyle=$	$\displaystyle\langle\!\langle d\bar{}\hskip 1.00006pt{\mathscr{Q}}\rangle\!\rangle=d\bar{}\hskip 1.00006pt{Q}^{\rm hk}+d\bar{}\hskip 1.00006pt{Q}^{\rm ex},$		(47)
	$\displaystyle d\bar{}\hskip 1.00006pt{W}$	$\displaystyle=$	$\displaystyle\langle\!\langle d\bar{}\hskip 1.00006pt{\mathscr{W}}\rangle\!\rangle=d\bar{}\hskip 1.00006pt{W}^{\rm hk}+d\bar{}\hskip 1.00006pt{W}^{\rm ex}.$		(48)

The system entropy is:

$\displaystyle S[p]=-\int_{\bm{x}}p\,\log p,$

(49)

whose differential can be calculated using the Fokker-Planck equation:

	$\displaystyle dS$	$\displaystyle=$	$\displaystyle-dt\int_{\bm{x}}\log p\,\mathscr{L}p$		(50)
		$\displaystyle=$	$\displaystyle-dt\int_{\bm{x}}\log p\,\partial_{i}\frac{T}{\gamma}(\partial_{i}+\partial_{i}U-\varphi_{i})p$
		$\displaystyle=$	$\displaystyle\frac{dt\,T}{\gamma}\int_{\bm{x}}\frac{1}{p}(\partial_{i}p)(\partial_{i}+\partial_{i}U-\varphi_{i})p,$

where in the last step we have integrated by parts.

The EP is defined as

$\displaystyle dS^{\rm tot}$

$\displaystyle\equiv$

$\displaystyle dS+dS^{\rm env}=dS-\beta\,d\bar{}\hskip 1.00006ptQ,$

(51)

where $dS^{\rm env}\equiv-\beta\,d\bar{}\hskip 1.00006ptQ$ is defined as the environmental entropy change. As explained previously, $dS^{\rm env}$ is only the part of environmental entropy change that can be captured by our theory of stochastic thermodynamics. This can be further decomposed into a housekeeping EP and an excess EP :

	$\displaystyle dS^{\rm tot}$	$\displaystyle=$	$\displaystyle dS^{\rm hk}+dS^{\rm ex},$		(52)
	$\displaystyle dS^{\rm hk}$	$\displaystyle=$	$\displaystyle-\beta d\bar{}\hskip 1.00006ptQ^{\rm hk}=\frac{T\,dt}{\gamma}\int_{\bm{x}}p(\varphi_{i})^{2}\geq 0,$		(53)
	$\displaystyle dS^{\rm ex}$	$\displaystyle=$	$\displaystyle dS-\beta\,d\bar{}\hskip 1.00006ptQ^{\rm ex}.$		(54)

In particular, in the NESS, the excess EP vanishes identically, whereas the housekeeping EP reduces to

$\displaystyle\frac{dS^{\rm hk}}{dt}=\frac{T\,dt}{\gamma}\int_{\bm{x}}e^{-U}(\varphi_{i})^{2}=\frac{\gamma}{T}\int_{\bm{x}}e^{U}\left(j_{i}^{\rm ss}\right)^{2},$

(55)

where $j_{i}^{\rm ss}$ is the NESS current given in Eq. (14).

Further using Eqs. (50) and (41), as well as the Gibbs gauge condition (15), we may rewrite the excess EP in the following apparently positive form:

$\displaystyle dS^{\rm ex}=\frac{T\,dt}{\gamma}\int_{\bm{x}}\frac{1}{p}\left[(\partial_{i}+\partial_{i}U)p\right]^{2}\geq 0.$

(56)

Hence EP is the sum of a positive housekeeping part and a positive excess part, a general feature of Markov systems with even variables and parameters that lack instantaneous detailed balance [17, 18, 19, 20].

Finally we may also define non-equilibrium free energy:

$\displaystyle F[p]\equiv\int_{\bm{x}}p(H+T\,\log p).$

(57)

It is then easy to verify the following differential forms:

$\displaystyle dF[p]=d\bar{}\hskip 1.00006ptW^{\rm ex}+d\bar{}\hskip 1.00006ptQ^{\rm ex}-T\,dS,$

(58)

which may be further rewritten as

$\displaystyle dS^{\rm ex}=dS-\beta\,d\bar{}\hskip 1.00006ptQ^{\rm ex}=\beta(d\bar{}\hskip 1.00006ptW^{\rm ex}-dF[p])\geq 0.$

(59)

To study fluctuation theorems, it is necessary to know the short-time transition probability of the Langevin process defined by Eq. (17a). Let ${\bm{x}},{\bm{x}}_{1}={\bm{x}}+d{\bm{x}}$ be respectively the initial position and the final position of an infinitesimal transition taking place during $dt$, and let ${\bm{x}}_{1/2}={\bm{x}}+d{\bm{x}}/2$ be the mid-point. A general expression for the short-time transition probability $p_{\bm{\varphi}}({\bm{x}}_{1}|{\bm{x}};dt)$ of the Langevin equation (6) was derived in Eqs. (A4) of Ref. [21], using the general result of time-slicing path integral in Ref. [23]. Specializing to the Langevin dynamics Eq. (17a), we find ²²2The notation may be unfortunately confusing since we need to integrate over $d{\bm{x}}$ to verify the normalization of this transition probability density function. To avoid this confusion, we may replace $d{\bm{x}},dt$ by $\Delta{\bm{x}},\Delta t$ and remember that they are infinitesimal quantities. (Note that the notations are slightly different here)

	$\displaystyle p_{\bm{\varphi}}({\bm{x}}_{1}|{\bm{x}};dt)=\frac{\gamma}{4\pi T\,dt}e^{-{\mathcal{A}}_{\bm{\varphi}}(d{\bm{x}};{\bm{x}}_{1/2},dt)},$				(60a)
where the action ${\mathcal{A}}_{\bm{\varphi}}(d{\bm{x}};{\bm{x}}_{1/2},dt)$ is given by
	$\displaystyle{\mathcal{A}}_{\bm{\varphi}}(d{\bm{x}};{\bm{x}}_{1/2},dt)$	$\displaystyle=$	$\displaystyle\frac{\gamma}{4T\,dt}(dx^{i}+\frac{T}{\gamma}(\partial_{i}U-\varphi_{i})dt)^{2}_{1/2}$
		$\displaystyle-$	$\displaystyle\frac{T\,dt}{2\gamma}(\partial_{i}^{2}U-\partial_{i}\varphi_{i})_{1/2}+o(dt),$
					(60b)

where the subscript $1/2$ in Eq. (60b) means that all functions inside the bracket are evaluated at ${\bm{x}}_{1/2}$. The action is expanded up to the first order in $dt$, which is sufficient to guarantees a correct continuum limit. In fact it is ok to evaluate the second term in the r.h.s. of Eq. (60b) at any point, the resulting error is of higher order than $dt$, and hence is negligible in the continuum limit. Note that we show explicitly the dependence of the action on the non-conservative force $\bm{\varphi}$.

Let us supply a heuristic explanation for Eqs. (60). First we note that the Wiener noises are infinitesimal Gaussian with basic properties (3). Using these we can readily construct their pdf:

$\displaystyle p(d\bm{W})=\frac{1}{{2\pi dt}}\exp\left(-\frac{dW_{x}^{2}+dW_{y}^{2}}{2dt}\right).$

(61)

Now given ${\bm{x}}$, the Langevin equation (17a) may be understood as a linear relation between $d\bm{W}$ and the infinitesimal displacement $d{\bm{x}}$. Hence we may obtain the pdf for $d{\bm{x}}$ directly from Eq. (61):

$\displaystyle p(d{\bm{x}})=\frac{\gamma}{4\pi T\,dt}\,e^{-\frac{\gamma}{4\pi T\,dt}(dx^{i}+\frac{T}{\gamma}(\partial_{i}U-\varphi_{i})dt)^{2}}.$

(62)

Note that the action appearing in the exponent above is formally identical to the first term in the action (60b). It is important to note however as a basic property of Ito-Langevin equation, the function $(\partial_{i}U-\varphi_{i})$ in (17a), which also appears in Eq. (62) is evaluated at the initial point ${\bm{x}}$. This should be contrasted with Eqs. (60), where the same function is evaluated at the mid-point ${\bm{x}}_{1/2}$. Because of the $dt$ appearing in the denominator of the actions, however, this difference is qualitatively important and is compensated by the second term in the action (60b).

Using Eq. (60a), we can also compute the backward transition probability from ${\bm{x}}_{1}$ to ${\bm{x}}$. All we need is to swap ${\bm{x}}_{1}$ and ${\bm{x}}$ in Eq. (60a). Note that ${\bm{x}}$ and ${\bm{x}}_{1}$ appear in Eq. (60a) only in the combinations $d{\bm{x}}$ and ${\bm{x}}_{1/2}$, which are respectively odd and even under the swap. Hence to obtain $p({\bm{x}}|{\bm{x}}_{1};dt)$ we only need to flip the sign of $d{\bm{x}}$. This leads to

	$\displaystyle p_{\bm{\varphi}}({\bm{x}}|{\bm{x}}_{1};dt)$	$\displaystyle=$	$\displaystyle\frac{\gamma}{4\pi T\,dt}e^{-{\mathcal{A}}_{\bm{\varphi}}(-d{\bm{x}};{\bm{x}}_{1/2},dt)},$		(63a)
	$\displaystyle{\mathcal{A}}_{\bm{\varphi}}(-d{\bm{x}};{\bm{x}}_{1/2},dt)$	$\displaystyle=$	$\displaystyle\frac{\gamma}{4T\,dt}(-dx^{i}+\frac{T}{\gamma}(\partial_{i}U-\varphi_{i})dt)^{2}_{1/2}$
		$\displaystyle-$	$\displaystyle\frac{T\,dt}{2\gamma}(\partial_{i}^{2}U-\partial_{i}\varphi_{i})_{1/2}+o(dt).$
					(63b)

Recall the adjoint process defined in Sec. II.2 is related to the original process by changing the sign of $\bm{\varphi}$. We can construct the corresponding transition probability for the adjoint process from Eqs. (60):

	$\displaystyle p_{-\bm{\varphi}}({\bm{x}}_{1}|{\bm{x}};dt)$	$\displaystyle=$	$\displaystyle\frac{\gamma}{4\pi T\,dt}e^{-{\mathcal{A}}_{-\bm{\varphi}}(d{\bm{x}};{\bm{x}}_{1/2},dt)},$		(64a)
	$\displaystyle{\mathcal{A}}_{-\bm{\varphi}}(d{\bm{x}};{\bm{x}}_{1/2},dt)$	$\displaystyle=$	$\displaystyle\frac{\gamma}{4T\,dt}(dx^{i}+\frac{T}{\gamma}(\partial_{i}U+\varphi_{i})dt)^{2}_{1/2}$
		$\displaystyle-$	$\displaystyle\frac{T\,dt}{2\gamma}(\partial_{i}^{2}U+\partial_{i}\varphi_{i})_{1/2}+o(dt),$
					(64b)

The backward transition probability of the adjoint process can be similarly obtained from Eqs. (63):

	$\displaystyle p_{-\bm{\varphi}}({\bm{x}}|{\bm{x}}_{1};dt)$	$\displaystyle=$	$\displaystyle\frac{\gamma}{4\pi T\,dt}e^{-{\mathcal{A}}_{-\bm{\varphi}}(-d{\bm{x}};{\bm{x}}_{1/2},dt)},$		(65a)
	$\displaystyle{\mathcal{A}}_{-\bm{\varphi}}(-d{\bm{x}};{\bm{x}}_{1/2},dt)$	$\displaystyle=$	$\displaystyle\frac{\gamma}{4T\,dt}(-dx^{i}+\frac{T}{\gamma}(\partial_{i}U+\varphi_{i})dt)^{2}_{1/2}$
		$\displaystyle-$	$\displaystyle\frac{T\,dt}{2\gamma}(\partial_{i}^{2}U+\partial_{i}\varphi_{i})_{1/2}+o(dt).$
					(65b)

Using Eqs. (60) and Eqs. (63)-(65), we readily obtain the following ratios:

	$\displaystyle\frac{p_{\bm{\varphi}}({\bm{x}}_{1}|{\bm{x}};dt)}{p_{\bm{\varphi}}({\bm{x}}|{\bm{x}}_{1};dt)}$	$\displaystyle=$	$\displaystyle e^{-(\partial_{i}U-\varphi_{i})\circ dx^{i}}.$		(66a)
Similarly, with the aid of Gibbs gauge condition Eq. (15), we may also prove
	$\displaystyle\frac{p_{\bm{\varphi}}({\bm{x}}_{1}|{\bm{x}};dt)}{p_{-\bm{\varphi}}({\bm{x}}_{1}|{\bm{x}};dt)}$	$\displaystyle=$	$\displaystyle\frac{p_{-\bm{\varphi}}({\bm{x}}|{\bm{x}}_{1};dt)}{p_{\bm{\varphi}}({\bm{x}}|{\bm{x}}_{1};dt)}=e^{\bm{\varphi}\circ d{\bm{x}}},$		(66b)
	$\displaystyle\frac{p_{\bm{\varphi}}({\bm{x}}_{1}|{\bm{x}};dt)}{p_{-\bm{\varphi}}({\bm{x}}|{\bm{x}}_{1};dt)}$	$\displaystyle=$	$\displaystyle\frac{p_{-\bm{\varphi}}({\bm{x}}_{1}|{\bm{x}};dt)}{p_{\bm{\varphi}}({\bm{x}}|{\bm{x}}_{1};dt)}=e^{-d_{\bm{x}}U}.\quad\quad$		(66c)

Equations (66) may be called the conditions of local detailed balance.

We keep the flow field fixed, and vary parameters $\lambda_{t}=\{K(t),{\bm{x}}_{0}(t)\}$ systematically, which fully determines the Langevin dynamics. We call $\{U({\bm{x}};\lambda_{t}),\bm{\varphi}({\bm{x}};\lambda_{t})\}$ the dynamic protocol in the Gibbs gauge, and $\{V({\bm{x}};\lambda_{t}),\bm{v}({\bm{x}})\}$ the dynamic protocol in the natural gauge. A dynamic process is determined by the initial pdf of $p({\bm{x}},t=0)$ together with the dynamic protocol either in the Gibbs gauge or in the natural gauge. Transformation between two dynamic protocols are given by Eqs. (12) and (5).

We define four processes as below, all of which start from $t=0$ and end at $t=\tau$:

1. 1.

   Forward process:   The dynamic protocol is

   	$\displaystyle U^{\rm F}$	$\displaystyle=$	$\displaystyle U({\bm{x}},\lambda_{t}),\quad$		(67a)
   	$\displaystyle\bm{\varphi}^{\rm F}$	$\displaystyle=$	$\displaystyle\bm{\varphi}({\bm{x}},\lambda_{t}).$		(67b)

   The initial pdf is $p^{\rm ss}({\bm{x}};\lambda_{0})$, defined in Eq. (30).

2. 2.

   Backward process:   The dynamic protocol is

   	$\displaystyle U^{\rm B}$	$\displaystyle=$	$\displaystyle U({\bm{x}},\lambda_{\tau-t}),\quad$		(68a)
   	$\displaystyle\bm{\varphi}^{\rm B}$	$\displaystyle=$	$\displaystyle\bm{\varphi}({\bm{x}},\lambda_{\tau-t}).$		(68b)

   The initial pdf is $p^{\rm ss}({\bm{x}};\lambda_{\tau})$.

3. 3.

   Adjoint process:   The dynamic protocol is

   	$\displaystyle U^{\rm Ad}$	$\displaystyle=$	$\displaystyle U({\bm{x}},\lambda_{t}),\quad$		(69a)
   	$\displaystyle\bm{\varphi}^{\rm Ad}$	$\displaystyle=$	$\displaystyle-\bm{\varphi}({\bm{x}},\lambda_{t}).$		(69b)

   The initial pdf is $p^{\rm ss}({\bm{x}};\lambda_{0})$.

4. 4.

   Adjoint backward process:   The protocol is

   	$\displaystyle U^{\rm AdB}$	$\displaystyle=$	$\displaystyle U({\bm{x}},\lambda_{\tau-t}),\quad$		(70a)
   	$\displaystyle\bm{\varphi}^{\rm AdB}$	$\displaystyle=$	$\displaystyle-\bm{\varphi}({\bm{x}},\lambda_{\tau-t}).$		(70b)

   The initial pdf is $p^{\rm ss}({\bm{x}};\lambda_{\tau})$.

Note that each of these processes starts from the NESS corresponding to the initial control parameter of the dynamic protocol. Such an initial state can be realized easily in experiments. Note also that in general, the system is not in a NESS at the end of any of these processes.