Statistical mechanics of a nonequilibrium steady-state classical particle system driven by a constan
本站小编 Free考研考试/2022-01-02
Jie Yao(姚婕)1,2, Yanting Wang(王延颋),1,21 CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of Sciences, 55 East Zhongguancun Road, PO Box 2735, Beijing, 100190 China 2 School of Physical Sciences, University of Chinese Academy of Sciences, 19A Yuquan Road, Beijing 100049, China
Abstract A classical particle system coupled with a thermostat driven by an external constant force reaches its steady state when the ensemble-averaged drift velocity does not vary with time. In this work, the statistical mechanics of such a system is derived solely based on the equiprobability and ergodicity principles, free from any conclusions drawn on equilibrium statistical mechanics or local equilibrium hypothesis. The momentum space distribution is determined by a random walk argument, and the position space distribution is determined by employing the equiprobability and ergodicity principles. The expressions for energy, entropy, free energy, and pressures are then deduced, and the relation among external force, drift velocity, and temperature is also established. Moreover, the relaxation towards its equilibrium is found to be an exponentially decaying process obeying the minimum entropy production theorem. Keywords：nonequilibrium steady state;constant external force;dissipation;statistical mechanics
PDF (386KB)MetadataMetricsRelated articlesExportEndNote|Ris|BibtexFavorite Cite this article Jie Yao(姚婕), Yanting Wang(王延颋). Statistical mechanics of a nonequilibrium steady-state classical particle system driven by a constant external force. Communications in Theoretical Physics[J], 2020, 72(11): 115601- doi:10.1088/1572-9494/abb7d2
1. Introduction
The theoretical framework for equilibrium statistical mechanics is well established, but that for nonequilibrium, except the linear response regime, is still in its early stage. Until now, the fluctuation theorems [1 –6], represented by the Jarzynski equality [6], are the only set of equalities that hold for nonequilibrium systems arbitrarily far from equilibrium. The lack of theoretical frameworks for nonequilibrium statistical mechanics hinders the investigations to real systems since most of the time they work at nonequilibrium states. A simpler case is the nonequilibrium steady state (NESS) whose thermodynamic properties do not change with time while a certain flux exists. For an NESS, a generalized force constantly input energy into the system, which is exactly counter-balanced by a certain dissipation process, and thus the ensemble-averaged flux does not vary with time. An applied NESS system can be the electrolyte in a battery driven by the electrostatic force between two electrodes (see, e.g. [7, 8]). Surprisingly, there are no general theories even for the case as simple as NESS. A direct consequence is that molecular dynamics simulations of such kind of systems still have to approximately use the equilibrium thermostats, such as the Berendsen thermostat [9] and the Nosé–Hoover thermostat [10, 11], despite the fact that it is known to introduce, sometimes significant, errors into simulation results [12].
In this paper, we establish the statistical mechanics for a classical particle system driven by a constant external force and coupled with a constant-temperature thermostat, whose NESS is arbitrarily far from equilibrium. None of the conclusions drawn on equilibrium systems or nonequilibrium systems with local equilibrium hypothesis are adopted. That is, similar to the equilibrium case, the deduction is solely based on the equiprobability and ergodicity principles. As shown in figure 1, the system contains N identical classical particles coupled to a thermostat with a constant temperature T driven by a constant force F along the x -axis. After the system reaches the steady state, the macroscopic drift velocity Vd along the x -axis is a constant. Because the system is still canonical (the external potential is stationary), the momentum and position spaces are independent of each other. Note that the choice of a cylindrical shape in figure 1 is for the convenience of calculations, and the specific shape does not make any differences when the system approaches the thermodynamic limit.
Figure 1.
New window|Download| PPT slide Figure 1.The schematic of the studied classical particle system in the NESS with N particles driven by a constant force F and coupled to a thermostat with a constant temperature T, resulting in a constant drift velocity Vd .
In the following sections, we will first determine the momentum space and position space distributions, followed by the calculations of system energy, entropy, free energy, and pressures. The relation among the external force F, the drift velocity Vd, and the temperature T is then established. Finally, the relaxation process towards its equilibrium is revealed to be an exponential decay obeying the minimum entropy production theorem.
2. Momentum space distributions
Since for a canonical classical system, the momentum space is independent of the position space, the first question to ask is what the distributions in the momentum space, i.e. the velocity distributions, are. Since the force is applied along the x -axis, obviously the velocities along the y and z directions ${v}_{y}$ and ${v}_{z}$ should both still obey the equilibrium Maxwell–Boltzmann distribution$ \begin{eqnarray}\begin{array}{l}f\left({v}_{y}\right){\rm{d}}{v}_{y}=\sqrt{\displaystyle \frac{m}{2\pi {k}_{{\rm{B}}}T}}\exp \left(-\displaystyle \frac{m{v}_{y}^{2}}{2{k}_{{\rm{B}}}T}\right){\rm{d}}{v}_{y},\\ f\left({v}_{z}\right){\rm{d}}{v}_{z}=\sqrt{\displaystyle \frac{m}{2\pi {k}_{{\rm{B}}}T}}\exp \left(-\displaystyle \frac{m{v}_{z}^{2}}{2{k}_{{\rm{B}}}T}\right){\rm{d}}{v}_{z},\end{array}\end{eqnarray}$ where ${k}_{{\rm{B}}}$ is the Boltzmann constant and m is the particle mass. For a random velocity along the x direction ${v}_{x},$ we assume that it is generated by n -times random walks with the magnitude of $\delta $ for each walk, and the direction can be either positive or negative. If the probability for taking the positive direction is p, then the probability of having k times positive walks in all n walks is described by the binomial distribution$ \begin{eqnarray}P\left(k| n\right)={{C}}_{n}^{k}{p}^{k}{\left(1-p\right)}^{n-k}.\end{eqnarray}$ When $n\to \infty ,$ the above approaches the Gaussian distribution$ \begin{eqnarray}f\left(k\right){\rm{d}}k=\displaystyle \frac{1}{\sqrt{2\pi }\sigma }\exp \left(-\displaystyle \frac{{\left(k-\mu \right)}^{2}}{2{\sigma }^{2}}\right){\rm{d}}k,\end{eqnarray}$ where $\mu =np$ and ${\sigma }^{2}=np\left(1-p\right).$ The corresponding velocity along the x -axis is ${v}_{x}=\left(2k-n\right)\delta ,$ which reaches its maximum value of ${V}_{\max }=n\delta $ when $k=n$ and its minimum value of $-{V}_{\max }=-n\delta $ when $k=0.$ Because $\left\langle k\right\rangle =\mu =np,$ the thermodynamic drift velocity$ \begin{eqnarray}{V}_{{\rm{d}}}=\left\langle {v}_{x}\right\rangle =\left(2p-1\right){V}_{\max }.\end{eqnarray}$ Since $k=\tfrac{{v}_{x}}{2\delta }+\tfrac{n}{2},$ equation (3 ) becomes$ \begin{eqnarray}\begin{array}{l}f\left({v}_{x}\right){\rm{d}}{v}_{x}=\displaystyle \frac{\sqrt{n}}{2{V}_{\max }\sqrt{2\pi }\sqrt{p(1-p)}}\exp \\ \,\times \left(-\displaystyle \frac{{\left({v}_{x}-\left(2p-1\right){V}_{\max }\right)}^{2}}{8p\left(1-p\right){V}_{\max }^{2}/n}\right){\rm{d}}{v}_{x}.\end{array}\end{eqnarray}$ By requiring that both the expectation value and the variance of equation (5 ) should be constants when $n\to \infty ,$ we obtain ${V}_{\max }\propto \sqrt{n}$ and $p-\tfrac{1}{2}\propto \tfrac{1}{\sqrt{n}},$ so$ \begin{eqnarray}f\left({v}_{x}\right){\rm{d}}{v}_{x}=\displaystyle \frac{\sqrt{n}}{{V}_{\max }\sqrt{2\pi }}\exp \left(-\displaystyle \frac{{\left({v}_{x}-{V}_{{\rm{d}}}\right)}^{2}}{2{V}_{\max }^{2}/n}\right){\rm{d}}{v}_{x}.\end{eqnarray}$ When ${V}_{{\rm{d}}}=0,$ equation (6 ) should go back to the equilibrium Maxwell–Boltzmann distribution, so we obtain $\tfrac{{V}_{\max }^{2}}{n}=\tfrac{{k}_{{\rm{B}}}T}{m}.$ The velocity distribution along the x -axis is thus$ \begin{eqnarray}f\left({v}_{x}\right){\rm{d}}{v}_{x}=\sqrt{\displaystyle \frac{m}{2\pi {k}_{{\rm{B}}}T}}\exp \left(-\displaystyle \frac{m{\left({v}_{x}-{V}_{{\rm{d}}}\right)}^{2}}{2{k}_{{\rm{B}}}T}\right){\rm{d}}{v}_{x}.\end{eqnarray}$
The same result can be reached by utilizing the local equilibrium hypothesis [13]. Nevertheless, our derivation indicates that this result is of a more general basis free of the local equilibrium hypothesis.
By utilizing the ${\chi }^{2}$ -distribution, it is easy to show that the velocity perpendicular to be drift velocity ${v}_{\perp }=\sqrt{{v}_{y}^{2}+{v}_{z}^{2}}$ obeys the probability distribution$ \begin{eqnarray}f({v}_{\perp })=\displaystyle \frac{m{v}_{\perp }}{{k}_{{\rm{B}}}T}\exp \left(-\displaystyle \frac{m{v}_{\perp }^{2}}{2{k}_{{\rm{B}}}T}\right),\,{v}_{\perp }\geqslant 0.\end{eqnarray}$
The derivation is as follows. We define two new variables obeying the standard normal distribution $N\left(0,1\right):$$ \begin{eqnarray}\begin{array}{c}\begin{array}{l}{\bar{v}}_{y}\equiv \sqrt{\displaystyle \frac{m}{{k}_{{\rm{B}}}T}}{v}_{y},\\ {\bar{v}}_{z}\equiv \sqrt{\displaystyle \frac{m}{{k}_{{\rm{B}}}T}}{v}_{z},\end{array}\end{array}\end{eqnarray}$ then the random variable $a\equiv {\bar{v}}_{\perp }^{2}={\bar{v}}_{y}^{2}+{\bar{v}}_{z}^{2}$ obeys the ${\chi }^{2}$ distribution:$ \begin{eqnarray}\begin{array}{c}g\left(a\right)={\chi }^{2}\left(k=2\right)={\left.\frac{{a}^{\displaystyle \frac{k}{2}-1}\exp \left(-\displaystyle \frac{a}{2}\right)}{{2}^{\displaystyle \frac{k}{2}}{\rm{\Gamma }}\left(\displaystyle \frac{k}{2}\right)}\right|}_{k=2}\\ \,=\,\frac{1}{2}\exp \left(-\frac{a}{2}\right),\end{array}\end{eqnarray}$ so we have$ \begin{eqnarray}f\left({\bar{v}}_{\perp }\right)=-\frac{{\rm{d}}{\bar{v}}_{\perp }^{2}}{{\rm{d}}{\bar{v}}_{\perp }}g\left({\bar{v}}_{\perp }^{2}\right)=-{\bar{v}}_{\perp }\exp \left(-\frac{{\bar{v}}_{\perp }^{2}}{2}\right).\end{eqnarray}$ Because ${v}_{\perp }=\sqrt{{v}_{y}^{2}+{v}_{z}^{2}}$ = $\sqrt{\tfrac{{k}_{{\rm{B}}}T}{m}}{\bar{v}}_{\perp }$ and $f\left({v}_{\perp }\right)=-\tfrac{{\rm{d}}{\bar{v}}_{\perp }}{{\rm{d}}{v}_{\perp }}f\left({\bar{v}}_{\perp }\right),$ we finally obtain equation (8 ).
3. Position space distribution
The derivation of the position space distribution is analogous to the equilibrium case [14]. The energy of the NESS system is denoted as ${U}_{{\rm{t}}},$ and the energy of the whole system composed of the thermostat and the NESS system is denoted as U . Both the whole system and the thermostat are treated as infinitely-large microcanonical ensemble systems, whose numbers of microstates are ${\rm{\Omega }}\left(U\right)$ and ${\rm{\Omega }}\left(U-{U}_{{\rm{t}}}\right),$ respectively. The NESS system is also infinitely large, but infinitely small compared to the thermostat, so the energy dissipates to the thermostat is negligible. Note that in the following derivation, the conclusions from equilibrium statistical mechanics are applied to the whole system and the thermostat, but not directly to the NESS. Given the entropy of a microcanonical system $S={k}_{{\rm{B}}}\,\mathrm{ln}$ Ω, the probability that the NESS system has the total energy ${U}_{{\rm{t}}}$ is$ \begin{eqnarray}\begin{array}{l}P\left({U}_{{\rm{t}}}\right)=\displaystyle \frac{1}{Z}\displaystyle \frac{{\rm{\Omega }}\left(U-{U}_{{\rm{t}}}\right)}{{\rm{\Omega }}\left(U\right)}=\displaystyle \frac{1}{Z}\exp \\ \,\,\,\,\times \,\left(\displaystyle \frac{S\left(U-{U}_{{\rm{t}}}\right)-S\left(U\right)}{{k}_{{\rm{B}}}}\right).\end{array}\end{eqnarray}$ By utilizing the Taylor expansion and $T={\left(\tfrac{\partial U}{\partial S}\right)}_{{U}_{{\rm{t}}}=0},$ we have$ \begin{eqnarray}P\left({U}_{{\rm{t}}}\right)=\displaystyle \frac{1}{Z}\exp \left(-\beta {U}_{{\rm{t}}}\right),\end{eqnarray}$ where the partition function, serving as the normalization factor, is written as$ \begin{eqnarray}Z=\displaystyle \frac{1}{N!{h}^{3N}}\displaystyle \int {\rm{d}}{{\boldsymbol{r}}}^{N}\displaystyle \int {\rm{d}}{{\boldsymbol{p}}}^{N}\exp \left(-\beta {U}_{{\rm{t}}}\left({{\boldsymbol{r}}}^{N},{{\boldsymbol{p}}}^{N}\right)\right),\end{eqnarray}$ in which h is the Plank constant, ${{\boldsymbol{r}}}^{N}$ and ${{\boldsymbol{p}}}^{N}$ are the collections of particle positions and momenta, respectively. For a classical canonical system, the momentum space is independent of the position space, so the partition function $Z=Q{Z}_{{\rm{id}}},$ where $Q=\displaystyle \int {\rm{d}}{{\boldsymbol{r}}}^{N}\exp \left(-\beta {U}_{{\rm{p}}}\left({{\boldsymbol{r}}}^{N}\right)\right)$ is the configuration integral, and ${Z}_{{\rm{id}}}$ is the partition function of the ideal-gas NESS system. Because ${v}_{x}\in \left(-\infty ,+\infty \right),$ according to equation (7 ), the expression of ${Z}_{{\rm{id}}}$ is exactly as the partition function of the equilibrium ideal gas system.
Although temperature is well defined for equilibrium systems, currently its rigorous definition for nonequilibrium systems is still lack [15]. In this work, it should be emphasized that the temperature T and associated relations are defined for the equilibrium microcanonical thermostat as well as the whole system, not directly defined for the NESS system.
4. Entropy and free energy
Because the system entropy expression $S=-{k}_{{\rm{B}}}\displaystyle {\sum }_{i}P\left({U}_{i}\right)\mathrm{ln}\,P\left({U}_{i}\right)$ can be deduced solely based on the equiprobability and ergodicity principles, it is also applicable to the NESS system which has a well-define probability distribution. Because ${v}_{x}\in \left(-\infty ,+\infty \right),$ according to equation (7 ), the expression of $P\left({U}_{i}\right)$ is identical for both equilibrium and NESS systems, and thus their entropy is the same. According to free energy $A={U}_{{\rm{t}}}-TS,$ the free energy of a NESS system is ${E}_{{\rm{d}}}\,\equiv \tfrac{N}{2}m{V}_{{\rm{d}}}^{2}$ larger than its equilibrium counterpart, exactly as the total energy difference. In other words, ideally the drift energy ${E}_{{\rm{d}}}$ can all be used to do external work.
5. Injection-dissipation relation
Next we determine the relation among the external force F, the drift velocity ${V}_{{\rm{d}}},$ and the temperature T according to the balance between energy injection and dissipation. During a time interval ${\rm{\Delta }}t,$ the system kinetic energy increment caused by F is$ \begin{eqnarray}{\rm{\Delta }}{E}_{{\rm{K}}}=N\displaystyle \int \displaystyle \frac{m}{2}\left[{\left({v}_{x}+\displaystyle \frac{F}{m}{\rm{\Delta }}t\right)}^{2}-{v}_{x}^{2}\right]f\left({v}_{x}\right){\rm{d}}{v}_{x},\end{eqnarray}$ so the injection power$ \begin{eqnarray}w=\mathop{\mathrm{lim}}\limits_{{\rm{\Delta }}t\to \infty }\displaystyle \frac{{\rm{\Delta }}{E}_{{\rm{k}}}}{{\rm{\Delta }}t}=NF{V}_{{\rm{d}}}.\end{eqnarray}$
On the other hand, the dissipation power can be calculated by considering the energy exchange between the system and the thermostat. Let us first consider the energy exchange by perfect elastic collisions on the boarder of the cylindrical container along the perpendicular direction (see figure 1 ). The energy change of an NESS particle in a single-time collision can be obtained as ${\rm{\Delta }}u=\tfrac{m}{2}\left({{\boldsymbol{v}}}_{1}^{2}-{{\boldsymbol{v}}}^{2}\right),$ where ${\boldsymbol{v}}$ and ${{\boldsymbol{v}}}_{1}$ are the velocities of the NESS and thermostat particles, respectively. However, the input energy cannot be sufficiently dissipated if the system particles only exchange energy on the boarder, because the external force F is applied to all system particles; when the system size goes to infinity, the input power grows with the volume but the dissipation power grows with the area. Therefore, we have to assume that all system particles can experience ‘virtual’ collisions with thermostat in the whole system volume, and the collision times ${n}_{{\rm{c}}}$ for each particle in a unit time interval is proportional to the particle velocity in the direction perpendicular to the force:$ \begin{eqnarray}{n}_{{\rm{c}}}\left({v}_{\perp }\right)={c}_{0}{v}_{\perp },\end{eqnarray}$ where ${c}_{0}$ is a coefficient reflecting the coupling strength between the system and the thermostat. Because the number of particles with a perpendicular velocity of ${v}_{\perp }$ is $Nf\left({v}_{\perp }\right),$ the total collision times for all particles with the perpendicular velocity of ${v}_{\perp }$ in a time interval ${\rm{\Delta }}t$ is$ \begin{eqnarray}{M}_{{\rm{c}}}\left({v}_{\perp }\right)=Nf\left({v}_{\perp }\right){n}_{{\rm{c}}}\left({v}_{\perp }\right){\rm{\Delta }}t,\end{eqnarray}$ where the probability distribution of the perpendicular velocity is given by equation (8 ). Therefore, the dissipation power, which has the exact value of w but an opposite sign, can be calculated as$ \begin{eqnarray}\begin{array}{l}-w=\displaystyle {\int }_{0}^{\infty }{\rm{d}}{v}_{\perp }\displaystyle \frac{M({v}_{\perp })}{{\rm{\Delta }}t}\displaystyle \int {\rm{d}}{\boldsymbol{v}}f({\boldsymbol{v}})\displaystyle \int {\rm{d}}{{\boldsymbol{v}}}_{1}f({{\boldsymbol{v}}}_{1}){\rm{\Delta }}u\\ \,\,\,\,=\,-N{C}_{0}\sqrt{T}{V}_{{\rm{d}}}^{2},\end{array}\end{eqnarray}$ where ${C}_{0}\equiv {c}_{0}\sqrt{\tfrac{\pi m{k}_{{\rm{B}}}}{8}}$ is a coefficient whose exact value may vary with the coupling strength between the NESS system and the thermostat. Therefore, the relation among F, ${V}_{{\rm{d}}},$ and T is$ \begin{eqnarray}F={C}_{0}{V}_{{\rm{d}}}\sqrt{T}.\end{eqnarray}$
If the NESS system is composed of charged particles driven by a constant electric force $F=qE,$ where $q$ is the particle charge and $E$ is the constant external electric field, then the system mobility can be calculated according to equation (20 ) as$ \begin{eqnarray}\mu \equiv \displaystyle \frac{{V}_{{\rm{d}}}}{E}=\displaystyle \frac{q}{{C}_{0}\sqrt{T}}.\end{eqnarray}$ Moreover, because the friction coefficient $\xi $ can be related to the diffusion coefficient $D$ by [16]$ \begin{eqnarray}\xi =\displaystyle \frac{{k}_{{\rm{B}}}T}{mD}.\end{eqnarray}$ According to the Einstein relation, the diffusion coefficient is connected to the mobility by a linear relation:$ \begin{eqnarray}\mu =\displaystyle \frac{D}{{k}_{{\rm{B}}}T}.\end{eqnarray}$ Equations (21 )–(23 ) together determines the friction coefficient to be$ \begin{eqnarray}\xi =\displaystyle \frac{1}{m\mu }=\displaystyle \frac{{C}_{0}\sqrt{T}}{mq}.\end{eqnarray}$
6. Pressures
We then calculate the pressures of the ideal-gas NESS system. For a more general case, the contribution from particle interactions can be added in the same way as for equilibrium systems. Under the perfect elastic collision assumption, the momentum change for a single collision of the NESS particle on the boarder is ${\rm{\Delta }}p\left({v}_{\perp }\right)=2m{v}_{\perp }.$ According to the ergodicity and equiprobability principles, the ensemble-averaged particle density in the spatial space is uniform, so the times for particles colliding the container surface in ${\rm{\Delta }}t$ is$ \begin{eqnarray}\begin{array}{lll}{M}_{{\rm{s}}}\left({v}_{\perp }\right) & = & Nf\left({v}_{\perp }\right)\displaystyle \frac{{\rm{\Delta }}V\left({v}_{\perp }\right)}{V}\\ & = & Nf\left({v}_{\perp }\right)\,\,\displaystyle \frac{\pi \left[{R}^{2}-\left(R-{v}_{\perp }^{2}{\rm{\Delta }}{t}^{2}\right)\right]L}{\pi {R}^{2}L}\\ \, & \approx & \,Nf\left({v}_{\perp }\right)\displaystyle \frac{2{v}_{\perp }{\rm{\Delta }}t}{R}=\rho S{v}_{\perp }{\rm{\Delta }}tf\left({v}_{\perp }\right),\end{array}\end{eqnarray}$ where the particle density $\rho \equiv \tfrac{N}{V},$R and L are the radius and side length of the cylinder, respectively, and $S=2\pi RL$ and $V=\pi {R}^{2}L$ are the side area and volume of the cylinder, respectively, as shown in figure 1 . Then the pressure perpendicular to the flux is$ \begin{eqnarray}{P}_{\perp }=\displaystyle {\int }_{0}^{\infty }{\rm{d}}{v}_{\perp }\displaystyle \frac{{M}_{{\rm{s}}}\left({v}_{\perp }\right)}{S{\rm{\Delta }}t}{\rm{\Delta }}p\left({v}_{\perp }\right)=4\rho {k}_{{\rm{B}}}T.\end{eqnarray}$
Note that the above expression includes the pressure along four directions, namely ${y}^{+},$${y}^{-},$${z}^{+},$ and ${z}^{-}.$ It satisfies the equilibrium ideal gas equation of state along each direction. Along the flux direction, the times for a particle colliding a cylindrical surface is $M\left({v}_{x}\right){\rm{d}}{v}_{x}=\rho \pi {R}^{2}{\rm{\Delta }}t{v}_{x}f\left({v}_{x}\right){\rm{d}}{v}_{x},$ and the momentum change for each collision is ${\rm{\Delta }}p\left({v}_{x}\right)=2m{v}_{x},$ so the pressure along the ${x}^{+}$ direction is$ \begin{eqnarray}\begin{array}{l}{P}_{{x}^{+}}=\displaystyle \frac{1}{\pi {R}^{2}{\rm{\Delta }}t}\displaystyle {\int }_{0}^{\infty }{\rm{d}}{v}_{x}M\left({v}_{x}\right)2m{v}_{x}\\ \,\,\,\,\,\,\,\,=\,2m\rho \displaystyle {\int }_{0}^{\infty }{v}_{x}^{2}\sqrt{\displaystyle \frac{m}{2\pi {k}_{{\rm{B}}}T}}\exp \left(-\displaystyle \frac{m{\left({v}_{x}-{V}_{{\rm{d}}}\right)}^{2}}{2{k}_{{\rm{B}}}T}\right){\rm{d}}{v}_{x}.\end{array}\end{eqnarray}$ Similarly, the pressure along the ${x}^{-}$ direction is$ \begin{eqnarray}\begin{array}{c}\begin{array}{l}{P}_{x-}=2m\rho \displaystyle {\int }_{-\infty }^{0}{v}_{x}^{2}\sqrt{\displaystyle \frac{m}{2\pi {k}_{{\rm{B}}}T}}\exp \left(-\displaystyle \frac{m{\left({v}_{x}-{V}_{{\rm{d}}}\right)}^{2}}{2{k}_{{\rm{B}}}T}\right){\rm{d}}{v}_{x}.\end{array}\end{array}\end{eqnarray}$
The integrals in equations (27 ) and (28 ) do not have analytical solutions and can be solved numerically.
7. Relaxation towards equilibrium
If the constant external force F is applied to the system when time $t\lt 0$ and suddenly removed at $t=0,$ the system will relax from NESS towards its equilibrium state by dissipating the drift energy $\tfrac{N}{2}m{V}_{{\rm{d}}}^{2}\left(t=0\right)$ into the thermostat gradually. During an infinitesimal time interval ${\rm{d}}t,$ the energy change equals to the energy dissipation:$ \begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\left(\displaystyle \frac{N}{2}m{V}_{{\rm{d}}}^{2}\left(t\right)\right)=-w=-N{C}_{0}\sqrt{T}{V}_{{\rm{d}}}^{2}\left(t\right).\end{eqnarray}$ The solution of the above equation is$ \begin{eqnarray}{V}_{{\rm{d}}}\left(t\right)={V}_{{\rm{d}}}\left(0\right)\exp \left(-\displaystyle \frac{{C}_{0}\sqrt{T}}{m}t\right),\end{eqnarray}$ where ${V}_{{\rm{d}}}\left(0\right)$ is given by equation (20 ).
As we have shown in section 4, for an NESS with an arbitrary ${V}_{{\rm{d}}}\left(t\right),$ its entropy equals to the corresponding equilibrium system (${V}_{{\rm{d}}}=0$ ). Therefore, we always have$ \begin{eqnarray}\displaystyle \frac{{\rm{d}}S\left(t\right)}{{\rm{d}}t}=\displaystyle \frac{{\rm{d}}{S}_{{\rm{exch}}}\left(t\right)}{{\rm{d}}t}+\displaystyle \frac{{\rm{d}}{S}_{\mathrm{int}}\left(t\right)}{{\rm{d}}t}=0,\end{eqnarray}$ at any time t during the decaying process, where ${\rm{d}}{S}_{{\rm{exch}}}\left(t\right)$ is the entropy change flowing from the system to the thermostat, ${\rm{d}}{S}_{\mathrm{int}}\left(t\right)$ is the internal entropy change inside the system. According to equation (16 ), the dissipation heat from the NESS system to the thermostat in an infinitesimal time interval ${\rm{d}}t$ is$ \begin{eqnarray}{\rm{d}}Q\left(t\right)=-w\left(t\right){\rm{d}}t=-N{C}_{0}\sqrt{T}{V}_{{\rm{d}}}^{2}\left(t\right){\rm{d}}t.\end{eqnarray}$ On the other hand, we have$ \begin{eqnarray}{\rm{d}}Q\left(t\right)=T{\rm{d}}{S}_{{\rm{exch}}}\left(t\right),\end{eqnarray}$ so$ \begin{eqnarray}\displaystyle \frac{{\rm{d}}{S}_{\mathrm{int}}\left(t\right)}{{\rm{d}}t}=-\displaystyle \frac{{\rm{d}}{S}_{{\rm{exch}}}\left(t\right)}{{\rm{d}}t}=\displaystyle \frac{N{C}_{0}}{\sqrt{T}}{V}_{{\rm{d}}}^{2}\left(t\right).\end{eqnarray}$ Substituting equation (30 ) into the above equation, we obtain the entropy production$ \begin{eqnarray}{P}_{S}\left(t\right)\equiv \displaystyle \frac{{\rm{d}}{S}_{\mathrm{int}}\left(t\right)}{{\rm{d}}t}=\displaystyle \frac{N{C}_{0}}{\sqrt{T}}{V}_{{\rm{d}}}^{2}\left(0\right)\exp \left(-2\displaystyle \frac{{C}_{0}\sqrt{T}}{m}t\right)\geqslant 0.\end{eqnarray}$ The change of the entropy production with time is thus$ \begin{eqnarray}\displaystyle \frac{{\rm{d}}{P}_{S}\left(t\right)}{{\rm{d}}t}=-2N\displaystyle \frac{{C}_{0}^{2}}{m}{V}_{{\rm{d}}}^{2}\left(0\right)\exp \left(-2\displaystyle \frac{{C}_{0}\sqrt{T}}{m}t\right)\leqslant 0.\end{eqnarray}$
The above two expressions indicate that the decaying process has its entropy production monotonically decreases and approaches its minimum value of 0 at $t\to \infty $ when the system approaches equilibrium, which agrees with the minimum entropy production theorem by Prigogine [17].
8. Conclusions
In this work, solely based on the ergodicity and equiprobability principles, we derive the statistical mechanics of an NESS system coupled to a thermostat with temperature T driven by a constant force F . A random-walk analysis determines that the momentum space distribution differs from its equilibrium counterpart only by a shift of the drift velocity Vd in the distribution of the velocity along the force direction. The position space distribution is determined in the same way as for an equilibrium system in the canonical ensemble and is found to have exactly the same expression as the equilibrium counterpart. The entropy of the NESS system is found to be exactly the same as the corresponding equilibrium system, which means that the drift energy equals to the free energy difference and can all be used to do external work. The relation among F, T, and Vd is established by considering the balance between the input power and dissipation power: Vd is proportional to F and T−1/2 . The pressure towards the force is larger than that along the directions perpendicular to the force, and the pressure along the force is smaller. Finally, the relaxation towards equilibrium is proved to be an exponential decay obeying the minimum entropy production theorem.
Although the studied NESS system is the simplest nonequilibrium thermodynamic system, it is remarkable that its statistical mechanics features far from equilibrium can be quantified without any approximations and only based on the ergodicity and equiprobability principles. The results of this work not only allow the development of a thermostat algorithm for molecular dynamics simulations of NESS systems, but also serve as a starting point for further investigations of more complex nonequilibrium systems more common in nature than equilibrium systems, such as biological and social systems.
Acknowledgments
This work was funded by the Strategic Priority Research Program of Chinese Academy of Sciences (Grant No. XDA17010504) and the National Natural Science Foundation of China (Nos. 11774357, 11947302).
,1,21 CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of Sciences, 55 East Zhongguancun Road, PO Box 2735, Beijing, 100190 
New window|Download| PPT slide
