Local-Equilibrium Approximation in Non-Equilibrium Thermodynamics of Diffusion

1. Introduction

The concept of local equilibrium in non-equilibrium thermodynamics comes up in analyses of transport processes such as heat- and mass transport. In so-called linear non-equilibrium thermodynamics, one assumes local equilibrium, meaning that the local version of the Gibbs-Duhem equation is assumed valid. The assumption is closely related to assuming that the local non-equilibrium entropy can be approximated by its equilibrium value in the entropy balance, which is the basis for finding the entropy production and coupled flux-force relations [1]. Whereas entropy is well understood and modeled for thermodynamic systems at equilibrium, the situation for systems out of equilibrium is less clear [2]. Usually, the assumption is introduced as an approximation and justified with successful consistency tests as "circumstantial evidence". This is not satisfactory, and the extremely large forces and fluxes used in computer simulations that are necessary for acceptable signal to noise ratios in the results, raise a need for a quantitative examination of the local-equilibrium approximation (LEA) [3].

An alternative approach is to assume the opposite, viz. local non-equilibrium, derive the corresponding theory whenever possible, and determine the deviation from local equilibrium. A powerful criterion for this deviation is the difference in equilibrium and non-equilibrium entropy. Several theories for the non-equilibrium entropy have been proposed for heat transport and diffusion, see e.g. [4,5], and references therein. Cimmelli et al. [4] have given a very thorough discussion of four different theories and compared their pros and cons. They discuss the theories’ relation to kinetic theory, valid for low-density gases and conclude that in general, continuum theories that are compatible with kinetic theory are less universal (e.g. limited to large Knudsen numbers, monatomic gas), but more predictable than other theories.

In the present work, we pursue developing a theory for transport processes that do not depend on the LEA. The theory is based on the Boltzmann equation and developed for diffusion in a binary gas mixture. It is then used to analyze results from molecular dynamics simulations of binary diffusion in a Lennard-Jones/spline (LJs) system. The two components are identical (a "color experiment"), which simplifies the analysis and enables comparison with self-diffusion data. This is one of the simplest cases one can imagine for a study of a transport process. Our intention is to force the diffusion process as much as the system allows in an attempt to provoke a deviation from local equilibrium. The results are compared with what we get from Fick’s law and with data from equilibrium simulations (mean square displacement).

2. Theory

2.1. Fick’s Law

In most common applications, diffusion is usually expressed in the form of Fick’s law, which states that the diffusion flux ${\mathbf{J}}_i$ of component i is proportional to the gradient in the concentration $n_i$ of the corresponding component in the mixture

$${\mathbf{J}}_i=-D\nabla n_i$$

(1)

where D is the diffusion coefficient1. Fick’s law for gas mixtures can be obtained ab initio from kinetic theory in the short mean free path limit, meaning that it holds for descriptions on length scales much longer than the mean free path. The divergence of Eq. (1) at constant temperature and mass density, combined with the mass conservation law

$${\partial }_t{n}_i+\nabla ·{\mathbf{J}}_i=0$$

(2)

yields the well-known diffusion equation, under the additional approximation that D is independent of $n_i$,

$${\partial }_t{n}_i=D{\nabla }^2{n}_i.$$

(3)

Fick’s law can also be obtained using the non-equilibrium thermodynamics methodology. Considering a binary mixture at isothermal conditions, the entropy differential is

$$\mathrm{d}s=-{\displaystyle \frac{{\mu }_1}{T}}\mathrm{d}{x}_1-{\displaystyle \frac{{\mu }_2}{T}}\mathrm{d}{x}_2=-{\displaystyle \frac{\mu }{T}}\mathrm{d}{x}_1$$

(4)

where s is the entropy per particle, T the temperature, and ${\mu }_i$ and $x_i$ are the chemical potential and mole fraction of component i, respectively. We have used that $\mathrm{d}{x}_1=-\mathrm{d}{x}_2$ and introduced $\mu ={\mu }_1-{\mu }_2$. In order to keep the analysis relatively simple, we will in the following neglect the Dufour effect in the special case we consider here, in other words neglect the heat flux accompanying diffusion. Using the mass conservation law, we obtain the entropy balance

$$n{\partial }_{t}s={\displaystyle \frac{\mu }{T}}\nabla ·{\mathbf{J}}_1=\nabla ·{\displaystyle \frac{\mu {\mathbf{J}}_1}{T}}-{\mathbf{J}}_1·\nabla {\displaystyle \frac{\mu }{T}}=-\nabla ·{\mathbf{J}}_s+\sigma$$

(5)

where we have identified the entropy flux ${\mathbf{J}}_s=-\mu {\mathbf{J}}_1/T$ and the entropy production rate $\sigma ={\mathbf{J}}_1·\left(-\nabla \mu /T\right)$. The latter is expressed on the standard form as a product of the flux ${\mathbf{J}}_1$ and its conjugate driving force $-\nabla \mu /T$. Assuming that this driving force is proportional to ${\mathbf{J}}_1$ by some proportionality factor $r>0$ guarantees the positivity of the entropy production rate, and gives the following relation at isothermal conditions:

$$r{\mathbf{J}}_1=-\nabla {\displaystyle \frac{\mu }{T}}\Rightarrow {\mathbf{J}}_1=-{\displaystyle \frac{1}{rT}}{\displaystyle \frac{\partial \mu }{\partial n_1}}\nabla n_1=-D\nabla n_1\Rightarrow r={\displaystyle \frac{\partial \mu /\partial n_1}{TD}}.$$

(6)

Thus, we arrive at Fick’s law by the general requirement that the flux is a linear combination of the thermodynamic driving forces in the system, where in this case $-\nabla \mu /T$ is the sole driving force. The identification of D from r clarifies the connection.

2.2. Generalized Diffusion Equation

While Fick’s law is adequate for describing diffusion processes when the system is everywhere in local equilibrium, significant departures from local equilibrium expose its unphysical features. For example, according to Fick’s law, a local change in the composition or flux anywhere in the system will immediately affect the evolution of the physical state everywhere else in the system. That is, Fick’s law propagates information across the spatial domain at an unphysical infinite speed, when it can easily be argued that this propagation rate must have a speed limit dictated by the molecular velocities involved (the "causality problem" [6]). We seek a more general evolution equation of diffusion from the Boltzmann equation, which is not constrained by assumptions of local equilibrium. Following the work of Snell, Aranow and Spangler, the general equation of motion for component i in an isothermal multicomponent mixture, as obtained from the Boltzmann equation, is [7]

$${\partial }_t{\rho }_i{\mathbf{u}}_i=\nabla ·\left[{\mathbf{A}}_i-{\mathbf{B}}_i-{\rho }_i\left({\mathbf{u}}_i\otimes \mathbf{u}+\mathbf{u}\otimes {\mathbf{u}}_i-\mathbf{u}\otimes \mathbf{u}\right)\right]+x_i\sum _k{\xi }_{ik}{x}_k\left({\mathbf{u}}_k-{\mathbf{u}}_i\right)-n_i\nabla {\mu }_i$$

(7)

where ${\rho }_i$ and ${\mathbf{u}}_i$ are the mass density and the mean velocity of component i, respectively, $\mathbf{u}$ is the local barycentric velocity of the mixture, and ${\xi }_{ik}$ is a friction coefficient describing the exchange of momentum between components i and k. The tensor ${\mathbf{A}}_i$ denotes viscous contributions,

$${\mathbf{A}}_i=x_i\sum _k\left[\left({\zeta }_{ik}-{\displaystyle \frac{2}{3}}{\eta }_{ik}\right)\mathbf{I}\nabla ·x_k{\mathbf{u}}_k+{\eta }_{ik}\left(\nabla x_k{\mathbf{u}}_k+{\left(\nabla x_k{\mathbf{u}}_k\right)}^{\mathrm{T}}\right)\right],$$

(8)

where ${\zeta }_{ik}$ and ${\eta }_{ik}$ are coefficients related to the bulk and shear viscosity, respectively. We refer the interested reader to [7] for their definitions. We do not go into those here, as they are not of consequence for the analysis in this work. The tensor ${\mathbf{B}}_i$ represents the traceless part of the diffusion momentum outer product (indicated by the symbol $\overline{\otimes }$):

$${\mathbf{B}}_i={\rho }_i\overline{\left({\mathbf{u}}_i-\mathbf{u}\right)\otimes \left({\mathbf{u}}_i-\mathbf{u}\right)}={\rho }_i\left[\left({\mathbf{u}}_i-\mathbf{u}\right)\otimes \left({\mathbf{u}}_i-\mathbf{u}\right)-{\displaystyle \frac{1}{3}}\mathbf{I}\mathrm{Tr}\left[\left({\mathbf{u}}_i-\mathbf{u}\right)\otimes \left({\mathbf{u}}_i-\mathbf{u}\right)\right]\right].$$

(9)

Consider a mixture with no barycentric flow, $\mathbf{u}=\mathbf{0}$. The equation of motion simplifies somewhat:

$${\partial }_t{\rho }_i{\mathbf{u}}_i+\nabla ·\left({\rho }_i\overline{{\mathbf{u}}_i\otimes {\mathbf{u}}_i}\right)=x_i\sum _k{\xi }_{ik}{x}_k\left({\mathbf{u}}_k-{\mathbf{u}}_i\right)-n_i\nabla {\mu }_i+\nabla ·{\mathbf{A}}_i.$$

(10)

We identify the fluxes ${\mathbf{J}}_i=n_i\left({\mathbf{u}}_i-\mathbf{u}\right)=n_i{\mathbf{u}}_i$ and consider a binary mixture. The difference between the equations of motion for the two components is

$$\begin{array}{cc}\hfill & {\partial }_t\left(m_1{\mathbf{J}}_1-m_2{\mathbf{J}}_2\right)+\nabla ·\left({\displaystyle \frac{m_1\overline{{\mathbf{J}}_1\otimes {\mathbf{J}}_1}}{n_1}}-{\displaystyle \frac{m_2\overline{{\mathbf{J}}_2\otimes {\mathbf{J}}_2}}{n_2}}\right)\hfill \\ \hfill & =2x_1{x}_2\xi \left({\displaystyle \frac{{\mathbf{J}}_2}{n_2}}-{\displaystyle \frac{{\mathbf{J}}_1}{n_1}}\right)-n_1\nabla {\mu }_1+n_2\nabla {\mu }_2+\nabla ·\left({\mathbf{A}}_1-{\mathbf{A}}_2\right).\hfill \end{array}$$

(11)

where $m_i$ is the molecular mass of component i and the friction coefficients are symmetric, ${\xi }_{ik}={\xi }_{ki}=\xi$. We now use that ${\mathbf{J}}_2=-{\mathbf{J}}_1$ and use ${\mu }_i={\mu }_i^{\circ }+\mathsf{\Gamma }{k}_{\mathrm{B}}Tln{x}_i$ where $\mathsf{\Gamma }$ is the thermodynamic factor. Eq. (11) can then be reduced to

$$\overline{m}{\partial }_t{\mathbf{J}}_1+{\displaystyle \frac{\overline{m}}{n}}\nabla ·\left[\overline{{\mathbf{J}}_1\otimes {\mathbf{J}}_1}\left({\displaystyle \frac{c_1}{x_1}}-{\displaystyle \frac{1-c_1}{1-x_1}}\right)\right]-{\displaystyle \frac{1}{2}}\nabla ·\left({\mathbf{A}}_1-{\mathbf{A}}_2\right)=-{\displaystyle \frac{\xi {\mathbf{J}}_1}{n}}-\mathsf{\Gamma }n{k}_{\mathrm{B}}T\nabla x_1$$

(12)

where $\overline{m}=(m_1+m_2)/2$ is the average molecular mass, and $c_1=m_1/(m_1+m_2)$. In the system we consider here, $m_1=m_2=m$, so that $c=1/2$, which is used in the following. We see that under conditions where the terms on the left-hand side of Eq. (12) are small, that is when temporal and spatial variations in the flux are small, we reproduce Fick’s law. This allows us to identify the friction coefficient $\xi$ in terms of the diffusion coefficient:

$${\mathbf{J}}_1=-{\displaystyle \frac{\mathsf{\Gamma }{n}^2{k}_{\mathrm{B}}T}{\xi }}\nabla x_1=-nD\nabla x_1\Rightarrow \xi ={\displaystyle \frac{\mathsf{\Gamma }n{k}_{\mathrm{B}}T}{D}}.$$

(13)

Using the product rule, we can split the term involving the divergence of the traceless kinetic energy tensor:

$$\begin{array}{cc}\hfill & \nabla ·\left[\overline{{\mathbf{J}}_1\otimes {\mathbf{J}}_1}\left({\displaystyle \frac{1}{x_1}}-{\displaystyle \frac{1}{1-x_1}}\right)\right]\hfill \\ \hfill & =\left({\displaystyle \frac{1}{x_1}}-{\displaystyle \frac{1}{1-x_1}}\right)\nabla ·\overline{{\mathbf{J}}_1\otimes {\mathbf{J}}_1}-\left({\displaystyle \frac{1}{x_1^2}}+{\displaystyle \frac{1}{{(1-x_1)}^2}}\right)\overline{{\mathbf{J}}_1\otimes {\mathbf{J}}_1}·\nabla x_1.\hfill \end{array}$$

(14)

Furthermore, when the two components are physically identical, the viscosity coefficients ${\zeta }_{ij}$ and ${\eta }_{ij}$ in Eq. (8) are equal for both components, so that all contributions to ${\mathbf{A}}_i$ are proportional to $\nabla \mathbf{u}$, which vanishes in the absence of barycentric motion. Thus, ${\mathbf{A}}_i=\mathbf{0}$ in this case. Inserting Eq. (14) into Eq. (12) and removing viscous terms, we find

$$\begin{array}{cc}\hfill & \tau \left[{\partial }_t{\mathbf{J}}_1+{\displaystyle \frac{1}{2n}}\left({\displaystyle \frac{1}{x_1}}-{\displaystyle \frac{1}{1-x_1}}\right)\nabla ·\overline{{\mathbf{J}}_1\otimes {\mathbf{J}}_1}\right]+{\mathbf{J}}_1\hfill \\ \hfill & =-D\left[1-{\displaystyle \frac{m}{2\mathsf{\Gamma }{n}^2{k}_{\mathrm{B}}T}}\left({\displaystyle \frac{1}{x_1^2}}+{\displaystyle \frac{1}{{(1-x_1)}^2}}\right)\overline{{\mathbf{J}}_1\otimes {\mathbf{J}}_1}·\right]\nabla n_1\hfill \end{array}$$

(15)

where we identified the characteristic time $\tau =mn/\xi =mD/\left(\mathsf{\Gamma }{k}_{\mathrm{B}}T\right)$. We note that when the terms involving the divergence of the traceless kinetic energy tensor can be neglected, this reduces to

$$\tau {\partial }_t{\mathbf{J}}_1+{\mathbf{J}}_1=-D\nabla n_1,$$

(16)

the divergence of which leads to the well-known telegrapher’s equation for $n_1$, which is a damped wave equation describing waves propagating with a characteristic speed $\sqrt{D/\tau }$ and attenuated over a characteristic time $\tau$. Eq. (1) corresponds to Eq. (16) in the limit $\tau \to 0$, which corresponds to a diverging wave speed that leads to the unphysically instantaneous propagation of information across the domain. The telegrapher’s equation takes into account the inertia of molecular motion that resists rapid changes to the diffusion flux, and its connection to non-equilibrium thermodynamics is thoroughly discussed in the work by Jou, Casas-Vázquez and Lebon [8]. There, the authors arrive at equations of the telegrapher type by means of the thirteen moment approximation, which is on the same order as Eq. (15) with respect to spatial and temporal correlations. We can thus expect that the relationship between $\tau$ and D at this level of approximation should be corrected for higher order fluxes in order to match real data, justifying an a priori unknown correction factor $\alpha$ such that the true relaxation time is $\alpha \tau$.

For sake of simplicity, we will in the following restrict our attention to self-diffusion along a single dimension in three-dimensional space. The case of two identical components (self-diffusion) is particularly simple because the viscous terms vanish. Introducing also the correction factor $\alpha$, Eq, (15) is replaced by the one-dimensional equation

$$\alpha \tau \left[{\partial }_t{J}_1+{\displaystyle \frac{1}{3n}}\left({\displaystyle \frac{1}{x_1}}-{\displaystyle \frac{1}{1-x_1}}\right)\nabla J_1^2\right]+J_1=-D\left[1-{\displaystyle \frac{\alpha m{J}_1^2}{3\mathsf{\Gamma }{n}^2{k}_{\mathrm{B}}T}}\left({\displaystyle \frac{1}{x_1^2}}+{\displaystyle \frac{1}{{(1-x_1)}^2}}\right)\right]\nabla n_1.$$

(17)

Rearranging Eq. (17) and using that $r{J}_1$ represents a generalized driving force (see Eq. (6)), we find

$$\begin{array}{cc}\hfill r{J}_1=& -\left[1-{\displaystyle \frac{\alpha m{J}_1^2}{3\mathsf{\Gamma }{n}^2{k}_{\mathrm{B}}T}}\left({\displaystyle \frac{1}{x_1^2}}+{\displaystyle \frac{1}{{(1-x_1)}^2}}\right)\right]{\displaystyle \frac{\nabla {\mu }_{1,\mathrm{eq}}}{T}}\hfill \\ \hfill -& {\displaystyle \frac{\alpha m}{n{x}_{1}T}}\left[{\partial }_t{J}_1+{\displaystyle \frac{1}{3n}}\left({\displaystyle \frac{1}{x_1}}-{\displaystyle \frac{1}{1-x_1}}\right)\nabla J_1^2\right].\hfill \end{array}$$

(18)

In Eq. (18), we have specified that the chemical potential ${\mu }_{1,\mathrm{eq}}$ is the equilibrium value, i.e. that obtained from the equation of state to distinguish it from the non-equilibrium value to be introduced below. We can now identify the right-hand side of Eq. (18) as a generalized driving force. The expected entropy production rate is therefore

$$\begin{array}{cc}\hfill \sigma =-\sum _i\{J_i& \left[1-{\displaystyle \frac{\alpha m{J}_i^2}{3\mathsf{\Gamma }{n}^2{k}_{\mathrm{B}}T}}\left[{\displaystyle \frac{1}{x_i^2}}+{\displaystyle \frac{1}{{(1-x_i)}^2}}\right]\right]\nabla {\displaystyle \frac{{\mu }_{i,\mathrm{eq}}}{T}}\hfill \\ \hfill +& J_i{\displaystyle \frac{\alpha m}{n{x}_{i}T}}\left[{\partial }_t{J}_i+{\displaystyle \frac{1}{3n}}\left({\displaystyle \frac{1}{x_i}}-{\displaystyle \frac{1}{1-x_i}}\right)\nabla J_i^2\right]\}.\hfill \end{array}$$

(19)

In order to connect the departure from Fickian diffusion to the deviation from local equilibrium, we consider a generalized entropy s, which reduces to the local equilibrium entropy $s_{\mathrm{eq}}$ when the system relaxes to local equilibrium. Taking the mole fractions and fluxes as independent variables, the differential is

$$\mathrm{d}s=-{\displaystyle \frac{{\mu }_1}{T}}\mathrm{d}{x}_1-{\displaystyle \frac{{\mu }_2}{T}}\mathrm{d}{x}_2+B_1\mathrm{d}{J}_1+B_2\mathrm{d}{J}_2$$

(20)

where ${\mu }_i$ is now a generalized chemical potential of component i. In order to have correspondence to the equilibrium entropy, we must have that ${lim}_{J_i\to 0}{\mu }_i={\mu }_{i,\mathrm{eq}}$. We can determine $B_i$ by assuming that the equation of motion arises from a generalized driving force that is proportional to $J_i$, again ensuring the positivity of the resulting entropy production rate. We note that the contribution from the time derivative of $J_i$ to $\sigma$ is

$$\sigma \propto -\sum _i{\displaystyle \frac{\alpha m{J}_i}{n{x}_{i}T}}{\partial }_t{J}_i,$$

(21)

which allows us to deduce that

$$n{\partial }_{t}s\propto -{\displaystyle \frac{\alpha m{J}_i}{n{x}_{i}T}}{\partial }_t{J}_i\Rightarrow \mathrm{d}s\propto -{\displaystyle \frac{\alpha m{J}_i}{n^2{x}_{i}T}}\mathrm{d}{J}_i\Rightarrow B_i=-{\displaystyle \frac{\alpha m{J}_i}{n^2{x}_{i}T}}.$$

(22)

Now, we use that $x_1+x_2=1$ to simplify Eq. (20) to

$$\mathrm{d}s=-{\displaystyle \frac{\mu }{T}}\mathrm{d}{x}_1+B_1\mathrm{d}{J}_1+B_2\mathrm{d}{J}_2$$

(23)

where $\mu ={\mu }_1-{\mu }_2$, equivalent to the definition in Eq. (4). We can then obtain the $J_i$-dependence of $\mu$ through symmetry of mixed derivatives

$$-{\displaystyle \frac{\partial \mu /T}{\partial J_i}}={\displaystyle \frac{{\partial }^{2}s}{\partial J_i\partial x_1}}={\displaystyle \frac{{\partial }^{2}s}{\partial x_1\partial J_i}}={\displaystyle \frac{\partial B_i}{\partial x_1}}=-{\displaystyle \frac{\alpha }{n^{2}T}}{\displaystyle \frac{\partial }{\partial x_1}}{\displaystyle \frac{m{J}_i}{x_i}}=\left\{\begin{array}{cc}{\displaystyle \frac{\alpha m{J}_1}{n^2{x}_1^{2}T}}\hfill & i=1\hfill \\ -{\displaystyle \frac{\alpha m{J}_2}{n^2{x}_2^{2}T}}\hfill & i=2\hfill \end{array}\right.$$

(24)

Combined with the correspondence requirement to $s_{\mathrm{eq}}$ and that $J_1+J_2=0$ we find

$$-{\displaystyle \frac{\mu }{T}}=-{\displaystyle \frac{{\mu }_{\mathrm{eq}}}{T}}+{\int }_0^{J_1}\mathrm{d}{J}^{\prime }{\displaystyle \frac{\alpha m{J}^{\prime }}{n^{2}T}}\left({\displaystyle \frac{1}{x_1^2}}+{\displaystyle \frac{1}{{(1-x_1)}^2}}\right)=-{\displaystyle \frac{{\mu }_{\mathrm{eq}}}{T}}+{\displaystyle \frac{\alpha m{J}_1^2}{2n^{2}T}}\left({\displaystyle \frac{1}{x_1^2}}+{\displaystyle \frac{1}{{(1-x_1)}^2}}\right)$$

(25)

We identify the second term on the right-hand side as the kinetic energy associated with species interdiffusion, multiplied by a factor $\alpha /T$. This implies that in the non-equilibrium situation, the equilibrium chemical potential is replaced by a generalized chemical potential where a term proportional to the kinetic energy of diffusion is subtracted. This agrees with the assessment by de Groot and Mazur, who arrived at a similar generalized driving force by considering the necessity of subtracting the kinetic energy of diffusion from the typically defined internal energy to obtain the truly internal energy [9]. We may quantify the departure from local equilibrium using the deviation of the local chemical potential from its equilibrium value:

$${\mu }_{\mathrm{eq}}-\mu ={\displaystyle \frac{\alpha m{J}_1^2}{2n^2}}\left({\displaystyle \frac{1}{x_1^2}}+{\displaystyle \frac{1}{{(1-x_1)}^2}}\right)=\alpha \left({\displaystyle \frac{m{u}_1^2}{2}}+{\displaystyle \frac{m{u}_2^2}{2}}\right),$$

(26)

which reproduces the de Groot and Mazur result in the special case $\alpha =1$. We restrict our attention to the steady-state equation of motion, where ${\partial }_t{J}_1=0$ and $\nabla J_1\propto {\partial }_t{x}_1=0$. This now reads

$$J_1=-D\left[1-{\displaystyle \frac{{\mu }_{\mathrm{eq}}-\mu }{\mathsf{\Gamma }{E}_k}}\right]\nabla n_1$$

(27)

where $E_k=3k_{\mathrm{B}}T/2$ is the average molecular translational kinetic energy. Fick’s law with constant D predicts a constant composition gradient in steady state, according to Eq. (3). According to Eq. (27), we expect to see a curved composition profile when the local deviation in chemical potential from the local equilibrium value is significant relative to $\mathsf{\Gamma }{E}_k$. The mechanism at play here is that of the diffusion process approaching its speed limit dictated by the local distribution of molecular velocities. Whereas a Fickian model would violate this speed limit by allowing the diffusion flux to increase beyond the rate of transport permissible by the ballistic motion of molecules at a given temperature, this generalized formulation properly enforces said speed limit in a continuous and physically consistent manner. Investigating the curvature of the steady-state composition profile resulting from a large imposed diffusion flux allows us to properly quantify the deviation from local equilibrium.

We note that for fixed $x_i$, we can use the result for $B_i$, Eq. (22), in Eq. (20) and integrate to give

$$s-s_{\mathrm{eq}}=-\sum _i{\displaystyle \frac{\alpha m{J}_i^2}{2n^2{x}_{i}T}}=-{\displaystyle \frac{\alpha m{J}_1^2}{2n^{2}T}}\left({\displaystyle \frac{1}{x_1}}+{\displaystyle \frac{1}{1-x_1}}\right)=-{\displaystyle \frac{\alpha m{J}_1^2}{2n^2{x}_1\left(1-x_1\right)T}},$$

(28)

which shows that the non-equilibrium entropy s is always smaller than the equilibrium value, which is consistent with the second law. The absolute deviation in the entropy from the local equilibrium value has a minimum at $x_1=0.5$.

3. MD Simulations

Non-equilibrium molecular dynamics (NEMD) simulations were carried out to provide data for use with the theory. The layout of the MD system is shown in Figure 1. The system consisted of $N=32,768$ Lennard-Jones/spline (LJs) particles [10] in a rectangular cuboid with aspect ratio $L_x:L_y:L_z=2:1:1$. The LJs potential is a smoothly truncated Lennard-Jones potential with cut-off distance $r_c\approx 1.7\sigma$. The MD box had periodic boundary conditions in all three directions. The box was divided into 64 layers of equal thickness in x-direction for recording local data of temperature, density, composition, etc.. The particles were of two types, 1 ("red") and 2 ("blue") with exactly the same mass and potential parameters.

Concentration gradients and diffusive mass fluxes were generated with an in-house NEMD code using the MEX algorithm [3]. The MEX algorithm works so that a randomly picked particle of type 2 is changed to type 1 in one of the control volumes at the box ends (see Figure 1) while a particle of type 1 is changed to type 2 in the central volume at the same time. The overall concentration was thus kept constant. The identity swapping occurred at regular intervals. By changing the swapping frequency, the diffusion flux was controlled from zero (equilibrium) to maximum (saturated control volumes, i.e. $x_1=1$ at the box ends and $x_1=0$ in the box center). If there was no candidate particles to swap in either control volume, the swapping was skipped.

The simulations started from a FCC lattice configuration with $N_1=N_2=N/2$. To prevent heating of the system due to the entropy production, the system temperature was controlled by thermostating each of the 64 layers to the same set temperature and the overall density was controlled by the total number of particles in the fixed box volume. Six cases were generated with the density and temperature as listed in Table 1. All values of physical properties reported here are in Lennard-Jones units. The mean free path, $\lambda$, was estimated from elementary kinetic theory as $\lambda =1/\left(\pi \sqrt{2}n\right)$ and compared with half box length, $L_x/2$. Case 1 is a very dilute (ideal) gas which is not well suited for MD simulations due to the fact that the particles are mostly in free flight, the mean free path is of the same order as $L_x/2$ (the distance between the swapping layers) in this case. Nevertheless, this low-density case was included here because of the kinetic-theory basis for the analysis. Case 6 is close to the LJs system’s dew point (at $n=0.026$ for $T=0.7$). When the mean free path is much larger than the range of the potential, particles may overlap at the end of a time step. To avoid strong accelerations if this happens, the particles were given a reflecting "shield" at $r_{ik}=0.8$.

A total of 4,000,000 time steps of length $\delta t=0.002$ in Lennard-Jones units were used to produce data for analysis. The final 2,000,000 time steps were used to represent steady-state with $J=0$.

4. Results and Discussion

According to Darken’s second equation [11], the mutual diffusion coefficient D in a binary mixture relates to the self-diffusion coefficients $D_1$ and $D_2$, as

$$D=(x_1{D}_1+x_2{D}_2)\mathsf{\Gamma }.$$

(29)

In our ideal case, $\mathsf{\Gamma }=1$ and $D_1=D_2$, so that the diffusion coefficient we compute by generating a diffusive flux in the two-component system is simply the self-diffusion coefficient. This enables a comparison between the binary diffusion coefficient, generated with extreme fluxes and forces, and self-diffusion coefficients generated at global equilibrium.

We shall now analyze the binary diffusion data using the theory developed in Section 2.2. Central in this analysis is the curvature of the mole-fraction profile and the deviation from local equilibrium. A typical mole fraction profile is shown in Figure 2. Three important features are: (i) the profile appears to be linear for $x/L_x\approx 0.25$ and $x/L_x\approx 0.75$, (ii) the gradient $\nabla x_1$ is slightly curved for $0.04<x/L_x<0.46$ and $0.54<x/L_x<0.96$, and (iii) there are significant jumps in $x_1$ between the swapping volumes and the bulk.

The jumps in $x_1$ are consequences of the large imposed mass flux. Similar jumps are seen for temperature profiles in simulations of heat tansport [12] and can be reduced by reducing the flux. The curvature in $x_1$ is very weak, meaning that the fraction $({\mu }_{\mathrm{eq}}-\mu )/\mathsf{\Gamma }{E}_k$ is small compared to 1 (D in Eq.(27) must be constant).

A central question in this work is to what extent the LEA is valid. Eq. (26) gives a direct measure of the deviation from the equilibrium chemical potential. We see that ${\mu }_{\mathrm{eq}}-\mu$ is positive. From Eq. (27) we see that Fick’s law will always underestimate the value of D. The difference can be quantified with the second term in the square bracket in Eq. (27). If the fraction $({\mu }_{\mathrm{eq}}-\mu )/\mathsf{\Gamma }{E}_k$ is small compared to 1, local equilibrium is a good approximation, reducing Eq. (27) to Fick’s law. Otherwise, the approximation is poor. Exactly what we mean by an acceptable approximation is a matter of choice, but the important point is that whatever choice is made, it can be quantified and communicated. Eq. (26) shows that the deviation from equilibrium depends on $x_1$ with the minimum value at $x_1=0.5$. The minimum value of $({\mu }_{\mathrm{eq}}-\mu )/\mathsf{\Gamma }{E}_k$ is then $8\alpha m{J}^2/\left(3n^{2}kT\right)$ with $\mathsf{\Gamma }=1$.

Lack of local equilibrium will depend on density, temperature, and the mass flux. Figure 3 shows the deviations from local equilibrium for the two extreme densities in this study ($n=0.001$ and $0.02$), $T=0.7$, and maximum imposed mass flux (saturated control volumes). The LEA appears to be good for the higher density for a wide range of mole fractions (approximately 2 % deviation from equilibrium) and very poor for the lower density (more than 25 % deviation). For $n=0.001$, molecular collisions are rare and mass transport is to a large extent ballistic. The divergence of $({\mu }_{\mathrm{eq}}-\mu )/\mathsf{\Gamma }{E}_k$ can be understood in terms of Eq. (26) considering the diffusive velocities of each component, which is $u_i=J_i/n_i$. For a constant flux $J_i$, this velocity diverges as $x_1\to 0$ or $x_1\to 1$.

Figure 4 shows the deviations from local equilibrium for three temperatures ($T=0.7,1.0$, and $2.0$), $n=0.01$, and maximum mass flux. The LEA is less than 4 % off for equimolar mixtures for the highest temperature and better for the lower temperatures.

The LEA will clearly be better for reduced mass flux. We simulated Case 4 for 6 additional imposed fluxes. The results are shown i Figure 5. According to Eqs. (26) and (27), the $({\mu }_{\mathrm{eq}}-\mu )/\mathsf{\Gamma }{E}_k$ is proportional with $J_1^2$, which is illustrated in Figure 5.

The parameter $\alpha$ in Eq. (27) was determined by requiring that the diffusion coefficient D was independent of $x_1$ and $\nabla x_1$. The curvature in $x_1$ shown in Figure 2, which balances the two components’ mean velocities in Eq. (26), was found to be small, meaning that the deviation from local equilibrium was small. The determined value of $\alpha$ is therefore uncertain, but sufficient to quantify the deviation from local equilibrium.

Numerical values for the diffusion coefficient obtained from Eq. (27) for the six cases are shown in Table 2. The coefficient obtained from Fick’s law is consistently smaller than D and the agreement between $D_{\mathrm{Fick}}$ and D is better for higher densities with fixed temperature and for higher temperature with fixed density. Both trends can be understood in terms of the system’s ability to dissipate the effect of the perturbation.

As mentioned in the Introduction, the main purpose of this work is to develop a tool for analysis of the LEA in non-equilibrium thermodynamics. We have considered a very special case, mutual diffusion of an ideal system in gas phase, and used kinetic theory as basis. As it stands, the method has limited value for higher densities, which requires an improved version of the kinetic theory. Another limitation is that we have approximated $\mathsf{\Gamma }$, ${\mu }_i$, and $s_{\mathrm{eq}}$ with ideal-gas values. These properties must be determined from a more realistic equation of state, which is available for the LJs model [13].

The theory was developed for binary diffusion at low density in general and simplified for a system with identical components. Application to binary systems with different components is therefore straightforward. Another interesting application is heat transport, where the same considerations of the combination of ballistic and diffusive transport mechanisms apply.

Eq. (28) gives us a tool to compute the difference between the non-equilibrium and equilibrium entropy. We found for $(s-s_{\mathrm{eq}})/s_{\mathrm{eq}}$:

$${\displaystyle \frac{s-s_{\mathrm{eq}}}{s_{\mathrm{eq}}}}=-{\displaystyle \frac{\alpha m}{2s_{\mathrm{ST}}{n}^2{x}_1(1-x_1)T}}{J}_1^2$$

(30)

where we have used the Sackur-Tetrode ideal gas entropy, $s_{\mathrm{ST}}$, as representing $s_{\mathrm{eq}}$. Using Case 2 as an example and assuming that the LJs model represents argon with the potential parameters, $ϵ/k_{\mathrm{B}}=120$ K and $\sigma =0.34$ nm, we find that the molar equilibrium entropy is 125 J ${\mathrm{mol}}^{-1}{\mathrm{K}}^{-1}$ and the difference between the molar non-equilibrium and equilibrium entropies is 0.5 J ${\mathrm{mol}}^{-1}{\mathrm{K}}^{-1}$ for $x_1=0.5$. The relative deviation is $(s-s_{\mathrm{eq}})/s_{\mathrm{eq}}=0.004$ for $x_1=0.5$ and 0.01 for $x_1=0.9$ (the extreme value in the bulk, cf Figure 2). These results are additional demonstrations of the excellence of the LEA in this diffusion case.

The importance of the kinetic energy of diffusion [9] can be illustrated by the ratio between the kinetic energy of diffusion and the thermal energy,

$$\mathrm{ratio}={\displaystyle \frac{m}{3k_{\mathrm{B}}T}}{\left({\displaystyle \frac{J_1}{n_1}}\right)}^2\left({\displaystyle \frac{1}{x_1^2}}+{\displaystyle \frac{1}{{(1-x_1)}^2}}\right).$$

(31)

Again using Case 2 as an example, we found that the ratio is 0.04 and 0.46 for $x_1=0.5$ and $x_1=0.9$, respectively. This result shows that even if the kinetic energy of diffusion is large compared with the thermal energy, the deviation from local equilibrium is small and LEA is good. A more detailed analysis of the velocities should include the full velocity distribution.

5. Conclusions

Using a theoretical framework based on the Boltzmann equation, we have quantitatively assessed the local equilibrium approximation (LEA) commonly used in non-equilibrium thermodynamics. Our analysis, which applies to both equilibrium and non-equilibrium conditions, was developed for an ideal diffusion process in a binary system composed of two physically identical components. This approach yields a nonlinear generalization of the telegrapher’s equation. In particular, the steady-state equation we derive, Eq. (27), provides a novel method for determining a corrected relaxation time from steady-state measurements, which are generally easier to analyze than transient phenomena. We then employ the relaxation time correction factor, $\alpha$, to quantify the deviation from equilibrium by comparing the equilibrium entropy with its non-equilibrium counterpart, Eq. (28). This framework offers a quantitative measure of the deviation from local equilibrium and enhances the assessment of the LEA in non-equilibrium thermodynamics.

Our results indicate that the LEA holds well for most conditions except at the lowest density examined (0.001 in Lennard-Jones units). At low density, a large mean free path and a high diffusion kinetic energy relative to the system’s thermal energy lead to significant deviations from local equilibrium. In contrast, for higher gas-phase densities, the LEA remains a robust approximation for the system considered. Notably, even under conditions of large fluxes and forces where, for argon, the mass fluxes and concentration gradients reach approximately ${10}^4$ mol ${\mathrm{m}}^{-1}{\mathrm{s}}^{-1}$ and ${10}^{10}$ mol ${\mathrm{m}}^{-4}$, respectively (Case 2), the estimated diffusion coefficient for argon at about 1 bar and 85 K is on the order of ${10}^{-6}{\mathrm{m}}^2{\mathrm{s}}^{-1}$, which agrees well with experimental data [14].