Perturbative Expansion of Irreversible Work in Fokker-Planck Equation a la Quantum Mechanics

TL;DR

This paper introduces a perturbative method inspired by quantum mechanics to systematically expand and calculate irreversible work in Fokker-Planck systems, simplifying computations and enabling optimization of control protocols.

Contribution

It develops a novel perturbative expansion approach for the Fokker-Planck equation, analogous to quantum adiabatic calculations, to efficiently compute irreversible work and optimize control protocols.

Findings

First order expansion yields exact results for harmonic potentials.

Method simplifies calculations by avoiding coupled differential equations.

Optimized protocols minimize irreversible work effectively.

Abstract

We discuss the systematic expansion of the solution of the Fokker-Planck equation with the help of the eigenfunctions of the time-dependent Fokker-Planck operator. The expansion parameter is the time derivative of the external parameter which controls the form of an external potential. Our expansion corresponds to the perturbative calculation of the adiabatic motion in quantum mechanics. With this method, we derive a new formula to calculate the irreversible work order by order, which is expressed as the expectation value with a pseudo density matrix. Applying this method to the case of the harmonic potential, we show that the first order term of the expansion gives the exact result. Because we do not need to solve the coupled differential equations of moments, our method simplifies the calculations of various functions such as the fluctuation of the irreversible work per unit time. We further investigate the exact optimized protocol to minimize the irreversible work by calculating its variation with respect to the control parameter itself.

Equations

$${\partial }_t\rho (x,t)=\left[\frac{1}{\nu \beta }{\partial }_x^2+\frac{1}{\nu }{\partial }_x{V}^{(1)}(x,a_t)\right]\rho (x,t)\equiv {\mathcal{L}}_t(x)\rho (x,t),$$

$${\mathcal{L}}_t(x){\rho }_n(x,a_t)$$

$${\mathcal{L}}_t^{†}(x){\tilde{\rho }}_n(x,a_t)=-{\bar{\lambda }}_n(t){\tilde{\rho }}_n(x,a_t),$$

$${\mathcal{L}}_t^{†}(x)=\left[\frac{1}{\nu \beta }{\partial }_x^2-\frac{1}{\nu }{V}^{(1)}(x,a_t){\partial }_x\right].$$

$$e^{G(x,a_t)/2}{\mathcal{L}}_t(x)e^{-G(x,a_t)/2}=e^{-G(x,a_t)/2}{\mathcal{L}}_t^{†}(x)e^{G(x,a_t)/2}=-{\mathcal{H}}_t,$$

$${\mathcal{H}}_t=-\frac{1}{\nu \beta }{\partial }_x^2-\frac{1}{2\nu }{V}^{(2)}(x,a_t)+\frac{\beta }{4\nu }(V^{(1)}(x,a_t){)}^2,$$

$$G(x,a_t)=\beta V(x,a_t).$$

$${\mathcal{H}}_t{u}_n(x,a_t)={\lambda }_n(t)u_n(x,a_t),$$

$$\int dx{u}_n(x,a_t)u_m(x,a_t)$$

$$\sum _n{u}_n(x,a_t)u_n(x^{\prime },a_t)$$

$${\rho }_n(x,a_t)$$

$${\tilde{\rho }}_n(x,a_t)$$

$${\bar{\lambda }}_n(t)={\lambda }_n(t).$$

$$\int dx{\tilde{\rho }}_n(x,a_t){\rho }_m(x,a_t)$$

$$\sum _n{\tilde{\rho }}_n(x,a_t){\rho }_n(x^{\prime },a_t)$$

$$0={\bar{\lambda }}_0(t)<{\bar{\lambda }}_1(t)<{\bar{\lambda }}_2(t)<\cdots ,$$

$${\rho }_0(x,a_t)$$

$${\tilde{\rho }}_0(x,a_t)$$

$$Z(a_t)=\int dx{e}^{-\beta V(x,a_t)}.$$

$${\mathcal{H}}_t=\frac{1}{\nu \beta }{B}^{†}B,$$

$$B$$

$$B^{†}$$

$${\bar{\lambda }}_n(t)={\lambda }_n(t)=\frac{1}{\nu \beta }\int dx|B{u}_n(x,a_t){|}^2\ge 0.$$

$$[{\mathcal{H}}_t,B]=-\frac{1}{\nu }{V}^{(2)}(x,a_t)B.$$

$${\mathcal{H}}_t(B{u}_0(x,a_t))=\left({\lambda }_0(t)-\frac{1}{\nu }{V}^{(2)}(x,a_t)\right)(B{u}_0(x,a_t)).$$

$$u_0(x,a_t)=\frac{1}{\sqrt{Z(a_t)}}{e}^{-\beta V(x,a_t)/2},$$

$$[{\mathcal{H}}_t,B^{†}]=\frac{1}{\nu }{B}^{†}{V}^{(2)}(x,a_t),$$

$${\rho }_n(x,a_t)$$

$${\tilde{\rho }}_n(x,a_t)$$

$$⟨x|{\hat{\mathcal{L}}}_t|n,a_t⟩$$

$$⟨n,a_t|{\hat{\mathcal{L}}}_t|x⟩$$

$$\int dx|x⟩⟨x|$$

$$⟨x|\hat{G}(a_t)|x^{\prime }⟩$$

$$⟨x|\hat{V}(a_t)|x^{\prime }⟩$$

$$⟨x|{\hat{\mathcal{L}}}_t|x^{\prime }⟩$$

$$⟨x|1|x^{\prime }⟩$$

$${\partial }_t|\rho (t)⟩={\hat{\mathcal{L}}}_t|\rho (t)⟩,$$

Full text

I Introduction

The accelerating development of experimental methods enables us to access individual thermal random processes at microscopic scales. Because the typical scale of the system is very small and then thermal fluctuations play important roles, thermodynamics is not directly applicable to describe this system. The establishment of the coarse-grained description of the small-fluctuating systems is an intriguing subject of statistical physics and nanophysics.

Such a system is often modeled by the Brownian motion hanngi-review ; broeck-review ; sekimoto-book ; sei-rev . Then the work which is exerted or absorbed by the system of the Brownian particle is represented as the change of the form of the external confinement potential. In practice, any protocol of this change is realized within a finite time period and thus we cannot avoid the loss of energy. However, we can still consider the optimized protocol which minimizes the irreversible work for a fixed time period sekimoto-book ; sekimoto-sasa ; schmiedl-seifert2007 ; schmiedl-seifert2008 ; gomez ; koning ; aurell ; geiger ; then ; dechant ; bona ; zulkow ; crook ; all . The investigation of the optimized protocol is important to construct more efficient nanomachines.

As is well-known, the distribution function of the Brownian particle is described by the Fokker-Planck equation and the irreversible work is calculated from it. In principle, the behavior can be solved numerically, but such an approach will not be promising for investigating the optimization problem. On the other hand, when the work is expressed as an analytic function of the control parameter, we can apply the variational scheme to find the optimized protocol sekimoto-sasa ; sekimoto-book . In fact, the optimization has been exclusively studied for the harmonic potential which can be solved exactly schmiedl-seifert2007 ; schmiedl-seifert2008 ; gomez ; speck ; mazon ; raybov ; imparato ; cohen ; kwon ; hol2 . 111The optimization in the logarithmic-harmonic potential is studied in Ref. hol .

In an exact calculation, the macroscopic quantities of the Brownian particle is expressed as the functions of the moments of the position of the Brownian particle schmiedl-seifert2007 ; schmiedl-seifert2008 ; gomez ; speck . Therefore, to calculate the fluctuations of, for example, the irreversible work, we need to solve highly coupled differential equations of the moments. In addition, all methods which have been employed to solve the exactly solvable models are not easily applicable to systems which have more general potentials blickle ; dieterich . To study the optimized protocol in more general situations, we thus need to develop a systematic method to find an analytic expression of the irreversible work.

In this work, we develop a perturbative expansion method to calculate the solution of the Fokker-Planck equation. The expansion basis is the eigenfunctions of the time-dependent Fokker-Planck operator, and the expansion parameter is the time derivative of the external parameter which controls the form of the external potential. Our expansion corresponds to the perturbative calculation of the adiabatic motion in quantum mechanics. Then all results are expressed in integral forms and we do not need to solve differential equations of the moments. This feature is an advantage of our approach compared to the moment method. To confirm the consistency of the expansion, we apply this method to the harmonic potential and derive the irreversible work and its fluctuation. We then show that the first and the second order correction terms of the expansion are sufficient to reproduce the exact results. We further calculate the variation of the irreversible work with respect to the control parameter itself, not to the moment (variance), and derive the exact equation to determine the optimized protocol. Because of the complex integro-differential equation, it is difficult to find the exact solution. Instead we discuss the approximated solution of the integro-differential equation.

This paper is organized as follows. In Sec. II, we derive the bi-orthogonal system with the eigenfunctions of the time-dependent Fokker-Planck operator. Using this, we develop the perturbative expansion of the solution of the Fokker-Planck equation in Sec. III, and derive the expansion formula to calculate the irreversible work in Sec. IV. In Sec. V, we apply the obtained result to the harmonic potential and investigate the exact the exact irreversible work, its fluctuation and the optimized control parameter. Section VI is devoted to the concluding remarks.

II Bi-orthogonal system

In this section, we discuss the definition of the expansion basis. A similar expansion basis is discussed, for example, in Ref. gardiner . In this section, we generalize the method to the case of the time-dependent Fokker-Planck operator namiki . Note that the eigenvalue theory of the time-periodic Fokker-Planck operator (the Kolmogorov operator) is discussed in Ref. caceres and the properties found below are consistent with the result.

We consider a Brownian particle which is confined in an external confinement potential $V$ and interacts with a thermal bath with a fixed temperature $T$. When the inertial term is negligible and the system is spatially one dimensional, the distribution function of the particle $\rho(x,t)$ is described by the Fokker-Planck equation,

[TABLE]

where $\nu$ is the constant friction coefficient and $\beta=1/(k_{B}T)$ with $k_{B}$ being the Boltzmann constant. Note that $V^{(n)}=\partial^{n}_{x}V$ with $V$ being the potential energy, and $a_{t}$ is a time dependent external parameter which controls the form of the potential. Differently from Ref. caceres , $a_{t}$ is an arbitrary function of time. In the following, we consider the case where the particle distribution function $\rho(x,t)$ vanishes quickly at infinite distance, ${\displaystyle\lim_{x\rightarrow\pm\infty}}\rho(x,t)=0$.

The eigenvalue and eigenfunction of the time-dependent Fokker-Planck operator ${\cal L}_{t}(x)$ are defined by

[TABLE]

Note that ${\cal L}_{t}(x)$ is not self-adjoint and the eigenfunctions do not form a complete set in general. However, we can further introduce the eigenfunctions defined by

[TABLE]

where

[TABLE]

Then $\rho_{n}(x,a_{t})$ and $\tilde{\rho}_{n}(x,a_{t})$ form a bi-orthogonal system as is shown soon later and the solution of the Fokker-Planck equation can be expanded by these eigenfunctions.

To construct $\rho_{n}(x,a_{t})$ and $\tilde{\rho}_{n}(x,a_{t})$, we should note that both of ${\cal L}_{t}$ and ${\cal L}^{\dagger}_{t}$ are characterized by ${\cal H}_{t}$, which is defined by

[TABLE]

where

[TABLE]

This quantity ${\cal H}_{t}$ can be regarded as the Hamiltonian operator which has the potential given by $-\frac{1}{2\nu}V^{(2)}(x,a_{t})+\frac{\beta}{4\nu}(V^{(1)}(x,a_{t}))^{2}$. Then we introduce the eigenfunctions of ${\cal H}_{t}$ as

[TABLE]

which form a complete orthonormal set,

[TABLE]

We should further notice that because there is no degeneracy for the eigenvalue of the one-dimensional Hamiltonian ($\lambda_{n}\neq\lambda_{m}$ for $n\neq m$), the eigenfunction $u_{n}$ is given by a real function ($u^{*}_{n}=u_{n}$).

Once the eigenvalues and eigenfunctions of ${\cal H}_{t}$ are found, owing to Eq. (4), the eigenfunctions of the time-dependent Fokker-Planck operators are given by

[TABLE]

and the corresponding eigenvalue is given by

[TABLE]

One can easily confirm that these eigenfunctions form the following bi-orthogonal system,

[TABLE]

We can further show the following properties as the universal natures of the present eigenvalue problem.

The smallest eigenvalue is given by zero and other eigenvalues are larger than zero. Because there is no degeneracy for the eigenvalues, we can set,

[TABLE]

without loss of generality. 2. 2.

The eigenfunctions $\rho_{0}(x,a_{t})$ and $\tilde{\rho}_{0}(x,a_{t})$ for $\bar{\lambda}_{0}(t)$ are given by

[TABLE]

respectively. Here we introduced

[TABLE]

These properties can be shown as follows. Note that the Hamiltonian operator can be reexpressed as

[TABLE]

where

[TABLE]

Therefore the eigenvalues are more than or equal to zero,

[TABLE]

From the absence of the degeneracy of the eigenfunctions, $\bar{\lambda}_{n}(t)\neq\bar{\lambda}_{m}(t)$ for $n\neq m$. Moreover, the Fokker-Planck equation has a stationary state, which can be interpreted as the eigenfunction of zero eigenvalue. Thus, without loss of generality, we can set as $0=\bar{\lambda}_{0}(t)<\bar{\lambda}_{1}(t)<\bar{\lambda}_{2}(t)<\cdots$.

From Eq. (17), it is seen that the normalizable eigenfunction for $\bar{\lambda}_{0}(t)$ is given by the solution of $Bu_{0}(x,a_{t})=0$. However, for the case of $V^{(2)}(x,a_{t})>0$, we can discuss as follows. The operator $B$ satisfies

[TABLE]

Therefore, we can show

[TABLE]

Because $u_{0}(x,a_{t})$ has the lowest eigenvalue of ${\cal H}$, it should be the solution of the equation $Bu_{0}(x,a_{t})=0$, leading to

[TABLE]

Then $\rho_{0}(x,a_{t})$ and $\tilde{\rho}_{0}(x,a_{t})$ are calculated using Eqs. (8) and (9).

There is an important remark. The above operators $B$ and $B^{\dagger}$, generally do not correspond to the lowering and raising operators in quantum mechanics. In fact, we can show that

[TABLE]

and $B^{\dagger}V^{(2)}(x,a_{t})\neq V^{(2)}(x,a_{t})B^{\dagger}$. Therefore we cannot obtain excited states by multiplying $B^{\dagger}$ to $u_{0}(x,a_{t})$. It is however not the case with the harmonic potential as is seen in Sec. V.1.

For the sake of simplicity, we introduce the following bra-ket notation as quantum mechanics,

[TABLE]

Here $|x\rangle$ ($\langle x|$) is an eigenfunction of the position operator $\hat{x}$, $\hat{x}|x\rangle=x|x\rangle$ ($\langle x|\hat{x}=\langle x|x$). Other operators are introduced as quantities satisfying $\langle x|\hat{A}|x^{\prime}\rangle=A(x)\delta(x-x^{\prime})$, for example,

[TABLE]

Then the Fokker-Planck equation is symbolically expressed as

[TABLE]

and the conditions satisfied for the eigenfunctions are summarized as

[TABLE]

In the following, we use this bi-orthogonal set to expand the solution of the Fokker-Planck equation. We should note the difference between our approach and quantum mechanics: the bra-ket vectors represent the probability amplitude in quantum mechanics, while we use it to expand the probability density. In addition, a bra (ket) vector is normalized by calculating the scalar product with the corresponding ket (bra) vector as is shown in Eq. (20), but the normalization of the probability density is determined by $\int dx\rho(x,t)=1$.

III Adiabatic expansion

There exist various expansion methods to find the solution of the Fokker-Planck equation. See Refs. hesam ; dreger ; talkner ; welles2014 and references therein. In our method, we expand the solution of the Fokker-Planck equation in terms of the eigenfunctions of the time-dependent Fokker-Planck operator which are obtained in the previous section. Expanding $|\rho(t)\rangle$ in terms of $|n,a_{t}\rangle$, we find

[TABLE]

where

[TABLE]

Here $t_{i}$ is an initial time and we introduced the notation

[TABLE]

Substituting this into the Fokker-Planck equation, we find the equation for the coefficient,

[TABLE]

On the other hand, operating the time derivative to Eq. (2), we find

[TABLE]

Solving this, for $m\neq n$, $\langle m,a_{t}|n,\dot{a}_{t}\rangle$ is reexpressed as

[TABLE]

Using this expression, we obtain the following differential equation

[TABLE]

which can be formally solved as

[TABLE]

From this equation, we can determine the coefficient iteratively.

This is similar to the derivation of the wave function in the adiabatic motion in quantum mechanics, and then $\theta_{n}(t)$ can be identified with the quantity corresponding to Berry’s geometrical phase berry .

It should be noted that the coefficient $C_{0}(t)$ is not affected by the second term on the right hand side of Eq. (27), because

[TABLE]

Here we used that $\langle 0,a_{t}|x\rangle$ is constant due to Eq. (15). Therefore $C_{0}(t)$ is given by the initial condition itself,

[TABLE]

Moreover, $\theta_{0}(t)$ is calculated as

[TABLE]

with $a_{i}\equiv a_{t_{i}}$.

III.1 Evolution from equilibrium state

In the following, we limit our discussion to the evolution from the equilibrium state defined by

[TABLE]

Then the initial coefficient is given by

[TABLE]

In this case, the expansion of $|\rho(t)\rangle$ is simplified as

[TABLE]

where the new coefficient $D_{n}(t)$ is determined by the following equation,

[TABLE]

Note that, because of Eq. (28),

[TABLE]

Then the expectation values can be expressed symbolically as

[TABLE]

where we introduced the pseudo density matrix defined by

[TABLE]

In this derivation, Eq. (15) is used.

One can easily confirm that the pseudo density matrix satisfies

[TABLE]

See also Eq. (42). Moreover the diagonal component is given by the solution of the Fokker-Planck equation,

[TABLE]

Differently from quantum mechanics, however, $\hat{\rho}(t)\neq\int dx\rho(x,t)|x\rangle\langle x|$, because, as indicated by Eq. (36), the pseudo density matrix does not describe the expectation value with the so-called mixed state in quantum mechanics.

The above representation of the expectation value with the pseudo density matrix is possible for the evolution from the equilibrium state. In more general cases, it should be expressed as

[TABLE]

IV Irreversible work

Our system is in equilibrium with the control parameter $a_{i}$ at the initial time $t_{i}$. Now we change this parameter so as to take the value $a_{f}\equiv a_{t_{f}}$ at the final time $t_{f}$. The work associated with this process is calculated through the expectation value of the change of the potential energy induced by its deformation. Then the mean work in this process is calculated as sekimoto-sasa ; sekimoto-book ; schmiedl-seifert2007

[TABLE]

On the other hand, by using Eq. (31), the expansion coefficient is iteratively expressed as

[TABLE]

Therefore the mean work is expanded in terms of $\dot{a}_{t}$.

The lowest order calculation of the mean work is given by substituting $D_{n}(t)$ with $\delta_{n,0}$ in the above equation, and then we find

[TABLE]

where $F(a_{t})$ is the Helmholtz free energy defined by

[TABLE]

In this derivation, we used

[TABLE]

Note that this simple relation is not satisfied for the ket vector, $|0,a_{i}\rangle\neq\sqrt{\frac{Z(a_{t})}{Z(a_{i})}}|0,a_{t}\rangle$.

This result corresponds to the work in the quasi-static process which is given by the difference of the Helmholtz free energy. Thus the irreversible work appears from the higher order terms in $D_{n}(t)$.

Calculating $D_{n}(t)$ up to the first order of $\dot{a}_{t}$, the irreversible work is expressed as

[TABLE]

where

[TABLE]

A similar but different formula to calculate the irreversible work is proposed in Refs.  sekimoto-book ; sekimoto-sasa , where the solution of the Fokker-Planck equation is expanded with the method which reminds us of the Chapman-Enskog expansion used for the Boltzmann equation. To reproduce the same result as Refs.  sekimoto-book ; sekimoto-sasa , we need to ignore the off-diagonal contributions in $\Lambda(t^{\prime},t)$, and replace the matrix element calculated from the derivative of the time-dependent Fokker-Planck operator ($-\langle n,a_{t}|(\partial_{a}\hat{\cal L}_{t})|0,a_{t}\rangle$) with that of the potential $(\langle n,a_{t}|(\partial_{a}\hat{V}(a_{t}))|0,a_{t}\rangle)$ in Eq. (44). Nevertheless, as is shown in Sec. V.3.2, our formula gives the same result as Ref. sekimoto-sasa when it is applied to the harmonic potential in the quasi-static limit.

IV.1 Sum rule

Besides Eqs. (35), (36) and (37), our pseudo density matrix $\hat{\rho}(t)$ satisfies another mathematical relation. When we consider the situation where the system reaches another equilibrium state at $t=t_{f}$ with $a_{f}$. Because then $D_{m}(t_{f})=0$ for $m\neq 0$ by definition, we find the following sum rule,

[TABLE]

In this derivation, we used the following mathematical property,

[TABLE]

It should be noticed that this relation is not the Jarzynski equality itself jar . As a matter of fact, the above expectation value $\langle~{}~{}\rangle$ is not the probability distribution of the work. However, if we consider a very short time evolution, the change of the particle position will be negligibly small and the expectation value with the work distribution might be identified with that of the (initial) distribution of the particle. Then, the above sum rule can be regarded as a special case of the Jarzynski equality.

V Application to harmonic potential

The irreversible work can be calculated exactly for the case of the harmonic potential as is shown in Ref. schmiedl-seifert2007 . To see the consistency of our expansion, we apply our result to this case.

V.1 Eigenfunctions

Then the time-dependent Fokker-Planck operator is given by

[TABLE]

This can be transformed as

[TABLE]

where

[TABLE]

We introduce the lowering and raising operators,

[TABLE]

and then the Hamiltonian operator is expressed as

[TABLE]

One can easily confirm that the following commutation relations are satisfied,

[TABLE]

That is, the eigenfunctions are constructed by operating the raising operators to the ground state as is the case of quantum mechanics. The eigenvalue of ${\cal H}_{t}$ is given by $\lambda_{n}(t)=a_{t}n/\nu$, and the corresponding eigenfunction is expressed using the Hermite polynomials as

[TABLE]

For the properties of the Hermite polynomials $H_{n}(x)$, see Appendix A. Therefore, from Eqs. (8) and (9), the bi-orthogonal system is constructed as

[TABLE]

with the eigenvalue,

[TABLE]

V.2 Irreversible work

Substituting the above bi-orthogonal system to our formula of the irreversible work (43), we obtain

[TABLE]

In this derivation, we used the following relations,

[TABLE]

The derivations are shown in Appendix A.

For the harmonic oscillator, higher order corrections do not contribute and Eq. (50) gives the exact expression of the irreversible work. Because of Eq. (53), the coefficient which can contribute to the calculation of $W_{irr}$ is only $D_{2}(t)$. The $N$-th order correction term in the calculation of $D_{2}(t)$ contains the following product of the matrix elements,

[TABLE]

On the other hand, from Eq. (51), the matrix $\langle m,a_{t}|\dot{\hat{\cal L}}_{t}|n,a_{t}\rangle$ has a finite contribution only when $m=n+2$ is satisfied, because $n=m$ is excluded in the sum. Therefore the contributions which are higher order than $N=1$ vanish.

V.3 Comparison with results from other works

V.3.1 irreversible work for harmonic potential

The exact irreversible work is calculated in Refs. schmiedl-seifert2007 ; schmiedl-seifert2008 ; gomez with the moment method, where the variance of the Brownian particle is introduced as

[TABLE]

and then the mean work is expressed as

[TABLE]

This is equivalent to our result as is shown below.

Note that, by using the partial integration formula, Eq. (50) can be reexpressed as

[TABLE]

On the other hand, as is shown in Ref. schmiedl-seifert2007 , the variance satisfies the differential equation,

[TABLE]

In short, solving this with the equilibrium initial condition where $\omega(t_{i})=1/(\beta a_{i})$, Eq. (50) can be cast into the following form,

[TABLE]

V.3.2 Small and large time-scale limits of external operation $\tau_{op}$

Let us investigate the behavior of the irreversible work when the control parameter $a_{t}$ is changed slowly. To clarify this limit, we introduce the adimensional time variable sekimoto-book ; sekimoto-sasa as

[TABLE]

where $\tau_{op}=t_{f}-t_{i}$ characterizes the time scale of the external operation. Then Eq. (50) is reexpressed as

[TABLE]

where $a_{t}$ is represented as a function of $\tau$, $\bar{a}_{\tau}=a_{\tau\tau_{op}+t_{i}}$.

For the consistency check of the later numerical calculations, it should be noted that, in the limit of the instantaneous jump where $\tau_{op}\rightarrow 0$, the irreversible work is simply given by

[TABLE]

This is the exact result independent of the choice of the protocol $\bar{a}_{\tau}$.

In the large limit of $\tau_{op}$, the integral for $\tau^{\prime}$ has the dominant contribution from $\tau^{\prime}\sim\tau$ because of the exponential factor, and thus Eq. (56) is approximately given by

[TABLE]

Here we used the approximation which is justified for the large $\tau_{op}$ limit,

[TABLE]

This approximated expression for the irreversible work is the same as that in Ref. sekimoto-sasa , which is obtained by using different formula for the irreversible work.

The above result (58) means that the irreversible work is proportional to $\tau_{op}^{-1}$ for the large $\tau_{op}$. It is worth emphasizing that this prediction is experimentally verified. See Ref. blickle2 for details.

V.4 Optimization

In the Fokker-Planck equation, we can show that the mean work is always larger than $\Delta F$ sekimoto-book ,

[TABLE]

Therefore there exists an optimized control protocol which minimizes the irreversible work for a given time period $\tau_{op}$.

To find the optimized protocol, we calculate the variation of the irreversible work for the following change of $a_{t}$,

[TABLE]

with the fixed initial and final values,

[TABLE]

After some of algebra, we find

[TABLE]

See Appendix B for more details.

It is difficult to determine the optimized protocol by solving this equation exactly. To find an approximated solution, we reexpress this as

[TABLE]

where

[TABLE]

Taking the large $\tau_{op}$ limit, we can employ the same approximation used to obtain the approximated irreversible work (58). Then, for $t_{i}<t<t_{f}$, Eq. (60) is reduced to

[TABLE]

and, near the initial and final times $t\sim t_{i},t_{f}$, we obtain

[TABLE]

The solution of Eq. (62) is

[TABLE]

which is the same result obtained in Refs. sekimoto-sasa ; schmiedl-seifert2007 for the large $\tau_{op}$ limit.

In Fig. 1, we plotted the irreversible works calculated using the control parameter given by Eq. (64). On the right and left panels, we consider the compression process choosing $\bar{a}_{1}=2\bar{a}_{0}$ and the expansion process choosing $\bar{a}_{0}=2\bar{a}_{1}$, respectively. We introduced ${\tau}_{*}=\nu/\sqrt{\bar{a}_{0}\bar{a}_{1}}$ to express the results in adimensional quantities. To see the applicability of our approximation, the two different irreversible works are plotted. The solid line represents the exact result (56), while the dashed line is the approximated one (58), which is expressed as

[TABLE]

One can observe that both irreversible works decrease as $\tau_{op}$ increases, because the large $\tau_{op}$ corresponds to the quasi-static limit and the contribution from the irreversible work disappears. For $\tau_{op}\gtrsim\tau_{*}$, our approximation is in good agreement with the exact one, while the approximation overestimates the irreversible work for the smaller $\tau_{op}$. In particular, the exact result shows that the magnitude of the irreversible work has a finite upper bound and does not diverge even in the limit of the instantaneous jump, $\tau_{op}\longrightarrow 0$. The same behavior is found also in Ref. schmiedl-seifert2007 . The above mentioned upper bound is given by Eq. (57).

In Ref. schmiedl-seifert2007 , the optimized protocol is obtained by calculating the variation of the irreversible work for the variance, not for the control parameter. The influence of the change of the variational variables is discussed in Sec. VI.

V.5 Fluctuation of irreversible work per unit $\tau$

The calculations of other correlation functions are simplified in our approach compared to the moment method because we do not need to solve the coupled differential equations of the moments. As an example, we calculate the fluctuation of the irreversible work per unit $\tau$, which is given by $\langle(d\hat{V}(\bar{a}_{\tau})/d\tau)^{2}\rangle$. Then we define the quantity which characterizes the fluctuation as

[TABLE]

where $\Delta\Omega_{re}(\tau)$ is the fluctuation in the quasi-static process, which is defined by

[TABLE]

This fluctuation can be calculated exactly for the harmonic oscillator. From the structure of the matrix element $\langle 0,\bar{a}_{0}|\hat{x}^{4}|0,\bar{a}_{\tau}\rangle$, the coefficients which contribute to the calculation are only $D_{2}(\tau)$ and $D_{4}(\tau)$. These are calculated exactly by the first and second order iterations of Eq. (31), respectively. Then we find the exact representation of the fluctuation of the irreversible work per unit $\tau$ as

[TABLE]

where

[TABLE]

and

[TABLE]

As is the case of the irreversible work, we can calculate the simple form of this fluctuation in the vanishing limit of $\tau_{op}$ as

[TABLE]

This is the exact result independent of the choice of $\bar{a}_{\tau}$.

We consider that the control parameter is again given by Eq. (64). Then the numerical results of the fluctuations are shown in Fig. 2. The left and right panels represent the compression process $\bar{a}_{1}=2\bar{a}_{0}$ and the expansion process $\bar{a}_{0}=2\bar{a}_{1}$, respectively. Each line represents the result for $\tau_{op}/\tau_{*}=0$, $1$, $5$, $10$, $20$ and $50$ from the most top line, respectively. Although it is not plotted, $\Delta\Omega_{irr}/\Delta\Omega_{re}$ disappears at $\tau=0$. One can see that the magnitudes of the fluctuations do not diverge and take finite values even for the instantaneous jump process.

Note that the time evolutions of the fluctuations are not monotonic and there exist peaks for lines of $\tau_{op}\gtrsim 5\tau_{*}$ for the compression. Such peaks, however, may disappear when we use the optimized parameter by exactly solving Eq. (59).

VI Concluding remarks

In this paper, we developed the framework for the systematic expansion of the solution of the Fokker-Planck equation with the help of the eigenfunctions of the time-dependent Fokker-Planck operator. The expansion parameter is the time derivative of the external parameter which controls the form of the external confinement potential. Our expansion corresponds to the perturbative calculation of the adiabatic motion in quantum mechanics. With this result, we derived a new formula to calculate irreversible work order by order, which is expressed as the expectation value with the pseudo density matrix which can describe only the expectation values with pure states.

By applying this to the harmonic potential, we confirmed that the first order calculation gives the exact irreversible work, which is the consistent result with that of the moment method schmiedl-seifert2007 ; schmiedl-seifert2008 ; gomez . By taking the large $\tau_{op}$ limit, we further verified that our formula reproduces the result of Refs. sekimoto-book ; sekimoto-sasa where the solution of the Fokker-Planck equation is expanded with the method which reminds us of the Chapman-Enskog expansion. However, the structure of the formula in Refs. sekimoto-book ; sekimoto-sasa is qualitatively different from ours and such a coincidence will not be seen when it is applied to other potentials.

Higher perturbative corrections can be calculated systematically in our approach and the accuracy of the prediction is improved up to arbitrary order, differently from the formula proposed in Refs. sekimoto-book ; sekimoto-sasa . Moreover, the calculations of correlation functions are simplified compared to the moment method because we do not need to solve the coupled differential equations of the moments. As an example, we showed the calculation of the fluctuation of the irreversible work per unit $\tau$. These are advantages of our approach.

Because the irreversible work is expressed as an analytic function of the control parameter $a_{t}$, we can apply the variational procedure to find the optimized protocol which minimizes the irreversible work. The derived equation, however, is a complex integro-differential equation and difficult to be solved. Instead, we discussed the procedure to find the approximated solution. As is shown in Fig. 1, this approximation reproduces the exact behavior of the irreversible work for the large $\tau_{op}$ region, and coincides with the results in Refs. sekimoto-sasa ; schmiedl-seifert2007 .

There is a comment for the optimization of the control parameter. In Ref. schmiedl-seifert2007 , the optimized parameter for the case of the harmonic potential is obtained by calculating the variation of the irreversible work with respect to the variance $\omega(t)$ of the position of the Brownian particle, instead of the control parameter $a_{t}$ itself. However, optimizations generally depend on the choice of the variational variables and we cannot change the variables without justification. Of course, in analytical mechanics, the variational procedure is known to be independent of the choice of coordinates, but it is because the variable transformation in such a case is given by the local function of time. That is, if $\omega(t)$ can be expressed as a function only of $a_{t}$, the variation for $\omega(t)$ leads to the same result as that for $a_{t}$. In the present case, however, $\omega(t)$ is the solution of the differential equation (54) and depends on the hysteresis of $a_{t}$. Then the variation for $\omega(t)$ does not necessarily coincide with that for $a_{t}$. For example, let us consider the modification of the optimized function by adding a term which is given by a function of $\omega(t_{f})$. Clearly the variation of this term with respect to $\omega(t)$ vanishes and the optimization is not affected. However, the added term induces another contribution in the variation with respect to $a_{t}$ because of the hysteresis of $a_{t}$ in $\omega(t)$. In the present optimization problem, therefore, we should consider the variation for $a_{t}$, not for $\omega(t)$.

For the sake of simplicity, we introduced the bra-ket notation then the Fokker-Planck equation (19) is expressed as if it is invariant for the choice of the representation. Then it might be possible to obtain a master equation by multiplying a discretized complete basis to Eq. (19), instead of $\langle x|$. Master equations are considered to be important in chemical reactions sekimoto-book and the corresponding eigenvalue problem for the case of time-periodic perturbations is discussed in Ref. caceres .

The present approach can be extended to more general potentials. Because it is generally difficult to find the analytic expressions of the eigenfunctions in non-linear potentials, we need to introduce another expansion to find the analytic forms of the expansion basis in our formula. For example, we often consider time-periodic protocols and then it will be useful to apply the method in Ref. shirley where the eigenvalue theory for a time-periodic Hamiltonian operator in quantum mechanics is developed. Moreover, when we consider systems in the higher spatial dimension, degenerated eigenvalues appear and then the present expansion method should be modified. These applications are left as future works.

In the present calculation, we have considered that the control parameter $a_{t}$ is a smooth deterministic function of time. However it is more realistic to consider the fluctuation of such a protocol and then we need the optimization with respect to stochastic variables. The stochastic variation has been discussed in the formulations of the quantum theory and hydrodynamics, but the applicability to the optimization of the irreversible work is still an open problem svm1 ; svm2 ; svm3 ; svm4 . It is also interesting whether we can apply the similar argument to relativistic kk-rel and quantum systems qse-koide .

The author acknowledges the ICE group of the Institute of Physics of URFJ for fruitful discussions and comments. This work is financially supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), project 307516/2015-6.

Appendix A Important formulae for calculations of harmonic potential

In the following we derive various formulae with the help of the Hermite polynomials which satisfy

[TABLE]

A.1 $\langle x|\dot{\hat{\cal L}}_{t}|n,a_{t}\rangle$

[TABLE]

where

[TABLE]

From the second to the third line, we used Eq. (75). Therefore

[TABLE]

A.2 $\langle n,a_{t}|n,\dot{a}_{t}\rangle$

[TABLE]

where we introduced

[TABLE]

A.3 $\langle 0,a_{i}|\frac{1}{2}\hat{x}^{2}|n,a_{t}\rangle$

[TABLE]

Appendix B Variation

The change of the irreversible work for the variation of $a_{t}$ is calculated as

[TABLE]

where

[TABLE]

Therefore the optimized protocol is described by the following equation,

[TABLE]