Optimization of stochastic thermodynamic machines

TL;DR

This paper optimizes the performance of stochastic thermodynamic machines by deriving optimal external potentials that minimize irreversible work and entropy production, providing explicit bounds and illustrative examples.

Contribution

It introduces a method to optimize stochastic thermodynamic machines using Fokker-Planck equations, variational techniques, and characteristics, deriving bounds on work, power, and efficiency.

Findings

Optimal external potentials minimize irreversible work and entropy production.

Explicit bounds for work output, power, and efficiency are derived.

Examples demonstrate the application of optimal protocols.

Abstract

The study of stochastic thermodynamic machines is one of the main topics in nonequilibrium thermodynamics. In this study, within the framework of Fokker-Planck equation, and using the method of characteristics of partial differential equation as well as the variational method, performance of stochastic thermodynamic machines is optimized according to the external potential, with the irreversible work $W_{irr}$, or the total entropy production $\Delta S_{\rm tot}$ equivalently, reaching its lower bound. Properties of the optimal thermodynamic machines are discussed, with explicit expressions of upper bounds of work output $W$, power $P$, and energy efficiency $\eta$ are presented. To illustrate the results obtained, typical examples with optimal protocols (external potentials) are also presented.

Equations

$${\partial }_{\tau }p(x,\tau )=-{\partial }_{x}J(x,\tau )=-{\partial }_x[p(x,\tau )u(x,\tau )],$$

$$Q_h={\int }_X\left({\int }_0^{t_h}{\partial }_{\tau }p(x,\tau )V(x,\tau )d\tau \right)dx.$$

$$V(x,\tau )={\int }_0^x\frac{[\xi {\int }_0^z{\partial }_{\tau }p(y,\tau )dy]-k_B{T}_h{\partial }_{z}p(z,\tau )}{p(z,\tau )}dz+C(\tau )$$

$$\begin{array}{rl}{Q}_h=& -{\int }_0^{x_0}\left({\int }_0^{t_h}{\partial }_{x}J(x,\tau )V(x,\tau )d\tau \right)dx\\ =& {\int }_0^{x_0}\left({\int }_0^{t_h}J(x,\tau ){\partial }_{x}V(x,\tau )d\tau \right)dx\\ \stackrel{(\ast )}{=}& -k_{B}T{\int }_0^{x_0}\left({\int }_0^{t_h}J(x,\tau ){\partial }_x[\ln p(x,\tau )]d\tau \right)dx\\ & -{\int }_0^{x_0}\left({\int }_0^{t_h}{\xi }_{h}J(x,\tau )u(x,\tau )d\tau \right)dx\\ =& k_{B}T{\int }_0^{x_0}\left({\int }_0^{t_h}{\partial }_{x}J(x,\tau )\ln p(x,\tau )d\tau \right)dx\\ & -{\int }_0^{x_0}\left({\int }_0^{t_h}{\xi }_{h}p(x,\tau )u^2(x,\tau )d\tau \right)dx\\ =:& -k_{B}T{\int }_0^{x_0}\left({\int }_0^{t_h}{\partial }_{\tau }p(x,\tau )\ln p(x,\tau )d\tau \right)dx-W_{irr}^h\\ \stackrel{(\#)}{=}& -k_{B}T{\int }_0^{x_0}\left({\int }_0^{t_h}{\partial }_{\tau }[p(x,\tau )\ln p(x,\tau )]d\tau \right)dx-W_{irr}^h\\ =:& T_h\Delta S-W_{irr}^h,\end{array}$$

$${\int }_0^{x_0}\left({\int }_0^{t_h}p(x,\tau ){\partial }_{\tau }[\ln p(x,\tau )]d\tau \right)dx={\int }_0^{x_0}\left({\int }_0^{t_h}{\partial }_{\tau }p(x,\tau )d\tau \right)dx=0.$$

$$W_{irr}^h:=$$

$$=$$

$${\partial }_{\tau }f(x,\tau )+u(x,\tau ){\partial }_{x}f(x,\tau )=0.$$

$$\frac{df(x,\tau )}{d\tau }$$

$$={\partial }_{\tau }f(x,\tau )+u(x,\tau ){\partial }_{x}f(x,\tau )$$

$$=0.$$

$$W=$$

$$=$$

$$\delta W_{irr}^h={\int }_0^{t_h}{\int }_0^{x_0}2{\xi }_h\left[\frac{{\partial }_{\tau }f}{{\partial }_{x}f}{\partial }_x\left(\frac{{\partial }_{\tau }f}{{\partial }_{x}f}\right)-{\partial }_{\tau }\left(\frac{{\partial }_{\tau }f}{{\partial }_{x}f}\right)\right]\delta fdxd\tau ,$$

$$\delta f(0,\tau )=\delta f(x_0,\tau )=\delta f(x,0)=\delta f(x,t_h)=0.$$

$$\left(\frac{{\partial }_{\tau }{f}^{\ast }}{{\partial }_x{f}^{\ast }}\right){\partial }_x\left(\frac{{\partial }_{\tau }{f}^{\ast }}{{\partial }_x{f}^{\ast }}\right)-{\partial }_{\tau }\left(\frac{{\partial }_{\tau }{f}^{\ast }}{{\partial }_x{f}^{\ast }}\right)=0.$$

$${\partial }_{\tau }{u}^{\ast }(x,\tau )+u^{\ast }(x,\tau ){\partial }_x{u}^{\ast }(x,\tau )=0.$$

$$\frac{d{u}^{\ast }(x,\tau )}{d\tau }={\partial }_{\tau }{u}^{\ast }(x,\tau )+\frac{dx}{d\tau }{\partial }_x{u}^{\ast }(x,\tau )=0.$$

$$\Gamma (x)=f_1^{-1}(f_0(x)).$$

$$x(z,\tau )=z+u^{\ast }(x,\tau )\tau =z+[\Gamma (z)-z]\tau /t_h.$$

$$W_{irr}^{h\ast }=$$

$$=$$

$$W_{irr}^{c\ast }=$$

$$=$$

$$S(t_h)=$$

$$=$$

$$\Delta S=S(t_h)-S(0)=k_B{\int }_0^{x_0}{p}_0(x)\ln {\Gamma }^{\prime }(x)dx.$$

$$Q_h^{\ast }=k_B{T}_h{\int }_0^{x_0}{p}_0(x)\ln {\Gamma }^{\prime }(x)dx-\frac{{\xi }_h}{t_h}{\int }_0^{x_0}{p}_0(x)[\Gamma (x)-x{]}^{2}dx,$$

$$W^{\ast }=$$

$$=$$

Full text

††journal: J. Stat. Phys.

1 Introduction

How we can improve the performance of thermodynamic machines is always an interesting but difficult problem in thermodynamics since the pioneer work of Carnot and Clausius Carnot1824 ; Clausius1856 . According to the second law of thermodynamics, the energy efficiency, defined as $\eta=W/Q_{h}$ with $W$ the output work and $Q_{h}$ the heat uptake from the hot heat bath, is not more than the Carnot efficiency $\eta_{C}=1-T_{c}/T_{h}$, where $T_{h}$ and $T_{c}$ are absolute temperatures of hot and cold heat baths, respectively. However, it is commonly thought the Carnot efficiency $\eta_{C}$ can only be achieved in quasistatic limit with output power vanished. So, in recent decades, most studies turned to discuss the efficiency at (or near) maximum power (EMP) $\eta_{\max P}$. For example, the Curzon-Ahlborn efficiency at maximum power $\eta_{CA}=1-\sqrt{T_{c}/T_{h}}$ is obtained in Yvon1955 ; Chambadal1957 ; Novikov1958 ; Curzon1975 ; Broeck2005 . The bounds of EMP, $\eta_{C}/2\leq\eta_{\max P}\leq\eta_{C}/(2-\eta_{C})$, is obtained in Schmiedl2008Efficiency1 ; Esposito2010Efficiency ; Izumida2011Efficiency . By methods of linear irreversible thermodynamics, general expressions for maximum power and maximum efficiency are obtained in Proesmans2016 . The upper bound of efficiency at arbitrary power is discussed in Holubec2016 ; Ryabov2016 ; Pietzonka2018 . For more studies which try to optimize the performance of thermodynamic machines, see Benenti2011 ; Izumida2011Efficiency ; Golubeva2012Efficiency ; Allahverdyan2013 ; Holubec2015 ; Calvo2015 ; Reyes2017 ; Polettini2017 . Or see Sekimoto2010 ; Seifert2012Stochastic ; Mart2016Colloidal ; Giuliano2017 for general reviews. Meanwhile, methods to operate heat engines infinitely close to the Carnot bound $\eta_{C}$ but with nonzero power are also suggested recently Benenti2011 ; Allahverdyan2013 ; Mart2016Brownian ; Lee2017Carnot ; Polettini2017 ; Holubec2018 , though there are also studies which show that it is impossible Hondou2000Unattainability ; Verley2014The .

One of the main difficulties in the study of stochastic thermodynamic machines is that, except for few special cases, no expressions for work $W$, power $P$ and efficiency $\eta$ can be obtained explicitly Schmiedl2007 ; Schmiedl2008Efficiency1 ; Holubec2014 ; Holubec2015 ; Zhang2019 . So, in references, they are usually optimized by choosing optimal work period $t$ to get maximum power, or just through extensive numerical calculations to get optimal values of thorough protocols Then2008 ; Horowitz2018 . In this study, using variational method, the external time-dependent potential, as well as the boundary distributions of isothermal steps in each work cycle, will be optimized generally to achieve the optimal performance of stochastic thermodynamic machines.

The same as in Schmiedl2008Efficiency1 ; Holubec2014 ; Zhang2019 , this study considers a thermodynamic machine with potential $V(x,\tau)$, which depends on spatial (state) variable $x$ and time variable $\tau$. For overdamped cases, the time evaluation of probability $p(x,\tau)$ to find the system at position $x$ at time $\tau$ is governed by the following Fokker-Planck equation

[TABLE]

with $J(x,\tau):=p(x,\tau)u(x,\tau)$ the flux of probability, $u(x,\tau)$ the instantaneous velocity given by $u(x,\tau):=-\partial_{x}[V(x,\tau)+k_{B}T\ln p(x,\tau)]/\xi$, and $\xi$ the friction coefficient, which satisfies Einstein relation $k_{B}T=\xi D$. Here $k_{B}$ is the Boltzmann constant, $T$ is the absolute temperature, and $D$ is the free diffusion constant.

The thermodynamic machine discussed in this study is assumed to work cyclicly, including two isothermal processes and two adiabatic transitions. In each work cycle, it performs sequently through the following four subprocesses Schmiedl2008Efficiency1 ; Zhang2019 . (1) Isothermal process with high temperature $T_{h}$ during time interval $0<\tau<t_{h}$. (2) Adiabatic transition (instantaneously) from high temperature $T_{h}$ to low temperature $T_{c}$ at time $\tau=t_{h}$. (3) Isothermal process with low temperature $T_{c}$ during time interval $t_{h}<\tau<t_{h}+t_{c}$. (4) Adiabatic transition from low temperature $T_{c}$ to high temperature $T_{h}$ at time $\tau=t_{h}+t_{c}$. The same as in Schmiedl2008Efficiency1 ; Holubec2014 ; Zhang2019 , adiabatic transitions are idealized as sudden jumps of potential, and assumed to occur instantaneously without heat exchange. Therefore, probability $p(x,\tau)$ and system entropy do not change during adiabatic transitions. For convenience, the work period is denoted as $t:=t_{h}+t_{c}$.

It is obvious that most of physical quantities of the thermodynamic machine, including output work $W$ in each work cycle, power $P=W/t$, and energy efficiency $\eta$, are functions of duration $t_{h}$ (or $t_{c}$), and period $t=t_{h}+t_{c}$, or even functions of temperatures $T_{h/c}$, friction coefficients $\xi_{h/c}$, and diffusion constants $D_{h/c}$, and so can be optimized according to them, as have been done previously Schmiedl2008Efficiency1 ; Esposito2010Efficiency ; Izumida2011Efficiency ; Holubec2016 ; Golubeva2012Efficiency ; Allahverdyan2013 ; Van2005Thermodynamic ; Cleuren2009Universality ; Van2012Efficiency ; Wang2013Efficiency ; Hooyberghs2013Efficiency . However, it is no doubt that $W,P$, and $\eta$ are also functionals of external potential $V(x,\tau)$ (or called protocols sometimes), probabilities $p(x,0)=p(x,t)$ and $p(x,t_{h})$. One main contribution of this study is that the performance of thermodynamic machine is optimized according to potential $V(x,\tau)$, as well as to probabilities $p(x,0)$ and $p(x,t_{h})$.

2 Methods of optimization

As stated above, this study assumes that the work cycle begins from the isothermal process with high temperature $T_{h}$. By definition, the heat uptake from the hot heat bath during one work cycle is

[TABLE]

In the following, we assume $X=[0,x_{0}]$. For other cases, such as $X=[0,+\infty)$, or $(-\infty,+\infty)$, results are similar, and corresponding examples will be presented at the end of this study. Here $x_{0}>0$ is a positive real number, which is used to describe the scale of the thermodynamic system, see Sec. 2.2 for some details.

From Eq. (1), and with the no flux boundary conditions at $x=0$ and $x_{0}$, it can be verified that potential $V(x,\tau)$ can be formally reformulated as

[TABLE]

for $0<\tau<t_{h}$ and $0\leq x\leq x_{0}$, with $C(\tau)$ an arbitrary function of time $\tau$. Substituting $V(x,\tau)$ into Eq. (2), and through routine calculations, one can get

[TABLE]

where the equality $(*)$ is from the definition of instantaneous velocity $u(x,\tau):=-\partial_{x}[V(x,\tau)+k_{B}T\ln p(x,\tau)]/\xi$, and the equality $(\#)$ is because

[TABLE]

The entropy production $\Delta S:=S(t_{h})-S(0)$ with system entropy $S(\tau):=-k_{B}\int_{0}^{x_{0}}p(x,\tau)\ln p(x,\tau)dx$, and the irreversible work $W^{h}_{irr}$ is defined by

[TABLE]

In Eq. (5), $f(x,\tau):=\int_{0}^{x}p(z,\tau)dz$ is the distribution function of the system at time $\tau$, which satisfies (see Eq. (1) and the no flux boundary condition at $x=0$)

[TABLE]

Eq. (6) implies that, along its characteristic curves, which are given by $dx/d\tau=u(x,\tau)$, distribution function $f(x,\tau)$ satisfies

[TABLE]

Which shows that $f(x,\tau)$ is constant along any characteristic curve of Eq. (6). For simplicity, this study assumes that, for any $0\leq\tau\leq t$, $f(x,\tau)$ is reversible as a function of $x$, or equivalently, probability $p(x,\tau)$ at any given time $\tau$ is positive at almost everywhere. Or, in other words, there is no interval $[x_{1},x_{2}]$ with $0\leq x_{1}<x_{2}\leq x_{0}$, such that $p(x,\tau)\equiv 0$ in it, i.e., $p(x,\tau)=0$ at the most for a zero measure set Loeb2016 . But all results obtained in this study hold true for any general cases (note that, more complicated details need to be added to discuss the general cases).

Similarly, heat exchange between thermodynamic machine and cold heat bath in one work cycle can be written as $Q_{c}=-T_{c}\Delta S-W^{c}_{irr}$, with irreversible work $W^{c}_{irr}$ defined similarly. By definition, the output work during one work cycle is

[TABLE]

Two strategies to improve work $W$ will be presented in this study. (I) For given probabilities $p_{0}(x):=p(x,0)$ and $p_{1}(x):=p(x,t_{h})$, i.e. for given $\Delta S$, find the minimum value of irreversible work $W^{h/c}_{irr}$ by optimization according to potential $V(x,\tau)$. (II) For given probability $p_{0}(x)$, find the maximum value of $\Delta S$ by optimization according to probability $p_{1}(x)$.

2.1 Optimization of distribution function $f(x,\tau)$

From Eq. (5), it can be shown that the variation of irreversible work $W^{h}_{irr}$ according to the distribution function $f(x,\tau)$ is as follows,

[TABLE]

where $\delta f$ is an arbitrary variation of $f(x,\tau)$, which satisfies

[TABLE]

The reason of these zero boundary conditions is that $f(x,\tau)$ satisfies $f(0,\tau)\equiv 0,f(x_{0},\tau)\equiv 1$, and $f(x,0)=\int_{0}^{x}p_{0}(z)dz,f(x,t_{h})=\int_{0}^{x}p_{1}(z)dz$.

Since $\delta W^{h}_{irr}=0$ at the minimum of $W^{h}_{irr}$ (as a functional of distribution function $f(x,\tau)$) for any variation $\delta f$, Eq. (8) indicates that, for the optimal distribution function $f(x,\tau)=f^{*}(x,\tau)$,

[TABLE]

Due to Eq. (6), $\partial_{\tau}f/\partial_{x}f=-u$. Therefore,

[TABLE]

Eq. (9) shows that, for the optimal distribution function $f^{*}(x,\tau)$, the slope $u^{*}(x,\tau)$ of its characteristic curves satisfies

[TABLE]

It means that, the slope $u^{*}(x,\tau)$ of any characteristic curves of the optimal distribution function $f^{*}(x,\tau)$ is constant. Therefore, characteristic curves of $f^{*}(x,\tau)$, which are given by $dx/d\tau=u(x,\tau)$, are all straight lines.

For convenience, we denote $f_{0}(x):=f(x,0)$, $f_{1}(x):=f(x,t_{h})$, and define a map $\Gamma:x\in[0,x_{0}]\rightarrow\Gamma(x)\in[0,x_{0}]$ which satisfies $f_{0}(x)=f_{1}(\Gamma(x))$, or equivalently

[TABLE]

Note, for extreme cases that at least one of functions $f_{0}(x)$ and $f_{1}(x)$ is irreversible, map $\Gamma(x)$ can also be well defined but with some complicated and tedious details, so will not be presented here. Obviously, map $\Gamma$ satisfies $\Gamma(0)=0,\Gamma(x_{0})=x_{0}$, and it increases monotonically with $x$, since $\Gamma^{\prime}(x)=p_{0}(x)/p_{1}(\Gamma(x))>0$.

The above analysis shows that, for optimal distribution function $f^{*}(x,\tau)$, the characteristic curve of Eq. (6), which is initiated from any $z\in[0,x_{0}]$ at time $\tau=0$, satisfies

[TABLE]

Which means that the characteristic curve of Eq. (6), which begins from any $z\in[0,x_{0}]$ at time $\tau=0$, is a straight line connecting point $(z,0)$ and point $(\Gamma(z),t_{h})$ in the 2-dimensional $x$-$\tau$ plane. For any given $z\in[0,x_{0}]$ at time $\tau=0$, the position $x(z,\tau)$ obtained by Eq. (10) satisfies $f^{*}(x(z,\tau),\tau)=f_{0}(z)\equiv f_{1}(\Gamma(z))$, see Eq. (6). For convenience, we denote the inverse function of $x(z,\tau)$ by $z(x,\tau)$. Then, for any given position $x\in[0,x_{0}]$ and time $0\leq\tau\leq t_{h}$, $z(x,\tau)$ satisfies $f^{*}(x,\tau)=f_{0}(z(x,\tau))\equiv f_{1}(\Gamma(z(x,\tau)))$. Finally, Eq. (10) is actually obtained by definition from the following ordinary differential equation, $dx(\tau)/d\tau=u^{*}(x,\tau)$, with $x(0)=z$. While $u^{*}(x,\tau)\equiv[\Gamma(z)-z]/t_{h}$ for the optimal distribution function $f^{*}(x,\tau)$.

Using function $x(z,\tau)$ and its inverse function $z(x,\tau)$, probability $p^{*}(x,\tau)$ corresponding to the optimal distribution function $f^{*}(x,\tau)$ can be obtained by $p^{*}(x,\tau)=\partial_{x}f^{*}(x,\tau)=\partial_{x}f_{0}(z(x,\tau))=p_{0}(z(x,\tau))\partial_{x}z(x,\tau)$. From Eq. (5), it can be shown that, the irreversible work $W^{h}_{irr}$ corresponding to the optimal distribution function $f^{*}(x,\tau)$ is

[TABLE]

Here, the third equality is obtained by change of variable $z=z(x,\tau)$.

Due to $f_{0}(y)=f_{1}(\Gamma(y))$, $p_{0}(y)=p_{1}(\Gamma(y))\Gamma^{\prime}(y)$ or $p_{1}(y)=p_{0}(\Gamma^{-1}(y))/\Gamma^{\prime}(\Gamma^{-1}(y))$. So, it can be shown that

[TABLE]

Here, the last equality is obtained by change of variable $y=\Gamma(z)$. Meanwhile, it can be verified that

[TABLE]

So, by definition,

[TABLE]

From Eqs. (4,11,13), with the optimal distribution function $f^{*}(x,\tau)$,

[TABLE]

and from Eqs. (7,11,12,13),

[TABLE]

where $\bar{t}_{h}:=t_{h}/t,\bar{t}_{c}:=t_{c}/t$. Note, $\bar{t}_{h}+\bar{t}_{c}\equiv 1$.

Generally, in one work cycle, the total entropy production of the system can be obtained by $\Delta S_{\rm tot}=-\int_{0}^{t}[\int_{X}\partial_{\tau}p(x,\tau)V(x,\tau)dx]/T(\tau)d\tau$ Brandner2015 ; Proesmans2015 . Within the framework of Fokker-Planck equation, $\Delta S_{\rm tot}$ can be reformulated as $\Delta S_{\rm tot}=\int_{0}^{t}[\int_{X}\xi p(x,\tau)u^{2}(x,\tau)dx]/T(\tau)d\tau$. For the particular case with optimal distribution function $f^{*}(x,\tau)$, it can be shown that the lower bound of total entropy production is

[TABLE]

One can easily show that, with optimal distribution function $f^{*}(x,\tau)$,

[TABLE]

where $\Sigma_{h,c}=\left.k_{B}\int_{0}^{x_{0}}p_{0}(x)[\Gamma(x)-x]^{2}dx\right/D_{h,c}$. Eq. (16) looks the same as the one assumed in discussions of the low dissipation cases Esposito2010Efficiency ; Izumida2011Efficiency ; Calvo2015 ; Gonzalezayala2016 ; Reyes2017 . However, Eq. (16) holds for any durations $t_{h}$ and $t_{c}$, and $\Sigma_{h,c}$ may be very large, though they are corresponding to the least values of entropy production. While in discussions of low dissipation cases, expressions like Eq. (16) are only the first order of their real values, and are reasonable for larger durations $t_{h}$ and $t_{c}$.

2.2 Properties of thermodynamic machine with optimal distribution function $f^{*}(x,\tau)$

To discuss the influence of the scale of $x_{0}$, normalized probabilities $\bar{p}_{0/1}(x)$, normalized distribution functions $\bar{f}_{0/1}(x)$, and the corresponding map $\bar{\Gamma}(x)$ are defined as follows,

[TABLE]

with $0\leq x\leq 1$. It can be easily verified that $\int_{0}^{1}\bar{p}_{0/1}(x)=1$, $\bar{\Gamma}(0)=0$, $\bar{\Gamma}(1)=1$, $\bar{f}_{0/1}(0)=0$, $\bar{f}_{0/1}(1)=1$, and $\bar{f}^{\prime}_{0/1}(x)=\bar{p}_{0/1}(x)$, $\bar{f}_{0}(x)=\bar{f}_{1}(\bar{\Gamma}(x))$.

From Eqs. (14,17),

[TABLE]

Similarly, from Eqs. (15,17),

[TABLE]

Note that $\bar{t}_{c}=1-\bar{t}_{h}$. Eq. (20) indicates that the value of work $W^{*}$ depends on parameters $t,\bar{t}_{h}$ (or $\bar{t}_{c}$), as well as functions $\bar{\Gamma}(x)$ and $\bar{p}_{0}(x)$, i.e. $W^{*}=W^{*}(\bar{t}_{h/c},t;\bar{\Gamma}(x),\bar{p}_{0}(x))$.

From Eq. (20), the stalling time of the machine, defined by $W^{*}(t^{\rm st})=0$ (see Schmiedl2008Efficiency1 ), is as follows,

[TABLE]

For given functions $\bar{\Gamma}(x)$, $\bar{p}_{0}(x)$, and period $t$, the work $W^{*}$ reaches its maximum when $\bar{t}_{h/c}=\sqrt{\xi_{h/c}}/(\sqrt{\xi_{h}}+\sqrt{\xi_{c}})$. This can be easily obtained by minimizing ${\xi_{h}}/{\bar{t}_{h}}+{\xi_{c}}/{\bar{t}_{c}}$ under constraint $\bar{t}_{h}+\bar{t}_{c}=1$. It can be shown that $\max_{0<\bar{t}_{h},\bar{t}_{c}<1;\bar{t}_{h}+\bar{t}_{c}=1}W^{*}=$

[TABLE]

It can also be shown that power $P(t):=W^{*}/t$ reaches its maximum when period $t=t^{*}:=2t^{\rm st}$ (see Appendix). The corresponding work for period $t=t^{*}$ is

[TABLE]

and the heat absorbed from hot heat bath is

[TABLE]

So the efficiency at maximum power is

[TABLE]

As has been found previously, $\eta_{C}/2<\eta^{*}_{\max P}<\eta_{C}/(2-\eta_{C})$ Schmiedl2008Efficiency1 ; Esposito2010Efficiency ; Izumida2011Efficiency ; Johal2017Heat ; Giuliano2017 ; Reyes2017 . It tends to lower bound $\eta_{C}/2$ if $(t_{c}\xi_{h})/(t_{h}\xi_{c})\to 0$, and tends to upper bound $\eta_{C}/(2-\eta_{C})$ if $(t_{h}\xi_{c})/(t_{c}\xi_{h})\to 0$. For particular cases with $\bar{t}_{h/c}=\sqrt{\xi_{h/c}}/(\sqrt{\xi_{h}}+\sqrt{\xi_{c}})$, with which work $W^{*}$ is maximized ($t^{*}$ and $t^{\rm st}$ are minimized), $\eta^{*}_{\max P}=\eta_{C}/[2-\eta_{C}\sqrt{\xi_{h}}/(\sqrt{\xi_{h}}+\sqrt{\xi_{c}})]$. If $\xi_{h}=\xi_{c}$, $\eta^{*}_{\max P}=\eta_{C}/[2-\eta_{C}t_{c}/(t_{c}+t_{h})]$. It can be easily known that, for $t_{h}=t_{c}=t/2$, the optimal probability $p^{*}(x,\tau)$ is mirror symmetrical to time $t/2$, i.e., $p^{*}(x,t/2-\tau)=p^{*}(x,t/2+\tau)$ for $0\leq\tau\leq t/2$. For such cases, $\eta^{*}_{\max P}=\eta_{C}/[2-\eta_{C}\xi_{h}/(\xi_{h}+\xi_{h})]$.

It can be verified that $\eta^{*}_{\max P}\leq\eta_{CA}$ if and only if $[(\xi_{c}/t_{c})/(\xi_{h}/t_{h})]^{2}\geq T_{c}/T_{h}$. This condition is equivalent to the one given in Schmiedl2008Efficiency1 , where the Curzon-Ahlborn efficiency $\eta_{CA}$ is recovered if $\alpha=\alpha_{CA}\equiv 1/(1+\sqrt{T_{c}/T_{h}})$, with parameter $\alpha=(\xi_{h}/t_{h})/(\xi_{h}/t_{h}+\xi_{c}/t_{c})$ in this study, see the above Eq. (25) and the Eq. (31) in Schmiedl2008Efficiency1 . For special cases with $\bar{t}_{h}=\bar{t}_{c}=1/2$, inequality $[(\xi_{c}/t_{c})/(\xi_{h}/t_{h})]^{2}\geq T_{c}/T_{h}$ reduces to $(\xi_{c}/\xi_{h})^{2}\geq T_{c}/T_{h}$. While for cases with $\bar{t}_{h/c}=\sqrt{\xi_{h/c}}/(\sqrt{\xi_{h}}+\sqrt{\xi_{c}})$, it reduces to $\xi_{c}/\xi_{h}\geq T_{c}/T_{h}$. So, with optimal distribution function $f^{*}(x,\tau)$ and optimal ratio of duration $t_{h}/t_{c}=\sqrt{\xi_{h}/\xi_{h}}$, and if diffusion constant $D_{h}=D_{c}$, which means $T_{c}/T_{h}=\xi_{c}/\xi_{h}$, then the efficiency at maximum power $\eta^{*}_{\max P}$ is exactly the Curzon-Ahlborn efficiency $\eta_{CA}$.

With optimal distribution function $f^{*}(x,\tau)$, the efficiency for any cyclic period $t>t^{\rm st}$ is, see Eqs. (18,20),

[TABLE]

where $\hat{t}_{h/c}:=x_{0}^{2}/(2D_{h/c})=\xi_{h/c}x_{0}^{2}/(2k_{B}T_{h/c})$ is the characteristic time of free diffusion, $\Omega(\bar{p}_{0},\bar{\Gamma})$ is a functional of $\bar{p}_{0}(x)$ and $\bar{\Gamma}(x)$,

[TABLE]

Eq. (26) shows that efficiency $\eta^{*}$ increases with both durations (fractions) $t_{h}(\bar{t}_{h})$ and $t_{c}(\bar{t}_{c})$, while decreases with both characteristic times $\hat{t}_{h}$ and $\hat{t}_{c}$. Meanwhile, $\eta^{*}$ increases with cyclic period $t$ but decreases with scale parameter $x_{0}$. Moreover, it can be shown from Eq. (26) that

[TABLE]

where the equality holds if $\Omega(\bar{p}_{0},\bar{\Gamma})=0$, or $t\to\infty$.

Denote $P^{*}:=W^{*}(t^{*})/t^{*}\equiv W^{*}_{\max P}/t^{*}$ as the maximum of power $P(t):=W^{*}(t)/t$, then for any $t\geq t^{\rm st}$,

[TABLE]

which equals to zero when $t=t^{\rm st}=t^{*}/2$ or $t\to\infty$, and reaches it maximum 1 when $t=t^{*}$, see Eq. (A2) in Appendix. Meanwhile, efficiency $\eta^{*}$ can be reformulated as (see Eq. (A3) in Appendix)

[TABLE]

It can be easily shown that $d\eta^{*}/dt>0$, so efficiency $\eta^{*}$ increases monotonically from 0 to $\eta_{C}$ as the cyclic period $t$ increases from stalling time $t^{\rm st}=t^{*}/2$ to infinity. For $t=t^{*}$, i.e., when power $P(t)$ reaches its maximum $P^{*}$, efficiency $\eta^{*}$ get its value $\eta_{\max P}^{*}$ as given in Eq. (25).

From Eqs. (27, 28), the ratio ${P(t)}/{P^{*}}$ can be reformulated as a function of efficiency $\eta^{*}$,

[TABLE]

which equals to zero if $\eta^{*}=0$ or $\eta^{*}=\eta_{C}$, and reaches its maximum 1 when $\eta^{*}=\eta_{\max P}^{*}$ as given in Eq. (25). It can be shown that ${P(t)}/{P^{*}}$ increases with $\eta^{*}$ for $0\leq\eta^{*}\leq\eta_{\max P}^{*}$, while decreases with $\eta^{*}$ for $\eta_{\max P}^{*}\leq\eta^{*}\leq\eta_{C}$.

It can also be shown that (see Appendix),

[TABLE]

Therefore, $W^{*}(t)/W^{*}_{\max P}$ increases from 0 to 2 as period $t$ increases from *stalling time * $t^{\rm st}=t^{*}/2$ to infinity, or equivalently as efficiency $\eta^{*}$ increases from 0 to $\eta_{C}$. $W^{*}(t)/W^{*}_{\max P}=1$, i.e., $W^{*}(t)=W^{*}_{\max P}$, if and only if $t=t^{*}$, or $\eta^{*}=\eta_{\max P}^{*}$ equivalently. As implied by Eq. (23), $W^{*}_{\max P}=W^{*}(t\to\infty)/2$.

Finally, from Eqs. (27, 30), it can be obtained that

[TABLE]

So, ratio $P(t)/P^{*}$ increases first and then decreases with ratio $W^{*}(t)/W^{*}_{\max P}$.

2.3 Optimal probability density $p_{1}(x)$

From Eq. (23), the variation of $W^{*}_{\max P}$ according to map $\bar{\Gamma}(x)$ is,

[TABLE]

where $\delta\bar{\Gamma}(x)$ is a (small) arbitrary variation of function $\bar{\Gamma}(x)$. Due to boundary conditions $\bar{\Gamma}(0)=0$ and $\bar{\Gamma}(1)=1$, variation $\delta\bar{\Gamma}(x)$ satisfies $\delta\bar{\Gamma}(0)=\delta\bar{\Gamma}(1)=0$.

Eq. (32) indicates that $W^{*}_{\max P}$ reaches its maximum when $\bar{\Gamma}(x)=\int_{0}^{x}\bar{p}_{0}(z)dz$, or equivalently when $\Gamma(x)=x_{0}\int_{0}^{x}p_{0}(z)dz$, see Eq. (17). Notice that, $\bar{\Gamma}(x)=\int_{0}^{x}\bar{p}_{0}(z)dz$ means $\bar{p}_{1}(x)\equiv 1$ and $\bar{f}_{1}(x)=x$, while $\Gamma(x)=x_{0}\int_{0}^{x}p_{0}(z)dz$ means $p^{*}_{1}(x)\equiv 1/x_{0}$ and $f^{*}_{1}(x)\equiv x/x_{0}$. In any way,

[TABLE]

2.4 Explanatory notes

As mentioned before, work $W$ is not only a function of period $t$ and duration $t_{h}$ (or $t_{c}$), but also a functional of the potential $V(x,\tau)$ and probability $p_{0}(x)$ (or $p_{1}(x)$ equivalently). So there are many ways to do optimization. Optimization according to distribution function $f(x,\tau)$ is actually equivalent to optimization according to potential $V(x,\tau)$, and optimization according to map $\bar{\Gamma}(x)$ is equivalent to optimization according to probability $p_{1}(x)$. In references, ratio $t_{h}/t_{c}$ is usually optimized to get the maximum of work $W$, while period $t$ is usually optimized to get the maximum of power $P$. However, in this study, for any given work period $t$, the potential $V(x,\tau)$ is optimized to get the minimum of entropy production $\Delta S_{\rm tot}$, or equivalently to get the maximum of work $W$.

From Eq. (17), $\int_{0}^{1}\bar{p}_{0}(x)\ln\bar{p}_{0}(x)dx=\ln x_{0}+\int_{0}^{x_{0}}p_{0}(x)\ln p_{0}(x)dx$. Therefore, $\max_{\bar{\Gamma}(x)}W^{*}_{\max P}$ increases (logarithmically) with scale parameter $x_{0}$ if the initial entropy $-k_{B}\int_{0}^{x_{0}}p_{0}(x)\ln p_{0}(x)dx$ is fixed, see Eq. (33). However, to keep power $P$ to be maximum, the optimal period $t^{*}=2t^{\rm st}$ should increase like $x_{0}^{2}$, see Eq. (21). Therefore, the corresponding value of maximum power $[\max_{\bar{\Gamma}(x)}W^{*}_{\max P}]/t^{*}$ actually decreases with scale parameter $x_{0}$.

For convenience, we denote the spatial coordinate of a point on the characteristic curve of Eq. (6), which begins from $z$ at time $\tau=0$, by $\Gamma(z,\tau)$, or $\Gamma_{\tau}(z)$ for simplicity. That is to say, $\Gamma(z,\tau)\equiv\Gamma_{\tau}(z):=x(z,\tau)$ with $x(z,\tau)$ satisfies Eq. (10). Particularly, by definition $\Gamma(x)=\Gamma(x,t_{h})$. In the following, methods to get the optimal probability $p^{*}(x,\tau)$ and optimal potential $V^{*}(x,\tau)$, which correspond to the optimal distribution function $f^{*}(x,\tau)$, will be presented for only $0\leq\tau\leq t_{h}$. For $t_{h}\leq\tau\leq t=t_{h}+t_{c}$, methods are almost the same.

By definition, $df(x,\tau)/d\tau=0$ along any characteristic curves of Eq. (6). This means that distribution function $f(x,\tau)$ satisfies $f_{0}(x)\equiv f(x,0)=f(\Gamma(x,\tau),\tau)=f(\Gamma(x,t_{h}),t_{h})\equiv f_{1}(\Gamma(x,t_{h}))$ for any $0\leq\tau\leq t_{h}$. Therefore, $f(x,\tau)=f_{0}(\Gamma_{\tau}^{-1}(x))=\int_{0}^{\Gamma_{\tau}^{-1}(x)}p_{0}(z)d\,z$. By taking derivatives on both side of $f_{0}(x)=f(\Gamma(x,\tau),\tau)$, one gets $p_{0}(x)=p(\Gamma(x,\tau),\tau)\partial_{x}\Gamma_{\tau}(x)$. So

[TABLE]

Particularly, for $\tau=t_{h}$, $p_{1}(x)=p_{0}(\Gamma^{-1}(x))/\Gamma^{\prime}(\Gamma^{-1}(x))$.

By definition $-\partial_{x}[V(x,\tau)+k_{B}T\ln p(x,\tau)]/\xi_{h}=u(x,\tau)$, so $V(x,\tau)=-\xi_{h}\int_{0}^{x}u(z,\tau)d\,z-k_{B}T\ln p(x,\tau)+C(\tau)$, with $C(\tau)$ an arbitrary function of time $\tau$. For cases with optimal distribution function $f^{*}(x,\tau)$, the instantaneous velocity $u^{*}(z,\tau)=[\Gamma(\Gamma_{\tau}^{-1}(z))-\Gamma_{\tau}^{-1}(z)]/t_{h}$. Therefore, the optimal potential is

[TABLE]

Here the last equality is obtained from Eqs. (10, 34) and $\Gamma(z,\tau)=x(z,\tau)$.

Since thermodynamic machines work cyclically, function $C(\tau)$ should satisfies $C(\tau+t)=C(\tau)$ for any time $\tau$. One can easily find that irreversible works $W_{irr}^{h/c*}$ are independent of the choice of function $C(\tau)$, so in all the following illustrative examples, we always set $C(\tau)\equiv 0$ simply.

3 Illustrative examples with optimal control

To illustrate the results obtained in this study, we present several examples of thermodynamic machine with the optimal distribution function $f^{*}(x,\tau)$, or equivalently with the optimal potential $V^{*}(x,\tau)$.

3.1 Examples for $0\leq x\leq x_{0}$ with $x_{0}$ a finite number

(I) Let

[TABLE]

one can show that

[TABLE]

From Eq. (10) and the definition of $\Gamma(x,\tau)$, we have

[TABLE]

Therefore,

[TABLE]

and potential $V^{*}(x,\tau)$ can be obtained by Eq. (35). For this particular case, the work output in one cyclic period (see Eq. (15)) is

[TABLE]

and the efficiency (see Eq. (26)) is

[TABLE]

with

[TABLE]

(II) Next, for

[TABLE]

it can be shown that $f_{0}(x)=\sqrt{x}$, $f_{1}(x)=x$, and

[TABLE]

So

[TABLE]

Similarly, $V^{*}(x,\tau)$ can be obtained by Eq. (35). For this particular case, the work (see Eq. (15)) is

[TABLE]

and the efficiency (see Eq. (26)) is

[TABLE]

with

[TABLE]

3.2 An example for $0\leq x<\infty$

We now give an example for $x_{0}=\infty$. Let $t_{h}=1$, and

[TABLE]

i.e., both $p_{0}(x)$ and $p_{1}(x)$ are Weibull distributions. Then for the optimal case, one can get $\Gamma(x)=\lambda_{1}x/\lambda_{0}$, and

[TABLE]

Consequently,

[TABLE]

with $\lambda(\tau)=(1-\tau)\lambda_{0}+\tau\lambda_{1}$. From Eq. (23), the work at maximum power is

[TABLE]

which decreases with $\lambda_{0}$ while increases with $\lambda_{1}$.

For $\lambda_{0}\to 0$, $p_{0}(x)$ tends to Dirac delta function $\delta(x)$. While for $\lambda_{1}\to\infty$, $p_{1}(x)$ tends to a constant. This is consistent with previous discussion for finite scale parameter $x_{0}$ cases, in which maximum of $W^{*}_{\max P}$ is reached when $p_{1}(x)\equiv 1/x_{0}$ is constant. As before, the optimal potential $V^{*}(x,\tau)$ can be obtained by Eq. (35). The work for this particular case is (see Eq. (15))

[TABLE]

with ${\bf\Gamma}(x)$ the Gamma function defined as

[TABLE]

For given value of parameter $\lambda_{0}$ in the expression of initial probability $p_{0}(x)$, there exists an optimal parameter

[TABLE]

with which work $W^{*}$ reaches its maximum. The efficiency for these particular cases (see Eq. (26) for definition) is also $\eta^{*}=\eta_{C}-(1-\eta_{C})\varepsilon$ but with

[TABLE]

One can easily show that $\varepsilon$ increases with $\lambda_{1}$ for $\lambda_{1}\geq\lambda_{0}$, and $\eta^{*}\to\eta_{C}$ when $\lambda_{1}\to\lambda_{0}$, but with which the work $W^{*}\to 0$.

3.3 An example for $-\infty\leq x<\infty$

Finally, we give an example in which state (spatial) variable $x$ is defined in the whole real domain, i.e., $x\in(-\infty$, $+\infty$). As before, for convenience, let $t_{h}=1$, but assume the following initial and final probabilities

[TABLE]

Usually, distributions with probability like these are called Laplace distributions. For this example, $\Gamma(x)=\beta_{1}(x-\alpha_{0})/\beta_{0}+\alpha_{1}$, and

[TABLE]

So the optimal distribution function

[TABLE]

with

[TABLE]

Here $\textrm{sgn}(x)$ is the sign function. It can be shown that the optimal probability is

[TABLE]

For these particular cases, the output work is

[TABLE]

Compared with the results for Weibull distribution with parameter $k=1$ (Note, $\mathfrak{F}(3)=2$), an extra term $(\xi_{h}/t_{h}+\xi_{c}/t_{c})(\alpha_{1}-\alpha_{0})^{2}$ is added to the irreversible work $W^{*}_{irr}$, see Eqs. (11, 12, 36, 37). This is due to the spatial translocation of the system from one mean position $\int_{-\infty}^{\infty}xp_{0}(x)dx=\alpha_{0}$ to another mean position $\int_{-\infty}^{\infty}xp_{1}(x)dx=\alpha_{1}$. If $\alpha_{0}=\alpha_{1}$, then properties of work $W^{*}$ and efficiency $\eta^{*}$ are the same as the ones discussed in the above subsection for Weibull distribution with parameter $k=1$.

4 Conclusions and Remarks

In this study, optimization methods for stochastic thermodynamic machines are discussed, and general properties of optimal machines are analyzed. To illustrate the results, several explicit examples with optimal protocol are also presented. The main idea of optimization is to find the optimal external potential, with which the heat dissipation into environment reaches its minimum. It should be noted that, within the theoretical framework of Langevin stochastic processes and using the method of Monge-Kantorovich optimal mass transport, similar results about the minimum of heat dissipation have been obtained generally in Aurell2011 ; Aurell2012 ; Bo2013 . But differently, in this study, the minimum of heat dissipation is obtained by the method of characteristics and variation, and within the framework of Fokker-Planck equation, which makes it easier to construct illustrating examples and to know more clearly what protocol (external potential) we should choose to reach the optimal one.

It seems that several results obtained in this study are similar as the ones have been presented in previous references, including those for efficiency and work at maximum power. But with the optimal distribution function $f^{*}(x,\tau)$, and consequently with the optimal potential $V^{*}(x,\tau)$, the period $t^{*}=2t^{\rm st}$ required to get the maximum power is less than those given in previous references. In fact, this study indicates that the period $t^{*}=2t^{\rm st}$, see Eq. (21), required to get the maximum power is the least one. In other words, the thermodynamic machines designed in this study, with optimal potential $V^{*}(x,\tau)$, can output work more fast and have the smallest stalling time $t^{\rm st}$ than any other ones. One of the main contributions of this study is that, the least value of irreversible work $W_{irr}^{*}$, or equivalently the least value of entropy production $\Delta S_{\rm tot}$, dissipated in one work cycle is obtained, and the optimal protocol which can lead to this lower bound is also presented.

Appendix: Details of analysis on maximum power

For convenience, in this appendix, the following notations will be employed,

[TABLE]

With these notations, and from Eqs. (14, 15, 21), one can easily show that stalling time $t^{\rm st}=A/B$, optimal period $t^{*}=2t^{\rm st}=2A/B$, work $W^{*}(t)=B-A/t$, heat $Q^{*}_{h}=B/\eta_{C}-AC/t$, power $P(t)=W^{*}(t)/t=B/t-A/t^{2}$, and the maximum power $P^{*}=P(t^{*})=W^{*}(t^{*})/t^{*}=B^{2}/4A$. Therefore,

[TABLE]

which equals to zero when $t=t^{\rm st}=t^{*}/2$ or $t\to\infty$, and reaches its maximum 1 for $t=t^{*}$. Meanwhile, with notations in Eq. (A1), the efficiency $\eta^{*}=W^{*}/Q_{h}^{*}$ can be written as $\eta^{*}=(B-A/t)/(B/\eta_{C}-AC/t)$, which gives that

[TABLE]

Since $C\eta_{C}<1$, it can be shown that efficiency $\eta^{*}$ increases monotonically from 0 to $\eta_{C}$ as the period $t$ increases from $t^{\rm st}$ to infinity. For $t=t^{*}$, i.e., when power $P(t)$ reaches its maximum $P^{*}=B^{2}/4A$, efficiency $\eta^{*}$ get its value $\eta_{\max P}^{*}=\eta_{C}/(2-C\eta_{C})$.

From Eqs. (A2, A3), the ratio ${P(t)}/{P^{*}}$ can be rewritten as follows

[TABLE]

which equals to zero if $\eta=0$ or $\eta=\eta_{C}$, and reaches its maximum 1 if $\eta=\eta_{C}/(2-C\eta_{C})$.

From Eqs. (23, A1), the work at maximal power $P^{*}$ is $W^{*}_{\max P}=B/2$. So $W^{*}(t)/W^{*}_{\max P}=2-t^{*}/t=2(1-C\eta_{C})\eta^{*}/[(1-C\eta^{*})\eta_{C}]$. Obviously, $W^{*}(t)/W^{*}_{\max P}$ increases monotonically from 0 to 2 as period $t$ increases from *stalling time * $t^{\rm st}=t^{*}/2$ to infinity, or equivalently as efficiency $\eta^{*}$ increases from 0 to Carnot efficiency $\eta_{C}$.