Theory of cross phenomena and their coefficients beyond Onsager theorem

Abstract

Cross phenomena, representing responses of a system to external stimuli, are ubiquitous from quantum to macro scales. The Onsager theorem is often used to describe them, stating that the coefficient matrix of cross phenomena connecting the driving forces and the fluxes of internal processes is symmetric. Here we show that this matrix is intrinsically diagonal when the driving forces are chosen from the gradients of potentials that drive the fluxes of their respective conjugate molar quantities in the combined law of thermodynamics including the contributions from internal processes. Various cross phenomena are discussed in terms of the present theory.

IMPACT STATEMENT

Through flux equations based on the combined law of thermodynamics, theory of cross phenomena is developed and applied to thermoelectricity, thermodiffusion, diffusion, electromigration, electrocaloric and electromechanical effects, and thermal expansion.

KEYWORDS:

• Thermoelectric
• thermodiffusion
• electromigration
• electrocaloric
• electromechanical

1. Introduction

In a recent overview article [1], the author discussed fundamental thermodynamics, thermodynamic modeling, and the applications of computational thermodynamics. Some kinetic aspects were briefly mentioned at the end of the overview including diffusion and Seebeck coefficients and general discussion on off-diagonal transport coefficients in connection with the Onsager reciprocal relationships (ORRs) [2]. The present paper aims to expand those discussions and present a broader perspective on thermodynamic basis of kinetic coefficients in the framework of the materials science and engineering (MSE) discipline. The hierarchical relationships among MSE discipline’s fundamental components, including chemistry and thermodynamics, processing and kinetics, structure and crystallography, and defect and property of materials, can be represented by the off-central circles shown in Figure 1 proposed by the author [3] coinciding with the time when he coined the term ‘Materials Genome®' [4]. Figure 1(i) and (ii) represent the engineering and scientific views of the MSE discipline, respectively, and demonstrate the core importance of thermodynamics and kinetics in the discipline. The off-center feature of the circles implies the author’s perspectives on both the broadness and coherency of material research in terms of the distinctive features of individual components and their integrations in shaping the strength of the MSE discipline.

Figure 1. Fundamental components of materials science and engineering with (i) engineering focus and (ii) science focus [3].

Following the discussions in the overview article [1], the present article focuses on the fundamentals of cross phenomena, such as transports of mass and electron due to heat conduction, and the predictions of their coefficients in connection with the ORRs [2]. Instead of the commonly used phenomenological kinetic equations proposed by Onsager and broadly used in the community, the present discussion starts from the combined first and second law of thermodynamics applicable to nonequilibrium systems, beyond its equilibrium form established by Gibbs [5], by including the entropy productions due to irreversible internal kinetic processes inside the system. It will be articulated that under the linear proportionality approximation of kinetic processes, the ORRs are naturally fulfilled with the off-diagonal coefficients being zero. Consequently, it follows that the coefficient of a cross phenomenon is the product of the kinetic coefficient of the transport process and the derivative of the potential driving the transport process with respect to the potential that triggers the cross phenomenon. The derivatives between potentials play a central role in materials functionalities and can be predicted from free energy containing independent internal state variables. In the present paper, free energy is used to describe all energies except the internal energy represented by the combined first and second law of thermodynamics, which is sometimes simplified as combined law of thermodynamics. While the Gibbs energy [6] and Helmholtz energy [7] are used without ‘free', as suggested by the International Union of Pure and Applied Chemistry.

The present article is organized as follows. In Section 2, the combined first and second law of thermodynamics applicable to both equilibrium and non-equilibrium systems is discussed by retaining the entropy productions due to irreversible internal processes in a system. In Section 3, the flux of an internal process is introduced for the change of a molar quantity based on the combined law that it is driven solely by the gradient of its conjugate potential, demonstrating that cross phenomena manifest through the derivatives between potentials. Various cross phenomena studied in the literature are systematically discussed in Section 4 in terms of the present theory, including thermoelectricity, thermodiffusion, uphill diffusion, electromigration, electrocaloric and electromechanical effects, and thermal expansion. In Section 5, our zentropy theory built on statistical mechanics and quantum mechanics in terms of first-principles calculations based on density functional theory (DFT) [8,9] is reviewed including fundamentals and applications in predicting cross phenomena in relation to critical points and associated property anomaly through Maxwell relations, including superconductivity. Finally, a summary and outlooks are presented.

2. Combined first and second law of thermodynamics for nonequilibrium systems

Thermodynamics is a science concerning the state of a system, whether it is stable, metastable, or unstable, when interacting with its surroundings. The first law of thermodynamics specifies that the exchange of heat, work, and chemical energies between a system and its surrounding must be represented by the change of the internal energy of the system, i.e. the energy is conserved. It does not prescribe whether the system is in internal equilibrium or not, though equilibrium is usually assumed in most textbooks. Furthermore, in most textbooks, the chemical energy is introduced much later, causing significant confusion on the concept of chemical potential. The second law represents an inequality, i.e. any irreversible internal processes (ip) in a system must generate entropy, resulting in positive entropy production written as $d_{ip}S>0$. Consequently, the total entropy change of a system can be obtained as follows [1,10,11]: (1) $dS=\frac{dQ}{T}+\sum _i{S}_{i}d{N}_i+d_{ip}S$(1) where $dQ$ and $d{N}_i$ are the exchanges of heat and moles of component $i$ from the surroundings to the system (negative when from the system to the surroundings), $T$ is the temperature, and $S_i$ is the partial entropy of component i defined as: (2) $S_i={\left(\frac{\mathrm{\partial }S}{\mathrm{\partial }{N}_i}\right)}_{dQ=0,d_{ip}S=0,d{N}_{j\ne i}=0}$(2)

The general form of the combined law of thermodynamics can thus be written as follows aided by Equation (1) [1,10,11]: (3) $\begin{array}{rl}dU& =dQ+dW+\sum _i{U}_{i}d{N}_i\\ & =TdS+dW+\sum _i{\mu }_{i}d{N}_i-T{d}_{ip}S\\ & =\sum _a{Y}^{a}d{X}^a-T{d}_{ip}S\end{array}$(3) where $dU$ is the change of the internal energy in the system, $dW$ is the exchange of work from the surroundings to the system, including mechanical, electric, and magnetic work, $U_i$ and ${\mu }_i=U_i-T{S}_i$ are the partial internal energy and chemical potential of component $i$. $Y^a$ and $X^a$ represent the pairs of conjugate variables with $Y^a$ for potentials, such as temperature, stress/pressure with negative sign [10,11], chemical potential, electrical and magnetic fields, and $X^a$ for molar quantities, such as entropy, strain/volume, moles of components, and electrical and magnetic displacements. It is important to note that the change of the internal energy does not depend on internal processes because this change concerns only the exchanges between the surroundings and the system as shown by the first part of Equation (3), i.e. for an isolated system with $dQ=dW=d{N}_i=0$, the internal energy of the system is constant with $dU=0$ independent of whether the system is at equilibrium or not. This can also be understood because $dS$ in Equation (3) contains the contribution from $d_{ip}S$ as shown by Equation (1).

All kinetics processes in a nonequilibrium system are irreversible and contribute to the total entropy change as shown by $d_{ip}S$ in Equation (1). Under the linear proportionality approximation, the last term in Equation (3) can be written as: (4) $T{d}_{ip}S=\sum _j{D}_j{d}_j$(4) where $d{\xi }_j$ represents the $j^{th}$ internal process for the change of the internal variable ${\xi }_j$, $D_j$ is the driving force for the $j^{th}$ internal process, and the summation goes over all irreversible internal processes inside the system. In reality, internal processes do not obey linear proportionality, but in principle, one can always select a small enough $d{\xi }_j$ so the higher-order terms are less important in defining Equation (4) and perform integrations along the pathways to obtain overall behaviors. It is also important to differentiate externally controlled variables, i.e. all $X^a$ in Equation (3), vs internal variables, i.e. all ${\xi }_j$ in Equation (4), as emphasized by Hillert [10]. Furthermore, to study the stability of an internal process or a system such as instability and critical points discussed in Sections 4.10 and 5.3, one would need to include higher order terms beyond linear proportionality.

As an example, let us consider the diffusion of component $i$ as an internal process, i.e. $d\xi =d{c}_i$ with $c_i=N_i/V$ being the concentration of component $i$ in terms of moles per volume, and the driving force is the decrease of the chemical potential of the component, i.e. $D=-\mathrm{\Delta }{\mu }_i$. This will be denoted by $D=-\mathrm{\Delta }{Y}_{\xi }$ or $D_j=-\mathrm{\Delta }{Y}_{{\xi }_j}$ with $Y_{{\xi }_j}$ being the conjugate potential of ${\xi }_j$ in the rest of the manuscript. When an internal process consists of changes of several molar quantities simultaneously, such as chemical reactions involving several components, the entropy production of the internal process can be written as the sum of entropy production of the change of each molar quantity as follows: (5) $Dd\xi =-\sum _i\mathrm{\Delta }{\mu }_{i}d{c}_i$(5)

It should be noted that there are internal variables that are related to microstructure such as grain size [12] and phase morphologies [13], which will not be discussed further in the present article because it is not directly related to cross phenomena. Another significant challenge is to quantify the entropy in a system because it contains the contributions from internal processes as shown in Equation (1) and is multiscale in nature [1,14], and these contributions are also important in equilibrium systems due to thermal fluctuations at finite temperatures above zero K [15–18] and will be discussed in Section 5.

For equilibrium systems, there are no irreversible internal processes, i.e. $T{d}_{ip}S=0$, and one obtains the combined law derived by Gibbs [5] as follows: (6) $dU=\sum _a{Y}^{a}d{X}^a$(6) All $d{X}^a$ in Equation (6) refer to exchanges from the surroundings to the system or negative exchanges from the system to the surroundings. The internal energy of the system is determined by the amount of each $X^a$ which is controlled from the surroundings, i.e. $X^a$ are independent variables of $U(X^a)$. All $X^a$ are called natural variables of the internal energy as they are defined from the combined law of thermodynamics [10]. All quantities in the system are also dependent on all $X^a$, including all internal variables, i.e. ${\xi }_j(X^a)$. The values of ${\xi }_j(X^a)$ can be determined by solving the following equilibrium equations for all conceivable internal processes in the system except thermal fluctuations, which will be further discussed later in the present paper: (7) $D_{j}d{\xi }_j=0$(7)

As mentioned above, for nonequilibrium systems, $dS$ in Equation (3) contains the contributions from irreversible internal processes as shown by Equation (1) with $T{d}_{ip}S>0$, and ${\xi }_j$ thus become independent variables of the system. All properties of the system are the functions of ${\xi }_j$, including the internal energy, i.e. $U(X^a,{\xi }_j)$, so are all the potentials and driving forces of internal processes, which are the partial derivative of the internal energy to the molar quantity as follows: (8) $\begin{array}{rl}{Y}^a(X^a,X^b,{\xi }_j)& ={\left(\frac{\mathrm{\partial }U}{\mathrm{\partial }{X}^a}\right)}_{X^b\ne X^a,{\xi }_j}\end{array}$(8) (9) $\begin{array}{rl}{D}_j(X^a,X^b,{\xi }_j,{\xi }_k)& ={\left(\frac{\mathrm{\partial }U}{\mathrm{\partial }{\xi }_j}\right)}_{X^a,X^b,{\xi }_{k\ne j}}\end{array}$(9) where $X^b$ denotes all other molar quantities except $X^a$, and ${\xi }_k$ represents all independent internal variables in the system except ${\xi }_j$. Equation (8) represents a critically important relationship that the values of $Y^a$ are not only affected by its conjugate molar quantity $X^a$, but also by all other non-conjugate molar quantities $X^b$ and all independent internal variables ${\xi }_j$. Equation (9) seems in contradiction with the statement just made that the internal energy is independent of internal variables, but it is not because when the entropy of the system is kept constant, there must be either exchange of heat or mass between the surroundings and the system as stipulated by Equation (1) if there are internal processes in the system, i.e. $d_{ip}S>0$, which results in the change of the internal energy of the system. In the rest of this paper, the variables kept constant for partial derivatives are those remaining natural variables and independent internal variables and will not be listed for the sake of brevity of formulas unless necessary.

Furthermore, the system can be constrained from the surroundings through the control of one or more potentials, $Y^c$, so that it becomes the independent variables of the system and the natural variable of a new free energy defined as follows along with its combined law: (10) $\begin{array}{rl}\mathrm{\Phi }& =U-\sum Y^c{X}^c\end{array}$(10) (11) $\begin{array}{rl}d\mathrm{\Phi }& =\sum _{a\ne c}{Y}^{a}d{X}^a-\sum _{c\ne a}{X}^{c}d{Y}^c-T{d}_{ip}S\end{array}$(11)

Equations (8) and (9) thus become: (12) $\begin{array}{rl}{Y}^a(X^a,X^b,Y^c,{\xi }_j)& ={\left(\frac{\mathrm{\partial }\mathrm{\Phi }}{\mathrm{\partial }{X}^a}\right)}_{X^b\ne (X^a,X^c),Y^c,{\xi }_j}\end{array}$(12) (13) $\begin{array}{rl}{D}_j(X^a,X^b,Y^c,{\xi }_j,{\xi }_k)& ={\left(\frac{\mathrm{\partial }\mathrm{\Phi }}{\mathrm{\partial }{\xi }_j}\right)}_{(X^a,X^b)\ne X^c,Y^c,{\xi }_{k\ne j}}\end{array}$(13)

This represents the core concept of cross phenomena such as thermoelectrics and thermodiffusion with $Y^c$ being temperature and pressure and $\mathrm{\Phi }$ being the Gibbs energy and will be discussed in the rest of the paper. It is noted that in both thermoelectrics and thermodiffusion, the systems are inhomogeneous and can be divided into sufficiently smaller regions or subsystems that may be considered to be homogeneous to apply Equations (12) and (13) [19]. Various constraints from the surroundings were presented in a recent publication [20], including the microcanonical, canonical, grand canonical, isothermal–isobaric (Gibbs), isoenthalpic-isobaric, and partial grand Gibbs ensembles, plus the isoentropic-isobaric ensemble, which is very hard to realize experimentally due to the contributions to entropy from internal processes as shown by Equation (1).

Let us introduce the hydrostatic, mechanical, and electrical works for the interest of the present work as follows: (14) $\begin{array}{rl}d{W}_h& =-PdV\end{array}$(14) (15) $\begin{array}{rl}d{W}_m& =\sigma d\epsilon \end{array}$(15) (16) $\begin{array}{rl}d{W}_e& =Ed\theta \end{array}$(16) where $P$ and $V$ are pressure and volume, $\epsilon$ and $\sigma$ are the strain and stress components of mechanical energy, and $\theta$ and $E$ are the components of the electrical displacement and electrical field, both with their tensor indexes omitted for the sake of brevity of formulas. $\theta$ is used here because $D_j$ is designated for the driving force of internal process. The negative sign in Equation (14) is because negative volume change implies the increase of the system energy, while the positive sign in Equations (15) and (16) is because both a compressive strain and an electrical displacement in the system are denoted by positive $\epsilon$ and $\sigma$ increase the system energy [21,22], though other conventions are also used in the literature [23].

3. Kinetic equations

3.1. Entropy production and flux through conjugate variables

From the above discussion, the rate of the entropy production per volume due to the $j^{th}$ internal process in a sufficiently small region with a thickness of $\mathrm{\Delta }z$ and area of $A$ can be written as: (17) $\frac{T{\dot{d}}_{ip}{S}_j}{V}=\frac{D_j\dot{d}{\xi }_j}{A\mathrm{\Delta }z}=-\frac{\dot{d}{\xi }_j}{A}\frac{\mathrm{\Delta }{Y}_{{\xi }_j}}{\mathrm{\Delta }z}=-\frac{\dot{d}{\xi }_j}{A}\mathrm{\nabla }{Y}_{{\xi }_j}$(17) with $D_j=-\mathrm{\Delta }{Y}_{{\xi }_j}$ and $Y_{{\xi }_j}$ being the conjugate potential of the internal variable of ${\xi }_j$ as discussed above and $\mathrm{\nabla }{Y}_{{\xi }_j}$ the gradient of $Y_{{\xi }_j}$. The flux of ${\xi }_j$, $J_{{\xi }_j}$, can then be defined as follows: (18) $J_{{\xi }_j}=\frac{\dot{d}{\xi }_j}{A}=-L_{{\xi }_j}\mathrm{\nabla }{Y}_{{\xi }_j}$(18) where $L_{{\xi }_j}$ is the kinetic coefficient for the change of ${\xi }_j$. Equation (17) can be re-written as: (19) $\frac{T{\dot{d}}_{ip}{S}_j}{V}=L_{{\xi }_j}(\mathrm{\nabla }{Y}_{{\xi }_j}{)}^2$(19)

Equation (18) can be further expanded to independent variables of $\mathrm{\nabla }{Y}_{{\xi }_j}$ as follows using Equation (12): (20) $\begin{array}{rl}{J}_{{\xi }_j}& =-L_{{\xi }_j}\mathrm{\nabla }{Y}_{{\xi }_j}\\ & =-L_{{\xi }_j}(\frac{\mathrm{\partial }{Y}_{{\xi }_j}}{\mathrm{\partial }{\xi }_j}\mathrm{\nabla }{\xi }_j+\sum _{{\xi }_k\ne {\xi }_j}\frac{\mathrm{\partial }{Y}_{{\xi }_j}}{\mathrm{\partial }{Y}_{{\xi }_k}}\mathrm{\nabla }{Y}_{{\xi }_k}+\sum _{{\xi }_l\ne {\xi }_j,{\xi }_k}\frac{\mathrm{\partial }{Y}_{{\xi }_j}}{\mathrm{\partial }{X}_{{\xi }_l}}\mathrm{\nabla }{X}_{{\xi }_l}\\ & +\sum _{{\xi }_n\ne {\xi }_j,{\xi }_k,{\xi }_l}\frac{\mathrm{\partial }{Y}_{{\xi }_j}}{\mathrm{\partial }{\xi }_n}\mathrm{\nabla }{\xi }_n)\end{array}$(20) where the 1^st term is the gradient of the conjugate internal molar quantity, the 2^nd and 3^rd terms denote the independent internal potentials and molar quantities controlled from the surroundings, respectively, and the 4^th term represents the contributions from the remaining independent internal molar quantities due to their redistributions in the system. It is noted that $L_{{\xi }_j}$ would also depend on all independent variables in Equation (20), i.e. ${\xi }_j$, $Y_{{\xi }_k}$, $X_{{\xi }_l}$, and ${\xi }_n$. The fluxes as a function of time are evaluated from both the initial and boundary conditions of the system depicted by Equation (12). The 2^nd and 3^rd terms represent the cross phenomena from externally controlled variables, while the 4^th term reflects the cross phenomena due to interactions among internal processes in the system.

The significance of the combined law of thermodynamics, Equation (3), and the flux equation, Equations (18) and (20), is as follows

1. The flow of a molar quantity is solely driven by the gradient of its conjugate potential with a characteristic kinetic coefficient under the linear proportionality approximation, derived from the combined law of thermodynamics without phenomenological considerations.

2. Both the potential gradient and the characteristic kinetic coefficient are functions of all independent variables in an internal process, resulting in the cross-phenomena shown by Equation (20) to be discussed in detail in the rest of the paper.

3. The product of the molar quantity and its conjugate potential results in the entropy production due to the internal process that contributes to the energy of the system as one term in the combined law of thermodynamics, enabling its applicability to non-equilibrium systems.

3.2. Onsager reciprocal relation, its limitations, and improvements

Onsager [2,24] proposed the phenomenological flux equation as follows (21) $J_{{\xi }_j}^{Onsager}=-\sum _{{\xi }_k}{L}_{{\xi }_j{\xi }_k}\mathrm{\nabla }{Y}_{{\xi }_k}$(21) where $L_{{\xi }_j{\xi }_k}$ is the kinetic coefficient for the flux of ${\xi }_j$ due to $\mathrm{\nabla }{Y}_{{\xi }_k}$ with the summation going through gradients of all potentials. Equation (21) implies that all driving forces contribute to the flux of ${\xi }_j$, in an apparent accordance with experimental observations of cross-phenomena. However, it is a postulation, neither from fundamental physics as pointed out by Hillert [10] nor the first and second laws of thermodynamics as mentioned by Balluffi et al. [25]. Following the phenomenological definition, i.e. Equation (21), the rate of entropy production is (22) $\frac{T{\dot{d}}_{ip}S}{V}=-J_{{\xi }_j}^{Onsager}\mathrm{\nabla }{Y}_{{\xi }_j}=\mathrm{\nabla }{Y}_{{\xi }_j}\sum _{{\xi }_k}{L}_{{\xi }_j{\xi }_k}\mathrm{\nabla }{Y}_{{\xi }_k}$(22)

Onsager’s fundamental theorem [2,10,24,26], commonly referred to as ORRs, states that the phenomenological coefficient matrix in Equation (21) is symmetric for independent fluxes and driving forces, i.e. (23) $L_{{\xi }_j{\xi }_k}=L_{{\xi }_k{\xi }_j}$(23)

However, since a symmetric matrix can be diagonalized through the standard principal component analysis, there must exist a set of unique and independent driving forces with $L_{{\xi }_j{\xi }_k}=0$ for ${\xi }_j\ne {\xi }_k$, i.e. eigenvectors of the coefficient matrix that are mutually orthogonal. This is exactly what Equation (18) represents. It is noted that there are two additional fundamental issues related to Equations (21) and (22):

1. When $\mathrm{\nabla }{Y}_{{\xi }_j}=0$, Equation (21) shows that $J_{{\xi }_j}^{Onsager}$ is not equal to zero unless all other $\mathrm{\nabla }{Y}_{{\xi }_k}$ or all off-diagonal $L_{{\xi }_j{\xi }_k}$ values are zero, with the latter equivalent to Equation (18). This means that even though the conjugate driving force for the flow of ${\xi }_j$ is zero, there could be a flux of ${\xi }_j$ due to non-conjugate driving forces based on Equation (21), which is incorrect for a steady-state system to be discussed in the present work.

2. Equation (22) indicates that there will be no entropy production if $\mathrm{\nabla }{Y}_{{\xi }_j}=0$ even when $J_{{\xi }_j}^{Onsager}$ is not zero. This means that the flow of ${\xi }_j$ does not result in entropy, while the second law of thermodynamics stipulates that any irreversible process must result in a positive entropy production, including the flow of ${\xi }_j$.

Both issues are resolved if the flux is represented by Equation (18) derived from the combined law of thermodynamics, while Equation (20) represents the cross phenomena that Onsager’s phenomenological Equation (21) was intended to describe. Equation (18) shows that when there is no gradient of a potential locally, its conjugate molar quantity will not change. While Equation (20) depicts that the vanishing potential gradient can be realized by a set of balanced gradients such as those in steady-state systems. It is self-evident that Equations (18) and (21) are equivalent with $L_{{\xi }_j{\xi }_k}=L_{{\xi }_k{\xi }_j}=0$ for ${\xi }_j\ne {\xi }_k$ so the ORR, Equation (23), is naturally fulfilled. It should be emphasized that $\mathrm{\nabla }{\xi }_j$ and $\mathrm{\nabla }{\xi }_k$ in Equation (20) are not the true driving force for internal processes based on Equations (3) and (4), one should not attempt to apply the ORR to Equation (20) as recently commented by the author [27].

It is noted that Coleman and Truesdell [28] also pointed out the self-contradictory in ORRs and mentioned that ‘unless “forces” and “fluxes” are defined by some property more specific than mere occurrence in the expression for production of entropy, there is no content in the statement that the matrix of phenomenological coefficients is or is not symmetric’. Ågren [29] recently revisited the ORRs and suggested that it can be interpreted as there is a frame of reference where all the transport processes are independent, through illustrations with isobarothermal diffusion and the Kirkendall effect, which is discussed in Section 4.6.

In the rest of the present paper, following cross phenomena will be discussed as examples to illustrate the applications of Equation (20):

1. Thermoelectricity for transport of electron or hole $(c_e\text{or}{c}_h)$, i.e. $\mathrm{\nabla }{\xi }_j=\mathrm{\nabla }{c}_e\text{or}\mathrm{\nabla }{c}_h$, with respect to the gradient of the non-conjugate externally controlled potential, temperature, i.e. $\mathrm{\nabla }{Y}_{{\xi }_j}=\mathrm{\nabla }T$, under the condition that there is no diffusion of other components.

2. Thermodiffusion for transport of a component $(c_i)$, i.e. $\mathrm{\nabla }{\xi }_j=\mathrm{\nabla }{c}_i$, with respect to the gradient of the non-conjugate externally controlled potential, temperature, i.e. $\mathrm{\nabla }{Y}_{{\xi }_j}=\mathrm{\nabla }T$, for a closed system. This situation can be complex if the diffusion of other components takes place at the same time, i.e. $\mathrm{\nabla }{c}_{k\ne i}\ne 0$, which is related to the next example.

3. Isothermal interdiffusion with $\mathrm{\nabla }{\xi }_j=\mathrm{\nabla }{c}_i$ for all components for a closed system without other externally controlled potentials. It is self-evident that the interdiffusion can happen simultaneously in thermodiffusion.

4. Electromigration for transport of component $(c_i)$ and electrocaloric effect for heat conduction $(S)$, i.e. $\mathrm{\nabla }{\xi }_j=\mathrm{\nabla }{c}_i$ and $\mathrm{\nabla }{\xi }_j=\mathrm{\nabla }S$ with respect to the gradient of the non-conjugate externally controlled potential, electrical field, i.e. $\mathrm{\nabla }{Y}_{{\xi }_j}=\mathrm{\nabla }E$.

5. Electromechanical effect for transport of electron $(c_e)$, i.e. $\mathrm{\nabla }{\xi }_j=\mathrm{\nabla }{c}_e$, with respect to the non-conjugate externally controlled potential, stress, i.e. $\mathrm{\nabla }{Y}_{{\xi }_j}=\mathrm{\nabla }\sigma$.

4. Kinetic coefficients and cross-phenomenon coefficients

There are four common types of transport phenomena in the MSE discipline: heat, electron, mass, and fluid, commonly represented by the Fourier’s, Ohm’s, Fick’s, and Darcy’s laws based on experimental observations, with ${\xi }_j=S$, $c_e$, $c_i$ and $V$, and $\mathrm{\nabla }{Y}_{{\xi }_j}=T$, $E$, ${\mu }_i$, and $-P$, respectively. In the Fourier’s, Ohm’s, and Darcy’s laws, Equation (18) is used with the gradient of conjugate potentials as the driving force, so their linear proportionality coefficients represent the kinetic coefficients, i.e. $L_{{\xi }_j}$ in Equation (18). On the other hand, the concentration gradients, i.e. $\mathrm{\nabla }{c}_i$, are used in Fick’s law in terms of Equation (20). Consequently, the linear proportionality coefficient is the product of the kinetic coefficient and thermodynamic factor, commonly referred as diffusion coefficient. It should be emphasized that when Equation (20) is used, all terms should be included, and Fick’s law for diffusion of component $i$ in a temperature gradient can be written in general as one of the following forms under fixed external pressure or volume, respectively, (24) $\begin{array}{rl}{J}_i& =-L_i\mathrm{\nabla }{\mu }_i\\ & =-L_i\left(\sum _j\frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }{c}_j}\mathrm{\nabla }{c}_j+\frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }V}\mathrm{\nabla }V+\frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }T}\mathrm{\nabla }T\right)\\ & =-L_i\left(\sum _j{\mathrm{\Phi }}_{ij}\mathrm{\nabla }{c}_j-P_i\mathrm{\nabla }V-S_i\mathrm{\nabla }T\right)\end{array}$(24) (25) $\begin{array}{rl}{J}_i& =-L_i\mathrm{\nabla }{\mu }_i\\ & =-L_i\left(\sum _j\frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }{c}_j}\mathrm{\nabla }{c}_j+\frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }P}\mathrm{\nabla }P+\frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }T}\mathrm{\nabla }T\right)\\ & =-L_i\left(\sum _j{\mathrm{\Phi }}_{ij}\mathrm{\nabla }{c}_j+V_i\mathrm{\nabla }P-S_i\mathrm{\nabla }T\right)\end{array}$(25) where the summation goes through all independent diffusion components including $i$, ${\mathrm{\Phi }}_{ij}=\frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }{c}_j}$ is the thermodynamic factor between components $i$ and ${\xi }_j$, and $P_i=\frac{\mathrm{\partial }\text{P}}{\mathrm{\partial }{c}_i}$, $S_i=\frac{\mathrm{\partial }\text{S}}{\mathrm{\partial }{c}_i}$, and $V_i=\frac{\mathrm{\partial }\text{V}}{\mathrm{\partial }{c}_i}$ are the partial pressure, partial entropy, and partial volume of component $i$, respectively. ${\mathrm{\Phi }}_{ij}$ can be further written as follows [10] (26) ${\mathrm{\Phi }}_{ij}=\frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }{c}_j}=\frac{\mathrm{\partial }({\mu }_i^0+RT\ln a_i)}{\mathrm{\partial }{c}_j}=\frac{RT}{a_i}\frac{\mathrm{\partial }{a}_i}{\mathrm{\partial }{c}_j}$(26) where ${\mu }_i^0$ and $a_i$ are the reference chemical potential and the activity of component $i$, and $R$ is the gas constant. $P_i$ and $V_i$ can be represented in terms of the Helmholtz energy, $F=U-TS$, and Gibbs energy, $G=F+PV$, under constant temperatures as follows: (27) $\begin{array}{rl}\frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }V}& =\frac{{\mathrm{\partial }}^{2}F}{\mathrm{\partial }V\mathrm{\partial }{c}_i}=-\frac{\mathrm{\partial }P}{\mathrm{\partial }{c}_i}=-P_i\end{array}$(27) (28) $\begin{array}{rl}\frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }P}& =\frac{{\mathrm{\partial }}^{2}G}{\mathrm{\partial }P\mathrm{\partial }{c}_i}=\frac{\mathrm{\partial }V}{\mathrm{\partial }{c}_i}=V_i\end{array}$(28)

The expressions for $S_i$ in Equations (24) and (25) need more discussion. Under the externally fixed pressure for Equation (24) or externally fixed volume for Equation (25), the system exhibits a volume gradient or a pressure gradient in the system, respectively. Therefore, $S_i$ in Equations (24) and (25) are represented in terms of the Helmholtz energy, $F=U-TS$, or Gibbs energy, $G=F+PV$, as follows: (29) $\begin{array}{rl}{\left(\frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }T}\right)}_P& ={\left(\frac{{\mathrm{\partial }}^{2}G}{\mathrm{\partial }T\mathrm{\partial }{c}_i}\right)}_P=-{\left(\frac{\mathrm{\partial }S}{\mathrm{\partial }{c}_i}\right)}_P=-(S_i{)}_P\end{array}$(29) (30) $\begin{array}{rl}{\left(\frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }T}\right)}_V& ={\left(\frac{{\mathrm{\partial }}^{2}F}{\mathrm{\partial }T\mathrm{\partial }{c}_i}\right)}_V=-{\left(\frac{\mathrm{\partial }S}{\mathrm{\partial }{c}_i}\right)}_V=-(S_i{)}_V\end{array}$(30)

Hillert [10] showed that for a stable system, the derivative of a potential to its conjugate molar quantity under constant volume is larger than its value under constant pressure. While no such conclusion for derivatives between two potentials, i.e. Equations (29) and (30), has been reported in the literature and is worth further investigations.

Let us have some general discussions using Equations (24) and (25). Since $S_i$ is positive, both equations indicate that the flux of each component is in the same direction of temperature gradient, i.e. migrating from a low temperature region to a high temperature region when all other gradients are zero. Similarly, since $V_i$ is positive, Equation (25) indicates that the flux of each component is in the opposite direction of pressure gradient, i.e. migrating from high pressure to low pressure regions when all other gradients are zero. Since the values of $S_i$ and $V_i$ are different for each component, their fluxes are different, resulting in inhomogeneous distributions of components along the gradients. Further complications occur as the migration introduces both the concentration and entropy gradients with the latter generating a temperature gradient. The final steady state is reached with the balance of contributions of $\mathrm{\nabla }{c}_j$, $\mathrm{\nabla }V$ or $\mathrm{\nabla }P$, and $\mathrm{\nabla }T$ to flux so the net flux is zero under the external constraints of the mass conservation and the constant temperature in the system. More detailed discussions on individual cases will be given in the next several sections.

From the above discussions, it is clear that only $L_{{\xi }_j}$: s are the kinetic coefficients, and the cross-effect coefficients are thermodynamic properties represented by derivatives in Equation (20) in general and Equations (24) and (25) for mass transports. The derivatives between two potentials are relatively easy to measure as temperature, pressure, and electrical field are controlled experimentally, while the derivative between molar quantities is straightforward to predict computationally [1,10,11]. This represents an integration of complimentary experimental and computational views of the same phenomenon connected by the Maxwell relation as follows with various derivatives between potentials further discussed in Section 4.10 (31) $\begin{array}{rl}\frac{\mathrm{\partial }{Y}^a}{\mathrm{\partial }{X}^b}& =\frac{{\mathrm{\partial }}^2\mathrm{\Phi }}{\mathrm{\partial }{X}^b\mathrm{\partial }{X}^a}=\frac{\mathrm{\partial }{Y}^b}{\mathrm{\partial }{X}^a}\end{array}$(31) (32) $\begin{array}{rl}\frac{\mathrm{\partial }{Y}^a}{\mathrm{\partial }{Y}^b}& =\frac{{\mathrm{\partial }}^2\mathrm{\Phi }}{\mathrm{\partial }{Y}^b\mathrm{\partial }{X}^a}=-\frac{\mathrm{\partial }{X}^b}{\mathrm{\partial }{X}^a}\end{array}$(32)

4.1. Kinetic coefficients, $L_{{\xi }_j}$

It is natural that theoretical calculations of $L_{{\xi }_j}$ are commonly approached from the kinetics of internal processes. Fundamentally, any internal processes introduce local short- or global long-range disturbances, resulting in a barrier against the internal process from its environment or the surroundings of the internal process. This is true even in Maxwell’s thought experiment where the Maxwell’s daemon must make some disturbances when opening or closing a small ‘massless' door between two chambers of gas. Consequently, $L_{X_j^a}$ may be generally represented by the following equation (33) $L_{{\xi }_j}=L_{{\xi }_j}^0{e}^{-\frac{Q_{{\xi }_j}}{k_{B}T}}$(33) where $L_{{\xi }_j}^0$ and $Q_{{\xi }_j}$ are the prefactor and the barrier for the internal process, respectively, both being functions of all independent variables, and $k_B$ is the Boltzmann constant. A schematical energy profile as a function of one internal process in Figure 2 depicts the transition from one state to another state through an energy barrier [1]. More complex energy landscape with multiple internal processes such as fractal free energy landscapes with simple basins, metabasins and fractal basins in structural glasses [30] can in principle be considered as being composed of many such individual energy profiles shown in Figure 2 at various time and spatial scales.

Figure 2. Schematic diagram of energy landscape as a function of one internal process [1].

Internal processes with $D_j>0$ may be broadly categorized into following cases in terms of the values of $Q_{{\xi }_j}$ with respect to $k_{B}T$:

1. $Q_{{\xi }_j}\le 0$. The system is unstable locally with respect to the internal process, resulting in a local stochastic evolution that is affected by interactions with its environment such as spinodal decompositions where self-assembled structures form [31] or in general the dissipated structures [32].

2. $Q_{{\xi }_j}\approx 0$. The internal process is near barrierless. A small disturbance from the surroundings will activate the internal process for it to continue perpetually, such as the electrical current in a superconductor.

3. $0<Q_{{\xi }_j}<k_{B}T$. The barrier is small so that intermediate states along the pathway of the internal process are not distinguishable due to the time and spatial resolutions of measurement techniques and can thus only be described by probability in space predicted in terms of wave functions such as migrations of electrons and phonons in electrical and thermal conductors.

4. $Q_{{\xi }_j}>k_{B}T$. The internal process is activated thermally with low probability, and the intermediate states along the pathway can be discretely investigated such as atomic diffusion.

5. $Q_{{\xi }_j}\gg k_{B}T$. The internal process is a rare event such as migrations of electrons and phonons in electrical and thermal insulators and transition of diamond to graphite under ambient conditions.

The present discussion concerns mainly internal processes in Case III and IV, and they all experience the scenario described by Case I during overcoming the barriers. It may be mentioned that we have developed a multiscale entropy approach (termed zentropy) to predict instability [14,20], which may have some implications in understanding the Case I and II internal processes in terms of their occurrences and driving forces and will be discussed in Section 5, including postulations on superconductivity in Section 5.5.

4.2. Boltzmann transport equation for electrical and thermal conductions

For electrical and thermal conductions in Case III, the Boltzmann transport equation (BTE) is commonly used [33–36]. The BTE describes the behavior of a nonequilibrium system in terms of a balance between scattering in and out of each possible state, with scalar scattering rates. In the widely used BoltzTraP program [35,37], the BTE is linearized under the relaxation time approximation (RTA) to describe the transport distribution function which is then used to calculate the moments of the generalized transport coefficients and the charge and heat currents. The kinetic coefficients are evaluated under two extreme cases, i.e. with only the electrical gradient or only the temperature gradient applied externally, respectively. The former is slightly different from the case without a temperature gradient in the system as discussed below. In the latter, there is usually an additional constraint, i.e. without electrical current in the system.

Let us first consider the transport of electrons in conductors or n-type semiconductors with electrons added to the conduction band. When an electrical field is applied externally to a system initially under equilibrium with homogeneous temperature and electron distributions, the initial internal process is the transport of electrons, i.e. $d\xi =d{c}_e$, resulting in an electrical current with the flux of electrons as follows: (34) $J_e=-L_e\mathrm{\nabla }{\mu }_e=L_{e}E$(34) where $L_e$ is the electrical conductivity, and ${\mu }_e$ the chemical potential of electrons. By evaluating $J_e$ under given $E$ applied externally to the system from the BTE simulations, $L_e$ can be obtained by Equation (34). It should be emphasized that both $L_e$ and ${\mu }_e$ depend on all independent internal variables and externally applied electrical field (see Equation (12)).

At the same time, the migrating electrons carry heat with them, inducing a heat current, $J_Q$, in the system as follows (35) $J_Q=T{S}_e{J}_e=-T{S}_e{L}_e\mathrm{\nabla }{\mu }_e$(35) where $S_e$ is the partial entropy of electrons, i.e. the heat carried by electrons. It is important to point out that this heat conduction is induced by the heat carried by the electrons in the charge current, not by the external temperature gradient. This results in the Peltier effect due to the fact that an electrical current is accompanied by a heat current in a homogeneous conductor even at constant temperature. The Peltier coefficient is thus evaluated by the division of heat current to the electrical current as follows, noting that the electrical current is in the opposite direction of that of electron flux (36) $\mathrm{\Pi }=\frac{J_Q}{-J_e}=-T{S}_e$(36)

In the second case, a temperature gradient is applied externally to the system and generates an internal process of heat current as follows (37) $J_Q=-L_Q\mathrm{\nabla }T$(37) where the thermal conductivity $L_Q$ can be evaluated by measuring $J_Q$ for given $\mathrm{\nabla }T$. Before the temperature gradient is applied, the electrons are homogeneously distributed in the system. When the temperature gradient is applied, the heat current results in another internal process, i.e. the migration of electrons, thus an electrical current as follows (38) $\begin{array}{rl}{J}_e& =-L_e\mathrm{\nabla }{\mu }_e=-L_e\left(\frac{\mathrm{\partial }{\mu }_e}{\mathrm{\partial }{c}_e}\mathrm{\nabla }{c}_e+\frac{\mathrm{\partial }{\mu }_e}{\mathrm{\partial }T}\mathrm{\nabla }T\right)\\ & =-L_e({\mathrm{\Phi }}_{ee}\mathrm{\nabla }{c}_e-S_e\mathrm{\nabla }T)\end{array}$(38) where ${\mathrm{\Phi }}_{ee}$ is the thermodynamic factor of electrons. As shown by Equation (35), this electrical current further contributes to heat conduction in the system. Furthermore, there can be a variation of volume along the way as shown by Equation (24), which is omitted here. With the initial condition of $\mathrm{\nabla }{c}_e=0$, $J_e=L_e{S}_e\mathrm{\nabla }T$, i.e. electrons migrate from the hot end to the cold end in the system instantaneously.

Under the condition that there is no electrical current across the system in a steady state, the inhomogeneous distribution of electrons results in an induced internal electrical field, i.e. a voltage between the two ends of the system due to electrostatic interactions. This voltage can be measured by forming an electrical circuit with an external voltage in the opposite direction, resulting in zero current in the circuit. By equating Equation (38) to zero, one obtains the following equation, noting the negative charge of electrons: (39) $J_e=-L_e(-\mathrm{\nabla }{V}_e-S_e\mathrm{\nabla }T)=0$(39)

The Seebeck coefficient is thus defined as follows: (40) $S_{T,e}=\frac{\mathrm{\nabla }{V}_e}{\mathrm{\nabla }T}=-S_e$(40)

Combining Equations (36) and (40), one obtains the Thomson relation, i.e. (41) $\mathrm{\Pi }=T{S}_{T,e}$(41)

The above discussion demonstrates that both Peltier and Seebeck coefficients are thermodynamic quantities related to $\mathrm{\partial }{\mu }_e/\mathrm{\partial }T$ as defined in Equation (24). This derivative between two potentials is related to the partial entropy of electron through the Maxwell relation from the combined law in terms of Gibbs energy used in Equation (32) with $T$ and $c_e$ as its natural variables as follows: (42) $S_{T,e}=\frac{\mathrm{\partial }{\mu }_e}{\mathrm{\partial }T}=\frac{{\mathrm{\partial }}^{2}G}{\mathrm{\partial }T\mathrm{\partial }{c}_e}=-\frac{\mathrm{\partial }S}{\mathrm{\partial }{c}_e}=-S_e$(42)

For p-type thermoelectric materials with positively charged holes added to the valence band, one has: (43) $\begin{array}{rl}{J}_h& =-L_h(\mathrm{\nabla }{V}_h-S_h\mathrm{\nabla }T)=0\end{array}$(43) (44) $\begin{array}{rl}{S}_{T,h}& =\frac{\mathrm{\nabla }{V}_h}{\mathrm{\nabla }T}=S_h\end{array}$(44) where $L_h$, $\mathrm{\nabla }{V}_h$, $S_h$, and $S_{T,h}$ are the kinetic coefficient, voltage, partial entropy, and Seebeck coefficient of holes in p-type thermoelectric materials, respectively. One thus has negative and positive Seebeck coefficients for n- and p-type thermoelectric materials, respectively. It needs to emphasize that this voltage is induced by the gradient of electron or hole concentrations in the system due to the heat conductions from the temperature gradient. As there is no external electrical field applied, the effect of temperature gradient on the chemical potential of electrons or holes is balanced by the internal voltage due to the concentration gradient of electrons or holes, resulting in homogeneous chemical potential of electrons or holes in the system and zero electrical current as shown by Equation (39) or Equation (43). Since both p- and n-type materials can be described by the same formula, only electrons will be used in the rest of the paper unless their differentiation is necessary.

It may be mentioned that the Green-Kubo formalism [38,39] based on statistical mechanics relates the kinetic coefficients to the integrated time correlation functions which involve the microscopic particle fluxes and heat flux. Recently, Green et al. derived a formula to describe the evolution of density excitations of the electron system when subject to a periodic potential with static or dynamical perturbations [40]. It is based on Green-Kubo formalism with a Mori memory function approach to include dissipation effects at all orders in the relaxation interaction. In their approach, the excitations of a given wave vector are kept as an active subspace, and those of different wave vectors are projected into a bath described by a memory function, while the slower variables within the active subspace are defined by the conservation laws. Their approach facilitates the approximation of the dissipation process into a single memory function that includes relaxation, dissipation, and quantum coherence. It is shown that this approach offers a practical method to evaluate the conductivity with existing electronic-structure codes and avoids the complications and the limitations of the Green-Kubo formula in the thermodynamic limit and the neglect of quantum coherence in the BTE approach. The potential importance of interband coherence on the transport process in systems with more than one band close to the Fermi level is emphasized beyond weak scattering captured by a relaxation time.

4.3. Atomic mobility by kinetic simulations and DFT-based calculations

For atomic diffusion in Case IV in Section 4.1, classic or ab-initio molecular dynamic (MD) [41–43] and kinetic Monte Carlo (KMC) [44–47] simulations directly follow the migration of diffusion component as a function of time and extract the self or tracer diffusivity ($D_i^{\ast }$) using Einstein’s relation in terms of the mean square displacement (MSD) [48]. The tracer-diffusivity can be further converted to atomic mobility ($M_i$), the kinetic coefficient discussed above ($L_i$), and intrinsic diffusivity (${}^i{D}_{ik}=L_i\frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }{c}_k}$) in the lattice-fixed frame of reference, and chemical diffusivity ($D_{ik}=\sum _j{L}_{ij}^{\mathrm{\prime }}\frac{\mathrm{\partial }{\mu }_j}{\mathrm{\partial }{c}_k}$) in the volume-fixed frame of reference with a complex relation between $L_{ij}^{\mathrm{\prime }}$ and $L_i$. Andersson and Ågren [49] presented an elegant discussion among all the diffusion coefficients and their computational implementation. The relationships among these diffusion coefficients and kinetic coefficients are shown in Figure 3.

Figure 3. Relationships among tracer diffusivity, atomic mobility, kinetic L parameters, and intrinsic and chemical diffusivities.

As indicated by Figure 1 and discussed above, thermodynamics and kinetics are closely related to each other, and it is thus natural that atomic diffusion process can be viewed from either thermodynamics or kinetics of the internal process. This can be accomplished by evaluating the free energy of vacancy formation and the individual microstates along the diffusion pathways. The free energy of vacancy formation can be obtained straightforwardly by the DFT-based first-principles calculations under the quasiharmonic phonon approximations, which also results in the vacancy concentration as a function of temperature [50,51]. In the classic transition state theory (TST) [52,53], the jump frequency for the migration of atom over the barrier is written in terms of enthalpy of barrier and an effective frequency, ${\nu }^{\ast }$, represented by the vibrational frequencies at the equilibrium and transition states with the imaginary vibrational modes at the transition state removed as follow: (45) ${\nu }^{\ast }=\frac{{\prod }_{i=1}^{3N-3}{\nu }_i}{{\prod }_{j=1}^{3N-4}{\nu }_j^{{}^{\mathrm{\prime }}}}$(45) where ${\nu }_i$ and ${\nu }_j^{\mathrm{\prime }}$ are the vibrational frequencies at the equilibrium and transition states, respectively. The enthalpy of migration and the vibrational frequencies can be obtained from the DFT-based first-principles calculations [54] in combination with the nudged elastic band (NEB) method [55] and the five frequency model [56] along with the predicted vacancy concentration. We were then able to predict the tracer diffusion coefficients for pure elements and dilute solutions for fcc [57–61], bcc [62], and hcp [62–65] structures, showing remarkable agreement with available experimental data. A similar approach was also developed for prediction of interstitial diffusion coefficients by Wimmer et al. [66]. It should be mentioned that this approach is predictive with all inputs from DFT-based first-principles calculations without empirical parameters. More recently, Allnatt et al.[67] developed a 13 frequency model to address diffusion kinetics for the full anisotropic three-dimensional hcp structure and verified the some of the results with Monte Carlo simulations.

In next several sections, some classic cross phenomena, i.e. thermoelectricity for electrical conduction and thermodiffusion for atomic diffusion, both induced by the temperature gradient, upfill diffusion for migration of fast diffusion component against its concentration gradient induced by inhomogeneity of slow diffusion components, electromigration for atomic diffusion and the electrocaloric effect for heat conduction, both induced by electrical field, and electromechanical effect for interactions between mechanical deformation and electrical field are discussed in detail by applying the theoretical framework presented above.

4.4. Thermoelectricity: a cross phenomenon between electrical conduction and temperature

The thermoelectricity case is discussed in Section 4.2 for a system under external electrical field or external temperature gradient, respectively. It is shown that the Seebeck coefficient is defined by the derivative of the chemical potential of electrons with respect to temperature as shown by Equations (42) and (44). At the same time, the Peltier coefficient and the Thomson relation are also obtained.

Based on this fundamental understanding, we predicted values of Seebeck coefficients for several n- and p-type thermoelectric materials using DFT-based first-principles calculations without empirical parameters [68,69]. In those calculations, the temperature dependence of the chemical potential of electrons are evaluated from the electron density of states (e-DOS) using the Mermin statistics [70,71] extended to DFT-based calculations at finite temperatures [50,72,73] with the electrons explicitly treated [74,75]. We start from the conservation equation: (46) $\int n(\epsilon ,V)fd\epsilon =N_e$(46) where $n(\epsilon ,V)$ is the e-DOS at the band energy $\epsilon$ and volume $V$, N_e the total number of electrons, and $f$ the Fermi-Dirac distribution expressed as follows: (47) $f=\frac{1}{exp\left[\frac{\epsilon +e\phi }{k_{B}T}\right]+1}$(47) where $e$ is the charge of electron, and $e\phi =-({\mu }_e-{\epsilon }_F)$ with ${\epsilon }_F$ being the Fermi energy. After a few rearrangements, we obtain the formulation to calculate the Seebeck coefficient under the isobaric condition as follows [68,69]: (48) $\begin{array}{rl}{S}_{T,e}& =-{\left(\frac{\mathrm{\partial }\phi }{\mathrm{\partial }T}\right)}_P=-\frac{k_{B}T{\alpha }_V}{e{c}_e}\int \frac{\mathrm{\partial }n(\epsilon ,V)}{\mathrm{\partial }V}fd\epsilon \\ & -\frac{1}{eT{c}_{e}V}\int n(\epsilon ,V)f(1-f)(\epsilon +e\phi )d\epsilon \end{array}$(48) (49) $\begin{array}{rl}{c}_e& =\frac{1}{V}\int n(\epsilon ,V)f(1-f)d\epsilon \end{array}$(49) where ${\alpha }_V$ is the volume thermal expansion coefficient, and $c_e$ the mobile charge carrier concentration. Therefore, the Seebeck coefficient can be computed from the readily accessible e-DOS without involving the electrical conductivity as in BTE, which is a much more challenging quantity to compute and requires the calculation of electron group velocity and the estimation of the relaxation time [76].

The integral $\int n(\epsilon ,V)f(1-f)d\epsilon$ in Equation (49) can be considered as the total number of mobile carriers in a thermoelectric solid, i.e. the number of electrons participating in the conduction process at finite temperatures. This can be understood as the pair product of $f(1-f)$ represents the probability that the electrons occupied ‘$f$ number' of electronic states with energy $\epsilon$, transmitted to (or the vice versa), the ‘$1-f$ number' of unoccupied electronic states with energy ε at finite temperature. Consequently, $n(\epsilon ,V)f(1-f)$ represents the e-DOS of charge carriers at finite temperatures.

For DFT-based computations, the issue related to the commonly underestimated band gap is resolved by following the approach developed by Singh [77] with the Engel-Vosko generalized gradient approximation (GGA) [78] plus spin-orbit coupling as implemented in the WIEN2k package for obtaining accurate electronic structures [79]. However, this approach has low accuracy for the total energy, which is instead computed by means of the Perdew-Burke-Ernzerhof exchange-correlational functional revised for solids (PBEsol) [80] as implemented in the Vienna ab initio simulation package (VASP, version 5.3) [74,75]. Furthermore, our mixed space approach [11,81–83] is used to account for the dipole-dipole interactions for phonon calculations of polar insulators in the real space using the supercell method. For simplicity, we assume the rigid band approximation, [34] i.e. the band structure is assumed to remain unchanged from doping or from the thermal electronic excitation at finite temperatures. Our approach demonstrated excellent agreement with experimental observations as shown in Figure 4 for PbTe [68], SnSe [68], and $L{a}_{3-x}T{e}_4$ [69]. In Figure 4(ii), the Seebeck coefficients reported in the literature [84] using the previous version of the BoltzTraP program [85] are superimposed, showing worse agreement with experimental data in comparison with our calculations without empirical parameters.

Figure 4. Calculated Seebeck coefficients for (i) PbTe for various p- and n-type doping levels [68]; (ii) p-type SnSe [68]; (iii) $L{a}_{2.75}T{e}_4$ [69].

It is important to emphasize that regardless of whether electrons are driven by external temperature gradient or external electrical field, it is the migration of electrons that governs the kinetic process of electrical conduction, i.e. the electrical conductivity and the chemical potential of electrons. Furthermore, both electrical conductivity and the chemical potential of electrons are functions of temperature and electron concentration, resulting in the cross-phenomenon of thermoelectricity. A commonly used parameter measuring the performance of thermoelectric materials is the figure of merit [86] defined as follows: (50) $z_T=\frac{L_e{S}_{T,e}^2}{L_Q}T$(50) which includes the two kinetic coefficients $L_e$ and $L_Q$ and the thermodynamic Seebeck coefficient $S_{T,e}$. The square of $S_{T,e}$ shows its importance in affecting the performance of thermoelectric materials and can be reliably predicted by our approach discussed here.

4.5. Thermodiffusion: a cross phenomenon between atomic diffusion and temperature

For thermodiffusion, the chemical potential of a diffusion component drives the transport of the component from a high chemical potential region to a low chemical potential region until its chemical potential becomes homogeneous in the system [19,87]. As in the case of thermoelectricity, when a temperature gradient is applied externally to a system initially at equilibrium with homogeneous compositions, the chemical potentials of component in the system become inhomogeneous, resulting in an internal process for the migration of diffusion component as shown by Equations (24) or (25). It is important to emphasize again that the fundamental process of thermodiffusion is the migration of the diffusion component, driven by the gradient of its conjugate chemical potential, and the chemical potential and kinetic coefficient are functions of temperature and compositions of all components plus all independent internal variables as shown by Equation (12).

Let us start with a system with homogeneous concentration, i.e. $\mathrm{\nabla }{c}_j=0$. When a temperature gradient is applied initially, two internal processes take place immediately in the system, i.e. the volume change and heat conduction, which will drive the initial diffusion flux as follows by re-writing Equation (24): (51) $J_i=-L_i(-P_i\mathrm{\nabla }V-S_i\mathrm{\nabla }T)$(51)

For a system with positive thermal expansion, $\mathrm{\nabla }V$ and $\mathrm{\nabla }T$ point to the same direction, so does $J_i$ for all diffusion components, though it should be noted that there are both liquid and solid phases that exhibit negative thermal expansion, i.e. $\mathrm{\nabla }V$ and $\mathrm{\nabla }T$ can point to the opposite directions in some systems [1,20,88].

As soon as the diffusion starts, the concentration gradients will be created and contribute to the flux as shown by Equation (24). As typical physical experiments are conducted for closed systems, the fluxes of all diffusion components must be zero at both ends of the system as the boundary condition. Therefore, when the system reaches a steady state with heat conduction as the only internal process in the system, $\mathrm{\nabla }V$, $\mathrm{\nabla }T$ and $\mathrm{\nabla }{c}_j$ balance each other to result in zero diffusion flux for each diffusion component, i.e.: (52) $\begin{array}{rl}{J}_i& =-L_i\mathrm{\nabla }{\mu }_i=-L_i\left(\sum _j{\mathrm{\Phi }}_{ij}\mathrm{\nabla }{c}_j-P_i\mathrm{\nabla }V-S_i\mathrm{\nabla }T\right)\\ & =0\end{array}$(52)

For non-ideal, non-dilute solutions, $L_i$, ${\mathrm{\Phi }}_{ij}$, $P_i$, and $S_i$ can strongly depend on temperature and compositions of the solutions, resulting in a diffusion component switching their segregation regions as a function of composition and temperature [89–96]. Equation (52) further indicates that at the steady state, $\mathrm{\nabla }{\mu }_i=0$ for each diffusion component, i.e. homogeneous chemical potential.

As most experiments reported in the literature are for binary systems, let us write Equation (52) for a binary A-B system under steady state as follows (53) $\begin{array}{rl}{J}_A& =-L_A\mathrm{\nabla }{\mu }_A=-L_A({\mathrm{\Phi }}_{AA}\mathrm{\nabla }{c}_A+{\mathrm{\Phi }}_{AB}\mathrm{\nabla }{c}_B\\ & -P_A\mathrm{\nabla }V-S_A\mathrm{\nabla }T)=0\end{array}$(53) (54) $\begin{array}{rl}{J}_B& =-L_B\mathrm{\nabla }{\mu }_B=-L_B({\mathrm{\Phi }}_{BA}\mathrm{\nabla }{c}_A+{\mathrm{\Phi }}_{BB}\mathrm{\nabla }{c}_B\\ & -P_B\mathrm{\nabla }V-S_B\mathrm{\nabla }T)=0\end{array}$(54)

Eliminating $\mathrm{\nabla }{c}_A$ from the equations results in: (55) $\begin{array}{rl}({\mathrm{\Phi }}_{BB}{\mathrm{\Phi }}_{AA}-{\mathrm{\Phi }}_{AB}{\mathrm{\Phi }}_{BA})\mathrm{\nabla }{c}_B& =(P_B{\mathrm{\Phi }}_{AA}-P_A{\mathrm{\Phi }}_{BA})\mathrm{\nabla }V\\ & +(S_B{\mathrm{\Phi }}_{AA}-S_A{\mathrm{\Phi }}_{BA})\mathrm{\nabla }T\end{array}$(55)

The Soret coefficient in thermodiffusion [97] is usually evaluated experimentally in the literature by the negative ratio of the ${\mathrm{\nabla }}_{C_B}$ with respect to ∇T as follows: (56) $\begin{array}{rl}{S}_{T,B}& =-\frac{\mathrm{\nabla }{c}_B}{c_B\mathrm{\nabla }T}\\ & =-\frac{1}{c_B}\frac{S_B{\mathrm{\Phi }}_{AA}-S_A{\mathrm{\Phi }}_{BA}}{{\mathrm{\Phi }}_{BB}{\mathrm{\Phi }}_{AA}-{\mathrm{\Phi }}_{AB}{\mathrm{\Phi }}_{BA}}\\ & \times \left(1+\frac{P_B{\mathrm{\Phi }}_{AA}-P_A{\mathrm{\Phi }}_{BA}}{S_B{\mathrm{\Phi }}_{AA}-S_A{\mathrm{\Phi }}_{BA}}\frac{\mathrm{\nabla }V}{\mathrm{\nabla }T}\right)\end{array}$(56)

Graphically, the Soret coefficient is thus the negative slope of the concentration in log scale plotted with respect to temperature (see e.g. the plot by Duhr and Braun [19]). Equation (56) shows that $S_{T,B}=0$ when the one of the following conditions is met: (57) $\begin{array}{rl}{S}_B{\mathrm{\Phi }}_{AA}-S_A{\mathrm{\Phi }}_{BA}& =0\end{array}$(57) (58) $\begin{array}{rl}{S}_B{\mathrm{\Phi }}_{AA}-S_A{\mathrm{\Phi }}_{BA}+(P_B{\mathrm{\Phi }}_{AA}-P_A{\mathrm{\Phi }}_{BA})\frac{\mathrm{\nabla }V}{\mathrm{\nabla }T}& =0\end{array}$(58)

For systems with negative thermal expansion at certain composition and temperature ranges [1,20,88], it may be more likely that the Soret coefficient changes its sign with respect to composition and temperature.

As pointed out by Duhr and Braun [19] and commonly agreed in the literature, a generally applied assumption for describing Soret effect is that it should be treated as a transport phenomenon. It is certainly true that Soret coefficients can be evaluated from kinetic processes as discussed in Section 4.2, such as the Green-Kubo formalism [38,39] for transport properties in liquid of binary systems using MD simulations [98–102]. Nevertheless, there were important efforts in the literature in addition to the present work to show that the Soret coefficient is a thermodynamic property of the system, rather than a transport property. Duhr and Braun [19] divided the nonequilibrium system into a succession of regions with small differences of local temperature and concentration and utilized a linearized Boltzmann distribution under local equilibria as follows: (59) $S_{T,B}dT=-\frac{d{c}_B}{c_B}=-\frac{d{\mu }_B}{k_{B}T}=\frac{S_B}{k_{B}T}dT$(59) where $c_B$, ${\mu }_B$, and $S_B$ are the concentration, chemical potential, and partial entropy of the particle $B$. In the publications by Duhr and Braun [19,103] $G$ was used in the place of ${\mu }_B$ in Equation (59), noting that ${\mu }_B={\mu }_B^0+k_{B}Tln{c}_B$ with ${\mu }_B^0$ being the reference state for ${\mu }_B$. They further correlated ${\mu }_B$ to the particle radius by considering that for solid particles, only the solvation energy at their surface can be temperature dependent. The Soret coefficient was thus obtained as follows: (60) $S_{T,B}=\frac{S_B}{k_{B}T}$(60)

It is noted that in their paper, there was a minus sign on the right-hand side of Equation (60), which is probably a typo.

It is evident that Equations (56) and (60) are substantially different. The issue is due to the relation defined by Equation (59). It is noted that ${\mu }_B$ is a function of $c_B$, $V$, and $T$ as shown in Equations (8), (12), (24), (53) and cannot be written as a function of either $c_B$ or $T$ as in Equation (59). The general differential form of ${\mu }_B$ in a binary system is to be written as follows: (61) $d{\mu }_B=\frac{\mathrm{\partial }{\mu }_B}{\mathrm{\partial }{c}_B}d{c}_B+\frac{\mathrm{\partial }{\mu }_B}{\mathrm{\partial }{c}_A}d{c}_A+\frac{\mathrm{\partial }{\mu }_B}{\mathrm{\partial }V}dV+\frac{\mathrm{\partial }{\mu }_B}{\mathrm{\partial }T}dT$(61)

Furthermore, Equation (60) implies that the Soret coefficient is always positive as $S_B$ is positive, which is in agreement with the experimental data by Duhr and Braun [19,103], but not true in general as reported in the literature [89–96]. While $d{c}_B+d{c}_A=0$ is commonly used in the literature, it is a good approximation for most liquid systems, but not true in general as the net flux of all diffusion components may not be zero, resulting in a net vacancy flow, i.e. the Kirkendall effect in solids to be discussed in Section 4.6.

Kocherginsky and Gruebele made an important contribution in a series of publications [104–106] through a systematic framework for chemical transport. Based on the continuity equation for diffusive transport and the Smoluchowski equation, they obtained the flux of a component $i$ as follows: (62) $J_i=-M_i{c}_i\mathrm{\nabla }{\mu }_i$(62) where $M_i$ (denoted by $U_i$ in their work), $c_i$ and ${\mu }_i$ (denoted by ${\mu }_{gi}$ in their work) are Einstein’s mobility, concentration, and physicochemical potential of component $i$, respectively. By comparing Equations (24) and (62), one obtains $L_i=M_i{c}_i$, which was shown by Andersson and Ågren [49] (see also Figure 3) for the lattice-fixed frame of reference. Even though the present author has some disagreements with Kocherginsky and Gruebele on how to extend the dependence of $\mathrm{\nabla }{\mu }_i$ to other gradients of independent variables (see Equations (20) and (24)) as discussed recently [27,107], Equation (62) derived by Kocherginsky and Gruebele [104–106] from kinetic considerations is significant as it demonstrates that all cross-effect coefficients for atomic diffusion are related to the derivatives of chemical potential with respect to other independent variables.

Köhler and Morozov [108] reviewed the experimental observations and theoretical work in the literature for the Soret effect in non-ionic molecular binary and ternary liquid mixtures, and Piazza and Parola [109] and Würger [110] reviewed the Soret effect in the ionic and surfacted colloids and micellar solutions. For binary systems, Hartmann et al. [95] defined the parameters for each pure substance and obtained the following equation for Soret coefficient defined by component 1 in a binary system as follows: (63) $S_{T,1}=\frac{f_1-f_2+C(\mathrm{\partial }{V}^E/\mathrm{\partial }{x}_1)}{1+\mathrm{\partial }ln{\gamma }_1/\mathrm{\partial }ln{x}_1}$(63) where ${\gamma }_1$ and $x_1$ are the chemical activity coefficient and mole fraction of component 1, $V^E$ is the non-ideal excess volume of mixing, $C$ is a constant, and $f_1$ and $f_2$ are model parameters related to Soret coefficients of pure components 1 and 2, respectively. They qualitatively correlated the model parameters to characteristic diameters of the solvent and solute, their interaction parameters, the solvent volume fraction, and the partial excess volume of mixing, thus implying that the Soret coefficient is a thermodynamic quantity. It can be seen that Equations (56) and (63) contain similar thermodynamic quantities such as thermodynamic factor and the effect of volume denoted by $\mathrm{\partial }ln{\gamma }_1/\mathrm{\partial }ln{x}_1$ and $V^E$ in Equation (63), and partial entropies in Equation (56) may be related to $f_1$ and $f_2$ in Equation (63) as shown below: (64) $S_i=S_i^0-R\ln x_i+\mathrm{\Delta }{S}_i^{ex}$(64) where $S_i^0$ is the partial entropy of the pure component $i$, and $\mathrm{\Delta }{S}_i^{ex}$ the excess non-ideal partial entropy of the component $i$ in the solution.

More recently, Schraml et al. [96] presented a detailed investigation of water ($H_{2}O$), ethanol ($ETH$), and triethylene glycol ($TEG$) ternary and the three constitutive binary systems and observed a characteristic sign change of the Soret coefficient as a function of temperature and/or composition. They correlated the decay of the Soret coefficient with respect to concentration to negative excess volumes of mixing as shown by Equation (63), i.e. the decrease of volume with the increase of concentration of the minor component. Since the minor component migrates towards the cold side, the higher concentration of the minor component at the low temperature region reduces its volume further and increases the ratio of $\mathrm{\nabla }V/\mathrm{\nabla }T$ (see Equations (56) and (58)). When the excess volume reaches the minimum, further increase of the concentration decreases the ratio $\mathrm{\nabla }V/\mathrm{\nabla }T$. This change along with the interplay of other properties in Equations (56) and (58) results in the sign change of the Soret coefficients in the systems studied.

Furthermore, Schraml et al. [96] showed the connection between the sign change of the Soret coefficients of the binaries and the signs of the Soret coefficients of the corresponding ternary mixtures as depicted in Figure 5. In particular, they observed that at least one ternary composition exists, where all three Soret coefficients vanish simultaneously, and no steady-state separation is observable. They postulated that this composition is on the path connecting the two compositions of zero Soret coefficients in the two binary systems. This can be understood by writing the flux equations for a ternary system as follows: (65) $\begin{array}{rl}{J}_A& =-L_A\mathrm{\nabla }{\mu }_A=-L_A({\mathrm{\Phi }}_{AA}\mathrm{\nabla }{c}_A+{\mathrm{\Phi }}_{AB}\mathrm{\nabla }{c}_B\\ & +{\mathrm{\Phi }}_{AC}\mathrm{\nabla }{c}_C-P_A\mathrm{\nabla }V-S_A\mathrm{\nabla }T)=0\end{array}$(65) (66) $\begin{array}{rl}{J}_B& =-L_B\mathrm{\nabla }{\mu }_B=-L_B({\mathrm{\Phi }}_{BA}\mathrm{\nabla }{c}_A+{\mathrm{\Phi }}_{BB}\mathrm{\nabla }{c}_B\\ & +{\mathrm{\Phi }}_{BC}\mathrm{\nabla }{c}_C-P_B\mathrm{\nabla }V-S_B\mathrm{\nabla }T)=0\end{array}$(66) (67) $\begin{array}{rl}{J}_C& =-L_C\mathrm{\nabla }{\mu }_C=-L_C({\mathrm{\Phi }}_{CA}\mathrm{\nabla }{c}_A+{\mathrm{\Phi }}_{CB}\mathrm{\nabla }{c}_B\\ & +{\mathrm{\Phi }}_{CC}\mathrm{\nabla }{c}_C-P_C\mathrm{\nabla }V-S_C\mathrm{\nabla }T)=0\end{array}$(67)

Figure 5. Signs of the Soret coefficients in the water (H2O), ethanol (ETH), and triethylene glycol (TEG) ternary system at 25 C [96]. The colored regions denote regions of negative Soret coefficients of the respective components. Point Z marks the intersection of the boundaries of the three colored regions, where all three Soret coefficients vanish simultaneously. The steady state optical signal vanishes along the dashed line.

Using these three equations, one can calculate the Soret coefficients of each component for given compositions in a ternary system. As observed by Schraml et al. [96], there is one composition for $S_{T,H_{2}O}=S_{T,ETH}=0$ in the $H_{2}O-ETH$ binary system and $S_{T,H_{2}O}=S_{T,TEG}=0$ in the $H_{2}O-TEG$ binary system, while in the $ETH-TEG$ binary system this happens near pure $ETH$. When the third component is added, Schraml et al. [96] discussed that there should be three composition pathways from each binary into the ternary with one Soret coefficient kept zero. They further articulated that these three composition pathways meet at one point with all three Soret coefficients being zero at the same time, i.e. $S_{T,H_{2}O}=S_{T,ETH}=S_{T,TEG}=0$. This can be understood that when two Soret coefficients are zero, the third one must be zero due to the mass conversation at the steady state with $\mathrm{\nabla }{c}_A+\mathrm{\nabla }{c}_B+\mathrm{\nabla }{c}_C=0$, though this approximation may bring some error if there is a net vacancy flow due to the Kirkendall effect as discussed in Section 4.6.

Based on the discussion above, it is clear that to predict the Soret coefficients, one needs to accurately model the thermodynamic properties of the solutions as a function of temperature, volume, and compositions. While if one would like to simulate thermodiffusion processes, the mobilities of diffusion components also need to be modeled. A very interesting work was reported by Höglund and Ågren [87] who simulated the thermodiffusion in an Fe-32%Ni-0.14%C (weight percent, wt%) alloy using an experimental temperature profile reported by Okafor et al. [111] and the Dictra software package [112,113] with thermodynamic and mobility databases developed by the CALPHAD method [1,114–118]. In the lattice-fixed frame of reference, $L_{NiC}=L_{CNi}=0$ were adopted in accordance with Equation (24). In the number/volume-fixed frame of reference, they derived the cross-coefficient as follows (see Figure 3): (68) $\frac{L_{kT}^{{}^{\mathrm{\prime }}}}{T^2}=\sum _{i=1}^n({\delta }_{ki}-u_k{a}_i)\frac{u_i}{V_S}{y}_{Va}{M}_{iVa}{Q}_i^{\ast }\frac{1}{T}$(68) where $k$ represents $Fe$, $Ni$, or $C$; the Kronecker delta ${\delta }_{kk}=1$ and ${\delta }_{ki}=0$ for $k\ne i$; $a_{Fe}=a_{Ni}=1$ for substitutional elements and $a_C$ = 0 for interstitial elements; $V_S$ is the molar volume per mole of substitutional atoms; $u_i$ is the mole fraction per mole of substitutional atoms; $y_{Va}$ is the mole fraction of vacant lattice sites; $M_{iVa}$ is the mobility of component $i$ through vacancy mechanism (see Figure 3); and $Q_i^{\ast }$ is the heat of transport for component $i$ commonly used in the literature.

As the mobilities of the substitutional atoms are many orders of magnitude lower than those of the interstitials, $Q_{Fe}^{\ast }=Q_{Ni}^{\ast }=0$ are assumed, and Equation (68) is approximated as follows for $C$: (69) $\frac{L_{CT}^{{}^{\mathrm{\prime }}}}{T^2}=\frac{u_C}{V_S}{y}_{Va}{M}_{CVa}{Q}_C^{\ast }\frac{1}{T}$(69)

The flux of $C$ is obtained as follows: (70) $J_C^{{}^{\mathrm{\prime }}}=-D_{CC}^{Fe{}^{\mathrm{\prime }}}\frac{1}{V_S}\left(\mathrm{\nabla }{u}_C+Q_C^{\ast }\frac{1}{Td{\mu }_C/d{u}_C}\mathrm{\nabla }T\right)$(70) where $D_{CC}^{F{e}^{\mathrm{\prime }}}$ is the interdiffusion coefficient of $C$ (see Figure 3), and the full derivative means that the effects of other components are included. Under the steady-state condition for a closed system with $J_C^{\mathrm{\prime }}=0$, $Q_C^{\ast }$ can be obtained as: (71) $Q_C^{\ast }=-T\frac{d{\mu }_C}{d{u}_C}\frac{\mathrm{\nabla }{u}_C}{\mathrm{\nabla }T}=-T\frac{d{\mu }_C}{d{u}_C}\frac{d{u}_C}{dT}=\frac{1}{T}\frac{d{\mu }_C}{d{u}_C}\frac{d{u}_C}{d(1/T)}$(71) where $\mathrm{\nabla }{u}_C/\mathrm{\nabla }T$ = $d{u}_C/dT$ is used though the physical meaning of the derivative is not very clear. Using Equation (56), one can obtain the following equation for the Soret coefficient: (72) $S_{T,C}=-\frac{\mathrm{\nabla }{u}_C}{\mathrm{\nabla }T}=\frac{Q_C^{\ast }}{T}/\frac{d{\mu }_C}{d{u}_C}$(72)

On the other hand, $Q_C^{\ast }$ can be obtained from the definition given in Equations (1) and (2) in the publication by Höglund and Ågren [87] in combination with Equation (24) in the present paper with $\mathrm{\nabla }V=0$ as follows: (73) $Q_C^{\ast }=T\frac{\mathrm{\partial }{\mu }_C}{\mathrm{\partial }T}=-T{S}_C$(73) Höglund and Ågren [87] assumed a constant $Q_C^{\ast }$ and found that $Q_C^{\ast }=-44\text{kJ/mol}$ can best reproduce the $C$ concentration profile after 102 hours reported by Okafor et al. [111], who assumed that the steady state were reached after 102 hours and obtained $Q_C^{\ast }=-12.3\text{ kJ/mol}$. The negative value of $Q_C^{\ast }$ is in agreement with Equation (73) as $S_C$ is positive. The simulated concentration profiles are shown in Figure 6 for 102 and 13,889 (marked infinite) hours superimposed with experimental data by Okafor et al. [111] after 102 hours. It is evident that the steady state was not reached after 102 hours. Furthermore, even with $Q_{Fe}^{\ast }=Q_{Ni}^{\ast }=0$, both Fe and Ni should diffuse due to the $C$ concentration profile. More detailed simulation results will be presented in a future publication.

Figure 6. Simulated $C$ concentration profiles for 102 and 13,889 hours [87] with experimental data from ref [111] after 102 hours superimposed.

4.6. Uphill diffusion: a cross phenomenon between diffusion components

Diffusion in multicomponent systems under the isothermal condition is important for materials processing and joining of dissimilar materials [119–121]. Following the discussion on thermodiffusion, let us consider an isolated single-phase multicomponent system with an initial inhomogeneous distribution of some diffusion components and a homogeneous temperature profile. The inhomogeneous composition results in variations of chemical potentials which drive the diffusion of individual components from high chemical potential regions to low chemical potential regions. One scenario worth further discussion is when one component diffuses much faster than other components and thus takes much short time to reach the chemical equilibrium than other slow diffusion components. In the case that the chemical potential of the fast diffusion component is significantly affected by a slow diffusion component, the concentration of the fast diffusion component can become more inhomogeneous than before, i.e. so-called uphill diffusion where the fast diffusion component migrates against its concentration gradient.

Uphill diffusion was postulated by Darken [122] who pointed out that ‘for a system of more than two components it is no longer necessarily true that a given element tends to diffuse toward a region of lower concentration even within a single phase region', due to non-ideal interactions in the solution phase that is so great that the concentration gradient and the chemical potential gradient may be of different signs [123]. Darken further performed experimental validation of his postulations through four weld-diffusion experiments with pairs of steel of virtually the same carbon content, but differing markedly in alloy contents [124]. In two diffusion couples with the single face-centered-cubic (fcc) austenite solid solution phase, i.e. (I) Fe-3.8Si-0.25Mn-0.49C/Fe-0.05Si-0.88Mn-0.45C and (II) Fe-3.8Si-0.25Mn-0.49C/Fe-0.14Si-6.45Mn-0.58C (all in wt%), Darken observed that C diffusion from the low concentration region to the high concentration region at 1050°C after about two weeks. The diffusion couple II has contributions from both Si and Mn, while the diffusion couple I is primarily due to the Si concentration difference.

For the sake of simplicity, let us consider the diffusion couple I in the Fe-Si-C ternary system and write the flux of C as follows by omitting the effects of partial volume and partial entropy: (74) $\begin{array}{rl}{J}_C& =-L_C\mathrm{\nabla }{\mu }_C\\ & =-L_C({\mathrm{\Phi }}_{CC}\mathrm{\nabla }{c}_C+{\mathrm{\Phi }}_{CSi}\mathrm{\nabla }{c}_{Si}+{\mathrm{\Phi }}_{CFe}\mathrm{\nabla }{c}_{Fe})\end{array}$(74)

In usual situations, the last two terms are less important when $\mathrm{\nabla }{c}_{Si}$ and $\mathrm{\nabla }{c}_{Fe}$ are small, and the flux is controlled by the first term, resulting in $J_C$ and $\mathrm{\nabla }{c}_C$ being in opposite directions, i.e. the diffusion of C from high concentration to low concentration regions. However, in dissimilar metal welds, the gradients of other components can be rather steep [125,126]. In the Fe-Si-C diffusion couple I, the Si contents in two parts of the diffusion couple are very different. The chemical potential of C at 1050°C is significantly increased by the Si content as shown in Figure 7(i) calculated using a CALPHAD thermodynamic database and Thermo-Calc software [112,113]. Since the mobility of Si is orders of magnitude lower than that of C, the second term ${\mathrm{\Phi }}_{CSi}\mathrm{\nabla }{c}_{Si}$ in Equation (74) will drive C from high Si content region to low Si content region in order to reduce the chemical potential gradient of C. With a CALPHAD mobility database and Dictra software [112,113], the diffusion process in the diffusion couple I can be simulated [112], and the C concentration profile after 13 days is plotted in Figure 7(ii) in comparison with the experimental data by Darken [124], showing excellent agreement. Furthermore, Darken drew a schematic diagram showing the change in composition of two points on opposite sides of the weld in ultimately approaching uniformity of composition (see figure 6 in ref. [124]). Such a diagram for the diffusion couple I after 13 days is plotted in Figure 7(ii) in qualitative agreement with Darken’s schematic diagram though much steeper.

Figure 7. (i) Chemical potential of C in Fe-Si-0.45C alloys plotted with respect to the Si content at 1050°C with the TCFE9 thermodynamic database [113]; (ii) C composition profile in diffusion couple I after 13 days; (iii) C and Si compositions in the diffusion couple I with the numbers used to calculate the distance from the high Si side with the formula shown in the diagram.

On the other hand, Smoluchowski mentioned a study of a diffusion couple between Fe-0.80C and Fe-4.0Co-0.80C alloys [127] without change in the carbon distribution after several days at 1000°C. It is indeed that the chemical potential of C in the Fe-0.80C alloy is increased only slightly by the Co content at 1000°C as shown in Figure 8, in agreement with Darken’s suggestion [124,128]. Darken also noted that the uphill diffusion in liquid was observed more than a decade earlier by Hartley [129] and continues to be a very important phenomenon in a wide range of applications [130,131], such as transport of ionic species in aqueous solutions and crossing of distillation boundaries that are normally forbidden in multicomponent mixtures with azeotropes.

Figure 8. Chemical potential of C in the Fe-0.80C alloy is increased only slightly by the Co content at 1000°C, calculated using a CALPHAD thermodynamic database and Thermo-Calc software [112,113].

Another commonly observed uphill diffusion is when the solution is unstable with respect to composition fluctuations, referred to as spinodal decomposition [13,130,132,133]. The limit of stability of a solution is located at ${\mathrm{\Phi }}_{ii}=0$, and inside a spinodal the solution is unstable with ${\mathrm{\Phi }}_{ii}<0$, and the flux and concentration gradient of the component can thus point to the same direction as shown in many equations above, e.g. Equation (52). This uphill diffusion thus irreversibly transfers an initially homogenous unstable solution to an inhomogeneous stable solution and ultimately a miscibility gap of the same phase with two distinct composition sets. The instability-driven spinodal decompositions exist in both binary and multicomponent systems and play a central role in the formation of patterns [13] and dissipative structures [32].

An important phenomenon should be discussed before the end of this section, i.e. the Kirkendall effect [134–136]. During the time of Darken’s work, it was commonly believed that atomic diffusion in solid took place through direct exchange of atoms or ring mechanism with equal diffusion coefficients of binary elements. During that time period, Huntington and Seitz [137] suggested the vacancy mechanism based on evaluation of the activation energy for self-diffusion in copper using the electron theory. Kirkendall measured the diffusion coefficients of Zn in the Cu-Zn fcc solution ($\alpha$-brass) using diffusion couples between $\beta$-brass with 56.41–58.63wt%Zn (an intermetallic phase) and pure Cu at three temperatures [134] and concluded that Zn diffuses faster than Cu in $\alpha$-brass [135]. Subsequently, using single-phase diffusion couples between Cu-30wt%Zn $\alpha$-brass and pure Cu with insoluble thin Mo wires inserted in the bonding interface, Smigelskas and Kirkendall [136] observed the shift of Mo wires towards the original $\alpha$-brass, i.e. a shrinkage of the original $\alpha$-brass, indicating that Zn diffuses out of the high Zn $\alpha$-brass faster than Cu diffuses into the region. The Kirkendall effect and vacancy diffusion mechanism were accepted by the community a couple of years later [138]. This effect can be written in general as follows: (75) $J_{Va}=-\sum _{i\in S}{J}_i$(75) where $J_{Va}$ is the flux of vacancy, and $i\in S$ denotes all substitutional components. It shows that the unbalanced fluxes of substitutional components on the right-hand side of Equation (75) can result in a net vacancy flow and may induce the formation of voids.

One additional effect of isothermal diffusion processes is that the diffusion of components can in principle result in a heat conduction due to the partial entropy differences of components similar to the discussion in Section 4.2 represented by Equation (35) and can be written as: (76) $J_Q=T\sum _i{S}_i{J}_i$(76)

It is noted that the heat flux in diffusion of atoms is probably much smaller than that in electron conduction.

4.7. Electromigration: a cross phenomenon between atomic diffusion and electrical field

Electromigration concerns the cross phenomenon of atomic diffusion under an external electrical field and is the most serious and persistent reliability issue in interconnect metallization and flip chip solder joints in electronic devices [139–141]. Early research activities on electromigration was reviewed by Black [142] as a critical issue in electronic devices due to the formation of voids at interfaces between metals and semiconductors and the formation of metallic whiskers resulting in electrical short circuit. Systematical investigations on electromigration were reported for Al and Al-alloys [143–145] and Cu and Cu-alloys [146–148]. Particularly, Blech [144] investigated the electromigration of Al as a function of temperature, electrical current, and sample size and observed significant diffusion of Al at high temperature, high electrical current, and large sample size. Electromigration related to Pb-free solder joints has also been extensively investigated including effects of thermodiffusion, electrical currents, and stress on the formation of Kirkendall voids as discussed in Section 4.6 above [149–152].

For electromigration under an externally applied electrical field, the initial internal process is the electrical current, followed by heat conduction and internal stress induced by electrical field and electrical current, that ultimately induces atomic diffusion. The flux for electromigration can thus be written as follows based on Equation (20): (77) $\begin{array}{rl}{J}_i& =-L_i\mathrm{\nabla }{\mu }_i\\ & =-L_i\left(\sum _j{\mathrm{\Phi }}_{ij}\mathrm{\nabla }{c}_j-S_i\mathrm{\nabla }T-{\epsilon }_i\mathrm{\nabla }\sigma -{\theta }_i\mathrm{\nabla }E\right)\end{array}$(77) where partial strain ${\epsilon }_i=\frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }{c}_i}$ and partial electrical displacement ${\theta }_i=\frac{\mathrm{\partial }\theta }{\mathrm{\partial }{c}_i}$ are introduced with the Gibbs energy as the characteristic function [10,21,22]. The summation in Equation (77) includes a term for electrons due to the electrical current. It is noted that $L_i$, ${\mathrm{\Phi }}_{ij}$, ${\epsilon }_i$, and ${\theta }_i$ are functions of all internal variables shown in Equation (77) plus additionally microstructural features such as dislocations, grain boundaries, and surfaces, which become more pronounced with the size reduction of electronic devices [141].

Kirchheim [153] and Basaran et al. [154] considered the contributions to flux from the concentration gradient of the diffusion component, the stress gradient, and the electrical current plus a source term for the non-equilibrium vacancy concentration with a characteristic vacancy generation or annihilation time when applying Fick’s second law. One may consider that the electrical current in their model includes both the contributions of electron density and externally imposed electrical field expressed in Equation (77). It is noted that the diffusion of vacancy is related to the diffusion of atomic components, and the vacancy concentration depends on all variables in Equation (77) and is included in $L_i$ unless it is controlled externally as an independent variable such as radiation. The net vacancy flux is denoted by Equation (75). It is thus not clear whether the vacancy needs to be separately considered as shown by Kirchheim [153] and Basaran et al. [154] as the mechanisms for non-equilibrium vacancy and the parameter for vacancy generation or annihilation time are not clearly defined.

For pure metals such as Al and Cu, ${\theta }_i\approx 0$, and one may also assume homogeneous heating due to the electrical current with $\mathrm{\nabla }{c}_e\approx 0$, $\mathrm{\nabla }T\approx 0$ and $\mathrm{\nabla }\sigma \approx 0$. Equation (77) can thus be approximated as follows: (78) $J_A=-L_A\mathrm{\nabla }{\mu }_A=-L_A{\mathrm{\Phi }}_{Ae}\mathrm{\nabla }{c}_e$(78)

Therefore, Al or Cu diffuses in the same direction of the decrease of electron concentration, i.e. the direction of electron flow from cathode to anode, and the vacancy diffuses in the opposite direction ($J_{Va}=-J_A$, see Equation (75)) resulting in the formation void on the cathode side [139]. The heat generated by the electrical current can be calculated by Equation (76) with the summation including both the element and electron. For multicomponent systems with a temperature gradient, thermal expansion or contraction, the equations derived in Sections 4.5 and 4.6 can be utilized along with a mechanical energy term [153,154] and the thermodynamic and atomic mobility databases by the CALPHAD method [1,114–118].

4.8. Electrocaloric effect: a cross phenomenon between heat conduction and electrical field

The electrocaloric effect (ECE) concerns the cross phenomenon of heat conduction under an external electrical field, promising for both heating and cooling devices [155–161], in the opposite direction of pyroelectricity with voltage generation due to heating or cooling. It may be mentioned that the caloric effect can be realized additionally by external magnetic and mechanical fields [160], which will be briefly discussed in Section 4.10. It is similar to electromigration with the electrons as the only diffusion component. When an external electrical field is applied to the system, an electrical current is generated with electrons carrying heat with them, resulting in a heat conduction as follows from the general form of Equation (20): (79) $\begin{array}{rl}{J}_Q& =-L_Q\mathrm{\nabla }T\\ & =-L_Q\left(\frac{\mathrm{\partial }T}{\mathrm{\partial }S}\mathrm{\nabla }S+\sum _i\frac{\mathrm{\partial }T}{\mathrm{\partial }{c}_i}\mathrm{\nabla }{c}_i+\frac{\mathrm{\partial }T}{\mathrm{\partial }\sigma }\mathrm{\nabla }\sigma +\frac{\mathrm{\partial }T}{\mathrm{\partial }E}\mathrm{\nabla }E\right)\end{array}$(79)

The derivatives in Equation (79) are represented below: (80) $\frac{\mathrm{\partial }T}{\mathrm{\partial }S}=\frac{T}{C_P}$(80) (81) $\frac{\mathrm{\partial }T}{\mathrm{\partial }{c}_i}=\frac{{\mathrm{\partial }}^2\mathrm{\Phi }}{\mathrm{\partial }S\mathrm{\partial }{c}_i}=\frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }S}=\sum _j\frac{\mathrm{\partial }{\mu }_i/\mathrm{\partial }{c}_j}{\mathrm{\partial }S/\mathrm{\partial }{c}_j}=\sum _j\frac{{\mathrm{\Phi }}_{ij}}{S_j}$(81) (82) $\frac{\mathrm{\partial }T}{\mathrm{\partial }\sigma }=\frac{{\mathrm{\partial }}^2\mathrm{\Phi }}{\mathrm{\partial }S\mathrm{\partial }\sigma }=-\frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }S}=-\frac{1}{S_{\epsilon }}$(82) (83) $\frac{\mathrm{\partial }T}{\mathrm{\partial }E}=\frac{{\mathrm{\partial }}^2\mathrm{\Phi }}{\mathrm{\partial }S\mathrm{\partial }E}=-\frac{\mathrm{\partial }\theta }{\mathrm{\partial }S}=-\frac{1}{S_{\theta }}$(83) where $C_P$, $S_{\epsilon }=\frac{\mathrm{\partial }S}{\mathrm{\partial }\epsilon }$, and $S_{\theta }=\frac{\mathrm{\partial }S}{\mathrm{\partial }\theta }$ are isobaric heat capacity under constant stress, and partial entropy with respect to strain and electrical displacement, respectively. The discussions on electrocaloric effects in the literature usually assume $\mathrm{\nabla }{c}_i=0$ and $\mathrm{\nabla }\sigma =0$, and Equation (79) becomes: (84) $J_Q=-L_Q\left(\frac{T}{C_P}\mathrm{\nabla }S-\frac{1}{S_{\theta }}\mathrm{\nabla }E\right)$(84)

More recently, Chen et al. [22] discussed the thermodynamics of the direct and indirect methods used in characterizing the ECE materials. In the direct method, the temperature increase under the adiabatic condition as a function of electrical field is measured, i.e. $(\mathrm{\Delta }T/\mathrm{\Delta }E{)}_{\mathrm{\Delta }Q=0}$. In the literature, the adiabatic condition with $\mathrm{\Delta }Q=0$ is often treated as identical to the isentropic condition with $\mathrm{\Delta }S=0$, which is true for a closed, equilibrium system, but not true for closed, non-equilibrium systems with internal processes as shown by Equation (1). The differential form of $(\mathrm{\Delta }T/\mathrm{\Delta }E{)}_S$ is shown by Equation (83) with the temperature change obtained through integration: (85) $\mathit{\Delta }{T}_{direct}=\underset{E_1}{\overset{E_2}{\int }}\frac{\mathrm{\partial }T}{\mathrm{\partial }E}dE=-\underset{E_1}{\overset{E_2}{\int }}\frac{1}{S_{\theta }}dE$(85)

In the indirect method, the heat production under the isothermal condition as a function of electrical field is measured, i.e. $(\mathrm{\Delta }Q/\mathrm{\Delta }E{)}_T$, where $\mathrm{\Delta }Q$ is again usually approximated to be the same as $\mathrm{\Delta }S$. Its differential form can be obtained as follows: (86) $\frac{\mathrm{\partial }S}{\mathrm{\partial }E}=-\frac{{\mathrm{\partial }}^{2}G}{\mathrm{\partial }E\mathrm{\partial }T}=\frac{\mathrm{\partial }\theta }{\mathrm{\partial }T}$(86)

The entropy change in the indirect method can be calculated through integration of Equation (86) and is then used to calculate the anticipated temperature increase using the heat capacity of the materials as follows: (87) $\mathit{\Delta }S=\underset{E_1}{\overset{E_2}{\int }}\frac{\mathrm{\partial }S}{\mathrm{\partial }E}dE=\underset{E_1}{\overset{E_2}{\int }}\frac{\mathrm{\partial }\theta }{\mathrm{\partial }T}dE=\underset{T_1}{\overset{T_{2,indirect}}{\int }}\frac{C_P}{T}dT$(87) Consequently, $\mathrm{\Delta }{T}_{indirect}=T_{2,indirect}-T_1$ can be obtained. While $\mathrm{\Delta }{T}_{direct}$ from Equation (85) and $\mathrm{\Delta }S$ and $\mathrm{\Delta }{T}_{indirect}$ from Equation (87) can all be used to characterize the performance of ECE materials, $\mathrm{\Delta }{T}_{direct}$ and $\mathrm{\Delta }{T}_{indirect}$ are likely different from each other since both methods start from the same state, i.e. $(T_1,S_1,E_1)$, but end at different states, i.e. $(T_{2,direct},E_2{)}_{S_1}$ vs $(S_{2,indirect},E_2{)}_{T_1}$. Therefore, it is in general that $\mathrm{\Delta }{T}_{indirect}\ne \mathrm{\Delta }{T}_{direct}$.

In both methods, the integration is important because the materials properties change nonlinearly with respect to the electrical field. This is particularly true near critical points or phase transition boundaries commonly referred as morphotropic phase boundaries (MPBs). Unfortunately, as shown by Chen et al. [22] (noting the sign difference of their Equation (1)), Equation (87) is usually approximated in the literature by the following equation (88) $T_1\mathrm{\Delta }S=C_P\mathrm{\Delta }{T}_{indirect}$(88)

This can thus introduce large errors due to the dramatic temperature dependence of $C_P$ near MPBs where ECE is mostly investigated. Chen et al. [22] also pointed out the importance of domains and domain interactions in microstructure as also discussed in the literature [158,162–165] and will be further examined in Section 5.

4.9. Electromechanical effect: a cross phenomenon between electrical current and stress

The electromechanical effect concerns the interactions between electrical field and elastic deformation usually under isothermal condition without atomic diffusion [166–170]. Giant electromechanical effects have been observed near MPBs [171–174]. In piezoelectricity, electrical charge is accumulated in response to externally applied elastic deformation. While in the converse piezoelectric effect, the externally applied electrical field creates elastic deformation in the system.

In piezoelectricity, the electrical current under an external applied strain or stress can be obtained from Equation (34), respectively. (89) $\begin{array}{rl}{J}_e& =-L_e\mathrm{\nabla }{\mu }_e=-L_e\left(\frac{\mathrm{\partial }{\mu }_e}{\mathrm{\partial }{c}_e}\mathrm{\nabla }{c}_e+\frac{\mathrm{\partial }{\mu }_e}{\mathrm{\partial }\epsilon }\mathrm{\nabla }\epsilon \right)\\ & =-L_e({\mathrm{\Phi }}_{ee}\mathrm{\nabla }{c}_e+{\sigma }_e\mathrm{\nabla }\epsilon )\end{array}$(89) (90) $\begin{array}{rl}{J}_e& =-L_e\mathrm{\nabla }{\mu }_e=-L_e({\mathrm{\Phi }}_{ee}\mathrm{\nabla }{c}_e+{\epsilon }_e\mathrm{\nabla }\sigma )\end{array}$(90) where ${\mathrm{\Phi }}_{ee}$ is the thermodynamic factor and can be calculated as shown in Section 4.4, and ${\sigma }_e=\frac{\mathrm{\partial }{\mu }_e}{\mathrm{\partial }\epsilon }=\frac{\mathrm{\partial }\sigma }{\mathrm{\partial }{c}_e}$ and ${\epsilon }_e=\frac{\mathrm{\partial }{\mu }_e}{\mathrm{\partial }\sigma }=-\frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }{c}_e}$ are the partial stress and strain with respect to electron concentration, respectively. Under the condition of zero electrical current, the voltage generated is as follows under an external applied strain or stress, similar to Equation (39), respectively. (91) $\begin{array}{rl}\frac{\mathrm{\nabla }{V}_e}{\mathrm{\nabla }\epsilon }& ={\sigma }_e\end{array}$(91) (92) $\begin{array}{rl}\frac{\mathrm{\nabla }{V}_e}{\mathrm{\nabla }\sigma }& ={\epsilon }_e\end{array}$(92)

In the converse piezoelectric effect, the elastic deformation is induced by an externally applied electrical field. The mechanical equilibrium can be defined by vanishing stress gradient as follows: (93) $\begin{array}{rl}\mathrm{\nabla }\sigma & =\frac{\mathrm{\partial }\sigma }{\mathrm{\partial }\epsilon }\mathrm{\nabla }\epsilon +\frac{\mathrm{\partial }\sigma }{\mathrm{\partial }{c}_e}\mathrm{\nabla }{c}_e+\frac{\mathrm{\partial }\sigma }{\mathrm{\partial }E}\mathrm{\nabla }E\\ & =c\mathrm{\nabla }\epsilon +{\sigma }_e\mathrm{\nabla }{c}_e-{\theta }_{\epsilon }\mathrm{\nabla }E=0\end{array}$(93) where $c$ represents the elastic coefficient tensor, and ${\theta }_{\epsilon }=\frac{\mathrm{\partial }\theta }{\mathrm{\partial }\epsilon }=-\frac{\mathrm{\partial }\sigma }{\mathrm{\partial }E}$ is the partial electrical displacement with respect to strain.

4.10. General discussion of cross phenomena

In the previous sections, the general flux equation denoted by Equation (20) is applied to cross phenomena on transport of electrons and atomic components induced by externally controlled temperature, electrical, stress, and strain fields. As shown by Equation (20) and the discussions above, the cross phenomena can be represented by the derivatives of one potential to its conjugate molar quantity, other potentials, and other non-conjugate molar quantities, which are the second derivatives of free energy to its natural variables. The choice of free energy depends on how the system is controlled from the surroundings with relevant terms subtracted from the right-hand side of the general free energy formula shown by Equation (10), i.e. the ensemble of the system as mentioned in Section 2. The derivatives between potential and molar quantities represent well known physical properties in the literature and are shown in Table 1 which were discussed in the recent overview by the author [1], except the last row which is for derivatives related to chemical reactivity of components due to external stimuli, $\mathrm{\partial }{N}_j/\mathrm{\partial }{Y}^a$, or other components for the last cell in the row, $\mathrm{\partial }{N}_j/\mathrm{\partial }{N}_i$. As discussed by the author [1], the diagonal properties in Table 1 are the derivatives of potentials with respect to their respective conjugate molar quantities $\left(\frac{\mathrm{\partial }{Y}^a}{\mathrm{\partial }{X}^a}\right)$ and are all positive for stable systems. The limit of stability is reached when they become zero $\left(\frac{\mathrm{\partial }{Y}^a}{\mathrm{\partial }{X}^a}=+0\right)$ and their inverses, i.e. the derivatives of molar quantities with respect to their respective conjugate potentials, diverge and become positive infinite, i.e. (94) $\frac{\mathrm{\partial }{X}^a}{\mathrm{\partial }{Y}^a}=+\mathrm{\infty }$(94)

Table 1. Physical quantities related to the first directives of molar quantities (first column) to potentials (first row), symmetric due to the Maxwell relations [1,11,345].

	T, Temperature	$\sigma$, Stress	$E$, Electrical field	$\mathcal{H}$, Magnetic field	${\mu }_i$, Chemical potential
S, Entropy	Heat capacity	Piezocaloric effect	Electrocaloric effect	Magnetocaloric effect	${\displaystyle \frac{\mathrm{\partial }S}{\mathrm{\partial }{\mu }_k}}$
$\epsilon$, Strain	Thermal expansion	Elastic compliance	Converse piezoelectricity	Piezomagnetic moduli	${\displaystyle \frac{\mathrm{\partial }{\epsilon }_{ij}}{\mathrm{\partial }{\mu }_k}}$
$\theta$, Electrical displacement	Pyroelectric coefficients	Piezoelectric moduli	Permittivity	Magnetoelectric coefficient	${\displaystyle \frac{\mathrm{\partial }{D}_i}{\mathrm{\partial }{\mu }_k}}$
$B$, Magnetic induction	Pyromagnetic coefficient	Piezomagnetic moduli	Magnetoelectric coefficient	Permeability	${\displaystyle \frac{\mathrm{\partial }{B}_i}{\mathrm{\partial }{\mu }_k}}$
$N_j$, Moles	Thermoreactivity	Stressoreactivity	Electroreactivity	Magnetoreactivity	${\displaystyle \frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }{\mu }_k}}$, Thermodynamic factor

On the other hand, the off-diagonal properties can be either positive or negative, so are their divergencies (95) ${\left(\frac{\mathrm{\partial }{X}^a}{\mathrm{\partial }{Y}^b}\right)}_{a\ne b}=\pm \mathrm{\infty }$(95)

This was discussed recently in terms of the zentropy theory [20] and will be further presented in Section 5. It should be mentioned that every derivative between potential and molar quantity in Table 1 denotes an internal process that represents how the system responds to the external stimulus and must result in entropy production for it to happen irreversibly as dictated by $T{d}_{ip}S>0$ in Equation (1). This includes thermal expansion, $\mathrm{\partial }V/\mathrm{\partial }T$ for the volumetric flow with respect to the external temperature change [20] that will be discussed in Section 5.3.

On the other hand, the derivatives between any two potentials are presented in Table 2 for all the cross-phenomenon coefficients. The last row of Table 2 denotes the derivatives between one chemical potential and other potentials and can represent either transport cross phenomena or chemical reactions or their combinations. For transport cross phenomena due to externally applied potentials, they include the thermodiffusion and electromigration discussed in the present work and stressmigration [175–177] and magnetomigration [178] discussed in the literature. The same quantities can be applied to electrons to result in thermoelectricity and electrochemical and electrocaloric effects discussed in the present work and the magnetochemical effect [179–182] in the literature. While for chemical reactions, they represent whether the formation of a component is enhanced or lessened by the external potentials as the increase or decrease of chemical potential of a components indicates the increase of decrease of its moles in a stable system. One example of chemical reactions is a quantum reaction for the formation of quasiparticles from electrons in quantum materials that will be discussed in Section 5.5.

Table 2. Cross phenomenon coefficients represented by derivatives between potentials [23].

	$T$, Temperature	$\sigma$, Stress	$E$, Electrical field	$\mathcal{H}$, Magnetic field	${\mu }_i$, Chemical potential
$T$	$1$	$-{\displaystyle \frac{\mathrm{\partial }S}{\mathrm{\partial }\epsilon }}$	$-{\displaystyle \frac{\mathrm{\partial }S}{\mathrm{\partial }\theta }}$	$-{\displaystyle \frac{\mathrm{\partial }S}{\mathrm{\partial }B}}$	$-{\displaystyle \frac{\mathrm{\partial }S}{\mathrm{\partial }{c}_i}}$ Partial entropy
$\sigma$	${\displaystyle \frac{\mathrm{\partial }\sigma }{\mathrm{\partial }T}}$	$1$	$-{\displaystyle \frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }\theta }}$	$-{\displaystyle \frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }B}}$	$-{\displaystyle \frac{\mathrm{\partial }\epsilon }{\mathrm{\partial }{c}_i}}$ Partial strain
$E$	${\displaystyle \frac{\mathrm{\partial }E}{\mathrm{\partial }T}}$	${\displaystyle \frac{\mathrm{\partial }E}{\mathrm{\partial }\sigma }}$	$1$	$-{\displaystyle \frac{\mathrm{\partial }\theta }{\mathrm{\partial }B}}$	$-{\displaystyle \frac{\mathrm{\partial }\theta }{\mathrm{\partial }{c}_i}}$ Partial electrical displacement
$\mathcal{H}$	${\displaystyle \frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }T}}$	${\displaystyle \frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }\sigma }}$	${\displaystyle \frac{\mathrm{\partial }\mathcal{H}}{\mathrm{\partial }E}}$	$1$	$-{\displaystyle \frac{\mathrm{\partial }B}{\mathrm{\partial }{c}_i}}$ Partial magnetic induction
${\mu }_i$	${\displaystyle \frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }T}}$ Thermodiffusion	${\displaystyle \frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }\sigma }}$ Stressmigration	${\displaystyle \frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }E}}$ Electromigration	${\displaystyle \frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }\mathcal{H}}}$ Magnetomigration	${\displaystyle \frac{\mathrm{\partial }{\mu }_i}{\mathrm{\partial }{\mu }_j}}=-{\displaystyle \frac{\mathrm{\partial }{c}_j}{\mathrm{\partial }{c}_i}}={\displaystyle \frac{{\mathrm{\Phi }}_{ii}}{{\mathrm{\Phi }}_{ji}}}$ Crossdiffusion

By combining Tables 1 and 2, one can define many cross phenomena using both derivatives with respect to potential, molar quantity, or a mixture of them, such as the electromechanical effects with respect to externally applied stress, strain, or a mixture of strain and stress [21]. It is noted that even though the derivatives between potentials can be understood straightforwardly through the combined law of thermodynamics as demonstrated in the present paper, they are usually not considered useful and rarely discussed in thermodynamic textbooks. This is probably due to the deeply rooted belief that thermodynamics is for equilibrium only as engraved by the combined law derived by Gibbs [5,183,184]. Since each potential is homogeneous in an equilibrium system, the derivatives between potentials are indeed not very meaningful when only equilibrium is considered. However, as discussed above, they play a central role in understanding and predicting internal processes and cross phenomena as the bridge between thermodynamics and kinetics.

The general flux equation, Equation (20), can be applied to all irreversible processes inside a system [185–187] driven by any or all five sets of potentials, i.e. temperature, mechanical, electrical, magnetic, and chemical with one or more of them controlled externally from the surroundings in experiments. Both Tables 1 and 2 are symmetrical due to the Maxwell relations through the second derivatives of free energy to its natural variables. This is particularly valuable for quantities in Table 2 as the derivatives between potentials, i.e. the lower left portion of the table, are relatively easy to measure experimentally, while the derivatives between molar quantities, i.e. the upper right portion of the table, are ideal for theoretically investigations, which are the focus of Section 5.

5. Zentropy theory for prediction of entropy and cross phenomenon coefficients

5.1. General discussion on entropy, fluctuation, and critical phenomena

The importance of temperature in all kinetic processes is reflected by the temperature dependence of kinetic coefficients shown by Equation (33). The derivative of a free energy to temperature gives one of the most important quantities in science, i.e. entropy, which is the core of the second law of thermodynamics shown by Equation (1) in the present paper and can be fundamentally understood by the Boltzmann distribution as discussed in Section 4.2. Entropy is intrinsically multi-scale from its presence in black holes [188–190], societies [191], forests [192], and quantum systems [17,193,194].

Thermodynamic modeling for metallic and some oxide and molten salt systems has been progressed significantly in last half century in terms of the semi-empirical CALPHAD method [1,114–118] as a function of temperature, compositions, and spontaneous magnetization, while the thermodynamic modeling with respect to ferroelectric polarization and stress/strain has been largely based on the Landau-Ginsburg-Devonshire (LGD) formalism [195,196] led by Cross and co-workers [197–199]. All those models are very useful for practical applications, particularly the extrapolations of the CALPHAD databases from binary and ternary systems to multicomponent systems for computational design of engineering materials [200–202] and the LGD formalism for ferroelectric materials [21]. However, those models are phenomenological, and their model parameters are evaluated from experimental data though DFT-based first-principles calculations have played more and more important roles recently [118,203].

While the DFT-based first-principles calculations can provide reliable free energy of a given configuration under the quasiharmonic approximation [50,51,204,205], including the recent stochastic self-consistent harmonic approximation [206], the prediction of anharmonicity remains challenging [207,208]. The state-of-the-art approaches are the ab initio molecular dynamic simulations (AIMD) for high temperatures [209] and the quantum Monte Carlo (QMC) simulations for low temperatures [210]. In addition to their highly computational intensity, they are not capable to systematically predict the critical point where a system changes from being stable to unstable, and the anharmonicity is at its maximum due to the divergency of molar quantities, including the thermodiffusion near critical point [211] and quantum critical point (QCP) [212] with the divergence of the quasiparticle effective mass (QEM) [213–216] which will be further discussed in Section 5.3. In Section 4.10, the concept of instability was introduced in terms of the divergence of the derivative of a molar quantity to its conjugate potential. A system usually changes from stable to metastable before reaching the limit of stability and becoming unstable, except at the critical point where the metastability region vanishes and the system switches from being stable to unstable directly.

The critical phenomena can be understood by the renormalization group (RG) theory in terms of statistical continuum limit with the fluctuations of many length or energy scales cascaded in a Hamiltonian that enters the partition function of the system [217–219], and the challenge is shifted to evaluating the Hamiltonians for various cases. Another challenge in implementing the RG theory is the explicit incorporation of multiscale interactions in the Hamiltonian. Consequently, predictive approaches are still lacking for critical phenomena at low temperatures for quantum critical points [210] and at high temperatures where Lennard-Jones potentials are usually used [220,221] or by fitted equation of state [222]. The more recently developed fluctuation theorem (FT) [15–17,223–226] also concerns the fluctuations in nonequilibrium steady-state systems, which may be considered as an extension of the equilibrium fluctuation often conveniently assumed in the RG theory. The development of FT started with the discussion that the 2^nd law of thermodynamics has a statistical probability to be ‘violated' in the individual microscopic trajectories of those fluctuations though not macroscopically with all the microscopic trajectories combined [15]. There have been a number of experimental and theoretical studies followed [225,227–230]. As discussed below, this is probably because an additional entropy contribution due to the probability of the microscopic trajectory is not accounted for [224]. This additional entropy contribution has been related to the information or internal statistical microstates of fluctuations [229,231,232] and will further be discussed in Section 5.4.

The lack of quantification of free energy of ferroelectric materials was also mentioned by Scott [158] who pointed out that there is no known analytic equation of state for most ferroelectric materials, and the theories based on the LGD formalism do not include ferroelectric domains even though much of the entropy change involves the creation and destruction of domains [162–165]. Scott further noted that in the LGD approximation, one typically assumes that all coefficients except one are independent of temperature, resulting in the entropy change at a first-order phase transition equal to the square of the polarization change timed by 2π. In other words, the Lorentz factor comes from a quadratic free energy expansion in polarization, which is rarely valid.

Even though the LGD formalism describes the ferroelectric (FE) and paraelectric (PE) behaviors well macroscopically, the modern microscopic theory of polarization (MTP) [233–236] is needed to understand thermodynamics of ferroelectric polarization. In the MTP framework, effective Hamiltonians have been constructed by representing the energy surface in terms of a Taylor expansion around the high-symmetry cubic structure along with several approximations for computation efficiency. Their model parameters were fitted to DFT-based first-principles calculations followed by Monte Carlo or molecular dynamics simulations [162,164,237–239]. Significant nonzero local polarization are predicted in the cubic PE state for BaTiO₃ [162,237] and PbTiO₃ [238] in agreement with the experimental observations in the literature [240–242] and the orthorhombic-like 90° domain structure in PbTiO₃ in agreements with experiments [164] and AIMD simulations [165]. These results further demonstrate the important contributions of fluctuations and domain walls to free energy of materials, which is discussed in the next section.

5.2. Zentropy theory for multiscale entropy

Entropy can be presented in terms of the integration of heat capacity (see Equation (87)), which is convenient in experiments, or Boltzmann distribution (see discussions in Sections 4.2 and 4.4), which is convenient in theory and computations. In the literature, one often differentiates the thermal electronic, phonon, and configurational entropies as they are at different scales, governed by different physics, and calculated by different methods [50,204,243]. Nevertheless, they are all configurational at their respective scales and distinctive physics. The distribution of electrons is represented by the Fermi-Dirac statistics as shown by Equation (47), while the phonon distribution is characterized by the Bose–Einstein statistics for indistinguishable particles. For classic distinguishable particles, the Maxwell-Boltzmann statistics applies.

Hierarchically, a macroscopic system of investigation (termed as macrostate in the present paper) can be considered as a statistical mixture of various microscopic configurations of the system (termed as microstates in the present paper), while each microstate can further be considered as a statistical mixture of various sub-microscopic configurations of the system. The scale higher than the macrostate is treated as the surroundings of the system [10,11], which dictates the statistical ensemble to be used to study the system as mentioned in Section 2 with a proper free energy represented by Equation (10). The configurational entropy of the macrostate can be represented by Gibbs-Boltzmann entropy as follows: (96) $S^{conf}=-k_B\sum _k{p}_{k}ln{p}_k$(96) where $p_k$ is the probability of microstate $k$ in the macrostate.

It is important to realize that in this hierarchical multiscale framework, each microstate and sub-microstate must fill the complete space of the system and experience the same constraints of the statistical ensemble of the system. It is thus evident that the total entropy of the macrostate needs to include the entropy of each microstate in addition to the configurational entropy among them as follows [14]: (97) $\begin{array}{rl}S& =\sum _k{p}_k{S}^k+S^{conf}=\sum _k{p}_k(S^k-k_{B}ln{p}_k)\\ & =S_{0K}+\underset{0}{\overset{T}{\int }}\frac{C_h}{T}dT\end{array}$(97) where $S^k$ is the entropy of microstate $k$ which can be further decomposed into its sub-microstates with the same formula as Equation (97), and $S_{0K}$ and $C_h$ are the entropy at $0\text{K}$ and the heat capacity of the macrostate. The superscript is used in $S^k$ as the subscript is already used for partial entropy in the present paper. The integration from the experimentally measured heat capacity in Equation (97) is usually obtained under the convention of the third law of thermodynamics, i.e. $S_{0K}=0$ at $0\text{K}$, which is also reflected by DFT stating that there is only one ground state at 0 K fully defined by its electron density distribution[8,9]. Nevertheless, contributions to $S_{0K}$ could be included if needed and available.

From Equation (97), the free energy and partition function of the macrostate and the probability of each microstate can be obtained as follows: (98) $\begin{array}{rl}\mathrm{\Phi }& =\sum _k{p}_k{\mathrm{\Phi }}^k+k_{B}T\sum _k{p}_k\ln p_k\end{array}$(98) (99) $\begin{array}{rl}Z& =e^{-\frac{\mathrm{\Phi }}{k_{B}T}}=\sum _k{e}^{-\frac{{\mathrm{\Phi }}^k}{k_{B}T}}=\sum _k{Z}_k\end{array}$(99) (100) $\begin{array}{rl}{p}_k& =\frac{Z_k}{Z}=\frac{Z_k}{{\sum }_k{Z}_k}=e^{-\frac{{\mathrm{\Phi }}^k-\mathrm{\Phi }}{k_{B}T}}\end{array}$(100)

The above equations represent a nested formula for many length or energy scales, resembling the RG theory discussed above though without explicit considerations of interactions between microstates, thus considerably simplifying the calculations. This is possible only if the microstates are ergodic, i.e. the interactions between microstates in the macrostate are represented internally in one or more microstates, similar to the ergodic dynamical systems considered in FT [16]. When the ergodicity is not met, an additional mean-field term needs to be added as shown by two approaches in predicting the critical point in $Ce$ [244,245]. It is noted that the free energy is used in the partition function in Equation (99) rather than the total energy used in the partition function shown by Landau and Lifshitz [246]. In their derivation, the quantum states are treated as the microstates with their entropy implicitly assumed to be zero, i.e. $S^k=0$, resulting in the configurational entropy among the microstates (Equation (96)) representing the total entropy of the system with the quasi-classical approximation. The detailed discussion on this comparison will be presented in a separate publication.

We have recently termed the nested formula, i.e. Equation (96) to Equation (100), as zentropy theory in an overview of our previous studies [20]. It is shown that the zentropy theory is capable of predicting critical points in magnetic materials and associated anomaly including the positive and negative divergency of volume with respect to temperature, i.e. positive and negative thermal expansions with all inputs from DFT-based quasiharmonic calculations without empirical parameters (see Section 5.3). The significances of the zentropy theory, which includes Equation (96) to Equation (100), are as follows:

1. The nested form of Equation (97) covers any scales of investigations through the combination of all lower scales and the scale of interest represented by $S^k$ and $S^{conf}$ and can thus in principle describe any systems of interest from the quantum scale to black holes.

2. In the MSE discipline, the macrostates of interest are microstructures consisting of individual phases and their temporal and spatial arrangements. To start with, one can take each individual phase as a macrostate of investigation and define its microstates in terms of atomic, magnetic, electrical, and defect configurations. The entropy of each microstate, $S^k$, includes the contributions from thermal electrons and phonons that can be predicted by DFT-based first-principles calculations using the Fermi-Dirac and the Bose–Einstein statistics, respectively [50,247].

3. $S^{conf}$ represents the statistical fluctuation of microstates and is the key to predict the critical phenomena and anharmonicity at high temperatures. It is partially represented by phenomenological models such as the LGD formalism [195,196] and the CALPHAD compound energy formalism [248].

4. The zentropy theory integrates the bottom-up approach based on quantum mechanics ($S^k$ in Equation (97)) and the top-down approach based on statistical mechanics ($S^{conf}$ in Equation (97)), i.e. (ii) and (iii) above, and connects the combination of the quantum mechanics and statistical mechanics to the macroscopic thermodynamics ($C_h$ and its integration in Equation (97)). Such characteristics of the zentropy theory preserve the quantum information from the DFT electron level to the macroscopic level of experimental investigations and enable the prediction of the nonlinear, emergent macroscopic behaviors of macroscopic functionalities at finite temperatures, including their ultimate extreme at critical points [20] and potential applicability to superconductivity to be discussed in Section 5.5.

5.3. Applications of zentropy theory to magnetic and ferroelectric transitions

It is mentioned above that thermal expansion is one type of cross phenomenon as it concerns the volumetric flow driven by temperature, though it is not commonly related to cross phenomena in the literature. The volumetric flow with respect to temperature is represented by the third term in Equation (20), i.e. $\frac{\mathrm{\partial }{Y}_j^a}{\mathrm{\partial }{X}_j^{b\ne a}}\mathrm{\nabla }{X}_j^{b\ne a}$ with $Y_j^a$ for pressure and $X_j^{b\ne a}$ for entropy. The inverse of the corresponding Maxwell relation is: (101) $\frac{\mathrm{\partial }S}{\mathrm{\partial }(-P)}=\frac{\mathrm{\partial }V}{\mathrm{\partial }T}$(101)

We applied the zentropy theory to magnetic transitions in anti-invar $Ce$ and invar $F{e}_{3}Pt$ with the microstates defined by magnetic spin configurations [20]. The thermal electronic and phonon entropic contributions for each magnetic microstate are predicted from the DFT-based first-principles calculations. Since the transition in $Ce$ is between nonmagnetic (NM) ground-state microstate and ferromagnetic (FM) high temperature non-ground-state microstate, we started with these two microstates, but had to add a mean-field magnetic flipping term in free energy as commonly done in the literature [244] in order to obtain the critical point.

It is known that the co-existence of the NM and FM microstates results in spin-flipping magnetic (SFM) microstates, generating domain walls. We subsequently added the antiferromagnetic (AFM) microstate and predicted the critical point and associated anomalies from the decrease of probability of the NM microstate and the increase of probabilities of the AFM and FM microstates without the mean-field term [1,245,249]. The predicted temperature-pressure potential phase diagram using the zentropy theory is show in Figure 9(i) with a critical point and the two-phase region represented by the line in agreement with the Gibbs phase rule [10,11]. Various experimental data are well reproduced with details discussed in the early publication [20,245].

Figure 9. Predicted phase diagrams of $Ce$ (i) temperature-pressure [245] with symbols for experimental data and (ii) temperature-volume [249] with purple diamond squares for thermal expansion anomaly and other symbols for experimental data.

To show the divergency of volume at the critical point, the pressure axis is replaced by its conjugate molar quantity, i.e. volume, resulting in a temperature-volume potential-molar phase diagram plotted in Figure 9(ii). It is noted that the two-phase region changes from a line in the temperature-pressure potential phase diagram to an area in the temperature-volume potential-molar phase diagram. This does not violate the Gibbs phase rule as the Gibbs phase rule concerns only the number of independent potentials without changing the number of phases in equilibrium and needs to be modified if applied to phase diagrams including molar quantities [10]. Several isobaric volume curves are plotted in Figure 9(ii) for $Ce$, depicting the thermal expansion anomaly marked by the pink diamond symbols. It can be seen that at the critical point marked by the green circle $\frac{\mathrm{\partial }T}{\mathrm{\partial }V}=0$ and $\frac{\mathrm{\partial }V}{\mathrm{\partial }T}=+\mathrm{\infty }$ in accordance with Equation (95). There is a temperature range on each isobaric volume curve that the thermal expansion shows abnormally larger than usual. The further is it away from the critical point, the less is the anomaly. In both Figure 9(i) and (ii), the agreement between experiments and predictions is remarkable with all inputs from DFT-based first-principles calculations and no empirical parameters.

While for $F{e}_{3}Pt$, the magnetic transition is between FM and paramagnetic (PM) macrostates, and many more SFM microstates are needed. We found that the 12 atom supercell with 9 $Fe$ atoms is large enough to predict experimentally observed anomalies [1,249,250]. This supercell results in $2^9=512$ magnetic microstates with 37 being unique due to symmetry. With their free energies predicted from DFT-based first-principles calculations including both thermal electronic and phonon contributions, the predicted temperature-pressure [250] is calculated using the zentropy theory and shown in Figure 10(i) with a critical point marked by the red circle in good agreement with experimental data. As in the case of $Ce$, the temperature-volume phase diagram of $F{e}_{3}Pt$ is plotted in Figure 10(ii) with several isobaric volume curves. The filled black circles on the ambient pressure (0 GPa) volume curve are experimental data with details discussed in the early publication [249], showing remarkable agreement with predicted values, including the negative thermal expansion temperature range. The purple diamonds mark the anomalous regions of negative thermal expansions, including the negative divergences at the critical point marked by the green circle, i.e. $\frac{\mathrm{\partial }T}{\mathrm{\partial }V}=0$ and $\frac{\mathrm{\partial }V}{\mathrm{\partial }T}=-\mathrm{\infty }$ in accordance with Equation (95). It was concluded that the negative thermal expansion originates from the decreasing probability of the ground-state FM microstate and the increasing probabilities of those non-ground-state SFM microstates with their volumes smaller than that of the FM microstate.

Figure 10. Predicted phases diagrams of $F{e}_{3}Pt$ (i) temperature-pressure [250] with experimental data superimposed, noting pressure decrease from left to right, and (ii) temperature-volume [249] with experimental data (black circles) on the ambient pressure (0 GPa) volume curve and the purple diamonds for the anomalous regions of negative thermal expansions.

We have applied the zentropy theory to the AFM-PM transitions in several rare-earth nicklates [251] and the FE-PE transitions in ferroelectric materials [165,252] and observed remarkable agreement between predicted and measured transition temperatures. Those results are being prepared for publications.

Another interesting experimental phenomenon related to the limit of stability is the divergence of QEM at a QCP [213–216] as mentioned in Section 5.1 with the examples of $YbR{h}_{2}S{i}_2$ [213], $YbR{h}_2(S{i}_{0.95}G{e}_{0.05}{)}_2$ [214], Al-doped CrAs [253], and $FeS{e}_{1-x}{S}_x$ [254,255] shown in Figure 11. Monthoux et al. [256] pointed out that quasiparticles reduce to bare electrons in the hypothetical absence of the electron-electron and electron-ion interactions. In conductors, a quasiparticle may be considered as an electron plus a co-moving screening cloud, while in superconductors, the quasiparticles can have effective masses two or more orders of magnitude greater than that of electrons and hence move relatively slowly [256]. According to itinerant spin fluctuation theory [212], the QEM divergence at a QCP is due to the abundance of low-lying and long-range spin fluctuations, mediating the interactions between the heavy quasiparticles, giving rise to pronounced deviations from Landau Fermi liquid behavior [213–216]. This is exactly how the zentropy theory works.

Figure 11. QCP phase diagrams: (i) temperature with respect to magnetic field of $YbR{h}_{2}S{i}_2$ with antiferromagnetic (AF), non-Fermi liquid (NFL), and Landau Fermi liquid (LFL) phase regions [213], (ii) temperature with respect to magnetic field of $YbR{h}_{2}S{i}_2$ and $YbR{h}_2(S{i}_{0.95}G{e}_{0.05}{)}_2$ with AF (blue left), NFL (yellow), and LFL (blue, right) phase regions [214], (iii) temperature with respect to pressure with contour map of electrical resistivity of Al-doped CrAs and pure CrAs [253], and (iv) temperature with respect to composition of $FeS{e}_{1-x}{S}_x$ with contour maps for electrical resistivity (ρ) and a parameter (A*) related to quasiparticle effective mass [255].

Let us discuss further the divergency of QEM at a QCP in terms of the limit of stability in Section 4.10 and the present section for $Ce$ and $F{e}_{3}Pt$ with more details presented in the literature [1,10,11]. One should keep in mind that even though a QCP is theoretically defined as a point at $0K$ where a microstate becomes unstable to new forms of microstates, experimental measurements are always performed at finite temperatures with thermal fluctuations. Let us consider a simple internal process by adding one electron to an existing quasiparticle represented by the following reaction: (102) $q_{n-1}+e=q_n$(102)

The combined law of thermodynamics can be re-written from Equation (3) as follows with $d\xi =dn$ and its driving force per quasiparticle as $D_n=-({\mu }_{q_n}-{\mu }_{q_{n-1}}-{\mu }_e)=-\mathrm{\Delta }{\mu }_{q_n}$: (103) $dU=Y^{a}d{X}^a+\mathrm{\Delta }{\mu }_{q_n}dn$(103) where ${\mu }_e$, ${\mu }_{q_n}$, ${\mu }_{q_{n-1}}$ are the chemical potentials of electron and quasiparticle with $n$ and $n-1$ electrons, respectively, and $Y^a$ and $X^a$ are the conjugate variables representing the external controlling variables, i.e. pressure, magnetic field, or composition (chemical potential) of dopants as shown in Figure 11. The limit of stability of the system in terms of Equation (95) can be written as: (104) $\frac{\mathrm{\partial }n}{\mathrm{\partial }{Y}^a}=\pm \mathrm{\infty }$(104) where the positive and negative signs represent the change of $Y^a$ or $X^a$ towards QCP or away from QCP as the driving force $D_n$ changes sign accordingly. To make the model more realistic, one needs to consider a distribution of quasiparticle sizes determined by thermodynamic equilibrium with a distribution of electrons as shown in Section 4.4. When $n$ diverges, its effective mass, $m_{q_n}^{\ast }\approx n{m}_e^0$ for relatively slow-moving quasiparticles [256] with $m_e^0$ being the resting mass of electron, also diverges as shown in Figure 11.

5.4. On microscopic ‘violation of second law of thermodynamics’

As mentioned above, the development of fluctuation theorem started with the discussion by Evans et al. [15] on probability of ‘violations of second law of thermodynamics' in a nonequilibrium steady state far from equilibrium. They derived an expression for the ratio of the probabilities to find a fluid with an induced shear stress in the direction of or opposite to the externally imposed shear rate on a phase space trajectory segment of a certain duration in a dynamical state. Through MD simulations, they observed statistically significant probability that the induced shear stress is in the opposite direction of the externally imposed shear rate, resulting in negative entropy production evaluated from the induced shear stress, thus an apparent ‘violation of the second law of thermodynamics' for those events or internal processes.

Independently, Jarzynski [224] derived an expression for the equilibrium free energy difference between two states of a system, $\mathrm{\Delta }F$, in terms of an ensemble of nonequilibrium (finite-time) measurements with the work, $W$, as follows: (105) $\mathrm{\Delta }F=-\frac{1}{k_{B}T}\ln \left(\overline{e^{-\frac{W}{k_{B}T}}}\right)$(105)

The work is computed for each trajectory in the ensemble with the quasiequilibrium statistical mechanics invoked. This is remarkable as $\mathrm{\Delta }F$ represents an equilibrium property independent of path, while $W$ is realized by trajectories with finite rate and (106) $\overline{W}\ge \mathrm{\Delta }F$(106)

The equality in Equation (106) holds for infinitely slow internal processes between two states, i.e. the conventionally defined reversible processes. Jarzynski [224] also showed the applicability of the above equations to the cases where the system of interest and its surroundings together constitute a larger, isolated Hamiltonian system.

The validity of Equations (105) and (106) does not prevent some trajectories with $W<\mathrm{\Delta }F$, and many theoretical, computational, and experimental investigations have demonstrated the existence of such trajectories [15,225,227–230]. Particularly, using a single-electron transistor two-level system with $\mathrm{\Delta }F=0$ and a single thermodynamic trajectory level coupled to a single heat bath, Maillet et al. [229] experimentally demonstrated that the amount of work up to large fractions of $k_{B}T$ can be extracted from the two carefully designed out-of-equilibrium driving cycles with probabilities substantially greater than 1/2. These observations result in an apparent ‘violation of the second law of thermodynamics' for those individual trajectories with the external intervention and energy consumption to drive and select those trajectories. By including the degrees of freedom of external intervention such as a feedback controller, Sagawa and Ueda [232] concluded the consistency with the conventional second law of thermodynamics because of the energy cost of the controller. However, this does not resolve the issue related to those trajectories with $W<\mathrm{\Delta }F$ as internal processes because second law of thermodynamics is not about the total entropy change, but the entropy production of each independent internal process.

Therefore, the interpretations of results from the theory, computations, and experiments in the framework of fluctuation theorem are worth another look, particularly with each individual trajectory as an internal process. We will approach this in terms of the relative magnitude of $Q_{{\xi }_j}$ and $k_{B}T$ discussed in Section 4.1. It is noted that Jarzynski [224] emphasized that $W$ must not be much greater than $k_{B}T$ in order to verify Equation (105) experimentally. If the fluctuations in $W$ between trajectories are much larger than $k_{B}T$ (Cases II and III discussed in Section 4.1), the ensemble average in Equation (105) will be dominated by trajectories with much lower $W$ (Case IV discussed in Section 4.1). The situation here is similar to the reversible Brownian motion in a macroscopically equilibrium system, which is used for the prediction of diffusivity in terms of DFT-based first-principles calculations discussed in Section 4.3 [57–64], with $\mathrm{\Delta }F=0$.

Let us follow each individual atom in a reversible Brownian motion and consider its every motion as an internal process or trajectory. When the atom sits still at its equilibrium position, such as at 0 K, both the work and entropy are zero. Next, when the atom starts to fluctuate from its equilibrium position at low temperatures, the entropy production of the internal process is represented by the phonon dispersion and the associated vibrational entropy calculated from the probability of each vibrational frequency that the atom can possess through integration over all frequencies and probabilities through the Bose-Einstein statistics. For conductors, thermal electronic entropy can be included through Fermi-Dirac statistics. When the fluctuation is large enough, the probability for the atom to move to the next vacant site becomes significant. This internal process now generates a new configurational entropy due to the probability for the atom to jump over to the previous vacant site or fall back to the new vacant site that it leaves behind. This new entropy contribution can be represented by Equation (96) as: (107) $S^{conf}=-k_B(p_{+}\ln p_{+}+p_{-}\ln p_{-})$(107) with $p_{+}$ and $p_{-}$ being the probability jumping forward and backward, respectively. This results in an additional positive entropy production with its maximum value as $k_B\ln 2$ with $p_{+}=p_{-}=0.5$, assuming that the two microstates are identical with the same entropy. The entropy production due to the Brownian motion can thus be written from Equation (97) as (108) ${\mathrm{\Delta }}_{ip}S=(S^{+}-S^{-})+k_B\ln 2=k_B\ln 2$(108) with $S^{+}$ and $S^{-}$ being the entropy of two microstates, respectively. Consequently, the second law of thermodynamics is maintained as the internal process of a single trajectory generates a positive entropy production.

The work along the trajectory for every successful jump of a Brownian motion is performed by the thermal energy. Part of the work could in principles be extracted by an external intervention, being a fraction of $T{S}^{conf}=k_{B}T\ln 2$. While in the experiments of the single-electron transistor two-level system by Maillet et al. [229], the electron is manipulated by external force, and they did point out the Shannon entropy difference between the equilibrium microstates before and after the quench, which is related to $S^{conf}$. The fluctuations can be considered in similar fashion for nonequilibrium steady state systems such as various cross phenomena discussed in Section 4 by dividing the nonequilibrium system into a succession of regions with small differences of local temperature and concentration [15]. Furthermore, the probability of a successful jump can be used to evaluate the diffusivity in the DFT-based first-principles calculations discussed in Section 4.3 [57–64]. It seems that in most discussions related to the ‘violation of second law of thermodynamics', $S^{conf}$ is not included, resulting in incomplete accounting of entropy production in a single trajectory.

Let us discuss further the entropy production of internal processes in general. In our previous work [14], it was shown that the entropy production of an internal process can be divided into four parts: (1) heat generation $(d_{ip}Q)$, (2) consumption of some nutrients $(d{N}_i^n)$, (3) production of some wastes $(d{N}_j^w),$ and (4) reorganization of its configurations $(d_{ip}{S}^{config})$, as follows: (109) $\begin{array}{rl}{d}_{ip}S& =\frac{d_{ip}Q}{T}-\sum _i{S}_i^{n}d{N}_i^n+\sum _j{S}_j^{w}d{N}_j^w+d_{ip}{S}^{config}\end{array}$(109) where $S_i^n$ and $S_j^w$ are the partial entropies of nutrient $i$ and waste $j$ of the internal process, respectively. For individual trajectories considered here, $d{N}_i^n=d{N}_j^w=0$ with no mass exchange, $d_{ip}Q=0$ without heat exchange between the internal process and its surroundings, and $d_{ip}{S}^{config}$ is obtained from Equation (97) as follows (110) $\begin{array}{rl}{d}_{ip}S& =d_{ip}{S}^{config}=\sum _{k}d(p_k{S}^k)-k_B\sum _{k}d(p_k\ln p_k)\end{array}$(110)

The first term in Equation (110) is due to the thermal electronic and vibrational contributions of each microstate. The second term makes a positive contribution to the entropy production of each trajectory due to introduction of new microstates and is the key to demonstrate the validity of the second law of thermodynamics in both microscopical and macroscopical scales. The fluctuations and their probabilities of individual trajectories are related to the relative magnitude of barrier and $k_{B}T$ and nested in our zentropy theory in terms of Equation (96) to Equation (100). It is further noted that the scales are relatives as mentioned at the beginning of Section 5.1, and a macrostate in one investigation can be a microstate in a larger scale of investigation [14].

5.5. Postulations on superconductivity as a cross phenomenon with criticality

Superconductivity is a phenomenon discovered in 1911 where electrical resistance of a material vanishes at temperatures below a critical temperature $T_C$, and correspondingly, its electrical conductivity diverges. It is a critical phenomenon with Schottky anomaly and a critical point between the 1st and 2nd order transitions [257]. The transition temperatures can be altered through composition, pressure, and magnetic field (see Figure 11). The Bardeen–Cooper–Schrieffer (BCS) theory [258,259] was developed in 1957 as the microscopic theory of superconductivity in terms of the condensation of Cooper pairs. This paired state of electrons is originated from the electron–phonon interaction and has a lower energy than the Fermi energy. Due to the weak pairing interaction $(\sim {10}^{-3}\text{eV})$, thermal energy $(k_{B}T$ with $k_B$ = $8.617{10}^{-5}\text{eV/K}$) can easily break the pairs, resulting in conventional superconductors only at low temperatures. The BCS theory can explain the conventional low $T_C$ superconductors, but fails for high $T_C$ superconductors [260], and the early predictions based on the phonon-induced attraction indicated a ceiling to $T_C$ below the boiling point of liquid nitrogen [256,261].

The current approaches used in the community to predict $T_C$ for conventional BCS superconductors is based on Eliashberg’s model for interactions between electrons and lattice vibrations [262]. McMillan [261] and Allen and Dynes [263] made further developments, with the latter pointing out that Eliashberg theory appears to place no bounds on achievable values of $T_C$ determined by two key parameters with one for the attractive effects of the electron-phonon coupling and another for the repulsive effects of direct Coulomb interactions between electrons. Ashcroft [264,265] predicted the superconductivity of metallic hydrogen either as a monatomic or paired metal using the following mean field formula in terms of the BCS theory: (111) $T_C=0.85{\mathrm{\Theta }}_D{e}^{-1/n({\epsilon }_F){\varphi }_{el-ph}}$(111) where ${\mathrm{\Theta }}_D$ is the Debye temperature derived from the highest-frequency vibrational mode in the system, $n({\epsilon }_F)$ the e-DOS at the Fermi level (see Equation (46)), and ${\varphi }_{el-ph}$ an effective electron-phonon attractive interaction [266,267]. This mean field formula is also used in predicting of $T_C$ of various hydrides [268–272] with strong anharmonicity included [273].

Stewart [274] compared the fundamental properties of 9 unconventional superconducting classes of materials including 4f-electron heavy fermions, cuprates, Fe-based superconductors, organic superconductors, and new emerging topological superconductors with the hope that common features would emerge to help theory explain and predict these phenomena. It was concluded that ‘a common understanding - in the manner of a BCS theory - is not possible' with the chief question being ‘what interaction is coupling the electrons?' Monthoux et al.[256] also pointed out the importance of understanding of effective charge-charge interaction between quasiparticles and its dependences on phonons or spins. Typical superconducting phase diagrams are presented in Figure 12 for (i) cooperates [275], (ii) $Ce{M}_2{X}_2$ [256], (iii) $Ba(F{e}_{1-x}{M}_x{)}_{2}A{s}_2$ [276], and (iv) carbon- and hydrogen-doped $H_{3}S$ [277], showing that the transitions between insulator, semiconductor, conductor, and superconductor of a material can be controlled through pressure and composition (doping).

Figure 12. Superconducting phase diagrams, (i) cuprates [275], (ii) $Ce{M}_2{X}_2$ [256], (iii) $Ba(F{e}_{1-x}{M}_x{)}_{2}A{s}_2$ [276], and (iv) carbon- and hydrogen-doped $H_{3}S$ [277].

The author’s interest on superconductivity started when working on processing of the $Mg{B}_2$ superconductor [278,279]. When the original concept of the zentropy theory was developed in 2008 [244], the author’s group explored the applicability of the zentropy theory to challenges facing the scientific community. After some literature search, five topics were identified: (1) thermoelectric materials, (2) multiferroic materials, (3) coexistence of conductivity and transparency, (4) superconductivity, and (5) metal-insulator transition. Five people collected relevant information in the literature, gave short presentations, and led the lively discussions during a retreat held on April 12, 2008. As discussed in Section 4.4, we have made good progress in predicting Seebeck coefficient [68,69] for topic 1. For topic 2, we are able to predict magnetic transitions [1,245,249,250] and are working on manuscripts applying the zentropy theory for more complex magnetic materials and ferroelectric materials [252]. We have not made any progresses on topics 3. We have been actively working on topics 4 and 5 as they are closely related as mentioned above.

Right after the retreat, the author wrote five postulations on superconductivity, which are re-organized as follows

1. Postulation One: Superconductivity at 0 K. A conductor could be made to be superconducting at 0 K through proper external constrains such as deformation or doping.

2. Postulation Two: Low $T_C$ superconductors. For crystals with high entropy due to phonons, the electron-phonon interactions will interfere with the migration of electron pairs at finite temperatures. When temperature is high enough, the above interactions will break the Cooper pairs, and superconductivity is lost. This group of materials exhibits low $T_C$ and follows the BCS theory, including $Mg{B}_2$.

3. Postulation Three: High $T_C$ superconductors. The layered structures in complex phases with or without magnetism significantly reduce the entropy due to phonons, thus preserving the Cooper pairs within the layers. There can be two scenarios here: (i) the entropy due to phonons becomes dominant with the increase of temperature, and the superconductivity is lost before magnetic disordering; (ii) as the magnetic disordering lifts the constraints on the entropy due to phonons, the superconductivity is lost at the same time of the magnetic disordering temperature.

4. Postulation Four: Doping Effects. Some possible effects of doping are: (a) increase (decrease) the number of electron pairs for higher (lower) $T_C$; (b) stiff (soften) the lattice and reduce (increase) the entropy due to phonons for higher (lower) $T_C$; (c) enhance (diminish) the magnetic strength for higher (lower) $T_C$.

5. Postulation Five: Transition between superconductor and non-superconductor. Similar to SFM microstates proposed for Invar materials such as $F{e}_{3}Pt$ [249,250] (which was submitted for publication in March 2008), the transition between superconductor and non-superconductor can be proposed as follows: (a) there are two extreme states in terms of superconductivity: superconducting state or non-superconducting state; (b) there can be many states between these two states; (c) the thermal populations of various states vary as a function of temperature, pressure, or magnetic field; (d) various mixed states can be represented through the Boltzmann distributions expressed in terms of partition functions (now called zentropy theory); (e) the transition between the superconductor and non-superconductor can be predicted using the zentropy theory in terms of Equation (98) to Equation (100); (f) the critical transition point with respect to magnetic field should also be reflected at 0 K. Therefore, one should be able to see a transition between the superconductor and non-superconductor on the temperature-pressure diagram at 0K and its extension to higher temperatures with respect to pressure (similar to Figures 9 and 10). Similarly, one should be able to see a transition in the temperature-magnetic field diagram at 0 K, which extends to higher temperatures with the decrease of the magnetic field. Therefore, 0 K stability should be examined first in detail; (g) The key challenge is to define the superconducting and various non-superconducting microstates.

Along these intuitive postulations, the author’s team applied the zentropy theory to predict the properties of $BaF{e}_{2}A{s}_2$ [280–283]. $BaF{e}_{2}A{s}_2$ is a parent compound of iron arsenide superconductors [284] and is orthorhombic and tetragonal at low and high temperatures, respectively. It was recently made into a superconductor through epitaxial superlattice heterostructures [285]. The ground state of $BaF{e}_{2}A{s}_2$ is a spin density wave (SDW) AFM microstate with a stripe-like Fe spin ordering pattern within the ab plane and antiparallel nearest-neighbor Fe spins along the c axis [280–283]. Figure 13(i) plots the predicted SDW ordering temperature ($T_{SDW}$) and the characteristic temperature ($T_{\ast }$) as a function of pressure using the zentropy theory with experimental measurements superimposed, showing remarkable agreement [282]. $T_{SDW}$ and $T_{\ast }$ are defined when the probability of the SDW-AFM microstate equals to 0.5 and 0.9999, respectively. The Fermi surfaces of the SDW and Stripe microstates are plotted in Figure 13(ii). The Stripe microstate is similar to the SDW microstate except that the nearest-neighbor Fe spins along the c axis are parallel [281,282]. The blue and red colors in Figure 13(ii) mark the single sheet hole-pocket and single sheet electron-pocket, respectively, with G, Z, Y, and T denoting the high symmetry points in the reciprocal space. The shape of the predicted Fermi surface near the Z points of the SDW microstate agree with the measurement reported in the literature [286].

Figure 13. Orthorhombic-$BaF{e}_{2}A{s}_2$ (i) SDW ordering temperature ($T_{SDW}$, red curve) and characteristic temperature ($T_{\ast }$, blue curve) plotted with respect to pressure with experimental data (symbols) [282], and (ii) Fermi surface, Isosurfaces: (a) Stripe and (b) SDW; Cuts: (c) Stripe and (d) SDW; Top views of isosurfaces: (e) Stripe and (f) SDW [281].

As demonstrated above, the essential component of the zentropy theory is to assemble the microstates at 0 K that are ergodic for the investigation and predict their free energies by DFT-based calculations. While we should continue to use our scientific intuitions to guide the design of important microstates, the author believes that the integration of domain knowledge and deep neural network machine learning approaches [287–289] has the potential to develop robust approaches to assemble the ergodic microstates for complex systems including superconductors aided by high throughput DFT-based calculations [290,291] and high throughput thermodynamic modeling [292–295]. The power of such integrated approach is demonstrated in our recent work in determining the crystal structure of $NdB{i}_2$ [296]. First, all the AB₂-type configurations (26,055) were identified within a large dataset of DFT-relaxed or experimental structures (∼1.3 millions at the time of extraction, now over 4 millions) [289]. Following the substitution with $NdB{i}_2$ into the microstates, all of the candidates were examined by our deep neural network machine learning models (SIPFENN: structure-informed prediction of formation energy using neural networks) [287,288] to identify configurations with low formation enthalpy, resulting in 500 AB₂-type configurations with low values up to 12 kJ/mol-atom variation. With initial DFT-based verifications, 20 configurations with the lowest values were selected for more detailed DFT-based calculations to obtain the energy versus volume equation of state and its stable crystal structure in accordance with the experimental structure data by X-ray diffractions [296].

The insulating and conducting microstates can be defined by their band gaps, which can be predicted by the DFT-based calculations with continuously improved accuracy and efficiency [68,297,298], while the superconducting microstate is often discussed in connection with the Fermi surface [299] (see Figure 13(ii)), but a succinct differentiation between superconducting and normal conducting microstates is lacking. Considering that the electronic properties in the DFT-based calculations are represented by electronic band structure, density of states, charge density, and Fermi surface [281,300–302], their combination may uniquely define the superconducting microstate of a superconductor at 0 K.

Among various superconductors studied, two are highlighted here, i.e. an optimally doped $HgB{a}_{2}Cu{O}_{4+d}$ (Hg-1201) [303] and $FeS{e}_{1-x}{S}_x$ [304] as shown in Figure 14. For Hg-1201 at lowest pressures, the behavior is metallic with a sharp superconducting transition and a well-defined tetragonal structure with a fluctuating 1D linear oxygen ordering along the two tetragonal axes. As pressure increases, $T_C$ increases slightly up to 5 GPa and then decreases with a widened superconducting transition, reaching 0 K at 11.5 GPa (see blue open circles Figure 14(i)). At higher pressures, it becomes an insulator. At 8 GPa, the superstructure diffraction spots start to appear and become more intense and diffuse at higher pressures (see green filled squares in Figure 14(i)). Since no extra diffuse lines develop in the a–b plane, the ordering is proposed to be along the c axis, resulting in a 3D ordering of the oxygen. Various phase regions are depicted in Figure 14(i) as a function of hole concentration in Hg-1201, including the superconducting phase (SC in the figure). These ordered structures can be considered as the microstates in Hg-1201 and will be used to study the transitions between superconductor, conductor and insulator in terms of the zentropy theory.

Figure 14. Phase diagrams of (i) $T_C$ and normalized intensity of the superstructure spots as a function of pressure in Hg-1201 [303], (ii) $T$ plotted against hole concentration in Hg-1201 with hole concentration for (i) marked by the vertical rectangle and various phase regions [303], (iii) $T$ plotted against pressure in $FeS{e}_{1-x}{S}_x$ [304] with various phase regions: SDW (see Figure 13(i)), Tetra (tetragonal), Ortho (orthorhombic), M (magnetic), Nematic (no magnetic long-range), and SC (superconducting).

A large number of experimental [305–315] and computational [316–318] investigations has been reported for $FeSe$ due to its unique properties as an exceptional member of the family of Fe-based superconductors [309] with complex phase diagrams involving magnetic, superconducting, and structural transitions. For example, the effects of $S$ [305] are shown in Figure 14(iii) with similar ones for $Te$ [307]. It is hoped that the zentropy theory will be able to help to develop better understanding of the superconducting transition through the quantum fluctuations of superconducting and non-superconducting microstates and aid the discovery of new superconductors.

5.6. Cross phenomena involving phase transformations with multiple internal processes

When commenting on Darken’s work [124], Bhadeshia [319] generalized Darken’s statement to ‘for a system of more than two components, it is no longer necessarily true that a given element tends to diffuse toward a region of lower chemical potential', which happens when some elements are trapped by a moving interface between two phases [320]. When a local equilibrium is assumed at the phase interface, the chemical potential of each element is the same in both phases at the interface, and the interface migration is controlled by diffusion in each phase from high chemical potential regions to low chemical potential regions in that phase. When the phase interfaces deviate from local equilibrium, the chemical potentials of some elements in the new phase increase, and those of other elements decrease with respect to their chemical potentials in the parent phase. At the extreme, some or all elements are completely trapped from the parent phase into the new phase and have the same compositions in both phases, resulting in their conjugate potentials being different in the two phases. The phase transformation is thus termed paraequilibrium for complete trapping of substitutional elements or partitionless for complete trapping of all elements, respectively [321,322]. It is important to point out that when independent variables are changed from potentials to conjugate molar quantities, the Gibbs phase rule needs to be modified as shown by Hillert [10]. If there are additional constraints to a system such as the coherency of phase interface, the Gibbs phase rule needs to be further revised [323].

The intermediate scenarios between local equilibrium, paraequilibrium, and partitionless transitions are determined by the interplay of thermodynamic and kinetic properties of the systems, resulting in a wide range of dissipation of energy and dissipative structures [13,32,322,324,325]. As depicted in Figure 15 for Gibbs energy diagram of an alloy system A-B, Ågren [326] showed that the chemical potential of element B in the growing α phase, ${\mu }_B^{\alpha /\beta }$, is higher than that in the parent β phase, ${\mu }_B^{\beta /\alpha }$, due to the diffusion over the interface between α and β that requires a driving force, and the chemical potential for local equilibrium is between these two values (not drawn). The opposite is true for the chemical potential of element A. This is called solute drag, a highly important topic for phase transformations and properties of multicomponent materials [322,326–333]. Therefore, it seems that the element B migrates from high chemical potential region to low chemical potential region across the interface, due to the overall reduction of free energy as pointed out by Bhadeshia [319].

Figure 15. Gibbs energy diagram of an alloy system A-B with the effect of a driving force for diffusion acting over an interface between α and β, and the tangents describing the chemical potentials on both sides of the interface rotated relative to each other [326].

This can be further related to the coupling of atomic diffusion when crossing the phase interface. When multiple internal processes occur simultaneously at the same location, i.e. the migration of interface and the diffusion across the interface in the present case, they are coupled, and a new driving force needs to be defined to describe the coupled outcome as shown by Equation (5). For example, in the case of paraequilibrium, since the diffusion of all substitutional elements is prevented due to the much faster migration of interface than the atomic diffusion rates, the substitutional elements are no longer independent, but act collectively. Their chemical potentials are combined into one chemical potential, ${\mu }_{\in S}$, as follows, resulting in a revised local equilibrium condition[321]: (112) ${\mu }_{\in S}=\sum _{i\in S}{c}_i{\mu }_i^{\alpha }=\sum _{i\in S}{c}_i{\mu }_i^{\beta }$(112)

Figure 16. Dissolution of cementite at 910°C in an Fe-2.06Cr-3.91C (at. pct) alloy, (i) Isothermal section [120], (ii) Cr concentration profile after 1000 s [120], (iii) isothermal section with C activity plotted on the x-axis [334], (iv) transmission electron microscope (TEM) bright field image of cementite exhibiting Widmanstatten plates of α-bcc phase (γ-fcc at dissolution temperature) [334], and (v) TEM bright field image of extraction replica showing lamellar $M_7{C}_3$ structure ($u_{Cr}=0.54$ at point A) and untransformed $M{C}_3$ (darker area) [334].

In a study on dissolution of cementite ($M{C}_3$ with $M$ representing the mixture of Fe and Cr) at 910°C in an Fe-2.06Cr-3.91C (at. pct) alloy, the author found that the newly formed austenite inherited the Cr content of cementite due to the low diffusivity of Cr, while the high diffusivity of C resulted in a nearly instantaneous dissolution of cementite with its volume fraction decreased from 14.1 to 12.4% [120] in a fraction of second, most likely by paraequilibrium transformation. After the C activity reached homogenous in the system, the further dissolution of cementite was controlled by the slower diffusion of Cr, resulting in the local equilibrium at the interface. This local equilibrium condition increased the Cr concentration in $M{C}_3$ at the interface to a very high value as shown by the isothermal section in Figure 16(i) [120]. Consequently, Cr diffuses from the interface towards the center of $M{C}_3$ particles as shown by the concentration profile in Figure 16(ii) after 1000 seconds holding at 910°C [120]. Computer simulations under the local equilibrium conditions showed good agreement with experimental measurements in Figure 16(ii).

Furthermore, Figure 16(i) shows that the compositions of $M{C}_3$ near the interface fall initially into the $\gamma +M_7{C}_3+M_{23}{C}_6$ three-phase region and move into the $\gamma +M_7{C}_3$ two-phare region when the activity of C increases with more $M{C}_3$ dissolved. This can be seen more clearly with the C composition replaced by its chemical activity in the isothermal section shown in Figure 16(iii) [334]. It is evident that both $M{C}_3$ compositions at the center and near the interface, denoted by the horizontal dashed line and the top dotted line, respectively, are within the $\gamma +M_7{C}_3$ stable two-phase regions. Indeed, the $M{C}_3(\text{low}Cr)\to \gamma +M{C}_3(\text{high }Cr)$ transformation was observed inside $M{C}_3$ particles (Figure 16(iv)) [334], along with the decomposition of $M{C}_3$ into the lamellar structure of $\gamma +M_7{C}_3$ (Figure 16(v)) [334].

Figure 17. Classification schemes of kinetic processes in relation to fundamental components of materials science and engineering.

This brings up an interesting topic on classification of phase transitions [335,336]. The cross phenomena discussed in the present article becomes more complex due to multiple internal processes with moving interfaces as mentioned above. The simplest cases are grain boundaries in polycrystals and domain walls in ferromagnetic and ferroelectric materials. In polycrystals, particularly nanograins, both the grain growth and redistributions of elements between the grain interior and grain boundary need to be considered as two independent internal processes simultaneously [12,337]. In ferroelectrics, the polarization domains need to be included explicitly in thermodynamic models in order to predict the FE-PE transitions as pointed out by Scott [158] and discussed in Section 5.1.

It is noted that phase transformations are usually classified by three schemes in terms of thermodynamics, mechanism, and morphology [335,336]. They are related to the three inner circles in Figure 1(ii), replotted in Figure 17 with the thermodynamic scheme being at the interface between the thermodynamics and kinetics circles, the mechanistic scheme within the circle of kinetics, and the morphological scheme between the crystallography and kinetics circles. Each scheme thus focuses on difference aspects of a phase transformation [336], and when combined together they give a complete understanding of phase transformations.

In the thermodynamic scheme, phase transformations are classified by the derivatives of free energy with respect to potentials, which are their conjugate molar quantities. The order of a phase transition is assigned when the highest order of derivatives between two phases is discontinuous. For example, in a first-order transformation, the molar quantities of the two phases are different, while for the second-order transformation, all quantities in Table 1 are different between the two phases. The first- and second-order transformations meet each other at a critical point as shown in Figures 9 and 10. Higher-order transformations are less discussed in the literature and can be important for smaller scale phase transitions that may include the superconducting transition mentioned above [338–341]. Furthermore, as pointed out by Ågren [336], one should use the molar quantity of the individual phases in defining the order of a transformation, not those of the system.

The microstructural scheme is divided into two major groups: homogenous and heterogeneous with the former taking place everywhere simultaneously and the latter starting from some regions and propagating through the system. While a heterogenous transformation is certainly first-order in the thermodynamic scheme, a homogeneous transformation can be either first- or second-order. Though simple to use, the microstructural scheme has less value as it depends on length- and time-scale of investigation and is qualitative in nature with limited fundamental understanding and practical usefulness.

The mechanistic scheme reflects the core activities in the MSE discipline, i.e. to understand the fundamental mechanisms of phase transformations so they can be controlled through chemistry and external stimuli to generate desired microstructures and properties. At the same time, it is a very complex one due to multiple, simultaneous internal processes with cross effects discussed in the present article. Furthermore, since phase transformations take place in three-dimensional (3D) space, while investigations are usually performed in two dimensions (2D), it is important to keep in mind that observations may be incomplete. For example, a pearlite colony in steels looks like that it consists of alternating lamellae of ferrite and cementite in 2D which were thought to be formed by repeated nucleation and growth, while 3D sectioning showed that a pearlite consists of only one ferrite crystal and one cementite crystal, and it is their cooperative growth with shortened diffusion distance that results in the 3D interwoven lamellar structure [342,343]. This demonstrates the importance to have a 3D perspective when investigating complex transformation mechanisms, which can be facilitated by 3D phase-field simulations [13,344].

6. Summary and outlooks

In the present article, the combined law of thermodynamics is presented without removing the contributions from irreversible internal processes and is thus applicable to both equilibrium and non-equilibrium systems. It is pointed out that the driving force for the flux of a molar quantity is the gradient of its conjugate potential, which enables the development of physics-based kinetic equations rather than the phenomenological ones proposed by Onsager. The kinetic coefficient matrix is thus diagonal and naturally fulfills the Onsager reciprocal relations for a symmetrical kinetic coefficient matrix. The dependence of the conjugate potential on other independent variables of the internal processes is manifested by the cross phenomena with the coupling coefficients being the derivatives of the potential to the other independent variables. The commonly studied cross phenomena are analyzed in terms of both fundamental understanding and DFT-based computational methodology through the Maxwell relations and our zentropy theory, including thermoelectricity, thermodiffusion, uphill diffusion, electromigration, electrocaloric and electromechanical effects, thermal expansion, and quantum criticality. The successful predictions of Seebeck coefficients in thermoelectric materials, critical phenomena with magnetic transitions and property divergence, and negative thermal expansion in Invar FePt without empirical parameters are reviewed. It is anticipated that this new theory of cross phenomena and the associated predictive computational approach will have the potential to enable more efficient discovery of materials with emergent behaviors such as superconductivity along with better understanding of their fundamentals.

Acknowledgements

The author feels privileged to work with all his current and former students and numerous collaborators over the years at Penn State and around the world as reflected in the references cited in this paper and previous overview paper [1]. He would like to particularly thank Yi Wang, Shun-Li Shang, and Jorge Sofo for simulating discussions on cross phenomena and superconductivity. The author would like to thank Johnn Ågren, Graeme Murch, and Irina Belova for many stimulating discussions on cross phenomena, John Mauro, Digby MacDonald, and Susan Sinnott for critical reading of the manuscript, and inspiring communications with Graeme Ackland, Alan Allnatt, Ping Ao, H.K.D.H. Bhadeshia, Dieter Braun, Nigel Goldenfeld, Martin Gruebele, Christopher Jarzynski, Werner Köhler, Manfred Martin, Piergiulio Tempesta, Anton Van Der Ven, George Verros, and Han-Ill Yoo.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The author is grateful for financial supports from many funding agencies in the United States as listed in the cited references and previous overview paper [1], plus some more recent ones including the National Science Foundation [NSF CMMI-2050069], the Department of Energy [DE-AR0001435, DE-NE0008945], and Office of Naval Research [N00014-21-1-2608]. The high-performance computing facilities include the Roar supercomputer at the Pennsylvania State University’s Institute for Computational and Data Sciences (ICDS), the National Energy Research Scientific Computing Center (NERSC) supported by the US Department of Energy Office of Science User Facility operated under [DE-AC02-05CH11231], and the Extreme Science and Engineering Discovery Environment (XSEDE) supported by NSF [ACI-1548562].