**5. Irreversible processes**

24 Will-be-set-by-IN-TECH

The temperature Green's function in Eq. (129) is very similar to Eq. (110), which includes the real single-particle energy spectrum. However, there is one important difference that the frequency is not a real frequency or energy, because *ωn* is a *discrete* variable. The computation

It is convenient to introduce a real-time Green's function corresponding to the real-time trace in Eq. (120) multiplied by the statistical grand distribution Eq. (121). The real-time Green's function is then written (Fetter & Walecka, 2003) in the momentum space representation as

Let us introduce a generic real-time Green's function of a complex variable *z* which satisfies

d*ω*� 2*π*

 ∞ −∞

<sup>1</sup> <sup>∓</sup> *<sup>e</sup>*−*βh*¯ *<sup>ω</sup> <sup>G</sup>*¯ *<sup>R</sup>* (**k**, *<sup>ω</sup>*) <sup>+</sup>

*<sup>i</sup>ω<sup>n</sup>* <sup>−</sup> *<sup>h</sup>*¯ <sup>−</sup><sup>1</sup>

*eiωn<sup>η</sup>*

− <sup>Σ</sup>*<sup>λ</sup>* (**k**, *<sup>ω</sup>n*)

<sup>1</sup> <sup>∓</sup> *<sup>e</sup>βh*¯ *<sup>ω</sup> <sup>G</sup>*¯ *<sup>A</sup>* (**k**, *<sup>ω</sup>*). (134)

*<sup>z</sup>* <sup>−</sup> *<sup>ω</sup>*� , (135)

*<sup>ω</sup>* <sup>−</sup> *<sup>ω</sup>*� <sup>+</sup> *<sup>i</sup><sup>η</sup>* , (136)

= 1 (138)

for |*ω*| → ∞. (139)

*<sup>i</sup>ω<sup>n</sup>* <sup>−</sup> *<sup>ω</sup>*� . (140)

. (137)

. (133)

*�*0 **<sup>k</sup>** − *μ* 

1

*A* (**k**, *ω*� )

) is the spectral weight function (Galitskii & Migdal, 1958), which yields the

d*ω*� 2*π*

d*ω*� 2*π*

*A* (**k**, *ω*� )

*A* (**k**, *ω*� )

*ω* − *ω*� − *iη*

 ∞ −∞

 ∞ −∞

 d3*k* (2*π*) 2 1 *<sup>β</sup>h*¯ <sup>∑</sup>*<sup>n</sup>*

of the energy spectra is available by using the temperature Green's function.

Γ¯ (**k**, *z*) =

*G*¯ *<sup>R</sup>* (**k**, *ω*) = Γ¯ (**k**, *ω* + *iη*) =

*<sup>G</sup>*¯ *<sup>A</sup>* (**k**, *<sup>ω</sup>*) <sup>=</sup> <sup>Γ</sup>¯ (**k**, *<sup>ω</sup>* <sup>−</sup> *<sup>i</sup>η*) <sup>=</sup>

 ∞ −∞ d*ω*� 2*π A* **k**, *ω*� 

for both bosons and fermions. With this sum rule, we have an asymptoctic behavior of the

 ∞ −∞

In other words, the function Γ¯ (**k**, *z*) in Eq. (135) determines the temperature Green's function <sup>G</sup> as well as the real-time Green's funcitons *<sup>G</sup>*¯ *<sup>R</sup>* and *<sup>G</sup>*¯ *<sup>A</sup>* by analytic continuation technique in the complex *ω* plane. The sum rule ensures the uniqueness of the analytic continuation (Baym & Mermin, 1947). Once we evaluate the temperature Green's function G only at the

d*ω*� 2*π*

*A* (**k**, *ω*� )

*<sup>G</sup>*¯ (**k**, *<sup>ω</sup>*) <sup>=</sup> <sup>1</sup>

integral representations of the *retarded* Green's function

The spectral weight function should satisfy the sum rule

<sup>Γ</sup>¯ (**k**, *<sup>ω</sup>*) <sup>∼</sup> <sup>1</sup>

*ω*

 ∞ −∞

The spectral function constructs also the temperature Green's function

G (**k**, *ωn*) =

d*ω*� 2*π A* **k**, *ω*� ∼ 1 *ω*

and the number of particles is

where *A* (**k**, *ω*�

Green's function

and the *advanced* Green's function

*N* (*T*, *V*, *μ*) = ∓*V* (2*s* + 1)

The Green's function method reviewed above describes the equilibrium thermodynamics of a matter composed of fermions with the interactions mediated by the corresponding bosons. The power of the Green's function method does not limited to the equilibrium thermodynamics, but does solve many irreversible thermodynamic problems for the situation not far from the equilbrium states. The information from the irreversible processes of not-so-far from the equilibrium is extremely important for building an equilibrium phase diagram.

The procedure for acquiring the thermodynamic information includes a certain perturbation to the system for equilibrium process. For example, one should apply the external pressure to the system for changing its volume, one should apply a magnetic field **H** to the system for measuring its magnetization, one should should add heat to the system for increasing its temperature, etc. Such realistic experimental procedure should be considered for comparing the theoretical and experimental thermodynamic data assessments. However, a certain kind of phases in the system does not stable under the given condition. Even the Gibbs free energies for some metastable/unstable structures are undefinable because they are dynamically unstatble. This kind of problems has been studies extensively in the CALPHAD society, see Kaufman (2001) for example.

Many difficulties in practical CALPHAD computations arise when the system has magnetic instabilities. The inclusion of magnetism in the CALPHAD method has been devoted for long time (Hillert & Jarl, 1978; Midownik, 1977). However, the current methods to build any phase diagram with magnetism are hardly thought to be completed, in terms of both the first-principles and the empirical one. The apparent reason for this failure is due to the nature of magnetism, which we have not understood its underlying physics yet.

Once magnetism involves, we are facing Invar effects (Kim, 1988; 1999), magnetostrictions (Lines, 1979), and magnetocaloric effects (de Oliveira & von Ranke, 2010). As stated in Sec. 3.2, there are spontaneous symmetry breaking related with elementary excitations, for example magnons or spin waves emerge when the spin rotation symmetry is broken. In a ferromagnetic system, the magnetic *susceptibility* which is defined by the magnetization with respect to the external magnetic field diverges below the Curie temperature *T*<sup>C</sup> without external magnetic fields. In general, the computation of the magnetic susceptibility should consider the correct underlying physical mechanism. The experimental procedure for measuring any susceptibilities is the irreversible processes. This is also in great importance to understand the phase transitions.

### **5.1 Linear response theory**

A theory for phase transition should be developed; the theory leaves the non-diverging error near the critical region, where a set of physical observables diverges. Kubo (1957) and Kubo et al. (1957) developed such a theory, called as the *linear reponse theory*, for the case of not-so-far from equilibrium, by the quantum mechanical interpretations on the Liouville' theorem in Eq. (17) and the Bloch equation in Eq. (122).

Let the Hamiltonian of the *natural motion* of a many-body system to be *H*ˆ . We would like to apply an external perturbation *H*ˆ � to the system, which yields now the total Hamiltonian

$$
\hat{H}\_{\text{total}} = \hat{H} + \hat{H}'.\tag{141}
$$

This differential equation is readily solved as

*ρ*ˆ

using Eqs. (144) and (150) to be

 *A*ˆ (*t*) <sup>=</sup> <sup>1</sup> *ih*¯ *t* −∞ d*t* � Tr

we have the linear response as

 *A*ˆ (*ω*) <sup>=</sup> <sup>±</sup> *<sup>i</sup> h*¯ ∞ 0

�

where *AI* (*t* − *t*

property of trace.

is then calculated as

� (*t*) <sup>=</sup> <sup>1</sup> *ih*¯ *e* <sup>−</sup>*i*K<sup>ˆ</sup> *<sup>t</sup>*/¯*<sup>h</sup>*

> <sup>=</sup> <sup>1</sup> *ih*¯ *t* −∞ d*t* � *e* <sup>−</sup>*i*Kˆ(*t*−*<sup>t</sup>* � )/¯*h H*ˆ � *t* � , *ρ*ˆ

quantities in the absence of the external perturbation.

*χμν* (*ω*) <sup>=</sup> *<sup>i</sup>*

*h*¯ ∞ 0 *e <sup>i</sup>ωτ*Tr *ρ*ˆ 

 *M*ˆ *<sup>μ</sup>* (*ω*) = *i h*¯ ∞ 0

*ρ*ˆ � *<sup>I</sup>* (*t*) <sup>=</sup> <sup>1</sup> *ih*¯ *t* −∞ *k t* � d*t*

where *c* is a constant. Let us assume that the system is in its equilibrium at *t* = −∞, *i.e.*, *<sup>H</sup>*<sup>ˆ</sup> � (*<sup>t</sup>* <sup>=</sup> <sup>−</sup>∞) <sup>=</sup> 0, and consequently *<sup>ρ</sup>*<sup>ˆ</sup> <sup>=</sup> *<sup>ρ</sup>*ˆ� <sup>=</sup> 0 for *<sup>t</sup>* <sup>=</sup> <sup>−</sup>∞. In this case, the constant *<sup>c</sup>* <sup>=</sup> 0.

Towards the Authentic *Ab Intio* Thermodynamics 569

Therefore, the linear response *A*ˆ (*t*) by the external perturbation *H*ˆ � is able to be computed by

When the external perturbation Hamiltonian in Eq. (142) can be written in an oscillatory form

*<sup>H</sup>*<sup>ˆ</sup> � <sup>=</sup> <sup>−</sup>*Be*<sup>ˆ</sup> *<sup>i</sup>ω<sup>t</sup>*

Since the external perturbation *H*ˆ � does not appear in Eq. (153) at all, we can treat *A*ˆ*<sup>I</sup>* as the Heisenberg representation with respect to the natural motion described by <sup>K</sup><sup>ˆ</sup> . Consequently, the linear response is determined from a correlation function of the *fluctuations* of the relevant

Computation of the dynamic magnetic susceptibility illustrates the usefulness of the linear response theory. When the oscillatory external field in Eq. (152) is applied in the *ν*-direction,

where we use the abbreviation *<sup>τ</sup>* <sup>=</sup> *<sup>t</sup>* <sup>−</sup> *<sup>t</sup>*� and replaces *<sup>B</sup>*<sup>ˆ</sup> by *<sup>M</sup>*<sup>ˆ</sup> *<sup>ν</sup>Hν*, where *<sup>H</sup><sup>ν</sup>* is the *<sup>ν</sup>*-the component of the external magnetic field **H**. Dynamic magnetic susceptibility, defined as

> *M*ˆ *<sup>μ</sup>* (*ω*)

> > *H<sup>ν</sup>*

the thermal expectation value of the magnetization in the *μ*-direction is obtained as

*χμν* (*ω*) ≡

*<sup>e</sup>iωτ*Tr *ρ*ˆ 

)/¯*h H*ˆ � *t* � , *ρ*ˆ ∓ *e <sup>i</sup>*Kˆ(*t*−*t*�

> *<sup>t</sup>* <sup>−</sup> *<sup>t</sup>* � , *B*ˆ ∓ d *t* − *t* �

) is in the interaction representation as in Eq. (148) and we use the cyclic

*M*ˆ *<sup>μ</sup>* (*τ*), *M*ˆ *<sup>ν</sup>*

*M*ˆ *<sup>μ</sup>* (*τ*), *M*ˆ *<sup>ν</sup>*

 − 

 − 

*<sup>e</sup>*−*i*Kˆ(*t*−*t*�

*eiω*(*t*−*<sup>t</sup>* � ) Tr *ρ*ˆ *A*ˆ*I*

Considering the interaction representation in Eq. (148), we arrive at the solution

 *t* −∞ d*t* � *k t* � *<sup>e</sup>i*K<sup>ˆ</sup> *<sup>t</sup>*/¯*<sup>h</sup>*

� + *c*, (149)

)/¯*hA*ˆ 

, (152)

)/¯*h*. (150)

. (151)

, (153)

*Hν*d*τ*, (154)

d*τ*. (156)

, (155)

<sup>∓</sup> *<sup>e</sup>i*Kˆ(*t*−*t*�

We consider the external perturbation being a time dependent driving force

$$
\hat{H}'(t) \sim -F(t) \,. \tag{142}
$$

When the system is in its natural motion, the statistical distribution function, or the *density matrix*, *ρ*ˆ, is calculated by Eq. (121). With the external perturbation *H*ˆ � , the *total* density matrix will be differ from its natural density matrix as

$$
\mathfrak{d}\_{\text{total}} = \mathfrak{d} + \mathfrak{d}'.\tag{143}
$$

For any physical observables *A*ˆ, the thermal expectation value of *A*ˆ is given by

$$
\langle \hat{A} \rangle = \text{Tr} \left( \not p' \left( t \right) \hat{A} \right) \equiv \left\langle A \left( t \right) \right\rangle. \tag{144}
$$

The initial ensemble which represents statistically the initial state of the system is specified by the density matrix *ρ*ˆ satisfying

$$\left[\left[\hat{\mathbf{H}},\hat{\boldsymbol{\rho}}\right]\_{\mp} = \left[\hat{\mathcal{K}},\hat{\boldsymbol{\rho}}\right]\_{\mp} = \mathbf{0},\tag{145}$$

where *A*ˆ, *B*ˆ <sup>∓</sup> is the Poisson bracket to be *<sup>A</sup>*ˆ*B*<sup>ˆ</sup> <sup>∓</sup> *<sup>B</sup>*<sup>ˆ</sup> *<sup>A</sup>*ˆ, for any quantum operators *<sup>A</sup>*<sup>ˆ</sup> and *<sup>B</sup>*ˆ. The symbol ∓ indicates the Poisson bracket for bosons and fermions, respectively. On the other hand, the motion of the ensemble under the perturbation Eq. (142) is represented by *ρ*ˆ� (*t*), which obeys the equation

$$\begin{split} i\hbar \frac{\partial}{\partial t} \left( \not{\rho} + \not{\rho}' \right) &= \left[ \hat{H} + \hat{H}', \not{\rho} + \not{\rho}' \right]\_{\mp} \\ &= \left[ \hat{\mathcal{K}} + \hat{H}', \not{\rho} + \not{\rho}' \right]\_{\mp} . \end{split} \tag{146}$$

By assuming that the external perturbation is sufficiently weak, we have a *linear* equation

$$i\hbar\frac{\partial}{\partial t}\not p' = \left[\hat{\mathbb{K}},\not p\right]\_{\mp} + \left[\hat{H}',\not p\right]\_{\mp'}\tag{147}$$

with Eq. (145) and the neglecting the second order term, *H*ˆ � , *ρ*ˆ� ∓. Let us solve Eq. (147) by introducing the interaction representation of the density matrix as

$$\boldsymbol{\beta}\_{I}^{\prime}\left(\boldsymbol{t}\right) = \boldsymbol{e}^{i\boldsymbol{\hat{\mathcal{K}}}\boldsymbol{t}/\hbar}\boldsymbol{\beta}^{\prime}\boldsymbol{e}^{-i\boldsymbol{\hat{\mathcal{K}}}\boldsymbol{t}/\hbar}.\tag{148}$$

Differntiation both sides of Eq. (148) with respect to *t* yields

$$i\hbar\frac{\partial}{\partial t}\rho\_I'(t) = -\left[\hat{\mathcal{K}}\_\prime\hat{\rho}'(t)\right]\_\mp + e^{i\hat{\mathcal{K}}t/\hbar}\left(i\hbar\frac{\partial}{\partial t}\rho'\right)e^{-i\hat{\mathcal{K}}t/\hbar}.$$

where the second term on the right hand side is replaced by Eq. (148) to obtain

$$i\hbar \frac{\partial}{\partial t} \rho\_I' \left( t \right) = e^{i\hat{\mathcal{K}}t/\hbar} \left[ \hat{H}' , \rho \right]\_\mp e^{-i\hat{\mathcal{K}}t/\hbar} \equiv k\left( t \right) \dots$$

26 Will-be-set-by-IN-TECH

Let the Hamiltonian of the *natural motion* of a many-body system to be *H*ˆ . We would like to apply an external perturbation *H*ˆ � to the system, which yields now the total Hamiltonian

*H*ˆtotal = *H*ˆ + *H*ˆ �

When the system is in its natural motion, the statistical distribution function, or the *density*

*ρ*ˆtotal = *ρ*ˆ + *ρ*ˆ

The initial ensemble which represents statistically the initial state of the system is specified by

symbol ∓ indicates the Poisson bracket for bosons and fermions, respectively. On the other hand, the motion of the ensemble under the perturbation Eq. (142) is represented by *ρ*ˆ� (*t*),

=

By assuming that the external perturbation is sufficiently weak, we have a *linear* equation

Let us solve Eq. (147) by introducing the interaction representation of the density matrix as

*<sup>i</sup>*K<sup>ˆ</sup> *<sup>t</sup>*/¯*hρ*<sup>ˆ</sup> �

> *H*ˆ � , *ρ*ˆ

*<sup>i</sup>*K<sup>ˆ</sup> *<sup>t</sup>*/¯*<sup>h</sup> ih*¯ *∂ ∂t ρ*ˆ � 

<sup>K</sup><sup>ˆ</sup> , *<sup>ρ</sup>*<sup>ˆ</sup> 

*H*ˆ + *H*ˆ �

<sup>K</sup><sup>ˆ</sup> <sup>+</sup> *<sup>H</sup>*<sup>ˆ</sup> �

<sup>∓</sup> is the Poisson bracket to be *<sup>A</sup>*ˆ*B*<sup>ˆ</sup> <sup>∓</sup> *<sup>B</sup>*<sup>ˆ</sup> *<sup>A</sup>*ˆ, for any quantum operators *<sup>A</sup>*<sup>ˆ</sup> and *<sup>B</sup>*ˆ. The

, *ρ*ˆ + *ρ*ˆ � ∓

, *ρ*ˆ + *ρ*ˆ � 

> *H*ˆ � , *ρ*ˆ� ∓.

<sup>∓</sup> *<sup>e</sup>*−*i*K<sup>ˆ</sup> *<sup>t</sup>*/¯*<sup>h</sup>* <sup>≡</sup> *<sup>k</sup>* (*t*).

�

We consider the external perturbation being a time dependent driving force

For any physical observables *A*ˆ, the thermal expectation value of *A*ˆ is given by

*matrix*, *ρ*ˆ, is calculated by Eq. (121). With the external perturbation *H*ˆ �

 *A*ˆ = Tr *ρ*ˆ � (*t*) *A*ˆ 

*ih*¯ *∂ ∂t ρ*ˆ + *ρ*ˆ � =

> *ih*¯ *∂ ∂t ρ*ˆ � = <sup>K</sup><sup>ˆ</sup> , *<sup>ρ</sup>*<sup>ˆ</sup> <sup>∓</sup> <sup>+</sup> *H*ˆ � , *ρ*ˆ

> > *ρ*ˆ � *<sup>I</sup>* (*t*) = *e*

> > > <sup>K</sup><sup>ˆ</sup> , *<sup>ρ</sup>*<sup>ˆ</sup> � (*t*) <sup>∓</sup> <sup>+</sup> *<sup>e</sup>*

*<sup>I</sup>* (*t*) <sup>=</sup> *<sup>e</sup>i*K<sup>ˆ</sup> *<sup>t</sup>*/¯*<sup>h</sup>*

where the second term on the right hand side is replaced by Eq. (148) to obtain

with Eq. (145) and the neglecting the second order term,

Differntiation both sides of Eq. (148) with respect to *t* yields

*<sup>I</sup>* (*t*) <sup>=</sup> <sup>−</sup>

*ih*¯ *<sup>∂</sup> ∂t ρ*ˆ �

*ih*¯ *∂ ∂t ρ*ˆ �  *H*ˆ , *ρ*ˆ <sup>∓</sup> <sup>=</sup>

will be differ from its natural density matrix as

the density matrix *ρ*ˆ satisfying

where

*A*ˆ, *B*ˆ 

which obeys the equation

. (141)

. (143)

≡ �*A* (*t*)�. (144)

<sup>∓</sup> <sup>=</sup> 0, (145)

<sup>∓</sup> . (146)

<sup>∓</sup> , (147)

*<sup>e</sup>*−*i*K<sup>ˆ</sup> *<sup>t</sup>*/¯*h*. (148)

*<sup>e</sup>*−*i*K<sup>ˆ</sup> *<sup>t</sup>*/¯*h*,

, the *total* density matrix

*<sup>H</sup>*<sup>ˆ</sup> � (*t*) ∼ −*<sup>F</sup>* (*t*). (142)

This differential equation is readily solved as

$$\beta\_I'(t) = \frac{1}{i\hbar} \int\_{-\infty}^t k\left(t'\right) \mathbf{d}t' + \mathbf{c},\tag{149}$$

where *c* is a constant. Let us assume that the system is in its equilibrium at *t* = −∞, *i.e.*, *<sup>H</sup>*<sup>ˆ</sup> � (*<sup>t</sup>* <sup>=</sup> <sup>−</sup>∞) <sup>=</sup> 0, and consequently *<sup>ρ</sup>*<sup>ˆ</sup> <sup>=</sup> *<sup>ρ</sup>*ˆ� <sup>=</sup> 0 for *<sup>t</sup>* <sup>=</sup> <sup>−</sup>∞. In this case, the constant *<sup>c</sup>* <sup>=</sup> 0. Considering the interaction representation in Eq. (148), we arrive at the solution

$$\begin{split} \hat{\rho}'(t) &= \frac{1}{i\hbar} e^{-i\hat{\mathcal{K}}t/\hbar} \int\_{-\infty}^{t} \mathbf{d}t' k \left(t'\right) e^{i\hat{\mathcal{K}}t/\hbar} \\ &= \frac{1}{i\hbar} \int\_{-\infty}^{t} \mathbf{d}t' e^{-i\hat{\mathcal{K}}(t-t')/\hbar} \left[\hat{H}'\left(t'\right), \rho\right]\_{\mp} e^{i\hat{\mathcal{K}}(t-t')/\hbar}. \end{split} \tag{150}$$

Therefore, the linear response *A*ˆ (*t*) by the external perturbation *H*ˆ � is able to be computed by using Eqs. (144) and (150) to be

$$\langle \hat{A}\left(t\right) \rangle = \frac{1}{i\hbar} \int\_{-\infty}^{t} \mathbf{d}t' \text{Tr}\left(e^{-i\hat{\mathcal{K}}(t-t')/\hbar} \left[\hat{H}'\left(t'\right), \hat{\rho}\right]\_{\mp} e^{i\hat{\mathcal{K}}(t-t')/\hbar} \hat{A}\right). \tag{151}$$

When the external perturbation Hamiltonian in Eq. (142) can be written in an oscillatory form

$$
\hat{H}' = -\hat{B}e^{i\omega t},
\tag{152}
$$

we have the linear response as

$$\langle \hat{A} \left( \omega \right) \rangle = \pm \frac{i}{\hbar} \int\_0^\infty e^{i\omega \left( t - t' \right)} \text{Tr} \left( \left[ \hat{A}\_I \left( \left( t - t' \right) \right), \hat{B} \right]\_{\mp} \right) \mathbf{d} \left( t - t' \right) \, \tag{153}$$

where *AI* (*t* − *t* � ) is in the interaction representation as in Eq. (148) and we use the cyclic property of trace.

Since the external perturbation *H*ˆ � does not appear in Eq. (153) at all, we can treat *A*ˆ*<sup>I</sup>* as the Heisenberg representation with respect to the natural motion described by <sup>K</sup><sup>ˆ</sup> . Consequently, the linear response is determined from a correlation function of the *fluctuations* of the relevant quantities in the absence of the external perturbation.

Computation of the dynamic magnetic susceptibility illustrates the usefulness of the linear response theory. When the oscillatory external field in Eq. (152) is applied in the *ν*-direction, the thermal expectation value of the magnetization in the *μ*-direction is obtained as

$$
\langle \hat{M}\_{\mu} \left( \omega \right) \rangle = \frac{i}{\hbar} \int\_{0}^{\infty} e^{i\omega \tau} \text{Tr} \left( \hat{\rho} \left[ \hat{M}\_{\mu} \left( \tau \right), \hat{M}\_{\nu} \right]\_{-} \right) H\_{\nu} \text{d} \tau,\tag{154}
$$

where we use the abbreviation *<sup>τ</sup>* <sup>=</sup> *<sup>t</sup>* <sup>−</sup> *<sup>t</sup>*� and replaces *<sup>B</sup>*<sup>ˆ</sup> by *<sup>M</sup>*<sup>ˆ</sup> *<sup>ν</sup>Hν*, where *<sup>H</sup><sup>ν</sup>* is the *<sup>ν</sup>*-the component of the external magnetic field **H**. Dynamic magnetic susceptibility, defined as

$$\chi\_{\mu\nu}\left(\omega\right) \equiv \frac{\left<\hat{M}\_{\mu}\left(\omega\right)\right>}{H\_{\nu}}\,,\tag{155}$$

is then calculated as

$$\chi\_{\mu\nu}\left(\omega\right) = \frac{i}{\hbar} \int\_0^\infty e^{i\omega\tau} \text{Tr}\left(\not p \left[\hat{M}\_\mu \left(\tau\right), \hat{M}\_\nu\right]\_-\right) \text{d}\tau. \tag{156}$$

The physical meaning of Eq. (166) is now clear: the correlation function of the thermal fluctuation determines the imaginary part of the linear response function, *i.e.*, the dissipation of the system. Therefore, we call Eqs. (165) and (166) as the *fluctuation-dissipation theorem*. One can easily find a similarity between the fluctuation-dissipation theorem in Eqs. (165)–(166) and Eq. (134). The spectral representation of the retarded real-time Green's function Eq. (136) is the momentum space Lehmann representation of the two-time retarded

Towards the Authentic *Ab Intio* Thermodynamics 571

In these definitions of the two-time Green's functions, *ψ*ˆ and *ψ*ˆ† are respectively the field

A two-time Green's function represents the correlation of thermal fluctuation of the process

two-time Green's function can be used to describe all the physical correlations of the system. In the linear response theory, the interaction representation of an operator is not different from the Heisenberg representation in the two-time Green's functions. Therefore, the linear response is nothing more than the retarded two-time Green's function. Note that the *causality* of any physical process stresses that the dissipation *χ* should be the imaginary part of the

Now, we have completed the formal theory for computing thermodynamic properties in terms of the quantum field theory. An experimental observation is indeed a measurement of a response of the system by a corresponding external perturbation. The response function is related to the retarded two-time Green's function in terms of the fluctuation-dissipation theory. The two-time Green's function provides the quasiparticle energy spectra. The two-time Green's function can be transformed to a temperature Green's function by introducing the spectral weight function. Once we have the temperature Green's function, we can calculate the thermodynamic information, such as the number of particles *N*, the internal energy *E*, and the grand potential Ω, directly. By performing derivatives on the grand potential, the remaining thermodynamic functions can be evaluated. When there is a phase transition, the response function diverges near the critical region, but the grand potential does not, because the grand potential is a integrated value of the temperature Green's function in

Kubo (1966) himself demonstrated the applications of the fluctuation-dissipation theorem to compute the density response, conduction and diffusion. The fluctuation-dissipation theorem (Kubo, 1957) is able to explain the irreversible processes related with fluctuations (Callen & Welton, 1951; Casimir, 1945; Nyquist, 1928; Onsager, 1931a;b). Further applications

functions are nothing more than the correlation functions of two operators *ψ*ˆ*<sup>α</sup>* and *ψ*ˆ†

that all the quantum operators are composed of the field operators *ψ*ˆ*<sup>α</sup>* and *ψ*ˆ†

*ψ*ˆ*<sup>α</sup>* (**x***t*), *ψ*ˆ†

*ψ*ˆ*<sup>α</sup>* (**x***t*), *ψ*ˆ†

*β* **x**� *t* � ∓ 

*β* **x**� *t* � ∓ 

*<sup>A</sup>*<sup>ˆ</sup> (**x***t*) <sup>≡</sup> *<sup>e</sup>i*K<sup>ˆ</sup> *<sup>t</sup>*/¯*hAe*<sup>ˆ</sup> <sup>−</sup>*i*K<sup>ˆ</sup> *<sup>t</sup>*/¯*h*. (169)

*t*� and it annihilated at **x***t*. In other words, the two-time Green's

. (167)

. (168)

*<sup>β</sup>*. Note

*<sup>β</sup>*, so that the

Green's function defined by

that a particle created at **x**�

retarded Green's function.

the complex *ω* plane.

are well summarized by Gammaitoni et al. (1998).

*iG*¯ *<sup>R</sup> αβ* **x***t*, **x**� *t* � = *θ t* − *t* � Tr *ρ*ˆ 

*iG*¯ *<sup>A</sup> αβ* **x***t*, **x**� *t* � = −*θ t* − *t* � Tr *ρ*ˆ 

The corresponding advanced Green's function is defined by

annihilation and creation operators in the Heisenberg representation

The extension to the case of the response to an magnetic field oscillating in space with the wave vector **q** is straightforward to be

$$\chi\_{\mu\nu}\left(\mathbf{q},\omega\right) = \frac{i}{\hbar} \int\_0^\infty e^{i\omega\tau} \text{Tr}\left(\not{p}\left[\hat{M}\_\mu\left(\mathbf{q},\tau\right),\not{M}\_\nu\left(-\mathbf{q}\right)\right]\_-\right) \text{d}\tau. \tag{157}$$

In general, a linear response function is given in the form

$$\chi\_{AB}\left(\mathbf{q},\omega\right) = \pm \frac{i}{\hbar} \int\_0^\infty e^{i\omega \tau} \text{Tr}\left(\hat{\rho}\left[\hat{A}\left(\mathbf{q},\tau\right),\hat{B}\left(-\mathbf{q}\right)\right]\_{\mp}\right) \text{d}\tau. \tag{158}$$

This is the virtue of the linear response theory.

#### **5.2 Fluctuation-dissipation theorem**

The Kubo formula of the linear response function, or the generalized susceptibilty, is given in Eq. (158). By extending the interval of integration to (−∞, ∞) in Eq. (158), we define a function *f* as

$$f\_{AB}\left(\mathbf{q},\omega\right) = \pm \frac{i}{\hbar} \int\_{-\infty}^{\infty} e^{i\omega \tau} \text{Tr}\left(\oint \left[\hat{A}\left(\mathbf{q},\tau\right),\hat{B}\left(-\mathbf{q}\right)\right]\_{\mp}\right) d\tau. \tag{159}$$

The integration is now devided into the negative and positive ranges as

$$f\_{AB}\left(\mathbf{q},\omega\right) = \pm \frac{i}{\hbar} \int\_0^\infty e^{i\omega \tau} \text{Tr}\left(\not{p}\left[\not{A}\left(\mathbf{q},\tau\right),\not{B}\left(-\mathbf{q}\right)\right]\_\mp\right) \text{d}\tau$$

$$\mp \frac{i}{\hbar} \int\_0^\infty e^{i\omega \tau} \text{Tr}\left(\not{p}\left[\not{A}\left(\mathbf{q}\right),\not{B}\left(-\mathbf{q},-\tau\right)\right]\_\mp\right) \text{d}\tau,\tag{160}$$

where we changed *τ* to −*τ* in the second term of the right hand side, and moved the factors exp (±*i*K*τ*/¯*h*) around *A* (**q**) to aound *B* (−**q**) using the cyclic property of trace. Since Eq. (158) has a symmetry

$$
\chi\_{AB}^\* \left( \mathbf{q}\_{\prime} \omega \right) = \chi\_{AB} \left( -\mathbf{q}\_{\prime} - \omega \right),
\tag{161}
$$

it is able to denote *fAB* as

$$f\_{AB} = \chi\_{AB} \left( \mathbf{q}\_{\prime} \omega \right) \mp \chi\_{BA}^\* \left( \mathbf{q}\_{\prime} \omega \right). \tag{162}$$

Eq. (159) is written as

$$f\_{AB}\left(\mathbf{q},\omega\right) = \left(1 \mp e^{-\beta\hbar\omega}\right) \frac{i}{\hbar} \int\_{-\infty}^{\infty} e^{i\omega\tau} \text{Tr}\left(\hat{\rho}\left[\hat{A}\left(\mathbf{q},\tau\right),\hat{B}\left(-\mathbf{q}\right)\right]\_{\mp}\right) d\tau,\tag{163}$$

because of the identity

$$\int\_{-\infty}^{\infty} e^{i\omega \tau} \text{Tr}\left(\hat{\rho}\left[\hat{A}\left(\tau\right), \hat{B}\right]\_{\mp}\right) \text{d}\tau = \pm e^{\mathcal{S}\hbar\omega} \int\_{-\infty}^{\infty} e^{i\omega \tau} \text{Tr}\left(\hat{\rho}\left[\hat{B}, \hat{A}\left(\tau\right)\right]\_{\mp}\right) \text{d}\tau. \tag{164}$$

The two representations of *fAB* in Eqs. (162) and (163) lead us a result

$$\int\_{-\infty}^{\infty} e^{i\omega t} \text{Tr}\left(\boldsymbol{\hat{\rho}} \left[\hat{A}\left(\mathbf{q}, t\right), \boldsymbol{\mathcal{B}}\left(-\mathbf{q}\right)\right]\_{\mp}\right) \text{d}t = \frac{i\hbar}{e^{-\beta\hbar\omega} \mp 1} \left\{ \chi\_{AB}\left(\mathbf{q}, \omega\right) \mp \chi\_{BA}^\*\left(\mathbf{q}, \omega\right) \right\}\,\mathrm{\,,}\tag{165}$$

where we changed the variable *τ* by *t*. Eq. (165) relates the linear response function and the ordinary correlation function. For a special case of *A* = *B*, Eq. (165) is reduced to

$$\int\_{-\infty}^{\infty} e^{i\omega t} \text{Tr}\left(\hat{\rho}\left[\hat{A}\left(\mathbf{q},t\right),\hat{A}\left(-\mathbf{q}\right)\right]\_{\mp}\right) \text{d}t = -\frac{\hbar}{e^{-\beta\hbar\omega} \mp 1} 2\hat{\c}\chi\_{AA}\left(\mathbf{q},\omega\right). \tag{166}$$

28 Will-be-set-by-IN-TECH

The extension to the case of the response to an magnetic field oscillating in space with the

The Kubo formula of the linear response function, or the generalized susceptibilty, is given in Eq. (158). By extending the interval of integration to (−∞, ∞) in Eq. (158), we define a

where we changed *τ* to −*τ* in the second term of the right hand side, and moved the factors exp (±*i*K*τ*/¯*h*) around *A* (**q**) to aound *B* (−**q**) using the cyclic property of trace. Since Eq.

> *eiωτ*Tr *ρ*ˆ

> > ∞ −∞ *e <sup>i</sup>ωτ*Tr *ρ*ˆ *B*ˆ, *A*ˆ (*τ*) ∓

<sup>d</sup>*<sup>t</sup>* <sup>=</sup> <sup>−</sup> *<sup>h</sup>*¯

*<sup>e</sup>*−*βh*¯ *<sup>ω</sup>* <sup>∓</sup> <sup>1</sup> {*χAB* (**q**, *<sup>ω</sup>*) <sup>∓</sup> *<sup>χ</sup>*<sup>∗</sup>

*<sup>e</sup>*−*βh*¯ *<sup>ω</sup>* <sup>∓</sup> <sup>1</sup>

<sup>d</sup>*<sup>τ</sup>* <sup>=</sup> <sup>±</sup>*eβh*¯ *<sup>ω</sup>*

<sup>d</sup>*<sup>t</sup>* <sup>=</sup> *ih*¯

where we changed the variable *τ* by *t*. Eq. (165) relates the linear response function and the

 ∓ 

*fAB* = *χAB* (**q**, *ω*) ∓ *χ*<sup>∗</sup>

 *i h*¯ ∞ −∞

*eiωτ*Tr *ρ*ˆ 

*<sup>M</sup>*<sup>ˆ</sup> *<sup>μ</sup>* (**q**, *<sup>τ</sup>*), *<sup>M</sup>*<sup>ˆ</sup> *<sup>ν</sup>* (−**q**)

*<sup>A</sup>*<sup>ˆ</sup> (**q**, *<sup>τ</sup>*), *<sup>B</sup>*<sup>ˆ</sup> (−**q**)

*<sup>A</sup>*<sup>ˆ</sup> (**q**, *<sup>τ</sup>*), *<sup>B</sup>*<sup>ˆ</sup> (−**q**)

*<sup>A</sup>*<sup>ˆ</sup> (**q**, *<sup>τ</sup>*), *<sup>B</sup>*<sup>ˆ</sup> (−**q**)

*<sup>A</sup>*<sup>ˆ</sup> (**q**), *<sup>B</sup>*<sup>ˆ</sup> (−**q**, <sup>−</sup>*τ*)

*AB* (**q**, *ω*) = *χAB* (−**q**, −*ω*), (161)

*<sup>A</sup>*<sup>ˆ</sup> (**q**, *<sup>τ</sup>*), *<sup>B</sup>*<sup>ˆ</sup> (−**q**)

 − 

 ∓ 

 ∓ 

 ∓ d*τ*

> ∓

*BA* (**q**, *ω*). (162)

 ∓ 

d*τ*. (157)

d*τ*. (158)

d*τ*. (159)

d*τ*, (160)

d*τ*, (163)

d*τ*. (164)

*BA* (**q**, *ω*)} , (165)

2�*χAA* (**q**, *ω*). (166)

wave vector **q** is straightforward to be

*χμν* (**q**, *<sup>ω</sup>*) <sup>=</sup> *<sup>i</sup>*

*h*¯ ∞ 0 *e <sup>i</sup>ωτ*Tr *ρ*ˆ 

> *h*¯ ∞ 0 *e <sup>i</sup>ωτ*Tr *ρ*ˆ

*h*¯ ∞ −∞

*h*¯ ∞ 0 *e <sup>i</sup>ωτ*Tr *ρ*ˆ 

∓ *i h*¯ ∞ 0 *e <sup>i</sup>ωτ*Tr *ρ*ˆ 

*χ*∗

<sup>1</sup> <sup>∓</sup> *<sup>e</sup>*−*βh*¯ *<sup>ω</sup>*

The two representations of *fAB* in Eqs. (162) and (163) lead us a result

 ∓ 

*<sup>A</sup>*<sup>ˆ</sup> (**q**,*t*), *<sup>A</sup>*<sup>ˆ</sup> (−**q**)

ordinary correlation function. For a special case of *A* = *B*, Eq. (165) is reduced to

The integration is now devided into the negative and positive ranges as

In general, a linear response function is given in the form

*<sup>χ</sup>AB* (**q**, *<sup>ω</sup>*) <sup>=</sup> <sup>±</sup> *<sup>i</sup>*

*fAB* (**q**, *<sup>ω</sup>*) <sup>=</sup> <sup>±</sup> *<sup>i</sup>*

*fAB* (**q**, *<sup>ω</sup>*) <sup>=</sup> <sup>±</sup> *<sup>i</sup>*

This is the virtue of the linear response theory.

**5.2 Fluctuation-dissipation theorem**

function *f* as

(158) has a symmetry

Eq. (159) is written as

because of the identity ∞ −∞

> ∞ −∞ *e iωt* Tr *ρ*ˆ

it is able to denote *fAB* as

*fAB* (**q**, *ω*) =

*eiωτ*Tr *ρ*ˆ *A*ˆ (*τ*), *B*ˆ ∓ 

 ∞ −∞ *e iωt* Tr *ρ*ˆ 

*<sup>A</sup>*<sup>ˆ</sup> (**q**, *<sup>t</sup>*), *<sup>B</sup>*<sup>ˆ</sup> (−**q**)

The physical meaning of Eq. (166) is now clear: the correlation function of the thermal fluctuation determines the imaginary part of the linear response function, *i.e.*, the dissipation of the system. Therefore, we call Eqs. (165) and (166) as the *fluctuation-dissipation theorem*. One can easily find a similarity between the fluctuation-dissipation theorem in Eqs. (165)–(166) and Eq. (134). The spectral representation of the retarded real-time Green's function Eq. (136) is the momentum space Lehmann representation of the two-time retarded Green's function defined by

$$i\bar{\mathcal{G}}\_{a\beta}^{R}\left(\mathbf{x}t,\mathbf{x}'t'\right) = \theta\left(t - t'\right)\text{Tr}\left(\hat{\rho}\left[\hat{\psi}\_{a}\left(\mathbf{x}t\right),\hat{\psi}\_{\beta}^{\dagger}\left(\mathbf{x}'t'\right)\right]\_{\mp}\right).\tag{167}$$

The corresponding advanced Green's function is defined by

$$i\bar{\mathbf{G}}\_{a\beta}^{A}\left(\mathbf{x}t,\mathbf{x}'t'\right)=-\theta\left(t-t'\right)\operatorname{Tr}\left(\boldsymbol{\upphi}\left[\boldsymbol{\upphi}\_{a}\left(\mathbf{x}t\right),\boldsymbol{\upphi}\_{\beta}^{\dagger}\left(\mathbf{x}'t'\right)\right]\_{\mp}\right).\tag{168}$$

In these definitions of the two-time Green's functions, *ψ*ˆ and *ψ*ˆ† are respectively the field annihilation and creation operators in the Heisenberg representation

$$
\hat{A}\left(\mathbf{x}t\right) \equiv e^{i\hat{\mathcal{K}}t/\hbar} \hat{A}e^{-i\hat{\mathcal{K}}t/\hbar}.\tag{169}
$$

A two-time Green's function represents the correlation of thermal fluctuation of the process that a particle created at **x**� *t*� and it annihilated at **x***t*. In other words, the two-time Green's functions are nothing more than the correlation functions of two operators *ψ*ˆ*<sup>α</sup>* and *ψ*ˆ† *<sup>β</sup>*. Note

that all the quantum operators are composed of the field operators *ψ*ˆ*<sup>α</sup>* and *ψ*ˆ† *<sup>β</sup>*, so that the two-time Green's function can be used to describe all the physical correlations of the system. In the linear response theory, the interaction representation of an operator is not different from the Heisenberg representation in the two-time Green's functions. Therefore, the linear response is nothing more than the retarded two-time Green's function. Note that the *causality* of any physical process stresses that the dissipation *χ* should be the imaginary part of the retarded Green's function.

Now, we have completed the formal theory for computing thermodynamic properties in terms of the quantum field theory. An experimental observation is indeed a measurement of a response of the system by a corresponding external perturbation. The response function is related to the retarded two-time Green's function in terms of the fluctuation-dissipation theory. The two-time Green's function provides the quasiparticle energy spectra. The two-time Green's function can be transformed to a temperature Green's function by introducing the spectral weight function. Once we have the temperature Green's function, we can calculate the thermodynamic information, such as the number of particles *N*, the internal energy *E*, and the grand potential Ω, directly. By performing derivatives on the grand potential, the remaining thermodynamic functions can be evaluated. When there is a phase transition, the response function diverges near the critical region, but the grand potential does not, because the grand potential is a integrated value of the temperature Green's function in the complex *ω* plane.

Kubo (1966) himself demonstrated the applications of the fluctuation-dissipation theorem to compute the density response, conduction and diffusion. The fluctuation-dissipation theorem (Kubo, 1957) is able to explain the irreversible processes related with fluctuations (Callen & Welton, 1951; Casimir, 1945; Nyquist, 1928; Onsager, 1931a;b). Further applications are well summarized by Gammaitoni et al. (1998).

The density of the system will contain a linear response to be

(Mermin, 1965) yields

condition can be written as

static case *ω* = 0.

where the density response function is

response function should satisfy the condition

<sup>−</sup> <sup>1</sup>

Then, the density fluctuation for constant *n*¯ is wrtten as

effective potential dependence on *T* and *q* as

*<sup>χ</sup>*−<sup>1</sup> (*q*, [*n*¯]) <sup>=</sup> <sup>1</sup>

*<sup>S</sup>* (*q*) <sup>=</sup> <sup>1</sup>

*n* (**x**) = *n*¯ +

*n* (**q**)

1 <sup>2</sup> ∑ **q**�=0

<sup>−</sup> *<sup>χ</sup>* (*q*, [*n*¯]) <sup>=</sup> <sup>1</sup>

<sup>−</sup> <sup>1</sup>

Towards the Authentic *Ab Intio* Thermodynamics 573

where *n*¯ is the mean density and *n* (**q**) is the density fluctuations. The Mermin condition

and *v*eff is the corresponding effective potential. For free energy minimum, the density

for the system is in an equilibrium state either absolutely stable or metastable. Therfore, <sup>−</sup>*χ*−<sup>1</sup> (*q*, [*n*¯]) <sup>=</sup> 0 determines the stability limit of the fluid (*qL*, *PL*, *TL*). The static stability

*<sup>β</sup><sup>S</sup>* (*q*) <sup>=</sup> <sup>1</sup>

where *S* (*q*) is the static mean density fluctuation. The condition Eq. (176) impiles that the mean density fluctuation *S* (*q*) with the particular wave number *qL* increases as this limit is approached. Near the stability limit only the region *q* � *qL* is important. We may set the

*<sup>β</sup>v*eff (*q*) <sup>=</sup> <sup>−</sup><sup>1</sup> <sup>+</sup> *<sup>A</sup>* (*<sup>T</sup>* <sup>−</sup> *TL*) <sup>+</sup> *<sup>B</sup>* (*<sup>q</sup>* <sup>−</sup> *qL*)

where *<sup>ξ</sup>*<sup>2</sup> <sup>=</sup> *<sup>B</sup>*/*<sup>A</sup>* (*<sup>T</sup>* <sup>−</sup> *TL*). As *<sup>q</sup>* approaches to *qL*, *<sup>S</sup>* (*q*) increases and the peak becomes sharper. The same effect is also seen in the dynamic mean density fluctuation *S* (*q*, *ω*) at the

This result is consistent with our everyday experiences. In a fluid state, there are random density fluctuations, which can be seen if there is an interface; a random wave distribution on a smooth lake surface is an example. Such density fluctuation is represented in the Fourier terms, by a linear combination of the wave vectors and frequencies. When the temperature is cooled down with a constant pressure, the fluid transforms to crystals, *e.g.*, ice. A crystal is by definition a regular arrangement of density, with a spatial period of a corresponding wave vector. Eq. (178) tells us that the static mean density fluctuation at the critical wave vector *qL* with the corresponding regularity will be sustained, even if there is no external

1/*ξ*<sup>2</sup> 1/*ξ*<sup>2</sup> + (*<sup>q</sup>* − *qL*)

*A* (*T* − *TL*)

*n* (**q**)*ei***q**·**<sup>x</sup>** + *c*.*c*.

*<sup>n</sup>*¯ <sup>=</sup> *<sup>χ</sup>* (*q*, [*n*¯]) *<sup>v</sup>*ext (**q**), (173)

*<sup>β</sup>*−<sup>1</sup> <sup>+</sup> *<sup>v</sup>*eff (*q*, [*n*¯]), (174)

*<sup>χ</sup>* (*q*, [*n*¯]) <sup>&</sup>gt; <sup>0</sup> (175)

*<sup>β</sup>* <sup>+</sup> *<sup>v</sup>*eff (*q*, [*n*¯]) <sup>=</sup> 0, (176)

<sup>2</sup> . (177)

<sup>2</sup> , (178)

, (172)

The remained work is mainly on the applications to the realistic matter.

### **6. Inhomogeneity is reality**

We have a brief review on the formalism of the quantum field theoretic approach for thermodynamics properties of the system in homogeneous (and isotropic) distributions. The homogeneous approximation of a matter has served as a model system of matter for investigating the nature of the complex phenomena occuring in the matter. The study on electron liquid is an example (Fetter & Walecka, 2003; Giuliani & Vignale, 2005; Kim, 1999; Pines, 1962; 1999). Phonon-mediated superconductivity can be solvable under the method after performing a suitable transformation for treating Cooper pairs (Schrieffer, 1988). A relativistic version of the method has been applied to plasmas (Melrose, 2008).

However, the realistic materials has internal (micro-)structures and even pure elemental materials are not fully homogeneous. Steels as representative structural materials have large variety of microstructures, which influences on the properties, *e.g.*, strength and elongations (Bhadeshia & Honeycombe, 2006). In metals and alloys, the internal structures of materials are governed by the equilibrium thermodynamics and the irreversible processes (Christian, 2002). It is necessary to deal with the inhomogeneity when we are dealing with thermodynamics and irreversible processes for materials sciences.

#### **6.1 Density functional theory**

The applications of the formal quantum theories, developed in a homogeneous model system, to the realistic inhomogeneous systems have long time history (Fermi, 1927; Slater, 1951; Thomas, 1938). An innovative formalism to treat the inhomogeneous electronic system was suggested by Hohenberg & Kohn (1964). They firstly proved that the ground state energy is a *unique* functional of density, *E* [*n* (**x**)]. Due to the uniqueness of the ground state energy, the discussions in the previous sections can be reformulated with density. This is the reason why the Hohenberg-Kohn theory is called the density functional theory (DFT). For a slowly vaying density profile, the density functional theory is a generalization of the Thomas-Fermi theory. They proved also that the essential density gradient in the inhomogeneous electron gas can be used for describing the energy functional.

An immediate extension of the Hohenberg-Kohn theory to the finite temperature was done by Mermin (1965). At a given temperature and chemical potential, no two local potential *v* (**x**) lead to the same equilibrium density. In other words, one can define a density functional *F* [*n* (**x**)] independent of *v* (**x**), such that the quantity

$$\Omega \equiv \int v\left(\mathbf{x}\right) n\left(\mathbf{x}\right) \mathbf{d}^3 \mathbf{x} + F\left[n\left(\mathbf{x}\right)\right],\tag{170}$$

is at a minimum and equal to the grand potential when *n* (**x**) is the equilibrium density in the grand canonical ensemble in the presence of *v* (**x**). Kohn & Sham (1965) themselves were able to derive the forms for the grand potential, specific heat, and spin susceptibility.

Schneider (1971) provided a general theory of the liquid-solid phase transition by employing the finite temperature density functional theory by virtue of the fluctuation-dissipation theorem on the density response. During the construction of the functional *F* [*n*], the external perturbation potential is assumed to be of form

$$v^{\rm ext}\left(\mathbf{x}\right) = \frac{1}{2} \sum\_{\mathbf{q}\neq 0} \left(v^{\rm ext}\left(\mathbf{q}\right)e^{i\mathbf{q}\cdot\mathbf{x}} + c.c.\right). \tag{171}$$

30 Will-be-set-by-IN-TECH

We have a brief review on the formalism of the quantum field theoretic approach for thermodynamics properties of the system in homogeneous (and isotropic) distributions. The homogeneous approximation of a matter has served as a model system of matter for investigating the nature of the complex phenomena occuring in the matter. The study on electron liquid is an example (Fetter & Walecka, 2003; Giuliani & Vignale, 2005; Kim, 1999; Pines, 1962; 1999). Phonon-mediated superconductivity can be solvable under the method after performing a suitable transformation for treating Cooper pairs (Schrieffer, 1988). A

However, the realistic materials has internal (micro-)structures and even pure elemental materials are not fully homogeneous. Steels as representative structural materials have large variety of microstructures, which influences on the properties, *e.g.*, strength and elongations (Bhadeshia & Honeycombe, 2006). In metals and alloys, the internal structures of materials are governed by the equilibrium thermodynamics and the irreversible processes (Christian, 2002). It is necessary to deal with the inhomogeneity when we are dealing with thermodynamics and

The applications of the formal quantum theories, developed in a homogeneous model system, to the realistic inhomogeneous systems have long time history (Fermi, 1927; Slater, 1951; Thomas, 1938). An innovative formalism to treat the inhomogeneous electronic system was suggested by Hohenberg & Kohn (1964). They firstly proved that the ground state energy is a *unique* functional of density, *E* [*n* (**x**)]. Due to the uniqueness of the ground state energy, the discussions in the previous sections can be reformulated with density. This is the reason why the Hohenberg-Kohn theory is called the density functional theory (DFT). For a slowly vaying density profile, the density functional theory is a generalization of the Thomas-Fermi theory. They proved also that the essential density gradient in the inhomogeneous electron gas can

An immediate extension of the Hohenberg-Kohn theory to the finite temperature was done by Mermin (1965). At a given temperature and chemical potential, no two local potential *v* (**x**) lead to the same equilibrium density. In other words, one can define a density functional

is at a minimum and equal to the grand potential when *n* (**x**) is the equilibrium density in the grand canonical ensemble in the presence of *v* (**x**). Kohn & Sham (1965) themselves were able

Schneider (1971) provided a general theory of the liquid-solid phase transition by employing the finite temperature density functional theory by virtue of the fluctuation-dissipation theorem on the density response. During the construction of the functional *F* [*n*], the external

to derive the forms for the grand potential, specific heat, and spin susceptibility.

<sup>2</sup> ∑ **q**�=0 *v*ext (**q**)*e*

*<sup>i</sup>***q**·**<sup>x</sup>** + *c*.*c*.

. (171)

*v* (**x**) *n* (**x**) d3**x** + *F* [*n* (**x**)] , (170)

The remained work is mainly on the applications to the realistic matter.

relativistic version of the method has been applied to plasmas (Melrose, 2008).

**6. Inhomogeneity is reality**

irreversible processes for materials sciences.

be used for describing the energy functional.

*F* [*n* (**x**)] independent of *v* (**x**), such that the quantity

perturbation potential is assumed to be of form

Ω ≡ 

*<sup>v</sup>*ext (**x**) <sup>=</sup> <sup>1</sup>

**6.1 Density functional theory**

The density of the system will contain a linear response to be

$$n\left(\mathbf{x}\right) = \overline{n} + \frac{1}{2} \sum\_{\mathbf{q}\neq 0} \left( n\left(\mathbf{q}\right) e^{i\mathbf{q}\cdot\mathbf{x}} + c.c.\right),\tag{172}$$

where *n*¯ is the mean density and *n* (**q**) is the density fluctuations. The Mermin condition (Mermin, 1965) yields

$$\frac{m\left(\mathbf{q}\right)}{\vec{n}} = \chi\left(\boldsymbol{q}, \left[\vec{n}\right]\right) v^{\text{ext}}\left(\mathbf{q}\right),\tag{173}$$

where the density response function is

$$-\chi\left(\boldsymbol{\eta}, \left[\vec{n}\right]\right) = \frac{1}{\beta^{-1} + \boldsymbol{\upsilon}^{\mathrm{eff}}\left(\boldsymbol{\eta}, \left[\vec{n}\right]\right)}\,\tag{174}$$

and *v*eff is the corresponding effective potential. For free energy minimum, the density response function should satisfy the condition

$$-\frac{1}{\chi\left(q\_{\prime}\left[\vec{n}\right]\right)} > 0\tag{175}$$

for the system is in an equilibrium state either absolutely stable or metastable. Therfore, <sup>−</sup>*χ*−<sup>1</sup> (*q*, [*n*¯]) <sup>=</sup> 0 determines the stability limit of the fluid (*qL*, *PL*, *TL*). The static stability condition can be written as

$$-\frac{1}{\chi^{-1}\left(q\_{\prime}\left[\bar{n}\right]\right)} = \frac{1}{\beta S\left(q\right)} = \frac{1}{\beta} + \boldsymbol{\upsilon}^{\mathrm{eff}}\left(q\_{\prime}\left[\bar{n}\right]\right) = 0,\tag{176}$$

where *S* (*q*) is the static mean density fluctuation. The condition Eq. (176) impiles that the mean density fluctuation *S* (*q*) with the particular wave number *qL* increases as this limit is approached. Near the stability limit only the region *q* � *qL* is important. We may set the effective potential dependence on *T* and *q* as

$$\beta \upsilon^{\rm eff}(q) = -1 + A \left( T - T\_L \right) + B \left( q - q\_L \right)^2. \tag{177}$$

Then, the density fluctuation for constant *n*¯ is wrtten as

$$S\left(\eta\right) = \frac{1}{A\left(T - T\_L\right)} \frac{1/\tilde{\xi}^2}{1/\tilde{\xi}^2 + \left(\eta - q\_L\right)^2},\tag{178}$$

where *<sup>ξ</sup>*<sup>2</sup> <sup>=</sup> *<sup>B</sup>*/*<sup>A</sup>* (*<sup>T</sup>* <sup>−</sup> *TL*). As *<sup>q</sup>* approaches to *qL*, *<sup>S</sup>* (*q*) increases and the peak becomes sharper. The same effect is also seen in the dynamic mean density fluctuation *S* (*q*, *ω*) at the static case *ω* = 0.

This result is consistent with our everyday experiences. In a fluid state, there are random density fluctuations, which can be seen if there is an interface; a random wave distribution on a smooth lake surface is an example. Such density fluctuation is represented in the Fourier terms, by a linear combination of the wave vectors and frequencies. When the temperature is cooled down with a constant pressure, the fluid transforms to crystals, *e.g.*, ice. A crystal is by definition a regular arrangement of density, with a spatial period of a corresponding wave vector. Eq. (178) tells us that the static mean density fluctuation at the critical wave vector *qL* with the corresponding regularity will be sustained, even if there is no external

where *T*<sup>0</sup> is the kinetic energy functional of a system of noninteracting electrons and *E*xc [*n*] is the exchange and correlation energy functional of an interacting system with density *n* (**x**).

Towards the Authentic *Ab Intio* Thermodynamics 575

the stationary condition of the energy functional Eq. (181) with respect to the density variation

where *v* (**x**) is the result of the functional derivatives of the static integration terms of Eq. (181) and *μ*xc is the result of the functional derivatives of the exchange-correlation functional *E*xc [*n*]. By assuming *μ*xc (*n* (**x**)) to be local, we can imagine a "ficticious" particle that moves by the

*ψ*KS

*<sup>i</sup>* (**x**) <sup>=</sup> *�*KS

*<sup>i</sup> <sup>ψ</sup>*KS

*<sup>i</sup>* (**x**), (184)

, (185)

*n* (**x**) *�<sup>c</sup>* (*n* (**x**)) d3**x**. (186)

d3**x**� = *�*KS

*<sup>i</sup> <sup>ψ</sup>*KS

. (188)

*<sup>i</sup>* (**x**), (187)

*n* (**x**) *�*xc (*n* (**x**)) d3**x**, (182)

*<sup>δ</sup><sup>n</sup>* (**x**) <sup>+</sup> *<sup>v</sup>* (**x**) <sup>+</sup> *<sup>μ</sup>*xc (*<sup>n</sup>* (**x**)) <sup>=</sup> 0, (183)

*E*xc [*n*] =

*δT*<sup>0</sup> [*n*]

<sup>∇</sup><sup>2</sup> <sup>+</sup> *<sup>ϕ</sup>* (**x**) <sup>+</sup> *<sup>μ</sup>*xc (*<sup>n</sup>* (**x**))

*n* (**x**) =

*E*xc [*n*] = *E*<sup>x</sup> [*n*] +

This assumption yields the Hartree-Fock-like second Kohn-Sham equation

*ψ*KS *<sup>i</sup>* (**x**) −

*n*1 **x**, **x**� = *N* ∑ *j*=1 *ψ*KS

*N* ∑ *i*=1 *<sup>ψ</sup>*KS *<sup>i</sup>* (**x**) 2

where *N* is the number of electrons. Since we have no knowledge on the density *n* (**x**), it should be obtained self-consistently starting from a guessed density profile *n* (**x**). We compute the Kohn-Sham orbitals *ψ*KS under the static potential *ϕ* (**x**) as well as the exchange-correlation potential *μ*xc (*n* (**x**)) and then evaluate the density from Eq. (185). Once we obtain the self-consistent density profile *n* (**x**), the energy of the system is obtained by Eq. (181), if we have an *explicit* functional form of *E*xc [*n*]. Indeed the knowledge of the exchange-correlation functional form is important to obtain the exchange-correlation potential during solving Eq. (184). The cost of computation Eq. (184) is only that of the Hartree equation times the number of self-consistent iterations. The cheap Kohn-Sham equation have get superiority in solving many-body problems as developing computer technology (Argaman & Makov, 2000). The second approximation by Kohn & Sham (1965) is based on the assumption that

*n*<sup>1</sup> (**x**, **x**�

*<sup>j</sup>* (**x**) *<sup>ψ</sup>*∗KS *j* **x**� 

This procedure is similar to the Hartree-Fock method with the local correlation effect, and its cost of computation is comparable to the computational labor intensive Hartree-Fock method due to its nonlocal density matrix operation on the Kohn-Sham orbitals *ψ<sup>i</sup>* (**x**) appeared in

) <sup>|</sup>**<sup>x</sup>** <sup>−</sup> **<sup>x</sup>**�| *<sup>ψ</sup>*KS *i* **x**� 

single-particle Schrödinger-like first *Kohn-Sham* equation

 −1 2

<sup>∇</sup><sup>2</sup> <sup>+</sup> *<sup>ϕ</sup>* (**x**) <sup>+</sup> *<sup>μ</sup>*<sup>c</sup> (*<sup>n</sup>* (**x**))

where *μc* (**x**) is the correlation potential and

By assuming

and setting

 −1 2

*δn* leads us to have an equation

density perturbation, at the given pressure and temperature. In the mean time, the remained fluctuation modes will be suppressed. This state with the density distribution of a specific spatial regularity, represented by *qL*, is the *crystal*. Hence, the denominator of Eq. (178) explains the critical properties of the solidification.

Once the liquid transformed to a solid, the description of the system state should be switched to a new language for solid states. Every crystalline solid is categorized by its three dimensional translation and rotation symmetries to be the one of the 230 space groups (Bravais, 1845). Bloch (1928) proved, by using the periodic boundary condition (Born & von Kármán, 1912; 1913) and the group representation theory (Wigner, 1927), that the electronic wave function in a crystalline solid should be in the form

$$
\psi\_{n\mathbf{k}}\left(\mathbf{x}\right) = e^{i\mathbf{k}\cdot\mathbf{R}} u\_{n\mathbf{k}}\left(\mathbf{x}\right) \,,\tag{179}
$$

where *n* is the band index and **R** is the Bravais lattice vector. The Bloch function has three consequencies (Aschcroft & Mermin, 1976): (1) The momentum eigenvalue ¯*h***k** of the momenum operator **p** = (*h*¯ /*i*) ∇ does not conserve, instead the crystal momentum ¯*h***k** + *ei***k**·**<sup>x</sup>** *<sup>h</sup>*¯ *<sup>i</sup>* ∇ conserves. (2) The wave vector **k** is confined to the first Brillouin zone, so that the index *n* appear in Eq. (179) represent the multiplicity of the **k** confinements. (3) The nonvanishing velocity of a particle in a level specified by band index *n* and wave vector **k** is given by

$$\mathbf{v}\_{\rm ll}\left(\mathbf{k}\right) = \frac{1}{\hbar} \nabla\_{\mathbf{k}} \varepsilon\_{\rm ll}\left(\mathbf{k}\right),\tag{180}$$

where ∇**<sup>k</sup>** is a gradient operator with respect to the wave vector **k**. Further consequencies can be found in many solid-state physics text books (Altmann, 1995; Aschcroft & Mermin, 1976; Callaway, 1974; Jones & March, 1973a;b; Kittel, 2005; Madelung, 1978).

On the other hand, the numerator, 1/*ξ*<sup>2</sup> in Eq. (178) explains another aspect of the solidificaiton. The numerator describes the width in the wave vector space of the divergying static mean density fluctuation *S* (*q*). As *q* → *qL*, the wave vectors in a narow range around *qL* also contributes to *S* (*q*). When two waves of near wave vector (or frequency) superposes in the same position, we observe a very well known phenomenon, *beats* (Landau & Lifshitz, 1978). The broadening in *S* (*q*) implies that there are many wave vectors around superposing in the space to exhibits a complex pattern of spatial beats during the crystallization process. In other words, the crystallization will not happen in whole the liquid simultaneously, but the crystal nucleates.

Note that the numerator 1/*ξ*<sup>2</sup> contains the details of the interaction of the liquid. As we seen in Sec. 5.2, the dissipation of the external perturbation or the system response is determined by the interal fluctuation of the system. The internal fluctuation is determined by the details of the interaction of the system. Hence, we need to have the detailed interaction information of the system also for understanding phase transition phenomena completely.

#### **6.2 Kohn-Sham equations**

Kohn & Sham (1965) suggested two approximation methods for treating an inhomogeneous system of interacting electrons. Let us describe the first approximation first. The ground-state energy of the system in a static potential *v* (**x**) can be written in the form, in Hartree atomic unit ¯*h* = *e*<sup>2</sup> = *me* = 1,

$$E = T\_0 \left[ n \right] + \int v \left( \mathbf{x} \right) n \left( \mathbf{x} \right) \mathbf{d}^3 \mathbf{x} + \frac{1}{2} \iint \frac{n \left( \mathbf{x} \right) n \left( \mathbf{x}' \right)}{\left| \mathbf{x} - \mathbf{x}' \right|} \mathbf{d}^3 \mathbf{x} \mathbf{d}^3 \mathbf{x}' + E\_{\mathbf{xc}} \left[ n \right], \tag{181}$$

where *T*<sup>0</sup> is the kinetic energy functional of a system of noninteracting electrons and *E*xc [*n*] is the exchange and correlation energy functional of an interacting system with density *n* (**x**). By assuming

$$E\_{\mathbf{x}\mathbf{c}}\left[n\right] = \int n\left(\mathbf{x}\right) \epsilon\_{\mathbf{x}\mathbf{c}}\left(n\left(\mathbf{x}\right)\right) \mathbf{d}^3 \mathbf{x},\tag{182}$$

the stationary condition of the energy functional Eq. (181) with respect to the density variation *δn* leads us to have an equation

$$\frac{\delta T\_0\left[n\right]}{\delta n\left(\mathbf{x}\right)} + v\left(\mathbf{x}\right) + \mu\_{\text{xc}}\left(n\left(\mathbf{x}\right)\right) = 0,\tag{183}$$

where *v* (**x**) is the result of the functional derivatives of the static integration terms of Eq. (181) and *μ*xc is the result of the functional derivatives of the exchange-correlation functional *E*xc [*n*]. By assuming *μ*xc (*n* (**x**)) to be local, we can imagine a "ficticious" particle that moves by the single-particle Schrödinger-like first *Kohn-Sham* equation

$$\left(-\frac{1}{2}\nabla^2 + \varrho\left(\mathbf{x}\right) + \mu\_{\text{xc}}\left(n\left(\mathbf{x}\right)\right)\right)\psi\_{\text{i}}^{\text{KS}}\left(\mathbf{x}\right) = \epsilon\_{\text{i}}^{\text{KS}}\psi\_{\text{i}}^{\text{KS}}\left(\mathbf{x}\right),\tag{184}$$

and setting

32 Will-be-set-by-IN-TECH

density perturbation, at the given pressure and temperature. In the mean time, the remained fluctuation modes will be suppressed. This state with the density distribution of a specific spatial regularity, represented by *qL*, is the *crystal*. Hence, the denominator of Eq. (178)

Once the liquid transformed to a solid, the description of the system state should be switched to a new language for solid states. Every crystalline solid is categorized by its three dimensional translation and rotation symmetries to be the one of the 230 space groups (Bravais, 1845). Bloch (1928) proved, by using the periodic boundary condition (Born & von Kármán, 1912; 1913) and the group representation theory (Wigner, 1927), that

where *n* is the band index and **R** is the Bravais lattice vector. The Bloch function has three consequencies (Aschcroft & Mermin, 1976): (1) The momentum eigenvalue ¯*h***k** of the momenum operator **p** = (*h*¯ /*i*) ∇ does not conserve, instead the crystal momentum ¯*h***k** +

*h*¯

where ∇**<sup>k</sup>** is a gradient operator with respect to the wave vector **k**. Further consequencies can be found in many solid-state physics text books (Altmann, 1995; Aschcroft & Mermin, 1976;

On the other hand, the numerator, 1/*ξ*<sup>2</sup> in Eq. (178) explains another aspect of the solidificaiton. The numerator describes the width in the wave vector space of the divergying static mean density fluctuation *S* (*q*). As *q* → *qL*, the wave vectors in a narow range around *qL* also contributes to *S* (*q*). When two waves of near wave vector (or frequency) superposes in the same position, we observe a very well known phenomenon, *beats* (Landau & Lifshitz, 1978). The broadening in *S* (*q*) implies that there are many wave vectors around superposing in the space to exhibits a complex pattern of spatial beats during the crystallization process. In other words, the crystallization will not happen in whole the liquid simultaneously, but the

Note that the numerator 1/*ξ*<sup>2</sup> contains the details of the interaction of the liquid. As we seen in Sec. 5.2, the dissipation of the external perturbation or the system response is determined by the interal fluctuation of the system. The internal fluctuation is determined by the details of the interaction of the system. Hence, we need to have the detailed interaction information

Kohn & Sham (1965) suggested two approximation methods for treating an inhomogeneous system of interacting electrons. Let us describe the first approximation first. The ground-state energy of the system in a static potential *v* (**x**) can be written in the form, in Hartree atomic

> 1 2

*n* (**x**) *n* (**x**�

)

<sup>|</sup>**<sup>x</sup>** <sup>−</sup> **<sup>x</sup>**�| <sup>d</sup>3**x**d3**x**� <sup>+</sup> *<sup>E</sup>*xc [*n*] , (181)

of the system also for understanding phase transition phenomena completely.

*v* (**x**) *n* (**x**) d3**x** +

*<sup>i</sup>* ∇ conserves. (2) The wave vector **k** is confined to the first Brillouin zone, so that the index *n* appear in Eq. (179) represent the multiplicity of the **k** confinements. (3) The nonvanishing velocity of a particle in a level specified by band index *n* and wave vector **k** is

*<sup>i</sup>***k**·**<sup>R</sup>***un***<sup>k</sup>** (**x**), (179)

∇**k***�<sup>n</sup>* (**k**), (180)

explains the critical properties of the solidification.

*ei***k**·**<sup>x</sup>** *<sup>h</sup>*¯

given by

crystal nucleates.

**6.2 Kohn-Sham equations**

*E* = *T*<sup>0</sup> [*n*] +

unit ¯*h* = *e*<sup>2</sup> = *me* = 1,

the electronic wave function in a crystalline solid should be in the form

*ψn***<sup>k</sup>** (**x**) = *e*

**<sup>v</sup>***<sup>n</sup>* (**k**) <sup>=</sup> <sup>1</sup>

Callaway, 1974; Jones & March, 1973a;b; Kittel, 2005; Madelung, 1978).

$$m\left(\mathbf{x}\right) = \sum\_{i=1}^{N} \left| \psi\_i^{\text{KS}}\left(\mathbf{x}\right) \right|^2,\tag{185}$$

where *N* is the number of electrons. Since we have no knowledge on the density *n* (**x**), it should be obtained self-consistently starting from a guessed density profile *n* (**x**). We compute the Kohn-Sham orbitals *ψ*KS under the static potential *ϕ* (**x**) as well as the exchange-correlation potential *μ*xc (*n* (**x**)) and then evaluate the density from Eq. (185). Once we obtain the self-consistent density profile *n* (**x**), the energy of the system is obtained by Eq. (181), if we have an *explicit* functional form of *E*xc [*n*]. Indeed the knowledge of the exchange-correlation functional form is important to obtain the exchange-correlation potential during solving Eq. (184). The cost of computation Eq. (184) is only that of the Hartree equation times the number of self-consistent iterations. The cheap Kohn-Sham equation have get superiority in solving many-body problems as developing computer technology (Argaman & Makov, 2000). The second approximation by Kohn & Sham (1965) is based on the assumption that

$$E\_{\mathbf{x}\mathbf{c}}\left[n\right] = E\_{\mathbf{x}}\left[n\right] + \int n\left(\mathbf{x}\right)\varepsilon\_{\mathbf{c}}\left(n\left(\mathbf{x}\right)\right)\mathbf{d}^3\mathbf{x}.\tag{186}$$

This assumption yields the Hartree-Fock-like second Kohn-Sham equation

$$\left(-\frac{1}{2}\nabla^2 + \boldsymbol{\varrho}\left(\mathbf{x}\right) + \boldsymbol{\mu}\_{\mathbb{C}}\left(\boldsymbol{n}\left(\mathbf{x}\right)\right)\right)\boldsymbol{\psi}\_{l}^{\text{KS}}\left(\mathbf{x}\right) - \int \frac{n\_{1}\left(\mathbf{x},\mathbf{x}'\right)}{\left|\mathbf{x}-\mathbf{x}'\right|}\boldsymbol{\psi}\_{l}^{\text{KS}}\left(\mathbf{x}'\right)\,\mathrm{d}^3\mathbf{x}' = \boldsymbol{\epsilon}\_{l}^{\text{KS}}\boldsymbol{\psi}\_{l}^{\text{KS}}\left(\mathbf{x}\right),\tag{187}$$

where *μc* (**x**) is the correlation potential and

$$m\_1\left(\mathbf{x}, \mathbf{x}'\right) = \sum\_{j=1}^{N} \boldsymbol{\psi}\_j^{\mathbf{KS}}\left(\mathbf{x}\right) \boldsymbol{\psi}\_j^{\ast \mathbf{KS}}\left(\mathbf{x}'\right). \tag{188}$$

This procedure is similar to the Hartree-Fock method with the local correlation effect, and its cost of computation is comparable to the computational labor intensive Hartree-Fock method due to its nonlocal density matrix operation on the Kohn-Sham orbitals *ψ<sup>i</sup>* (**x**) appeared in

we will have a Lehmann representation of the real-time Green's function for real *ω* such as

Towards the Authentic *Ab Intio* Thermodynamics 577

3 *δ* 

*<sup>ω</sup>* <sup>−</sup> *<sup>h</sup>*¯ <sup>−</sup><sup>1</sup> (*Kn* <sup>−</sup> *Km*)

<sup>−</sup>*βKm* (2*π*)

 ∞ −∞

> 3 *δ*

*<sup>ω</sup>* <sup>−</sup> *<sup>h</sup>*¯ <sup>−</sup><sup>1</sup> (*Kn* <sup>−</sup> *Km*)

Is it able to calculate the eigenvalues of the interacting system at a finite temperature? Let us come back to the idea of Born & Oppenheimer (1927) adiabatic separation of the motions of electrons and ions. In the Born-Oppenheimer scheme, we first calculate the motions of electrons in a given ionic background and then calculate the motions of ions at the given electronic one. It is well konw that the given symmetry and many-body effects will emerges as elementary excitations (Pines, 1999), which can be extended even to a *macroscopic* levels (Anderson, 1997). Some of the elementary excitations may originated from either electrons or ions. For example, plasoms are the elementary excitations of collective motions of electrons, while phonons are the elementary excitation of the collective motions of ions. One can also interesting to more complex excitations, which are composites of those elementary excitations. Anyhow, all such properties are computed from the electronic structures of the

<sup>−</sup>*βKm* (2*π*)

*<sup>ω</sup>* <sup>−</sup> *<sup>h</sup>*¯ <sup>−</sup><sup>1</sup> (*Kn* <sup>−</sup> *Km*)

where *s* is the spin degeneracy and ℘ denotes a principal value. The imaginary part of Eq.

3 *δ* 

d*ω*� *π*

Considering the spectral representation form in Eq. (135) we have the detailed spectral weight

Since the spectral weight function *A* (**k**, *ω*) relates the real-time Green's function and the temperature Green's function, we can, in principle, solve all the thermodynamic problems with detailed interaction information of the system: The thermal Green's function describes the equation of states, while the imaginary part of the retarded Green's function from the real-time Green's funciton yields the fluctuations to describe the phase transtion, including the transport properties. The remained difficulty is to calculate the eigenvalues *Km* and **P***m* in

<sup>1</sup> <sup>±</sup> *<sup>e</sup>*

�*G*¯ (**k**, *<sup>ω</sup>*�

*ω* − *ω*�

**<sup>k</sup>** <sup>−</sup> *<sup>h</sup>*¯ <sup>−</sup><sup>1</sup> (**P***<sup>n</sup>* <sup>−</sup> **<sup>P</sup>***m*)

**<sup>k</sup>** <sup>−</sup> *<sup>h</sup>*¯ <sup>−</sup><sup>1</sup> (**P***<sup>n</sup>* <sup>−</sup> **<sup>P</sup>***m*)

<sup>1</sup> <sup>∓</sup> *<sup>e</sup>*−*β*(*Kn*−*Km*)

<sup>1</sup> <sup>±</sup> *<sup>e</sup>*

**<sup>k</sup>** <sup>−</sup> *<sup>h</sup>*¯ <sup>−</sup><sup>1</sup> (**P***<sup>n</sup>* <sup>−</sup> **<sup>P</sup>***m*)

)

<sup>1</sup> <sup>∓</sup> *<sup>e</sup>*−*β*(*Kn*−*Km*)

−*β*(*Kn*−*Km*)

−*β*(*Kn*−*Km*)

<sup>1</sup> <sup>∓</sup> *<sup>e</sup>*−*βω*�

 �*m*<sup>|</sup> *<sup>ψ</sup>*<sup>ˆ</sup>

�*m*<sup>|</sup> *<sup>ψ</sup>*ˆ*<sup>α</sup>* <sup>|</sup>*n*�

<sup>1</sup> <sup>±</sup> *<sup>e</sup>*−*βω*� . (195)

*<sup>α</sup>* |*n*� 

<sup>2</sup> . (196)

�*m*<sup>|</sup> *<sup>ψ</sup>*ˆ*<sup>α</sup>* <sup>|</sup>*n*�

 2

, (193)

 2

. (194)

<sup>−</sup>*βKm* (2*π*)

1 *<sup>ω</sup>* <sup>−</sup> *<sup>h</sup>*¯ <sup>−</sup><sup>1</sup> (*Kn* <sup>−</sup> *Km*)

*<sup>G</sup>*¯ (**k**, *<sup>ω</sup>*) <sup>=</sup> *<sup>e</sup>β*<sup>Ω</sup>

�*G*¯ (**k**, *<sup>ω</sup>*) <sup>=</sup> <sup>−</sup> *<sup>π</sup>eβ*<sup>Ω</sup>

Landau (1958) proved a relation

*<sup>A</sup>* (**k**, *<sup>ω</sup>*) <sup>=</sup> *<sup>e</sup>β*<sup>Ω</sup>

(193) is written as

function of form

system.

<sup>2</sup>*<sup>s</sup>* <sup>+</sup> <sup>1</sup> ∑ *m*,*n e*

<sup>−</sup>*iπδ*

<sup>2</sup>*<sup>s</sup>* <sup>+</sup> <sup>1</sup> ∑ *m*,*n e*

�*G*¯ (**k**, *<sup>ω</sup>*) <sup>=</sup> <sup>−</sup><sup>℘</sup>

×*δ* 

<sup>2</sup>*<sup>s</sup>* <sup>+</sup> <sup>1</sup> ∑ *m*,*n e*

Eqs. (189)–(190) in the case of interacting systems.

<sup>×</sup>2*πδ*

× ℘

Eq. (187). Hence, the application of the Kohn-Sham equation of the type in Eq. (187) has been hesitated. Although the practical success of the first Kohn-Sham equation in Eq. (184), the approximations are remained uncontrollable, because there is no known explicit exchange-correlation functional. Recent developments of computer technology, many approaches has been tried to utilize the second type Kohn-Sham equation (Kümmel & Kronik, 2008). This kind of approaches has solved many problems in the first-type Kohn-Sham equation.

However, there is a critical issue on the solutions of the Kohn-Sham equations. Koopmans (1934) proved that the *i*th occupied eigenvalues of the Hartree-Fock equation is identified as the energy required to remove an electron from that orbital without perturbing the rest of the system. This theorem is not valid in the density functional theory, in general. Only the highest occupied eigenvalue of the Kohn-Sham equations has a rigorous physical meaning that the partial derivative is exactly the chemical potential of the system (Almbladh & von Barth, 1985; Levy et al., 1984). In other words, the Kohn-Sham orbitals *ψ*KS *<sup>i</sup>* (**x**) except the highest occupied orbital are physically meaningless, but they are useful only for building the density in Eq. (185). Indeed, the Kohn-Sham charge density generated by Eq. (185) differs hardly to that of the Hartree-Fock charge density (Stowasser & Hoffmann, 1999). A quantitative linear scailing interpretation of Kohn-Sham eigenvalues with respect to corresponding Hartree-Fock eigenvalues was suggested (Stowasser & Hoffmann, 1999).

#### **6.3 Spectral weight function and quasiparticle spectra**

As seen in Sec. 6.1, the detailed interaction information of the system is equally important to understand the phase transition as well as the equilibrium properties of the system. In the density functional theoretical framework, the Kohn-Sham equations have a critial problem of the unknown physical meanings of their eigenvalues. The resultant Kohn-Sham orbitals ensure the ground state energy, the density *n* (**x**), and the chemical potential. The other information obtained by using the Kohn-Sham orbitals are *fictitious*. It is fortunate or unfortunate that the Kohn-Sham orbitals are close to the quasiparticle orbitals, which have the concrete physical meaning through the Koopmans' theorem (Koopmans, 1934). Since there is no further proof on the physical meaning of the Kohn-Sham orbitals, a promising way to build the physically seamless theoretical framework for phase diagrams is starting from the quantum field theory described in Sec. 4, armed with the Kubo's framework in Sec. 5.1–5.2.

The Hamiltonian operator *H*ˆ , the momentum operator **P**ˆ, and the number opertator *N*ˆ commute each other. By symmetry, those operators have a set of eigenstates {|*m*�} simultaneously (Tinkam, 1964), such that

$$
\hat{\mathbf{P}}\left|m\right> = \mathbf{P}\_m\left|m\right>\,\tag{189}
$$

$$\left(\hat{H} - \mu \hat{N}\right)|m\rangle = \left(E\_m - \mu N\_m\right)|m\rangle. \tag{190}$$

Note that Eq. (190) is understood as

$$
\hat{\mathcal{K}}\left|m\right> = \mathcal{K}\_m\left|m\right>\,. \tag{191}
$$

With the aid of the Sokhotsky-Weierstrass theorem (Byron & Fuller, 1992),

$$\frac{1}{\omega \pm i\eta} = \wp \frac{1}{\omega} \mp i\pi\delta\left(\omega\right),\tag{192}$$

34 Will-be-set-by-IN-TECH

Eq. (187). Hence, the application of the Kohn-Sham equation of the type in Eq. (187) has been hesitated. Although the practical success of the first Kohn-Sham equation in Eq. (184), the approximations are remained uncontrollable, because there is no known explicit exchange-correlation functional. Recent developments of computer technology, many approaches has been tried to utilize the second type Kohn-Sham equation (Kümmel & Kronik, 2008). This kind of approaches has solved many problems in the first-type Kohn-Sham

However, there is a critical issue on the solutions of the Kohn-Sham equations. Koopmans (1934) proved that the *i*th occupied eigenvalues of the Hartree-Fock equation is identified as the energy required to remove an electron from that orbital without perturbing the rest of the system. This theorem is not valid in the density functional theory, in general. Only the highest occupied eigenvalue of the Kohn-Sham equations has a rigorous physical meaning that the partial derivative is exactly the chemical potential of the system (Almbladh & von Barth,

occupied orbital are physically meaningless, but they are useful only for building the density in Eq. (185). Indeed, the Kohn-Sham charge density generated by Eq. (185) differs hardly to that of the Hartree-Fock charge density (Stowasser & Hoffmann, 1999). A quantitative linear scailing interpretation of Kohn-Sham eigenvalues with respect to corresponding Hartree-Fock

As seen in Sec. 6.1, the detailed interaction information of the system is equally important to understand the phase transition as well as the equilibrium properties of the system. In the density functional theoretical framework, the Kohn-Sham equations have a critial problem of the unknown physical meanings of their eigenvalues. The resultant Kohn-Sham orbitals ensure the ground state energy, the density *n* (**x**), and the chemical potential. The other information obtained by using the Kohn-Sham orbitals are *fictitious*. It is fortunate or unfortunate that the Kohn-Sham orbitals are close to the quasiparticle orbitals, which have the concrete physical meaning through the Koopmans' theorem (Koopmans, 1934). Since there is no further proof on the physical meaning of the Kohn-Sham orbitals, a promising way to build the physically seamless theoretical framework for phase diagrams is starting from the quantum field theory described in Sec. 4, armed with the Kubo's framework in Sec. 5.1–5.2. The Hamiltonian operator *H*ˆ , the momentum operator **P**ˆ, and the number opertator *N*ˆ commute each other. By symmetry, those operators have a set of eigenstates {|*m*�}

**<sup>P</sup>**<sup>ˆ</sup> <sup>|</sup>*m*� <sup>=</sup> **<sup>P</sup>***<sup>m</sup>* <sup>|</sup>*m*�, (189)


<sup>K</sup><sup>ˆ</sup> <sup>|</sup>*m*� <sup>=</sup> *Km* <sup>|</sup>*m*�. (191)

∓ *iπδ* (*ω*), (192)

*<sup>i</sup>* (**x**) except the highest

1985; Levy et al., 1984). In other words, the Kohn-Sham orbitals *ψ*KS

eigenvalues was suggested (Stowasser & Hoffmann, 1999).

 *<sup>H</sup>*<sup>ˆ</sup> <sup>−</sup> *<sup>μ</sup>N*<sup>ˆ</sup> 

With the aid of the Sokhotsky-Weierstrass theorem (Byron & Fuller, 1992), 1 *<sup>ω</sup>* <sup>±</sup> *<sup>i</sup><sup>η</sup>* <sup>=</sup> <sup>℘</sup>

1 *ω*

**6.3 Spectral weight function and quasiparticle spectra**

simultaneously (Tinkam, 1964), such that

Note that Eq. (190) is understood as

equation.

we will have a Lehmann representation of the real-time Green's function for real *ω* such as

$$\tilde{G}\left(\mathbf{k},\omega\right) = \frac{e^{\beta\Omega}}{2s+1} \sum\_{m,n} e^{-\beta K\_m} \left(2\pi\right)^3 \delta\left[\mathbf{k} - \hbar^{-1} \left(\mathbf{P}\_n - \mathbf{P}\_m\right)\right] \left|\left\right|^2$$

$$\times \left\{\wp \frac{1}{\omega - \hbar^{-1} \left(K\_n - K\_m\right)} \left(1 \mp e^{-\beta\left(K\_n - K\_m\right)}\right)\right.$$

$$-i\pi\delta \left[\omega - \hbar^{-1} \left(K\_n - K\_m\right)\right] \left(1 \pm e^{-\beta\left(K\_n - K\_m\right)}\right)\right\},\tag{193}$$

where *s* is the spin degeneracy and ℘ denotes a principal value. The imaginary part of Eq. (193) is written as

$$\mathbb{S}\mathbb{G}\left(\mathbf{k},\omega\right) = -\frac{\pi e^{\beta\Omega}}{2s+1} \sum\_{m,n} e^{-\beta\mathbf{K}\_m} \left(2\pi\right)^3 \delta\left[\mathbf{k} - \hbar^{-1} \left(\mathbf{P}\_n - \mathbf{P}\_m\right)\right] \left|\left\right|^2$$

$$\times \delta\left[\omega - \hbar^{-1} \left(\mathbf{K}\_n - \mathbf{K}\_m\right)\right] \left(1 \pm e^{-\beta(\mathbf{K}\_n - \mathbf{K}\_m)}\right). \tag{194}$$

Landau (1958) proved a relation

$$\mathfrak{R}\ddot{\mathbf{G}}\left(\mathbf{k},\omega\right) = -\wp \int\_{-\infty}^{\infty} \frac{\mathbf{d}\omega'}{\pi} \frac{\mathfrak{J}\ddot{\mathbf{G}}\left(\mathbf{k},\omega'\right)}{\omega - \omega'} \frac{1 \mp e^{-\beta\omega'}}{1 \pm e^{-\beta\omega'}}.\tag{195}$$

Considering the spectral representation form in Eq. (135) we have the detailed spectral weight function of form

$$A\left(\mathbf{k},\omega\right) = \frac{e^{\beta\Omega}}{2s+1} \sum\_{m,n} e^{-\beta\mathbf{K}\_m} \left(2\pi\right)^3 \delta\left[\mathbf{k} - \hbar^{-1} \left(\mathbf{P}\_n - \mathbf{P}\_m\right)\right]$$

$$\times 2\pi\delta \left[\omega - \hbar^{-1} \left(\mathbf{K}\_n - \mathbf{K}\_m\right)\right] \left(1 \mp e^{-\beta(\mathbf{K}\_n - \mathbf{K}\_m)}\right) \left|\left\right|^2\right.\tag{196}$$

Since the spectral weight function *A* (**k**, *ω*) relates the real-time Green's function and the temperature Green's function, we can, in principle, solve all the thermodynamic problems with detailed interaction information of the system: The thermal Green's function describes the equation of states, while the imaginary part of the retarded Green's function from the real-time Green's funciton yields the fluctuations to describe the phase transtion, including the transport properties. The remained difficulty is to calculate the eigenvalues *Km* and **P***m* in Eqs. (189)–(190) in the case of interacting systems.

Is it able to calculate the eigenvalues of the interacting system at a finite temperature?

Let us come back to the idea of Born & Oppenheimer (1927) adiabatic separation of the motions of electrons and ions. In the Born-Oppenheimer scheme, we first calculate the motions of electrons in a given ionic background and then calculate the motions of ions at the given electronic one. It is well konw that the given symmetry and many-body effects will emerges as elementary excitations (Pines, 1999), which can be extended even to a *macroscopic* levels (Anderson, 1997). Some of the elementary excitations may originated from either electrons or ions. For example, plasoms are the elementary excitations of collective motions of electrons, while phonons are the elementary excitation of the collective motions of ions. One can also interesting to more complex excitations, which are composites of those elementary excitations. Anyhow, all such properties are computed from the electronic structures of the system.

Equating this relation to result of the operation of the gradient operator on Eq. (202) yields

Towards the Authentic *Ab Intio* Thermodynamics 579

This means that the motion of the quasiparticle of mass *m* may be thought as that of a free particle of mass *m*∗, of which all the many-body effects are corrected. Now we can call *m*∗ as the *effective mass* of the quasiparticle. The noninteracting quasiparticle Green's function may

where *Z* is called the *renormalization factor*, which contains the effective mass effect in Eq. (203).

In other words, the results of a photoemission spectroscopy are just the renormalization factor

The idea of renormalization was utilized by Landau (1957a;b; 1959), known as the phenomenological *Landau Fermi-liquid theory*. Landau classified the quasiparticles as the levels of a Fermi liquid corresponding to the classification of the levels of noninteracting particles. The interaction is slowly *turned on* to the particles gradually become quasiparticles. The methods described above can be directly applied for the quasiparticles. Using this concept, Landau (1957a;b; 1959) was able to calculate the macroscopic quantities, such as specific heat, entroy, magnetic susceptibility, kinetic motions, thermal conductivity, viscocity, dispersion and absorption of sound, and fluctuations in the distribution function. The relation with the quantum field theory of the Landau Fermi liquid theory was investigated by Nozières & Luttinger (1962) and Luttinger & Nozières (1962), but it is not clear yet. Good introductory reviews on the Ladau Fermi-liquid theory are given by Abrikosov & Khalatnikov

Bohm & Pines (1952) suggested a different point of view with a canonical transformation to treat a particle and its consequent collective excitations simultaneously. In this formalism, the canonical transformation stresses the elementary ingredients of the Hamiltonian are the particle-hole pairs, instead of individual particles. Bohm & Pines (1952) approximated the coupling between the excitations corresponding to different momentum transfers *q* and *q*�

This ignorance is originated from an intuition that the momentum transfer phase *q* − *q*� will be random to be enough for cancelling themselves. Hence, the method is known as the *random phase approximation (RPA)*. The RPA was applied to the electronic system with plasmons and refined (Sawada, 1957; Sawada et al., 1957). The RPA picture as the similarity with that of (Goldstone, 1957) and it deals with the correlation function corresponding to the collective

These features of the RPA have been utilized in susceptibility investigations. For instance, Kim et al. (1973) calculated a spin polarized susceptibility<sup>9</sup> based on the fluctuation-dissipation theorem and the free energy by considering higher order self-energy correction (Kim, 1974) withn the RPA. Kim (1982; 1989; 1999) used the RPA technique to demonstrate, within the homogeneous jeullium model, that all the thermodynamic properties, including the ferromagnetic phase transition, of the metallic ferromagnetism can be obtained by considering the spin-polarized electrons with the collective excitations of phonons and

*<sup>G</sup>*¯ (**k**, *<sup>ω</sup>*) <sup>=</sup> *<sup>Z</sup>*

Considering the spectral function relation with the Green's function, we have

of the system. The photoemission experiments gives the quasiparticle spectra.

(1959); Giuliani & Vignale (2005); v. Löhneysen et al. (2006).

<sup>9</sup> A practical self-consistent version was given by Lee et al. (1989).

*h*¯ ∇**k**Σ (**k**) *�*0 **k**

*<sup>ω</sup>* <sup>−</sup> *<sup>h</sup>*¯ <sup>−</sup><sup>1</sup> (*�***<sup>k</sup>** <sup>−</sup> *<sup>μ</sup>*)

*A* (**k**, *ω*) = *Zδ* (*ω* + *μ* − *h*¯ **k**). (205)

. (203)

, (204)

.

*m <sup>m</sup>*<sup>∗</sup> <sup>=</sup> <sup>1</sup> <sup>+</sup>

be in the form

excitations.

It is instructive to consider the formal form in Eq. (129) of the Dyson's equation. Since the Dyson's equation expresses the difference between the noninteracting and interacting temperature Green's functions in terms of the self-energy, Σ (**k**, *ωn*); all the many-body effects are contained in the self-energy. Considering the real space Dyson's equation form in Eq. (107), we may have a Hartree-Fock-like self-consistent equation (Fetter & Walecka, 2003)

$$
\left[-\frac{\hbar}{2m}\nabla\_1^2 + v\left(\mathbf{x}\_1\right)\right]\psi\_j\left(\mathbf{x}\_1\right) + \hbar \int \mathbf{d}^3 \mathbf{x}\_2 \Sigma\left(\mathbf{x}\_1, \mathbf{x}\_2\right)\psi\_j\left(\mathbf{x}\_2\right) = \epsilon\_j \psi\_j\left(\mathbf{x}\_1\right) \tag{197}
$$

where *v* (**x**1) is the static one-body potential to describe the *inhomogeneity* as in the density functional theory, <sup>∇</sup><sup>2</sup> <sup>1</sup> is the Laplacian operator with respect to the coordinate **x**1, and *μ* is the chemical potential obtained by

$$N\left(T\_{\prime}V\_{\prime}\mu\right) = \left(2s+1\right)\sum\_{j}n\_{j\prime}\tag{198}$$

with the number distribution function

$$m\_{\!\!\!\!-1} = \frac{1}{\mathcal{e}^{\beta\left(\epsilon\_{\!\!\!-\!\!\!-\!\!\!-\!\!\!-\!\!\!+\!\!\!\!\/)}}}.\tag{199}\tag{199}$$

The chemical potential *μ* (*T*, *V*, *N*) is, in principle, able to be inverted from Eq. (198) if *N* is fixed. In this derivation, we assumed that all the necessary frequency summed are already done for making Σ being depending explicitly on *T* and *μ*, but being independent to the frequency *ωn* of the Hartree-Fock form and the temperature Green's function

$$\mathcal{G}\left(\mathbf{x}\_{1},\mathbf{x}\_{2},\omega\_{n}\right) = \sum\_{\mathbf{j}} \frac{\psi\_{\mathbf{j}}\left(\mathbf{x}\_{1}\right)\psi\_{\mathbf{j}}^{\*}\left(\mathbf{x}\_{2}\right)}{\mathbf{i}\omega\_{n} - \hbar^{-1}\left(\boldsymbol{\varepsilon}\_{\boldsymbol{\hat{j}}} - \boldsymbol{\mu}\right)},\tag{200}$$

where *ψ<sup>j</sup>* is a single-particle wave function with energy *�<sup>j</sup>* for forming a complete orthonormal set, forms the same form of the noninteracting temperature Green's function <sup>G</sup><sup>0</sup> (**x**1, **<sup>x</sup>**2, *<sup>ω</sup>n*) with *ψ*<sup>0</sup> *<sup>j</sup>* , the solution of a differential equation

$$\left[i\hbar\omega\_n + \frac{\hbar}{2m}\nabla\_1^2 + \mu - v\left(\mathbf{x}\_1\right)\right]\psi\_j^0\left(\mathbf{x}\_1\right) = \epsilon\_i^0 \psi\_j^0\left(\mathbf{x}\_1\right). \tag{201}$$

In this self-consistent formalism, the quasiparticle energy spectrum becomes

$$
\boldsymbol{\varepsilon}\_{\mathbf{k}} = \boldsymbol{\varepsilon}\_{\mathbf{k}}^{0} + \hbar \boldsymbol{\Sigma} \left( \mathbf{k} \right),
\tag{202}
$$

where we recover the momentum space representation.

This equation tells us that the noninteracting quasiparticle moving with the energy spectrum *�*0 **<sup>k</sup>** *is dressed* by the self-energy Σ (**k**) by the sophisticated many-body interactions. The energy dispersion of a free particle of mass *m*∗ will have a dispersion

$$
\nabla\_{\mathbf{k}} \epsilon\_{\mathbf{k}} = \frac{\mathbf{k}}{m^\*}.
$$

36 Will-be-set-by-IN-TECH

It is instructive to consider the formal form in Eq. (129) of the Dyson's equation. Since the Dyson's equation expresses the difference between the noninteracting and interacting temperature Green's functions in terms of the self-energy, Σ (**k**, *ωn*); all the many-body effects are contained in the self-energy. Considering the real space Dyson's equation form in Eq. (107), we may have a Hartree-Fock-like self-consistent equation (Fetter & Walecka, 2003)

where *v* (**x**1) is the static one-body potential to describe the *inhomogeneity* as in the density

*N* (*T*, *V*, *μ*) = (2*s* + 1)∑

*nj* <sup>=</sup> <sup>1</sup>

frequency *ωn* of the Hartree-Fock form and the temperature Green's function

G (**x**1, **<sup>x</sup>**2, *<sup>ω</sup>n*) = ∑

*<sup>j</sup>* , the solution of a differential equation

where we recover the momentum space representation.

dispersion of a free particle of mass *m*∗ will have a dispersion

*h*¯ 2*m* ∇2

 *ih*¯ *ω<sup>n</sup>* + *<sup>e</sup>β*(*�j*−*μ*) <sup>∓</sup> <sup>1</sup>

*ψ<sup>j</sup>* (**x**1) *ψ*<sup>∗</sup>

*<sup>i</sup>ω<sup>n</sup>* <sup>−</sup> *<sup>h</sup>*¯ <sup>−</sup><sup>1</sup>

 *ψ*0 *<sup>j</sup>* (**x**2)

 *�<sup>j</sup>* − *μ*

*<sup>j</sup>* (**x**1) <sup>=</sup> *�*<sup>0</sup>

*i ψ*0

**<sup>k</sup>** + *h*¯Σ (**k**), (202)

The chemical potential *μ* (*T*, *V*, *N*) is, in principle, able to be inverted from Eq. (198) if *N* is fixed. In this derivation, we assumed that all the necessary frequency summed are already done for making Σ being depending explicitly on *T* and *μ*, but being independent to the

*j*

<sup>1</sup> + *μ* − *v* (**x**1)

In this self-consistent formalism, the quasiparticle energy spectrum becomes

*�***<sup>k</sup>** = *�*<sup>0</sup>

where *ψ<sup>j</sup>* is a single-particle wave function with energy *�<sup>j</sup>* for forming a complete orthonormal set, forms the same form of the noninteracting temperature Green's function <sup>G</sup><sup>0</sup> (**x**1, **<sup>x</sup>**2, *<sup>ω</sup>n*)

This equation tells us that the noninteracting quasiparticle moving with the energy spectrum

**<sup>k</sup>** *is dressed* by the self-energy Σ (**k**) by the sophisticated many-body interactions. The energy

<sup>∇</sup>**k***�***<sup>k</sup>** <sup>=</sup> **<sup>k</sup>**

*m*∗ .

<sup>1</sup> is the Laplacian operator with respect to the coordinate **x**1, and *μ* is the

*j*

d3**x**2Σ (**x**1, **x**2) *ψ<sup>j</sup>* (**x**2) = *�jψ<sup>j</sup>* (**x**1), (197)

*nj*, (198)

. (199)

, (200)

*<sup>j</sup>* (**x**1). (201)

 − *h*¯ 2*m* ∇2

chemical potential obtained by

with the number distribution function

functional theory, <sup>∇</sup><sup>2</sup>

with *ψ*<sup>0</sup>

*�*0

<sup>1</sup> + *v* (**x**1)

*ψ<sup>j</sup>* (**x**1) + *h*¯

Equating this relation to result of the operation of the gradient operator on Eq. (202) yields

$$\frac{m}{m^\*} = 1 + \frac{\hbar \nabla\_\mathbf{k} \Sigma\left(\mathbf{k}\right)}{\epsilon\_\mathbf{k}^0}.\tag{203}$$

This means that the motion of the quasiparticle of mass *m* may be thought as that of a free particle of mass *m*∗, of which all the many-body effects are corrected. Now we can call *m*∗ as the *effective mass* of the quasiparticle. The noninteracting quasiparticle Green's function may be in the form

$$\bar{G}\left(\mathbf{k},\omega\right) = \frac{Z}{\omega - \hbar^{-1}\left(\varepsilon\_{\mathbf{k}} - \mu\right)},\tag{204}$$

where *Z* is called the *renormalization factor*, which contains the effective mass effect in Eq. (203). Considering the spectral function relation with the Green's function, we have

$$A\left(\mathbf{k},\omega\right) = Z\delta\left(\omega + \mu - \hbar \mathbf{k}\right). \tag{205}$$

In other words, the results of a photoemission spectroscopy are just the renormalization factor of the system. The photoemission experiments gives the quasiparticle spectra.

The idea of renormalization was utilized by Landau (1957a;b; 1959), known as the phenomenological *Landau Fermi-liquid theory*. Landau classified the quasiparticles as the levels of a Fermi liquid corresponding to the classification of the levels of noninteracting particles. The interaction is slowly *turned on* to the particles gradually become quasiparticles. The methods described above can be directly applied for the quasiparticles. Using this concept, Landau (1957a;b; 1959) was able to calculate the macroscopic quantities, such as specific heat, entroy, magnetic susceptibility, kinetic motions, thermal conductivity, viscocity, dispersion and absorption of sound, and fluctuations in the distribution function. The relation with the quantum field theory of the Landau Fermi liquid theory was investigated by Nozières & Luttinger (1962) and Luttinger & Nozières (1962), but it is not clear yet. Good introductory reviews on the Ladau Fermi-liquid theory are given by Abrikosov & Khalatnikov (1959); Giuliani & Vignale (2005); v. Löhneysen et al. (2006).

Bohm & Pines (1952) suggested a different point of view with a canonical transformation to treat a particle and its consequent collective excitations simultaneously. In this formalism, the canonical transformation stresses the elementary ingredients of the Hamiltonian are the particle-hole pairs, instead of individual particles. Bohm & Pines (1952) approximated the coupling between the excitations corresponding to different momentum transfers *q* and *q*� . This ignorance is originated from an intuition that the momentum transfer phase *q* − *q*� will be random to be enough for cancelling themselves. Hence, the method is known as the *random phase approximation (RPA)*. The RPA was applied to the electronic system with plasmons and refined (Sawada, 1957; Sawada et al., 1957). The RPA picture as the similarity with that of (Goldstone, 1957) and it deals with the correlation function corresponding to the collective excitations.

These features of the RPA have been utilized in susceptibility investigations. For instance, Kim et al. (1973) calculated a spin polarized susceptibility<sup>9</sup> based on the fluctuation-dissipation theorem and the free energy by considering higher order self-energy correction (Kim, 1974) withn the RPA. Kim (1982; 1989; 1999) used the RPA technique to demonstrate, within the homogeneous jeullium model, that all the thermodynamic properties, including the ferromagnetic phase transition, of the metallic ferromagnetism can be obtained by considering the spin-polarized electrons with the collective excitations of phonons and

<sup>9</sup> A practical self-consistent version was given by Lee et al. (1989).

called *πn* (*R*, *x*). The fundamental group is called the *first homotopy group* and *πn* is then called the *n*th homotopy group. For example, *π*<sup>2</sup> (*S*2, *x*) is called when the maps *f* at the base point *x* is on the sphere *S*, *i.e.*, say, *u* and *v* into *R*, of three dimensional order-parameter space. For an integer *Z*, if *π*<sup>1</sup> (*S*1) = *Z*, the group structure does not depend on the choice of base point *x* and the map is topologically the same as the circle. Any mapping of a loop into a sphere can be shrunk to a point, so that the fundamental group of the surface of a sphere consists only of the identity, *π*<sup>1</sup> (*S*1) = 0. Trivial point defects are associated with mappings that can be

Towards the Authentic *Ab Intio* Thermodynamics 581

Let us say the group *G* containing translations as well as proper rotations, *i.e.* the proper part of the full Euclidean group. Let *H* be the subgroup of *G* containing those rigid body rotations that leave the reference system invariant. Naively the order-parameter space *R* is the coset space *G*/*H*. If we ignore of the rotational symmetries of crystals, the full proper Euclidean group can be replaced by the subgroup *T* (3) of translations. On the other hand the isotropy subgroup *H* becomes the subgroup of *T* (3) consisting of translations through Bravais lattice vectors. Since *T* (3) is parameterized by all of Euclidean three-space, it is connected and simply connected. Since *H* is discrete, we identifies *π*<sup>1</sup> (*G*/*H*) with *H* itself. Thus the line defects, in our case dislocations, are characterized by Bravais lattice vectors. The vector is usually known as the Burgers vectors. Because *H* is discrete *π*<sup>2</sup> (*G*/*H*) = 0. In this case, the homotopy theory says notopologically stable *point* defects. However, we know thermodynamically stable point defects, *e.g.*, vacancies or voids. The physical means of removing a void are reminiscent of the nonlocal means available for the elimination of

Now let us recover the operations of the rigid body rotations. The fundamental group is identified with the double group. The three dimensional rotational parts of the group operations are lifted from SO(3) rotation to SU(2). Discrete operations in the space group with no translational part correspond to line defects, in which the local crystal structure rotates through an angle associated with a point group operation as the line is encircled. Such defects are known as disclinations. Although the *π*<sup>3</sup> (*R*) has the classification problem (Mermin, 1979), the domain walls (grain boundaries) are thought as the string bounded topological defects by the broken discrete symmetry (Vilenkin, 1985), hence so that textures evolve.

Once the order parameter class is identified, the partition functions are divided by the corresponding topologically invariant number. Indeed, Burakovski et al. (2000) rederived the free energy density as a function of dislocation density depending on temperature, starting from the partition function computation for the indistinguishable dislocation loops. The resultant free energy density function is very similar to that out of the constitutional approach (Cotterill, 1980). The free-energy is in the form of the two-dimensional *xy* model, which should undergoes a two-dimensional Kosterlitz-Thouless phase transtion (Kosterlitz & Thouless, 1973). The theory break down point, *i.e.*, the melting temperature is also derived. They also showed that the dislocation density *ρ* cannot be a continuous function

0, *T* < *Tm*,

where *Tm* is the melting temperature. The dislocation-mediated melting mechanism was observed in a bulk colloidal crystals (Alsayed et al., 2005). A good review on two-dimensional

*<sup>ρ</sup>* (*Tm*) <sup>∼</sup> *<sup>b</sup>*−2, *<sup>T</sup>* <sup>=</sup> *Tm*, (207)

deformed to the constant map.

topologically stable defects.

**7.2 Renormalization group and gauge theory**

of temperature, but it is a gap function

*ρ* (*T*) =

melting phenomena is given by Strandburg (1986; 1989).

magnons. All those quasiparticles should interact each other, in contrast to the typical considerations that quasiparticles hardly interacts each other. Lee et al. (2006) confirmed that there are interactions between phonon and magnetism, and Bodraykov (2007) in his phenomenological model reproduced the invar effect of the temperature dependent thermal expansion coeffient by considering the results of the interacting phonon-magnon (Kim, 1982).
