**Towards the Authentic** *Ab Intio* **Thermodynamics**

In Gee Kim

*Graduate Institute of Ferrous Technology, Pohang University of Science and Technology, Pohang Republic of Korea*

#### **1. Introduction**

542 Thermodynamics – Interaction Studies – Solids, Liquids and Gases

[31] He, G. Z.; Pan, G.; Zhang, M. Y.; Wu, Z. Y., *J. Phys. Chem. C* 2009, 113 (39), 17076-17081.

[32] Zhang, M. Y.; He, G. Z.; Pan, G., *J. Colloid Interface Sci.* 2009, 338 (1), 284-286.

A phase diagram is considered as a starting point to design new materials. Let us quote the statements by DeHoff (1993):

A phase diagram is a map that presents the domains of stability of phases and their combiations. A point in this space, which represents a state of the system that is of interest in a particular application, lies within a specific domain on the map.

In practice, for example to calculate the lattice stability, the construction of the phase diagram is to find the phase equilibria based on the comparison of the Gibbs free energies among the possible phases. Hence, the most important factor is the accuracy and precesion of the given Gibbs free energy values, which are usually acquired by the experimental assessments. Once the required thermodynamic data are obtained, the phase diagram construction becomes rather straightforward with modern computation techniques, so called CALPHAD (CALculation of PHAse Diagrams) (Spencer, 2007). Hence, the required information for constructing a phase diagram is the reliable Gibbs free energy information. The Gibbs free energy *G* is defined by

$$G = E + PV - TS\_{\prime} \tag{1}$$

where *E* is the internal energy, *P* is the pressure, *V* is the volume of the system, *T* is the temperature and *S* is the entropy. The state which provides the minimum of the free energy under given external conditions at constant *P* and *T* is the equilibrium state. However, there is a critical issue to apply the conventional CALPHAD method in general materials design. Most thermodynamic information is relied on the experimental assessments, which do not available occasionally to be obtained, but necessary. For example, the direct thermodynamic information of silicon solubility in cementite had not been available for long time (Ghosh & Olson, 2002; Kozeschnik & Bhadeshia, 2008), because the extremely low silicon solubility which requires the information at very high temperature over the melting point of cementite. The direct thermodynamic information was available recently by an *ab initio* method (Jang et al., 2009). However, the current technology of *ab initio* approaches is usually limited to zero temperature, due to the theoretical foundation; the density functional theory (Hohenberg & Kohn, 1964) guarrentees the unique total energy of the ground states only. The example demonstrates the necessity of a systematic assessment method from first principles. In order to obtain the Gibbs free energy from first principles, it is convenient to use the equilibrium statistical mechanics for grand canonical ensemble by introducing the *grand*

In principle, we can calculate any macroscopic thermodynamic states if we have the complete knowledge of the (grand) partition function, which is abled to be constructed from first principles. However, it is impractical to calculates the partition function of a given system

Towards the Authentic *Ab Intio* Thermodynamics 545

Struggles have been devoted to calculate the summation of all accessible states. The number of all accessible states is evaluated by the constitutents of the system and the types of interaction among the constituents. The general procedure in statistical mechanics is nothing more than the calculation of the probability of a specific number of dice with the enormous number of repititions of the dice tosses. The fundamental principles of statistical mechanics of a mechanical system of the degrees of freedom *s* is well summarized by Landau & Lifshitz (1980). The state of a mechanical system is described a point of the *phase space* represented by the generalized coordinates *qi* and the corresponding generalized momenta *pi*, where the index *i* runs from 1 to *s*. The time evolution of the system is represented by the trajectory in the phase space. Let us consider a closed large mechanical system and a part of the entire system, called *subsystem*, which is also large enough, and is interacting with the rest part of the closed system. An exact solution for the behavior of the subsystem can be obtained only

Let us assume that the subsystem is in the small phase volume Δ*p*Δ*q* for short intervals. The

where *D* is the long time interval in which the short interval Δ*t* is included. Defining the

d*p*d*q* = d*p*1d*p*<sup>2</sup> ...d*ps*d*q*1d*q*<sup>2</sup> ...d*qs*,

where *ρ* is a function of all coordinates and momenta in writing for brevity *ρ* (*p*, *q*). This function *ρ* represents the density of the probability distribution in phase space, called

One should note that the statistical distribution of a given subsystem does not depend on the initial state of any other subsystems of the entire system, due to the entirely outweighed

A physical quantity *f* = *f*(*p*, *q*) depending on the states of the subsystem of the solved entire system is able to be evaluated, in the sense of the statistical average, by the distribution

By definition Eq. (12) of the probability, the statistical averaging is exactly equivalent to a time

1 *D*  *D* 0

(*statistical*) *distribution function*. Obviously, the distribution function is normalized as

¯ *f* = 

> ¯ *f* = lim *D*→∞

Δ*t*

d*w* = *ρ* (*p*1, *p*2,..., *ps*, *q*1, *q*2,..., *qs*) d*p*d*q*, (13)

*ρ* (*p*, *q*) d*p*d*q* = 1. (14)

*f*(*p*, *q*)*ρ*(*p*, *q*)d*p*d*q*. (15)

*f* (*t*) d*t*. (16)

*<sup>D</sup>* , (12)

*w* = lim *D*→∞

probability *w* for the subsystem stays in the Δ*p*Δ*q* during the short interval Δ*t* is

because the number of all accessible microstates, indexed by *ζ*, is enormously large.

by solving the mechanical problem for the entire closed system.

probability d*w* of states represented in the phase volum,

effects of the initial state over a sufficiently long time.

averaging, which is established as

may be written

function as

*partition function*

$$\Xi\left(T, V, \left\{\mu\_{i}\right\}\right) = \sum\_{N\_{i}} \sum\_{\tilde{\zeta}} \exp\left(-\beta \left(E\_{\tilde{\zeta}}\left(V\right) - \sum\_{i} \mu\_{i} N\_{i}\right)\right),\tag{2}$$

where *β* is the inverse temperature (*k*B*T*) <sup>−</sup><sup>1</sup> with the Boltzmann's constant *<sup>k</sup>*B, *<sup>μ</sup><sup>i</sup>* is the chemical potential of the *i*th component, *Ni* is the number of atoms. The sum of *ζ* runs over all accessible microstates of the system; the microstates include the electronic, magnetic, vibrational and configurational degrees of freedom. The corresponding *grand potential* Ω is found by

$$
\Omega\left(T, V, \{\mu\_i\}\right) = -\beta^{-1} \ln \Xi. \tag{3}
$$

The Legendre transformation relates the grand potential Ω and the *Helmholtz free energy F* as

$$
\Omega\left(T, V, \{\mu\_i\}\right) = F - \sum\_i \mu\_i N\_i = E - TS - \sum\_i \mu\_i N\_i. \tag{4}
$$

It is noticeable to find that the Helmholtz free energy *F* is able to be obtained by the relation

$$F\left(T, V, N\right) = -\beta^{-1} \ln Z\_{\prime} \tag{5}$$

where *Z* is the *partition function* of the canonical ensemble defined as

$$Z\left(T, V, N\right) = \sum\_{\zeta} \exp\left(-\beta E\_{\zeta}\left(V, N\right)\right). \tag{6}$$

Finally, there is a further Legendre transformation relationship between the Helmholtz free energy and the Gibbs free energy as

$$
\mathcal{G} = \mathcal{F} + PV.\tag{7}
$$

Let us go back to the grand potential in Eq. (4). The total differential of the grand potential is

$$\mathbf{d}\Omega = -\mathbf{Sd}T - P\mathbf{d}V - \sum\_{i} \mathbf{N}\_{i}\mathbf{d}\mu\_{i\nu} \tag{8}$$

with the coefficients

$$S = -\left(\frac{\partial \Omega}{\partial T}\right)\_{V\mu}{}^{\prime} \quad P = -\left(\frac{\partial \Omega}{\partial V}\right)\_{T\mu}{}^{\prime} \quad N\_{\text{i}} = -\left(\frac{\partial \Omega}{\partial \mu\_{\text{i}}}\right)\_{TV} \tag{9}$$

The Gibss-Duhem relation,

$$E = TS - PV + \sum\_{i} \mu\_{i} N\_{i\prime} \tag{10}$$

yields the thermodynamic functions as

$$F = -PV + \sum\_{i} \mu\_{i} N\_{i\nu} \quad G = \sum\_{i} \mu\_{i} N\_{i\nu} \quad \Omega = -PV. \tag{11}$$

Since the thermodynamic properties of a system at equilibrium are specified by Ω and derivatives thereof, one of the tasks will be to develop methods to calculate the grand potential Ω.

2 Will-be-set-by-IN-TECH

chemical potential of the *i*th component, *Ni* is the number of atoms. The sum of *ζ* runs over all accessible microstates of the system; the microstates include the electronic, magnetic, vibrational and configurational degrees of freedom. The corresponding *grand potential* Ω is

The Legendre transformation relates the grand potential Ω and the *Helmholtz free energy F* as

*i*

*ζ*

exp

Finally, there is a further Legendre transformation relationship between the Helmholtz free

Let us go back to the grand potential in Eq. (4). The total differential of the grand potential is

 *∂*Ω *∂V* 

*<sup>E</sup>* = *TS* − *PV* + ∑

*μiNi*, *G* = ∑

Since the thermodynamic properties of a system at equilibrium are specified by Ω and derivatives thereof, one of the tasks will be to develop methods to calculate the grand potential

<sup>d</sup><sup>Ω</sup> = −*S*d*<sup>T</sup>* − *<sup>P</sup>*d*<sup>V</sup>* − ∑

, *P* = −

It is noticeable to find that the Helmholtz free energy *F* is able to be obtained by the relation

*<sup>E</sup><sup>ζ</sup>* (*V*) − ∑

*i μiNi*

<sup>Ω</sup> (*T*, *<sup>V</sup>*, {*μi*}) <sup>=</sup> <sup>−</sup>*β*−<sup>1</sup> ln <sup>Ξ</sup>. (3)

*i*

*<sup>F</sup>* (*T*, *<sup>V</sup>*, *<sup>N</sup>*) <sup>=</sup> <sup>−</sup>*β*−<sup>1</sup> ln *<sup>Z</sup>*, (5)

*G* = *F* + *PV*. (7)

 *∂*Ω *∂μ<sup>i</sup>* 

*<sup>μ</sup>iNi* = *<sup>E</sup>* − *TS* − ∑

−*βE<sup>ζ</sup>* (*V*, *N*)

*i*

, *Ni* = −

*Tμ*

*i*

*i*

<sup>−</sup><sup>1</sup> with the Boltzmann's constant *<sup>k</sup>*B, *<sup>μ</sup><sup>i</sup>* is the

, (2)

*μiNi*. (4)

. (6)

*Ni*d*μi*, (8)

*TV*

*μiNi*, (10)

*μiNi*, Ω = −*PV*. (11)

. (9)

<sup>Ξ</sup> (*T*, *<sup>V</sup>*, {*μi*}) = ∑

where *β* is the inverse temperature (*k*B*T*)

energy and the Gibbs free energy as

*S* = −

yields the thermodynamic functions as

*∂*Ω *∂T* 

*Vμ*

*<sup>F</sup>* = −*PV* + ∑

*i*

with the coefficients

Ω.

The Gibss-Duhem relation,

*Ni* ∑ *ζ* exp −*β* 

<sup>Ω</sup> (*T*, *<sup>V</sup>*, {*μi*}) = *<sup>F</sup>* − ∑

where *Z* is the *partition function* of the canonical ensemble defined as

*Z* (*T*, *V*, *N*) = ∑

*partition function*

found by

In principle, we can calculate any macroscopic thermodynamic states if we have the complete knowledge of the (grand) partition function, which is abled to be constructed from first principles. However, it is impractical to calculates the partition function of a given system because the number of all accessible microstates, indexed by *ζ*, is enormously large.

Struggles have been devoted to calculate the summation of all accessible states. The number of all accessible states is evaluated by the constitutents of the system and the types of interaction among the constituents. The general procedure in statistical mechanics is nothing more than the calculation of the probability of a specific number of dice with the enormous number of repititions of the dice tosses. The fundamental principles of statistical mechanics of a mechanical system of the degrees of freedom *s* is well summarized by Landau & Lifshitz (1980). The state of a mechanical system is described a point of the *phase space* represented by the generalized coordinates *qi* and the corresponding generalized momenta *pi*, where the index *i* runs from 1 to *s*. The time evolution of the system is represented by the trajectory in the phase space. Let us consider a closed large mechanical system and a part of the entire system, called *subsystem*, which is also large enough, and is interacting with the rest part of the closed system. An exact solution for the behavior of the subsystem can be obtained only by solving the mechanical problem for the entire closed system.

Let us assume that the subsystem is in the small phase volume Δ*p*Δ*q* for short intervals. The probability *w* for the subsystem stays in the Δ*p*Δ*q* during the short interval Δ*t* is

$$w = \lim\_{D \to \infty} \frac{\Delta t}{D} \tag{12}$$

where *D* is the long time interval in which the short interval Δ*t* is included. Defining the probability d*w* of states represented in the phase volum,

$$\mathbf{d}p\mathbf{d}q = \mathbf{d}p\_1\mathbf{d}p\_2\dots \mathbf{d}p\_s\mathbf{d}q\_1\mathbf{d}q\_2\dots \mathbf{d}q\_{s'}$$

may be written

$$\mathbf{dw} = \rho \begin{pmatrix} p\_1, p\_2, \dots, p\_s, q\_1, q\_2, \dots, q\_s \end{pmatrix} \mathbf{dp} \mathbf{d} \boldsymbol{\eta}, \tag{13}$$

where *ρ* is a function of all coordinates and momenta in writing for brevity *ρ* (*p*, *q*). This function *ρ* represents the density of the probability distribution in phase space, called (*statistical*) *distribution function*. Obviously, the distribution function is normalized as

$$
\int \rho \left( p, q \right) \mathrm{d}p \mathrm{d}q = 1. \tag{14}
$$

One should note that the statistical distribution of a given subsystem does not depend on the initial state of any other subsystems of the entire system, due to the entirely outweighed effects of the initial state over a sufficiently long time.

A physical quantity *f* = *f*(*p*, *q*) depending on the states of the subsystem of the solved entire system is able to be evaluated, in the sense of the statistical average, by the distribution function as

$$f = \int f(p, q)\rho(p, q)\mathbf{d}p\mathbf{d}q.\tag{15}$$

By definition Eq. (12) of the probability, the statistical averaging is exactly equivalent to a time averaging, which is established as

$$\overline{f} = \lim\_{D \to \infty} \frac{1}{D} \int\_0^D f\left(t\right) dt. \tag{16}$$

for the minima points *A* and *B*. For the case of the *second-order phase transition*, it is required

Towards the Authentic *Ab Intio* Thermodynamics 547

The second derivative must vanish because the curve changes from concave to convex and the third derivative must vanish to ensure that the critical point is a minimum. It is convenient to reduce the variables in the vicinity of the critical point *t* ≡ *T* − *T*<sup>C</sup> and *h* ≡ *H* − *Hc* = *H*, where *T*<sup>C</sup> is the Curie temperature and *Hc* is the critical external field, yielding the Landau

Enforcing the inversion symmetry, L (*m*, *H*, *T*) = L (−*m*, −*H*, *T*), the Landau function will be

<sup>L</sup> (*m*, *<sup>h</sup>*, *<sup>t</sup>*) <sup>=</sup> *<sup>d</sup>*2*tm*<sup>2</sup> <sup>+</sup> *<sup>b</sup>*4*m*4. In order to see the dependency to the external field *H*, we add an arbitrary *H* field coupling

> 1 2

Let us consider the second-order phase transition with *H* = 0. For *T* > *T*C, the minimum of L is at *m* = 0. For *T* = *T*C, the Landau function has zero curvature at *m* = 0, where the point is still the global minimum. For *T* < *T*C, the Landau function Eq. (26) has two degenerate

<sup>−</sup>*at*

When *H* �= 0, the differentiation of L with respect to *m* gives the magnetic equation of state

The isothermal magnetic susceptibility is obtained by differentiating Eq. (28) with respect to

where *m* (*H*) is the solution of Eq. (28). Let us consider the case of *H* = 0. For *t* > 0, *m* = 0 and *<sup>χ</sup><sup>T</sup>* <sup>=</sup> 1/ (2*at*), while *<sup>m</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup>*at*/*<sup>b</sup>* and *<sup>χ</sup><sup>T</sup>* <sup>=</sup> <sup>−</sup>1/ (4*at*). As the system is cooled down, the nonmagnetized system, *m* = 0 for *t* > 0, occurs a spontaneous magnetization of (−*at*/*b*)

below the critical temperature *t* < 0, while the isothermal magnetic susceptibility *χ<sup>T</sup>* diverges

For the first-order phase transition, we need to consider Eq. (25) with *c*<sup>1</sup> = 0 and changing the

1 2 <sup>=</sup> <sup>1</sup> 2 

*at* + 3*b* (*m* (*H*))<sup>2</sup>

*atm* <sup>+</sup> *bm*<sup>3</sup> <sup>=</sup> <sup>1</sup>

 *T*

*<sup>∂</sup>m*<sup>3</sup> <sup>=</sup> 0, *<sup>∂</sup>*4<sup>L</sup>

*an* (*H*, *T*) �→ *an* (*h*, *t*) = *bn* + *cnh* + *dnt*, (24)

2 *b*:

*bm*<sup>4</sup> <sup>−</sup> *Hm*. (26)

*<sup>b</sup>* , for *<sup>t</sup>* <sup>&</sup>lt; 0. (27)

*<sup>m</sup>*<sup>4</sup> <sup>+</sup> *Cm*<sup>3</sup> <sup>−</sup> *Hm*. (30)

<sup>2</sup> *<sup>H</sup>*. (28)

, (29)

1 2

<sup>L</sup> (*m*, *<sup>h</sup>*, *<sup>t</sup>*) <sup>=</sup> *<sup>c</sup>*1*hm* <sup>+</sup> *<sup>d</sup>*2*tm*<sup>2</sup> <sup>+</sup> *<sup>c</sup>*3*hm*<sup>3</sup> <sup>+</sup> *<sup>b</sup>*4*m*4, *<sup>d</sup>*<sup>2</sup> <sup>&</sup>gt; 0, *<sup>b</sup>*<sup>4</sup> <sup>&</sup>gt; 0. (25)

*<sup>∂</sup>m*<sup>4</sup> <sup>&</sup>gt; 0. (23)

*<sup>∂</sup>m*<sup>2</sup> <sup>=</sup> *<sup>∂</sup>*3<sup>L</sup>

*∂*L *<sup>∂</sup><sup>m</sup>* <sup>=</sup> *<sup>∂</sup>*2<sup>L</sup>

and then the Ladau function near the critical point is

minima at *ms* = *ms* (*T*), which is explicitly

term and change the symbols of the coefficients *d*<sup>2</sup> to *a* and *b*<sup>4</sup> to <sup>1</sup>

<sup>L</sup> <sup>=</sup> *atm*<sup>2</sup> <sup>+</sup>

*ms* (*t*) = ±

*<sup>χ</sup><sup>T</sup>* (*H*) <sup>≡</sup> *<sup>∂</sup><sup>m</sup>* (*H*)

as 1/*t* for *t* → 0 both for the regions of *t* > 0 and *t* < 0.

*∂H*

<sup>L</sup> <sup>=</sup> *atm*<sup>2</sup> <sup>+</sup>

that

coefficient

for small *m* as

coefficient symbols to yield

*H*:

In addition, the Liouville's theorem

$$\frac{\mathbf{d}\rho}{\mathbf{d}t} = \sum\_{i=1}^{s} \left( \frac{\partial \rho}{\partial q\_i} \dot{q}\_i + \frac{\partial \rho}{\partial p\_i} \dot{p}\_i \right) = 0 \tag{17}$$

tells us that the distribution function is constant along the phase trajectories of the subsystem. Our interesting systems are (quantum) mechanical objects, so that the counting the number of accessible states is equivalent to the estimation of the relevant phase space volume.

#### **2. Phenomenological Landau theory**

A ferromagnet in which the magnetization is the order parameter is served for illustrative purpose. Landau & Lifshitz (1980) suggested a phenomenological description of phase transitions by introducing a concept of *order parameter*. Suppose that the interaction Hamiltonian of the magnetic system to be

$$\sum\_{i,j} f\_{ij} \mathbf{S}\_i \cdot \mathbf{S}\_j \tag{18}$$

where **S***<sup>i</sup>* is a localized Heisenberg-type spin at an atomic site *i* and *Jij* is the interaction parameter between the spins **S***<sup>i</sup>* and **S***j*.

In the ferromagnet, the total magnetization **M** is defined as the thermodynamic average of the spins

$$\mathbf{M} = \left\langle \sum\_{i} \mathbf{S}\_{i} \right\rangle\_{\prime} \tag{19}$$

and the magnetization **m** denotes the magnetization per spin

$$\mathbf{m} = \left\langle \frac{1}{N} \sum\_{i} \mathbf{S}\_{i} \right\rangle\_{\prime} \tag{20}$$

where *N* is the number of atomic sites. The physical order is the alignment of the microscopic spins.

Let us consider a situation that an external magnetic field **H** is applied to the system. Landau's idea1 is to introduce a function, <sup>L</sup> (**m**, **<sup>H</sup>**, *<sup>T</sup>*), known as the Landau function, which describes the "thermodynamics" of the system as function of **m**, **H**, and *T*. The minimum of L indicates the system phase at the given variable values. To see more details, let us expand the Ladau function with respect to the order parameter **m**:

$$\mathcal{L}\left(m, H, T\right) = \sum\_{n}^{4} a\_{n}\left(H, T\right)m^{n}\,. \tag{21}$$

where we assumed that both the magnetization **m** and the external magnetic field **H** are aligned in a specific direction, say ˆ**z**. When the system undergoes a *first-order phase transition*, the Landau function should have the properties

$$
\left.\frac{\partial \mathcal{L}}{\partial m}\right|\_{m\_A} = \left.\frac{\partial \mathcal{L}}{\partial m}\right|\_{m\_B} = 0, \qquad \mathcal{L}\left(m\_A\right) = \mathcal{L}\left(m\_B\right), \tag{22}
$$

<sup>1</sup> The description in this section is following Negele & Orland (1988) and Goldenfeld (1992).

4 Will-be-set-by-IN-TECH

tells us that the distribution function is constant along the phase trajectories of the subsystem. Our interesting systems are (quantum) mechanical objects, so that the counting the number of

A ferromagnet in which the magnetization is the order parameter is served for illustrative purpose. Landau & Lifshitz (1980) suggested a phenomenological description of phase transitions by introducing a concept of *order parameter*. Suppose that the interaction

where **S***<sup>i</sup>* is a localized Heisenberg-type spin at an atomic site *i* and *Jij* is the interaction

In the ferromagnet, the total magnetization **M** is defined as the thermodynamic average of the

where *N* is the number of atomic sites. The physical order is the alignment of the microscopic

Let us consider a situation that an external magnetic field **H** is applied to the system. Landau's idea1 is to introduce a function, <sup>L</sup> (**m**, **<sup>H</sup>**, *<sup>T</sup>*), known as the Landau function, which describes the "thermodynamics" of the system as function of **m**, **H**, and *T*. The minimum of L indicates the system phase at the given variable values. To see more details, let us expand the Ladau

> 4 ∑*n*

where we assumed that both the magnetization **m** and the external magnetic field **H** are aligned in a specific direction, say ˆ**z**. When the system undergoes a *first-order phase transition*,

 ∑ *i* **S***i* 

*<sup>q</sup>*˙*<sup>i</sup>* <sup>+</sup> *∂ρ ∂pi p*˙*i* 

= 0 (17)

*Jij***S***<sup>i</sup>* · **S***j*, (18)

, (19)

, (20)

*an* (*H*, *T*) *mn*, (21)

= 0, L (*mA*) = L (*mB*), (22)

In addition, the Liouville's theorem

**2. Phenomenological Landau theory**

Hamiltonian of the magnetic system to be

parameter between the spins **S***<sup>i</sup>* and **S***j*.

spins

spins.

d*ρ* <sup>d</sup>*<sup>t</sup>* <sup>=</sup>

*s* ∑ *i*=1  *∂ρ ∂qi*

accessible states is equivalent to the estimation of the relevant phase space volume.

∑ *i*,*j*

**M** =

**m** = 1 *<sup>N</sup>* ∑ *i* **S***i* 

L (*m*, *H*, *T*) =

<sup>1</sup> The description in this section is following Negele & Orland (1988) and Goldenfeld (1992).

<sup>=</sup> *<sup>∂</sup>*<sup>L</sup> *∂m mB*

and the magnetization **m** denotes the magnetization per spin

function with respect to the order parameter **m**:

the Landau function should have the properties

*∂*L *∂m mA* for the minima points *A* and *B*. For the case of the *second-order phase transition*, it is required that

$$
\frac{\partial \mathcal{L}}{\partial m} = \frac{\partial^2 \mathcal{L}}{\partial m^2} = \frac{\partial^3 \mathcal{L}}{\partial m^3} = 0, \qquad \frac{\partial^4 \mathcal{L}}{\partial m^4} > 0. \tag{23}
$$

The second derivative must vanish because the curve changes from concave to convex and the third derivative must vanish to ensure that the critical point is a minimum. It is convenient to reduce the variables in the vicinity of the critical point *t* ≡ *T* − *T*<sup>C</sup> and *h* ≡ *H* − *Hc* = *H*, where *T*<sup>C</sup> is the Curie temperature and *Hc* is the critical external field, yielding the Landau coefficient

$$a\_{\rm nl} \left( H, T \right) \mapsto a\_{\rm nl} \left( h, t \right) = b\_{\rm nl} + c\_{\rm nl} h + d\_{\rm nl} t \,, \tag{24}$$

and then the Ladau function near the critical point is

$$\mathcal{L}\left(m,h,t\right) = c\_1 lm + d\_2 t m^2 + c\_3 lm^3 + b\_4 m^4, \quad d\_2 > 0, \quad b\_4 > 0. \tag{25}$$

Enforcing the inversion symmetry, L (*m*, *H*, *T*) = L (−*m*, −*H*, *T*), the Landau function will be

$$\mathcal{L}\left(m,h,t\right) = d\_2 t m^2 + b\_4 m^4.$$

In order to see the dependency to the external field *H*, we add an arbitrary *H* field coupling term and change the symbols of the coefficients *d*<sup>2</sup> to *a* and *b*<sup>4</sup> to <sup>1</sup> 2 *b*:

$$
\mathcal{L} = atm^2 + \frac{1}{2}bm^4 - Hm.\tag{26}
$$

Let us consider the second-order phase transition with *H* = 0. For *T* > *T*C, the minimum of L is at *m* = 0. For *T* = *T*C, the Landau function has zero curvature at *m* = 0, where the point is still the global minimum. For *T* < *T*C, the Landau function Eq. (26) has two degenerate minima at *ms* = *ms* (*T*), which is explicitly

$$m\_s\left(t\right) = \pm \sqrt{\frac{-at}{b}}, \quad \text{for } t < 0. \tag{27}$$

When *H* �= 0, the differentiation of L with respect to *m* gives the magnetic equation of state for small *m* as

$$
tau + bm^3 = \frac{1}{2}H.\tag{28}$$

The isothermal magnetic susceptibility is obtained by differentiating Eq. (28) with respect to *H*:

$$\chi\_T(H) \equiv \left. \frac{\partial m}{\partial H} \right|\_T = \frac{1}{2 \left\{ at + 3b \left( m \left( H \right) \right)^2 \right\}} \,\tag{29}$$

where *m* (*H*) is the solution of Eq. (28). Let us consider the case of *H* = 0. For *t* > 0, *m* = 0 and *<sup>χ</sup><sup>T</sup>* <sup>=</sup> 1/ (2*at*), while *<sup>m</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup>*at*/*<sup>b</sup>* and *<sup>χ</sup><sup>T</sup>* <sup>=</sup> <sup>−</sup>1/ (4*at*). As the system is cooled down, the nonmagnetized system, *m* = 0 for *t* > 0, occurs a spontaneous magnetization of (−*at*/*b*) 1 2 below the critical temperature *t* < 0, while the isothermal magnetic susceptibility *χ<sup>T</sup>* diverges as 1/*t* for *t* → 0 both for the regions of *t* > 0 and *t* < 0.

For the first-order phase transition, we need to consider Eq. (25) with *c*<sup>1</sup> = 0 and changing the coefficient symbols to yield

$$\mathcal{L} = atm^2 + \frac{1}{2}m^4 + \mathbb{C}m^3 - Hm.\tag{30}$$

estimated as

theory is self-consistent.

*ξ* for *t* → 0. A simple arithematics yields *m* ∼ |*t*|

ferromagnet. This result leaves us the error

significantly the phase transitions of the system.

nucleus. It is customary to call valence electrons as electrons.

reduces the burdens of calculation of the distribution function of electrons.

**3. Matters as noninteracting gases**

**3.1 Electronic subsystem as Fermi gas**

(*a*/*R*)

*ε<sup>R</sup>* �

 *T T*C

where *a* is the lattice constant and *d* is the dimensionality of the interaction. In Eq. (36),

Towards the Authentic *Ab Intio* Thermodynamics 549

On the other hand, the correlation length grows toward infinity near the critical point; *R* �

which tends to infinity as *t* → 0. Hence, the Landau theory based on the mean-field approximation has error which diverges as the system approaches to the critical point. Mathematically, the Landau theory expands the Landau function in terms of the order parameter. The landau expansion itself is mathematically non-sense *near* the critical point for dimensions less than four. Therefore, the Landau theory is not a good tool to investigate

Materials are basically made of atoms; an atom is composed of a nucleus and the surrounding electrons. However, it is convenient to distinguish two types of electrons; the *valence electrons* are responsible for chemical reactions and the *core electrons* are tightly bound around the nucleus to form an *ion* for screening the strongly divergent Coulomb potential from the

The decomposition into electrons and ions provides us at least two advantages in treating materials with first-principles. First of all, the motions of electrons can be decoupled adiabatically from the those of ions, since electrons reach their equilibrium almost immediately by their light mass compared to those factors of ions. The decoupling of the motions of electrons from those of ions is accomplished by the Born-Oppenheimer adiabatic approximation (Born & Oppenheimer, 1927), which decouples the motions of electrons approximately begin independent adiabatically from those of ions. In practice, the motions of electrons are computed under the external potential influenced by the ions at their *static* equilibrium positions, before the motions of ions are computed under the external potential influenced by the electronic distribution. Hence, the fundamental information for thermodynamics of a material is its electronic structures. Secondly, the decoupled electrons of spin half are identical particles following the Fermi-Dirac statistics (Dirac, 1926; Fermi, 1926). Hence, the statistical distribution function of electrons is a closed fixed form. This feature

The consequence of the decoupling electrons from ions allows us to treat the distribution functions of distinguishable atoms, for example, an iron atom is distinguished from a carbon atom, can be treated as the source of external potential to the electronic subsystem. Modelling of electronic subsystem was suggested firstly by Drude (1900), before the birth of quantum mechanics. He assumes that a metal is composed of electrons wandering on the positive homogeneous ionic background. The interaction between electrons are cancelled to allow us

*<sup>ε</sup><sup>R</sup>* <sup>∼</sup> <sup>1</sup> |*t*| 2*β a R d*

 *a R d*

<sup>−</sup>*<sup>d</sup>* is essentially the corrdination number *<sup>z</sup>* > 1, so that *<sup>ε</sup><sup>R</sup>* < 1 and the mean field

, (36)

, (37)

<sup>2</sup> for a

*<sup>β</sup>*, where a *critial exponent β* is <sup>1</sup>

For *H* = 0, the equilibrium value of *m* is obtained as

$$
\hbar m = 0, \quad m = -c \pm \sqrt{c^2 - at/b}, \tag{31}
$$

where *c* = 3*C*/4*b*. The nonzero solution is valid only for *t* < *t*∗, by defining *t* <sup>∗</sup> <sup>≡</sup> *bc*2/*a*. Let *Tc* is the temperature where the coefficient of the term quadratic in *m* vanishes. Suppose *t*<sup>1</sup> is the temperature where the value of L at the secondary minimum is equal to the value at *m* = 0. Since *t*∗ is positive, this occurs at a temperature greater than *Tc*. For *t* < *t* ∗, a secondary minimum and maximum have developed, in addition to the minimum at *m* = 0. For *t* < *t*1, the secondary minimum is now the global minimum, and the value of the order parameter which minimizes L jumps *discontinuously* from *m* = 0 to a non-zero value. This is a first-order transition. Note that at the first-order transition, *m* (*t*1) is not arbitrarily small as *t* → *t* − <sup>1</sup> . In other words, the Landau theory is *not* valid. Hence, the first-order phase transition is arosen by introducing the cubic term in *m*.

Since the Landau theory is fully phenomenological, there is no strong limit in selecting order parameter and the corresponding conjugate field. For example, the magnetization is the order parameter of a ferromagnet with the external magnetic field as the conjugate coupling field, the polarization is the order parameter of a ferroelectric with the external electric field as the conjugate coupling field, and the electron pair amplitude is the order parameter of a superconductor with the electron pair source as the conjugate coupling field. When a system undergoes a phase transition, the Landau theory is usually utilized to understand the phase transition.

The Landau theory is motivated by the observation that we could replace the interaction Hamiltonian Eq. (18)

$$\sum\_{i,j} I\_{\bar{i}\bar{j}} S\_{\bar{i}} S\_{\bar{j}} = \sum\_{\bar{i}} S\_{\bar{i}} \sum\_{\bar{j}} I\_{\bar{i}\bar{j}} \left[ \langle S\_{\bar{j}} \rangle + \left\{ S\_{\bar{j}} - \langle S\_{\bar{j}} \rangle \right\} \right] \tag{32}$$

by ∑*ij SiJij*�*Sj*�. If we can replace *SiSj* by *Si*�*Sj*�, it is also possible to replace �*SiSj*� by �*Si*��*Sj*� on average if we assume the translational invariance. The fractional error implicit in this replacement can be evaluated by

$$\varepsilon\_{ij} = \frac{\left| \langle \mathbb{S}\_i \mathbb{S}\_j \rangle - \langle \mathbb{S}\_i \rangle \langle \mathbb{S}\_j \rangle \right|}{\langle \mathbb{S}\_i \rangle \langle \mathbb{S}\_j \rangle},\tag{33}$$

where all quantities are measured for *T* < *T*<sup>C</sup> under the Landau theory. The numerator is just a correlation function *C* and the interaction range **r***<sup>i</sup>* − **r***<sup>j</sup>* <sup>∼</sup> *<sup>R</sup>* will allow us to rewrite *<sup>ε</sup>ij* as

$$\varepsilon\_R = \frac{|\mathbb{C}\left(R\right)|}{m\_s^2},\tag{34}$$

where we assume the correlation function being written as

$$\mathbb{C}\left(R\right) = \mathbb{g}f\left(\frac{R}{\tilde{\xi}}\right),\tag{35}$$

where *f* is a function of the correlation length *ξ*. For *T* � *T*C, the correlation length *ξ* ∼ *R*, and the order parameter *m* is saturated at the low temperature value. The error is roughly estimated as

6 Will-be-set-by-IN-TECH

*Tc* is the temperature where the coefficient of the term quadratic in *m* vanishes. Suppose *t*<sup>1</sup> is the temperature where the value of L at the secondary minimum is equal to the value at

minimum and maximum have developed, in addition to the minimum at *m* = 0. For *t* < *t*1, the secondary minimum is now the global minimum, and the value of the order parameter which minimizes L jumps *discontinuously* from *m* = 0 to a non-zero value. This is a first-order transition. Note that at the first-order transition, *m* (*t*1) is not arbitrarily small as *t* → *t*

other words, the Landau theory is *not* valid. Hence, the first-order phase transition is arosen

Since the Landau theory is fully phenomenological, there is no strong limit in selecting order parameter and the corresponding conjugate field. For example, the magnetization is the order parameter of a ferromagnet with the external magnetic field as the conjugate coupling field, the polarization is the order parameter of a ferroelectric with the external electric field as the conjugate coupling field, and the electron pair amplitude is the order parameter of a superconductor with the electron pair source as the conjugate coupling field. When a system undergoes a phase transition, the Landau theory is usually utilized to understand the phase

The Landau theory is motivated by the observation that we could replace the interaction

by ∑*ij SiJij*�*Sj*�. If we can replace *SiSj* by *Si*�*Sj*�, it is also possible to replace �*SiSj*� by �*Si*��*Sj*� on average if we assume the translational invariance. The fractional error implicit in this

�*SiSj*�−�*Si*��*Sj*�

 **r***<sup>i</sup>* − **r***<sup>j</sup>* 

 *R ξ* 

where all quantities are measured for *T* < *T*<sup>C</sup> under the Landau theory. The numerator is just

*<sup>ε</sup><sup>R</sup>* <sup>=</sup> <sup>|</sup>*<sup>C</sup>* (*R*)<sup>|</sup> *m*2 *s*

where *f* is a function of the correlation length *ξ*. For *T* � *T*C, the correlation length *ξ* ∼ *R*, and the order parameter *m* is saturated at the low temperature value. The error is roughly

*C* (*R*) = *g f*

*Sj* − �*Sj*�

  �*Si*��*Sj*� , (33)

<sup>∼</sup> *<sup>R</sup>* will allow us to rewrite *<sup>ε</sup>ij* as

, (34)

, (35)

*<sup>c</sup>*<sup>2</sup> − *at*/*b*, (31)

<sup>∗</sup> <sup>≡</sup> *bc*2/*a*. Let

∗, a secondary

− <sup>1</sup> . In

(32)

*<sup>m</sup>* <sup>=</sup> 0, *<sup>m</sup>* <sup>=</sup> <sup>−</sup>*<sup>c</sup>* <sup>±</sup>

where *c* = 3*C*/4*b*. The nonzero solution is valid only for *t* < *t*∗, by defining *t*

*m* = 0. Since *t*∗ is positive, this occurs at a temperature greater than *Tc*. For *t* < *t*

For *H* = 0, the equilibrium value of *m* is obtained as

∑ *i*,*j*

a correlation function *C* and the interaction range

where we assume the correlation function being written as

*JijSiSj* = ∑

*i Si* ∑ *j Jij* �*Sj*� + 

*εij* =

 

by introducing the cubic term in *m*.

transition.

Hamiltonian Eq. (18)

replacement can be evaluated by

$$
\varepsilon\_R \simeq \left(\frac{T}{T\_\mathbb{C}}\right) \left(\frac{a}{R}\right)^d. \tag{36}
$$

where *a* is the lattice constant and *d* is the dimensionality of the interaction. In Eq. (36), (*a*/*R*) <sup>−</sup>*<sup>d</sup>* is essentially the corrdination number *<sup>z</sup>* > 1, so that *<sup>ε</sup><sup>R</sup>* < 1 and the mean field theory is self-consistent.

On the other hand, the correlation length grows toward infinity near the critical point; *R* � *ξ* for *t* → 0. A simple arithematics yields *m* ∼ |*t*| *<sup>β</sup>*, where a *critial exponent β* is <sup>1</sup> <sup>2</sup> for a ferromagnet. This result leaves us the error

$$
\varepsilon\_R \sim \frac{1}{|t|^{2\beta}} \left(\frac{a}{R}\right)^d,\tag{37}
$$

which tends to infinity as *t* → 0. Hence, the Landau theory based on the mean-field approximation has error which diverges as the system approaches to the critical point. Mathematically, the Landau theory expands the Landau function in terms of the order parameter. The landau expansion itself is mathematically non-sense *near* the critical point for dimensions less than four. Therefore, the Landau theory is not a good tool to investigate significantly the phase transitions of the system.

#### **3. Matters as noninteracting gases**

Materials are basically made of atoms; an atom is composed of a nucleus and the surrounding electrons. However, it is convenient to distinguish two types of electrons; the *valence electrons* are responsible for chemical reactions and the *core electrons* are tightly bound around the nucleus to form an *ion* for screening the strongly divergent Coulomb potential from the nucleus. It is customary to call valence electrons as electrons.

The decomposition into electrons and ions provides us at least two advantages in treating materials with first-principles. First of all, the motions of electrons can be decoupled adiabatically from the those of ions, since electrons reach their equilibrium almost immediately by their light mass compared to those factors of ions. The decoupling of the motions of electrons from those of ions is accomplished by the Born-Oppenheimer adiabatic approximation (Born & Oppenheimer, 1927), which decouples the motions of electrons approximately begin independent adiabatically from those of ions. In practice, the motions of electrons are computed under the external potential influenced by the ions at their *static* equilibrium positions, before the motions of ions are computed under the external potential influenced by the electronic distribution. Hence, the fundamental information for thermodynamics of a material is its electronic structures. Secondly, the decoupled electrons of spin half are identical particles following the Fermi-Dirac statistics (Dirac, 1926; Fermi, 1926). Hence, the statistical distribution function of electrons is a closed fixed form. This feature reduces the burdens of calculation of the distribution function of electrons.

#### **3.1 Electronic subsystem as Fermi gas**

The consequence of the decoupling electrons from ions allows us to treat the distribution functions of distinguishable atoms, for example, an iron atom is distinguished from a carbon atom, can be treated as the source of external potential to the electronic subsystem. Modelling of electronic subsystem was suggested firstly by Drude (1900), before the birth of quantum mechanics. He assumes that a metal is composed of electrons wandering on the positive homogeneous ionic background. The interaction between electrons are cancelled to allow us

and the chemical potential from the relation *N* = (*∂* (*PV*) /*∂μ*)*TV* as

 <sup>1</sup> <sup>−</sup> *<sup>π</sup>*<sup>2</sup>

*Vμ*

It is thus the heat capacity of the noninteracting homogeneous electron gas to be

 *∂S ∂T* 

 *μ* 0 *�* 3 <sup>2</sup> <sup>d</sup>*�* <sup>=</sup> *<sup>g</sup>* 4*π*<sup>2</sup>

<sup>12</sup> <sup>1</sup> *β�*<sup>F</sup>

Towards the Authentic *Ab Intio* Thermodynamics 551

<sup>=</sup> *gV* 4*π*<sup>2</sup>

*VN*

The internal energy is simply calculated by a summation of the microstate energy of all the

All the necessary thermodynamic information of the homogeneous noninteracting electron

For the case of ions, the treatment is rather complex. One can immediately raise the same treatment of the homogeneous noninteracting ionic gas model as we did for the electronic subsystem. Ignoring the nuclear spins, any kinds of ions are composed of fully occupied electronic shells to yield the effective zero spin; ions are massive bosons. It seems, if the system has single elemental atoms, that the ionic subsystem can be treated as an indistinguishable homogeneous noninteracting bosonic gas, following the Bose-Einstein statistics (Bose, 1926; Einstein, 1924; 1925). However, the ionic subsystem is hardly treated as a boson gas. Real materials are not elemental ones, but they are composed of many different kinds of elements; it is possible to distinguish the atoms. They are partially distinguishable each other, so that a combinatorial analysis is required for calculating thermodynamic properties (Ruban & Abrikosov, 2008; Turchi et al., 2007). It is obvious that the ions in a material are approximately distributed in the space isotropically and homogeneously. Such phases are usually called *fluids*. As temperature goes down, the material in our interests usually crystalizes where the homogeneous and isotropic symmetries are broken spontaneously and

In quantum field theoretical language, there is a massless boson, called Goldstone boson, if the Lagrangian of the system possesses a continuous symmetry group under which the the ground or vacuum state is not invariant (Goldstone, 1961; Goldstone et al., 1962). For example, phonons are emerged by the violation of translational and rotational symmetry of the solid crystal; a longitudinal phonon is emerged by the violation of the gauge invariance in liquid helium; spin waves, or magnons, are emerged by the violation of spin rotation symmetry (Anderson, 1963). These quasi-particles, or elementary excitations, have known in many-body theory for solids (Madelung, 1978; Pines, 1962; 1999). One has to note two facts: (i) the elementary exciations are not necessarily to be a Goldstone boson and (ii) they are not

<sup>2</sup>

2*m h*¯

<sup>=</sup> *<sup>π</sup>*<sup>2</sup> <sup>2</sup> *Nk*<sup>B</sup>

 3 <sup>2</sup> 2 3 2*π*<sup>2</sup> 4 *k*B *<sup>β</sup>* <sup>+</sup> ···

2*m h*¯ 2

1

3 <sup>2</sup> 2 5 *μ* 5

<sup>+</sup> ···

. (45)

*�*F*<sup>β</sup>* . (47)

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

. (46)

*μ* = *�*<sup>F</sup>

 *∂* (*PV*) *∂T*

*CV* = *T*

The low temperature limit entropy *S* is calculated as

*E <sup>V</sup>* <sup>=</sup> *<sup>g</sup>* 4*π*<sup>2</sup>

**3.2 Elementary excitation as massive boson gas**

individual atoms all occupy nearly fixed positions.

*S* (*β*, *V*, *μ*) =

occupied states to yield

subsystem is acquired.

for treating the electrons as a noninteracting gas. Albeit the Drude model oversimplifies the real situation, it contains many useful features of the fundamental properties of the electronic subsystem (Aschcroft & Mermin, 1976; Fetter & Walecka, 2003; Giuliani & Vignale, 2005). As microstates is indexed as *i* of the electron subsystem, the Fermi-Dirac distribution function is written in terms of occupation number of the state *i*,

$$m\_i^0 = \frac{1}{e^{\beta(\epsilon\_i - \mu)} + 1},\tag{38}$$

where *�<sup>i</sup>* is the energy of the electronic microstate *i* and *μ* is the chemical potential of the electron gas. At zero-temperature, the Fermi-Dirac distribution function becomes

$$\frac{1}{e^{\beta(\varepsilon-\mu)}+1} = \theta \left(\mu-\varepsilon\right) \tag{39}$$

and the chemical potential becomes the Fermi energy *�*F. In the high-temperature limit, the Fermi-Dirac distribution function recudes to

$$n^0 = e^{\beta(\varepsilon-\mu)},\tag{40}$$

the Maxwell-Boltzmann distribution function. With the nonrelativistic energy spectrum

$$
\epsilon\_\mathbf{p} = \frac{\mathbf{p}^2}{2m} = \frac{\hbar^2 k^2}{2m} = \epsilon\_\mathbf{k} \tag{41}
$$

where **p** is the single-particle momentum, **k** is the corresponding wave vector, the grand potential in Eq. (3) is calculated in a continuum limit2 as

$$-\beta \Omega\_0 = \beta PV = \frac{2}{3} \frac{gV}{3\pi^2} \left(\frac{2m}{\hbar^2}\right)^{\frac{2}{2}} \int\_0^\infty \text{d}\varepsilon \frac{\varepsilon^{\frac{3}{2}}}{e^{\beta(\varepsilon-\mu)}+1} \tag{42}$$

and the number density is written as

$$\frac{N}{V} = \frac{g}{4\pi^2} \left(\frac{2m}{\hbar^2}\right)^{\frac{3}{2}} \int\_0^\infty \mathbf{d}\epsilon \frac{\varepsilon^{\frac{1}{2}}}{e^{\beta(\varepsilon-\mu)}+1},\tag{43}$$

where *g* is 2, the degeneracy factor of an electron. After math (Fetter & Walecka, 2003), we can obtain the low-temperature limit (*T* → 0) of the grand potential of the noninteractic homogeneous electron gas as

$$PV = \frac{2}{3} \frac{gV}{4\pi^2} \left(\frac{2m}{\hbar}\right)^{\frac{3}{2}} \left[\frac{2}{5}\mu^{\frac{5}{2}} + \beta^{-2} \frac{\pi^2}{4}\mu^{\frac{1}{2}} + \dotsb\right] \tag{44}$$

<sup>2</sup> It is convenient to convert a summation over single-particle spectra to an integral over wavenumbers according to ∑*<sup>i</sup>* → *g* d3*n* = *gV* (2*π*) <sup>−</sup><sup>3</sup> d3*k* for a very large periodic system, hence a continuum case. If we have knowledge of the single-particle energy dispersion relation, the wavenumber integral is also replaced by an integral over energy as *gV* (2*π*) <sup>−</sup><sup>3</sup> <sup>d</sup>3*kF* (*�***k**) <sup>→</sup> *<sup>g</sup>* <sup>∞</sup> <sup>−</sup><sup>∞</sup> <sup>d</sup>*�*<sup>D</sup> (*�*) *<sup>F</sup>* (*�*), where <sup>D</sup> (*�*) is the density of states.

8 Will-be-set-by-IN-TECH

for treating the electrons as a noninteracting gas. Albeit the Drude model oversimplifies the real situation, it contains many useful features of the fundamental properties of the electronic subsystem (Aschcroft & Mermin, 1976; Fetter & Walecka, 2003; Giuliani & Vignale, 2005). As microstates is indexed as *i* of the electron subsystem, the Fermi-Dirac distribution function

*<sup>i</sup>* <sup>=</sup> <sup>1</sup>

*eβ*(*�i*−*μ*) + 1

where *�<sup>i</sup>* is the energy of the electronic microstate *i* and *μ* is the chemical potential of the

and the chemical potential becomes the Fermi energy *�*F. In the high-temperature limit, the

*n*<sup>0</sup> = *eβ*(*�*−*μ*)

<sup>2</sup>*<sup>m</sup>* <sup>=</sup> *<sup>h</sup>*¯ <sup>2</sup>*k*<sup>2</sup>

where **p** is the single-particle momentum, **k** is the corresponding wave vector, the grand

2*m h*¯ 2

 3 <sup>2</sup> ∞ 0 d*�*

> *�* 1 2 *eβ*(*�*−*μ*) + 1

> > 4 *μ* 1

<sup>2</sup> + ···

<sup>−</sup><sup>3</sup> d3*k* for a very large periodic system, hence a continuum case.

<sup>∞</sup>

<sup>2</sup> <sup>+</sup> *<sup>β</sup>*−<sup>2</sup> *<sup>π</sup>*<sup>2</sup>

<sup>−</sup><sup>3</sup> <sup>d</sup>3*kF* (*�***k**) <sup>→</sup> *<sup>g</sup>*

the Maxwell-Boltzmann distribution function. With the nonrelativistic energy spectrum

*�***<sup>p</sup>** <sup>=</sup> **<sup>p</sup>**<sup>2</sup>

3 *gV* 3*π*<sup>2</sup>

2*m h*¯ 2

2*m h*¯

3 <sup>2</sup> 2 5 *μ* 5

 3 <sup>2</sup> ∞ 0 d*�*

where *g* is 2, the degeneracy factor of an electron. After math (Fetter & Walecka, 2003), we can obtain the low-temperature limit (*T* → 0) of the grand potential of the noninteractic

<sup>2</sup> It is convenient to convert a summation over single-particle spectra to an integral over wavenumbers

If we have knowledge of the single-particle energy dispersion relation, the wavenumber integral is also

, (38)

, (40)

<sup>2</sup>*<sup>m</sup>* <sup>=</sup> *�***k**, (41)

*<sup>e</sup>β*(*�*−*μ*) <sup>+</sup> <sup>1</sup> (42)

, (43)

<sup>−</sup><sup>∞</sup> <sup>d</sup>*�*<sup>D</sup> (*�*) *<sup>F</sup>* (*�*), where <sup>D</sup> (*�*) is

(44)

*�* 3 2

*<sup>e</sup>β*(*�*−*μ*) <sup>+</sup> <sup>1</sup> <sup>=</sup> *<sup>θ</sup>* (*<sup>μ</sup>* <sup>−</sup> *�*) (39)

*n*0

electron gas. At zero-temperature, the Fermi-Dirac distribution function becomes 1

is written in terms of occupation number of the state *i*,

Fermi-Dirac distribution function recudes to

and the number density is written as

homogeneous electron gas as

according to ∑*<sup>i</sup>* → *g*

the density of states.

potential in Eq. (3) is calculated in a continuum limit2 as

<sup>−</sup> *<sup>β</sup>*Ω<sup>0</sup> <sup>=</sup> *<sup>β</sup>PV* <sup>=</sup> <sup>2</sup>

*N <sup>V</sup>* <sup>=</sup> *<sup>g</sup>* 4*π*<sup>2</sup>

*PV* <sup>=</sup> <sup>2</sup> 3 *gV* 4*π*<sup>2</sup>

d3*n* = *gV* (2*π*)

replaced by an integral over energy as *gV* (2*π*)

and the chemical potential from the relation *N* = (*∂* (*PV*) /*∂μ*)*TV* as

$$
\mu = \epsilon\_{\rm F} \left[ 1 - \frac{\pi^2}{12} \left( \frac{1}{\beta \epsilon\_{\rm F}} \right)^2 + \dotsb \right]. \tag{45}
$$

The low temperature limit entropy *S* is calculated as

$$S\left(\beta, V, \mu\right) = \left(\frac{\partial \left(PV\right)}{\partial T}\right)\_{V\mu} = \frac{gV}{4\pi^2} \left(\frac{2m}{\hbar}\right)^{\frac{3}{2}} \frac{2}{3} \left[\frac{2\pi^2}{4} \frac{k\_B}{\beta} + \cdots\right].\tag{46}$$

It is thus the heat capacity of the noninteracting homogeneous electron gas to be

$$\mathbf{C}\_{V} = T \left( \frac{\partial \mathbf{S}}{\partial T} \right)\_{VN} = \frac{\pi^2}{2} N k\_{\mathbf{B}} \frac{1}{\epsilon\_{\mathbf{F}} \beta}. \tag{47}$$

The internal energy is simply calculated by a summation of the microstate energy of all the occupied states to yield

$$\frac{E}{V} = \frac{\mathcal{g}}{4\pi^2} \int\_0^\mu \epsilon^{\frac{3}{2}} \text{d}\epsilon = \frac{\mathcal{g}}{4\pi^2} \left(\frac{2m}{\hbar^2}\right)^{\frac{3}{2}} \frac{2}{5} \mu^{\frac{5}{2}}.\tag{48}$$

All the necessary thermodynamic information of the homogeneous noninteracting electron subsystem is acquired.

#### **3.2 Elementary excitation as massive boson gas**

For the case of ions, the treatment is rather complex. One can immediately raise the same treatment of the homogeneous noninteracting ionic gas model as we did for the electronic subsystem. Ignoring the nuclear spins, any kinds of ions are composed of fully occupied electronic shells to yield the effective zero spin; ions are massive bosons. It seems, if the system has single elemental atoms, that the ionic subsystem can be treated as an indistinguishable homogeneous noninteracting bosonic gas, following the Bose-Einstein statistics (Bose, 1926; Einstein, 1924; 1925). However, the ionic subsystem is hardly treated as a boson gas. Real materials are not elemental ones, but they are composed of many different kinds of elements; it is possible to distinguish the atoms. They are partially distinguishable each other, so that a combinatorial analysis is required for calculating thermodynamic properties (Ruban & Abrikosov, 2008; Turchi et al., 2007). It is obvious that the ions in a material are approximately distributed in the space isotropically and homogeneously. Such phases are usually called *fluids*. As temperature goes down, the material in our interests usually crystalizes where the homogeneous and isotropic symmetries are broken spontaneously and individual atoms all occupy nearly fixed positions.

In quantum field theoretical language, there is a massless boson, called Goldstone boson, if the Lagrangian of the system possesses a continuous symmetry group under which the the ground or vacuum state is not invariant (Goldstone, 1961; Goldstone et al., 1962). For example, phonons are emerged by the violation of translational and rotational symmetry of the solid crystal; a longitudinal phonon is emerged by the violation of the gauge invariance in liquid helium; spin waves, or magnons, are emerged by the violation of spin rotation symmetry (Anderson, 1963). These quasi-particles, or elementary excitations, have known in many-body theory for solids (Madelung, 1978; Pines, 1962; 1999). One has to note two facts: (i) the elementary exciations are not necessarily to be a Goldstone boson and (ii) they are not

The classical chemical potential *μc* is now calculated as

1 *β*0 <sup>=</sup> *<sup>h</sup>*¯ <sup>2</sup> 2*m*

counted, using the Bose-Einstein distribution function Eq. (49), by

of the Bose particles with energies *�* > 0 is computed by Eq. (53) to be

*N�*=<sup>0</sup> *<sup>V</sup>* <sup>=</sup> *<sup>N</sup> V*

massive boson gas for *β* > *β*<sup>0</sup> is then computed (Fowler & Jones, 1938) as

� 2*m βh*¯ <sup>2</sup>

while the number density at the ground state is evaluated to be

*E <sup>V</sup>* <sup>=</sup> *<sup>g</sup>* 4*π*<sup>2</sup>

*N�*><sup>0</sup> *<sup>V</sup>* <sup>=</sup> *<sup>N</sup> V*

⎛ ⎝

*g*Γ �<sup>3</sup> 2 � *ζ* � <sup>3</sup> 2 �

*N* = ∑ *i*

with *μ* = 0 to be

integrand relative to its value at *β*0.

fact that *�<sup>i</sup>* = 0 vanishes the denominator *�*

*βμ<sup>c</sup>* = ln

⎡ ⎣ *N gV* �

2*πh*¯ <sup>2</sup> *m*

As *β* increases at fixed density, *βμc* passes through zero and becomes positive, diverging to infinity at *β* → ∞. This contradicts to the requirement Eq. (54). The critical temperature *β*0, where the chemical potential of an ideal boson gas vanishes, is calculated by using Eq. (53)

Towards the Authentic *Ab Intio* Thermodynamics 553

4*π*<sup>2</sup>

where Γ and *ζ* are Gamma function and zeta function, respectively. For *μ* = 0 and *β* > *β*0, the integral in Eq. (53) is less than *N*/*V* because these conditions increase the denominator of the

The breakdown of the theory was noticed by Einstein (1925) and was traced origin of the breakdown was the converting the conversion of the summation to the integral of the occupation number counting in Eq. (53). The total number of the ideal massive Bose gas is

It is obvious that the bosons tends to occupy the ground state for the low temperature range *β* > *β*0, due to the lack of the limitation of the occupation number of bosons. As temperature goes down, the contribution of the ground state occupation to the number summation increases. However, the first term of Eq. (60) is omitted, in the Bose-Einstein distribution, during the conversion to the integral Eq. (53) as *μ* → 0<sup>−</sup> for *β* > *β*0, because the

1

� 1 − � *β β*0

with the chemical potential *μ* = 0− for *β* > *β*0. The internal energy density of the degenerate

�3 <sup>2</sup> 1 *β* Γ �5 2 � *ζ* �5 2 �

1 *<sup>e</sup>β*(*�i*−*μ*) <sup>−</sup> <sup>1</sup>

> � *β β*0

�<sup>−</sup> <sup>3</sup> 2

> �<sup>−</sup> <sup>3</sup> 2 �

�3 2 *β* 3 2 ⎤

⎞ ⎠ 2 <sup>3</sup> � *<sup>N</sup> V* �2 3

⎦ . (58)

, (59)

. (60)

, (61)

, (62)

. (63)

<sup>2</sup> of the integrand in Eq. (53). The number density

necessarily limited to the ionic subsystem, but also electronic one. If the elementary excitations are fermionic, thermodynamics are basically calculable as we did for the non-interacting electrons gas model, in the beginning of this section. If the elementary excitations are (Goldstone) bosonic, such as phonons or magnons, a thermodynamics calculation requires special care. In order to illustrative purpose, let us see the thermodynamic information of a system of homogeneous noninteracting massive bosons.

The Bose-Einstein distribution function gives the mean occupation number in the *i*th state as

$$m\_i^0 = \frac{1}{e^{\beta(\epsilon\_i - \mu)} - 1}.\tag{49}$$

Since the chemical potential of a bosonic system vanishes at a certain temperature *T*0, a special care is necessary during the thermodynamic property calculations (Cornell & Wieman, 2002; Einstein, 1925; Fetter & Walecka, 2003). The grand potential of an ideal massive boson gas, where the energy spectrum is also calculated as in Eq. (41), is

$$-\beta \Omega\_0 = \beta PV = -\frac{gV}{4\pi^2} \left(\frac{2m}{\hbar^2}\right)^{\frac{3}{2}} \int\_0^\infty \mathbf{d} \epsilon \,\epsilon^{\frac{1}{2}} \ln\left(1 - e^{\beta(\mu - \varepsilon)}\right). \tag{50}$$

The integration by part yields

$$PV = \frac{gV}{4\pi^2} \left(\frac{2m}{\hbar^2}\right)^{\frac{3}{2}} \frac{2}{3} \int\_0^\infty \mathrm{d}\varepsilon \frac{e^{\frac{3}{2}}}{e^{\beta(\varepsilon-\mu)}-1}.\tag{51}$$

The internal energy is calculated to be

$$E = \sum\_{i} n\_i^0 \varepsilon\_i = \frac{3}{2} PV = \frac{gV}{4\pi^2} \left(\frac{2m}{\hbar^2}\right)^{\frac{3}{2}} \int\_0^\infty \mathrm{d}\varepsilon \frac{\varepsilon^{\frac{3}{2}}}{e^{\beta(\varepsilon-\mu)} - 1},\tag{52}$$

and the number density is calculated to be

$$\frac{N}{V} = \frac{g}{4\pi^2} \left(\frac{2m}{\hbar^2}\right)^{\frac{3}{2}} \int\_0^\infty \mathrm{d}\epsilon \frac{\varepsilon^{\frac{1}{2}}}{e^{\beta(\mu-\epsilon)} - 1}.\tag{53}$$

A care is necessary in treating Eq. (53), because it is meaningful only if *�* − *μ* ≥ 0, or

$$
\mu \le 0 \tag{54}
$$

with the consideration of the fact *�* ≥ 0. In the classical limit *T* → ∞, or *β* → 0, for fixed *N*, we have

$$
\beta \mu \to -\infty.\tag{55}
$$

Recall that the classical limit yields the Maxwell-Boltzmann distribution

$$m\_i^0 = e^{-\beta(\epsilon\_i - \mu)}\tag{56}$$

for both fermions and bosons, and the corresponding grand potential becomes

$$
\Omega\_0 = -PV = -\frac{1}{\beta} \sum\_i e^{\beta(\mu - \varepsilon\_i)}.\tag{57}
$$

10 Will-be-set-by-IN-TECH

necessarily limited to the ionic subsystem, but also electronic one. If the elementary excitations are fermionic, thermodynamics are basically calculable as we did for the non-interacting electrons gas model, in the beginning of this section. If the elementary excitations are (Goldstone) bosonic, such as phonons or magnons, a thermodynamics calculation requires special care. In order to illustrative purpose, let us see the thermodynamic information of a

The Bose-Einstein distribution function gives the mean occupation number in the *i*th state as

Since the chemical potential of a bosonic system vanishes at a certain temperature *T*0, a special care is necessary during the thermodynamic property calculations (Cornell & Wieman, 2002; Einstein, 1925; Fetter & Walecka, 2003). The grand potential of an ideal massive boson gas,

> 3 <sup>2</sup> ∞ 0

2*m h*¯ 2

 3 <sup>2</sup> ∞ 0 d*�*

<sup>d</sup>*� �* <sup>1</sup> <sup>2</sup> ln 

> *�* 3 2 *<sup>e</sup>β*(*�*−*μ*) <sup>−</sup> <sup>1</sup>

*�* 1 2 *<sup>e</sup>β*(*μ*−*�*) <sup>−</sup> <sup>1</sup>

*<sup>e</sup>β*(*�i*−*μ*) <sup>−</sup> <sup>1</sup>

. (49)

. (50)

. (51)

. (53)

. (57)

, (52)

<sup>1</sup> <sup>−</sup> *<sup>e</sup>β*(*μ*−*�*)

*�* 3 2 *<sup>e</sup>β*(*�*−*μ*) <sup>−</sup> <sup>1</sup>

*μ* ≤ 0 (54)

*βμ* → −∞. (55)

*<sup>i</sup>* <sup>=</sup> *<sup>e</sup>*−*β*(*�i*−*μ*) (56)

*<sup>i</sup>* <sup>=</sup> <sup>1</sup>

2*m h*¯ 2

> 3 <sup>2</sup> 2 3 ∞ 0 d*�*

*n*0

4*π*<sup>2</sup>

2*m h*¯ 2

*PV* <sup>=</sup> *gV* 4*π*<sup>2</sup>

> 2*m h*¯ 2

A care is necessary in treating Eq. (53), because it is meaningful only if *�* − *μ* ≥ 0, or

 3 <sup>2</sup> ∞ 0 d*�*

system of homogeneous noninteracting massive bosons.

where the energy spectrum is also calculated as in Eq. (41), is

*PV* <sup>=</sup> *gV* 4*π*<sup>2</sup>

> *N <sup>V</sup>* <sup>=</sup> *<sup>g</sup>* 4*π*<sup>2</sup>

In the classical limit *T* → ∞, or *β* → 0, for fixed *N*, we have

Recall that the classical limit yields the Maxwell-Boltzmann distribution

*n*0

<sup>Ω</sup><sup>0</sup> <sup>=</sup> <sup>−</sup>*PV* <sup>=</sup> <sup>−</sup> <sup>1</sup>

*<sup>β</sup>* ∑ *i*

*eβ*(*μ*−*�i*)

for both fermions and bosons, and the corresponding grand potential becomes

<sup>−</sup> *<sup>β</sup>*Ω<sup>0</sup> <sup>=</sup> *<sup>β</sup>PV* <sup>=</sup> <sup>−</sup> *gV*

The integration by part yields

The internal energy is calculated to be

*E* = ∑ *i n*0 *<sup>i</sup> �<sup>i</sup>* <sup>=</sup> <sup>3</sup> 2

and the number density is calculated to be

with the consideration of the fact *�* ≥ 0.

The classical chemical potential *μc* is now calculated as

$$
\beta \mu\_c = \ln \left[ \frac{N}{gV} \left( \frac{2\pi \hbar^2}{m} \right)^{\frac{3}{2}} \beta^{\frac{3}{2}} \right]. \tag{58}
$$

As *β* increases at fixed density, *βμc* passes through zero and becomes positive, diverging to infinity at *β* → ∞. This contradicts to the requirement Eq. (54). The critical temperature *β*0, where the chemical potential of an ideal boson gas vanishes, is calculated by using Eq. (53) with *μ* = 0 to be

$$\frac{1}{g\_0} = \frac{\hbar^2}{2m} \left( \frac{4\pi^2}{g^\Gamma \left(\frac{3}{2}\right) \zeta\left(\frac{3}{2}\right)} \right)^{\frac{2}{3}} \left( \frac{N}{V} \right)^{\frac{2}{3}}.\tag{59}$$

where Γ and *ζ* are Gamma function and zeta function, respectively. For *μ* = 0 and *β* > *β*0, the integral in Eq. (53) is less than *N*/*V* because these conditions increase the denominator of the integrand relative to its value at *β*0.

The breakdown of the theory was noticed by Einstein (1925) and was traced origin of the breakdown was the converting the conversion of the summation to the integral of the occupation number counting in Eq. (53). The total number of the ideal massive Bose gas is counted, using the Bose-Einstein distribution function Eq. (49), by

$$N = \sum\_{i} \frac{1}{e^{\beta(\varepsilon\_i - \mu)} - 1}. \tag{60}$$

It is obvious that the bosons tends to occupy the ground state for the low temperature range *β* > *β*0, due to the lack of the limitation of the occupation number of bosons. As temperature goes down, the contribution of the ground state occupation to the number summation increases. However, the first term of Eq. (60) is omitted, in the Bose-Einstein distribution, during the conversion to the integral Eq. (53) as *μ* → 0<sup>−</sup> for *β* > *β*0, because the fact that *�<sup>i</sup>* = 0 vanishes the denominator *�* 1 <sup>2</sup> of the integrand in Eq. (53). The number density of the Bose particles with energies *�* > 0 is computed by Eq. (53) to be

$$\frac{N\_{\varepsilon>0}}{V} = \frac{N}{V} \left(\frac{\beta}{\beta\_0}\right)^{-\frac{3}{2}},\tag{61}$$

while the number density at the ground state is evaluated to be

$$\frac{N\_{\varepsilon=0}}{V} = \frac{N}{V} \left[ 1 - \left( \frac{\beta}{\beta\_0} \right)^{-\frac{3}{2}} \right],\tag{62}$$

with the chemical potential *μ* = 0− for *β* > *β*0. The internal energy density of the degenerate massive boson gas for *β* > *β*<sup>0</sup> is then computed (Fowler & Jones, 1938) as

$$\frac{E}{V} = \frac{g}{4\pi^2} \left(\frac{2m}{\beta\hbar^2}\right)^{\frac{3}{2}} \frac{1}{\beta} \Gamma\left(\frac{5}{2}\right) \zeta\left(\frac{5}{2}\right). \tag{63}$$

In these conditions, the mean occupation number is following the Planck distribution function

Towards the Authentic *Ab Intio* Thermodynamics 555

which was originally suggested for describing the distribution function problem of black-body radiations. Note that the Planck distribution function is a special case of the

Considering the relation Eq. (4) and the condition Eq. (69), the grand potential of a massless Bose gas subsystem becomes the same as the Helmholtz free energy. The details of the thermodynamic properties are depending on the dispersion relation of the bosons. Photons are quantized radiations based on the fact of the *linearity* of electrodynamics,<sup>3</sup> so that photons

where *k* includes the definite polarizations. The photon Helmholtz free energy is calculated

<sup>d</sup>*ωω*<sup>2</sup> ln

*k*4 B <sup>1</sup> <sup>−</sup> *<sup>e</sup>*−*βh*¯ *<sup>ω</sup>*

 ∞ 0

*<sup>F</sup>* <sup>=</sup> <sup>−</sup><sup>4</sup> 3 *σV c*

*<sup>σ</sup>* <sup>=</sup> *<sup>π</sup>*<sup>2</sup> 60

> *<sup>∂</sup><sup>T</sup>* <sup>=</sup> <sup>16</sup> 3 *σV c*

which propotional to the fourth power of the temperature; *Boltzmann's law*. The constant

<sup>=</sup> <sup>16</sup>*σ<sup>V</sup> c*

> <sup>=</sup> <sup>4</sup>*<sup>σ</sup>* 3*c*

*<sup>S</sup>* <sup>=</sup> <sup>−</sup> *<sup>∂</sup><sup>F</sup>*

*<sup>E</sup>* <sup>=</sup> <sup>4</sup>*σ<sup>V</sup> c*

> *∂E ∂T V*

> > *∂F ∂V T*

*PV* <sup>=</sup> <sup>1</sup> 3

<sup>3</sup> The nonlinear character appeared in the qunatum electrodynamics will not be discussed here.

*CV* =

*P* = −

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

, (70)

*ω* = *ck*, (71)

*T*4, (73)

*<sup>h</sup>*¯ <sup>3</sup>*c*<sup>2</sup> . (74)

*<sup>T</sup>*<sup>4</sup> <sup>=</sup> <sup>−</sup>3*F*, (76)

*T*3. (75)

*T*3, (77)

*T*4. (78)

*E*. (79)

. (72)

**<sup>k</sup>** <sup>=</sup> <sup>1</sup>

*n*0

Bose-Einstein distribution function with zero chemical potential.

*<sup>F</sup>*<sup>0</sup> <sup>=</sup> *<sup>T</sup> <sup>V</sup> π*2*c*<sup>3</sup>

where *σ* is called the *Stefan-Boltzmann constant* defined as

The standard integration method yields

The total radiation energy *E* = *F* + *TS* is

volume heat capacity of the radiation is

Hence, the equation of states of the photon gas is

do not interact with one another. The photon dispersion relation is linear,

(Planck, 1901)

as

The entropy is

and the pressure is

The constant-volume heat capacity for *β* > *β*<sup>0</sup> becomes

$$\mathbf{C}\_{V} = \frac{5}{2} \left[ \frac{\Gamma\left(\frac{5}{2}\right) \zeta\left(\frac{5}{2}\right)}{\Gamma\left(\frac{3}{2}\right) \zeta\left(\frac{3}{2}\right)} Nk\_{\rm B} \left(\frac{\beta}{\beta\_{0}}\right)^{-\frac{3}{2}} \right],\tag{64}$$

and the pressure for *β* > *β*<sup>0</sup> becomes

$$P = \frac{2}{3} \frac{2\sqrt{2}}{4\pi^2} \Gamma\left(\frac{5}{2}\right) \zeta\left(\frac{5}{2}\right) \frac{m^{\frac{3}{2}}g}{\hbar^3} \beta^{-\frac{5}{2}}.\tag{65}$$

It is interesting to find that the pressure approaches to zero as temperature goes to zero, *i.e. β* → ∞. In other words, the Bose gas exerts no force on the walls of the container at *T* = 0, because all of the particles condensate in the zero-momentum state. The pressure is independent of the number density *N*/*V*, depending only on temperature *β*. The two different summations above and below the temperature *β*<sup>0</sup> lead us that the heat capacity as a function of temperature has discontinuity in its slope (Landau & Lifshitz, 1980) as

$$
\Delta \left( \frac{\partial \mathbb{C}\_V}{\partial V} \right)\_{\beta\_0} = -\frac{27}{4} \left[ \frac{\Gamma \left( \frac{3}{2} \right) \zeta \left( \frac{3}{2} \right)}{\pi} \right]^2 N k\_\text{B}^2 \beta\_0^{-1}. \tag{66}
$$

This implies that a homogeneous massive ideal Bose gas system exhibits a *phase transition* at *β*<sup>0</sup> *without* interaction. This phenomenon is known as the *Bose-Einstein condensation*. A good review for the realization of the Bose-Einstein condensation is provided by Cornell & Wieman (2002).

#### **3.3 Elementary excitations as massless boson gas**

As we stated previously, some elementary excitations emerged by the spontaneous symmetry breaking are *massless* bosons (Goldstone, 1961; Goldstone et al., 1962) as well as gauge bosons, which are elemental particles, *e.g.* photons, arosen from the fundamental interactions, *electromagnetic fields* for photons, with gauge degrees of freedom. Whether a boson is a Goldstone boson or a gauge one, the procedure described above is not appliable to the massless character of the boson, because its energy spectrum is not in the form of Eq. (41).

One has to remind that a masselss boson does not carry mass, but it carries momentum and energy. The energy spectrum of a massless boson is given by

$$
\epsilon = \hbar \omega \tag{67}
$$

and the frequency *ω* is obtained by the corresponding momentum **p** = *h*¯ **k** through a dispersion relation

$$
\omega = \omega \left( \mathbf{k} \right). \tag{68}
$$

The number of bosons *N* in the massless boson gas is a variable, and not a given constants as in an ordinary gas. Therefore, *N* itself must be determined from the thermal equilibrium condition, the (Helmholtz) free energy minimum (*∂F*/*∂N*)*T*,*<sup>V</sup>* = 0. Since (*∂F*/*∂N*)*T*,*<sup>V</sup>* = *μ*, this gives

$$
\mu = 0.\tag{69}
$$

12 Will-be-set-by-IN-TECH

� *β β*0

� *m*<sup>3</sup> <sup>2</sup> *g <sup>h</sup>*¯ <sup>3</sup> *<sup>β</sup>*<sup>−</sup> <sup>5</sup>

�<sup>−</sup> <sup>3</sup> 2 ⎤

⎦ , (64)

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

<sup>0</sup> . (66)

The constant-volume heat capacity for *β* > *β*<sup>0</sup> becomes

and the pressure for *β* > *β*<sup>0</sup> becomes

Δ � *∂CV ∂V*

**3.3 Elementary excitations as massless boson gas**

energy. The energy spectrum of a massless boson is given by

(2002).

dispersion relation

this gives

*CV* <sup>=</sup> <sup>5</sup> 2

> *<sup>P</sup>* <sup>=</sup> <sup>2</sup> 3 2 √2 <sup>4</sup>*π*<sup>2</sup> <sup>Γ</sup>

> > �

*β*0

⎡ ⎣ Γ � <sup>5</sup> 2 � *ζ* �<sup>5</sup> 2 �

Γ � <sup>3</sup> 2 � *ζ* �<sup>3</sup> 2 � *Nk*<sup>B</sup>

> �5 2 � *ζ* �5 2

function of temperature has discontinuity in its slope (Landau & Lifshitz, 1980) as

⎡ ⎣ Γ � <sup>3</sup> 2 � *ζ* �3 2 �

This implies that a homogeneous massive ideal Bose gas system exhibits a *phase transition* at *β*<sup>0</sup> *without* interaction. This phenomenon is known as the *Bose-Einstein condensation*. A good review for the realization of the Bose-Einstein condensation is provided by Cornell & Wieman

As we stated previously, some elementary excitations emerged by the spontaneous symmetry breaking are *massless* bosons (Goldstone, 1961; Goldstone et al., 1962) as well as gauge bosons, which are elemental particles, *e.g.* photons, arosen from the fundamental interactions, *electromagnetic fields* for photons, with gauge degrees of freedom. Whether a boson is a Goldstone boson or a gauge one, the procedure described above is not appliable to the massless character of the boson, because its energy spectrum is not in the form of Eq. (41). One has to remind that a masselss boson does not carry mass, but it carries momentum and

and the frequency *ω* is obtained by the corresponding momentum **p** = *h*¯ **k** through a

The number of bosons *N* in the massless boson gas is a variable, and not a given constants as in an ordinary gas. Therefore, *N* itself must be determined from the thermal equilibrium condition, the (Helmholtz) free energy minimum (*∂F*/*∂N*)*T*,*<sup>V</sup>* = 0. Since (*∂F*/*∂N*)*T*,*<sup>V</sup>* = *μ*,

*π*

⎤ ⎦ 2

*Nk*<sup>2</sup> B*β*−<sup>1</sup>

*�* = *h*¯ *ω* (67)

*ω* = *ω* (**k**). (68)

*μ* = 0. (69)

<sup>=</sup> <sup>−</sup><sup>27</sup> 4

It is interesting to find that the pressure approaches to zero as temperature goes to zero, *i.e. β* → ∞. In other words, the Bose gas exerts no force on the walls of the container at *T* = 0, because all of the particles condensate in the zero-momentum state. The pressure is independent of the number density *N*/*V*, depending only on temperature *β*. The two different summations above and below the temperature *β*<sup>0</sup> lead us that the heat capacity as a

In these conditions, the mean occupation number is following the Planck distribution function (Planck, 1901)

$$m\_{\mathbf{k}}^{0} = \frac{1}{e^{\beta \hbar \omega\_{\mathbf{k}}} - 1},\tag{70}$$

which was originally suggested for describing the distribution function problem of black-body radiations. Note that the Planck distribution function is a special case of the Bose-Einstein distribution function with zero chemical potential.

Considering the relation Eq. (4) and the condition Eq. (69), the grand potential of a massless Bose gas subsystem becomes the same as the Helmholtz free energy. The details of the thermodynamic properties are depending on the dispersion relation of the bosons. Photons are quantized radiations based on the fact of the *linearity* of electrodynamics,<sup>3</sup> so that photons do not interact with one another. The photon dispersion relation is linear,

$$
\omega = ck,\tag{71}
$$

where *k* includes the definite polarizations. The photon Helmholtz free energy is calculated as

$$F\_0 = T \frac{V}{\pi^2 c^3} \int\_0^\infty \mathbf{d}\omega \omega^2 \ln\left(1 - e^{-\beta \hbar \omega}\right) \,. \tag{72}$$

The standard integration method yields

$$F = -\frac{4}{3}\frac{\sigma V}{c}T^4,\tag{73}$$

where *σ* is called the *Stefan-Boltzmann constant* defined as

$$
\sigma = \frac{\pi^2}{60} \frac{k\_\text{B}^4}{\hbar^3 c^2}. \tag{74}
$$

The entropy is

$$S = -\frac{\partial F}{\partial T} = \frac{16}{3} \frac{\sigma V}{c} T^3. \tag{75}$$

The total radiation energy *E* = *F* + *TS* is

$$E = \frac{4\sigma V}{c}T^4 = -3F\_\prime \tag{76}$$

which propotional to the fourth power of the temperature; *Boltzmann's law*. The constant volume heat capacity of the radiation is

$$\mathcal{C}\_V = \left(\frac{\partial E}{\partial T}\right)\_V = \frac{16\sigma V}{c} T^3 \,\text{.}\tag{77}$$

and the pressure is

$$P = -\left(\frac{\partial F}{\partial V}\right)\_T = \frac{4\sigma}{\mathfrak{D}c}T^4. \tag{78}$$

Hence, the equation of states of the photon gas is

$$PV = \frac{1}{3}E.\tag{79}$$

<sup>3</sup> The nonlinear character appeared in the qunatum electrodynamics will not be discussed here.

to treat high temperature phonon thermodynamics. The thermal energy of the noninteracting

Towards the Authentic *Ab Intio* Thermodynamics 557

= *Nk*<sup>B</sup> (*βh*¯ *ω*)

On the other hand, Debye (1912) approximated that the velocity of sound is taken as constant for each phonon mode with the dispersion relation *ω* ∼ *q*. In this case, the method used for the photon gas is directly applied. At low temperature limit, the Debye extracted the *T*<sup>3</sup> law

<sup>5</sup> *nk*<sup>B</sup>

The discussions in Sec. 3 has succeed to describe the thermodynamic properties of materials in many aspects, and hence such descriptions were treated in many textbooks. However, the oversimplified model fails the many important features on the material properties. One of the important origin of such failures is due to the ignorance of the *electromagnetic* interaction among the constituent particles; electrons and ions, which carry electric charges. However, the inclusion of interactions among the particles is enormously difficult to treat. To date the quantum field theory (QFT) is known as the standard method in dealing with the interacting particles. There are many good textbooks on the QFT (Berestetskii et al., 1994; Bjorken & Drell, 1965; Doniach & Sondheimer, 1982; Fetter & Walecka, 2003; Fradkin, 1991; Gross, 1999; Itzykson & Zuber, 1980; Mahan, 2000; Negele & Orland, 1988; Parisi, 1988; Peskin & Schroeder, 1995; Zinn-Justin, 1997) in treating the interacting particles systematically in various aspects. In this article, the idea of the treatments will be reviewed briefly, instead

The idea of noninteracting particles inspires an idea to deal with the electronic subsystem as a sum of independent particles under a given potential field (Hartree, 1928), with the consideration of the effect of Pauli exclusion principle (Fock, 1930), which it is known as the exchange effect. This idea, known as the Hartree-Fock method, was mathematically formulated by introducing the Slater determinant (Slater, 1951) for the many-body electronic wave function. The individual wave function of an electron can be obtained by solving either Schrödinger equation (Schrödinger, 1926a;b;c;d) for the nonrelativistic cases or Dirac equation

Since an electron carries a fundamental electric charge *e* in its motion, it is necessary to deal with electromagnetic waves or their quanta photons. Immediate necessity was arosen in order to deal with both electrons and photons in a single quantum theoretical framework in consideration of the Einstein's special theory of relativity (Einstein, 1905). Jordan & Pauli (1928) and Heisenberg & Pauli (1929) suggested that a new formalism to treat both the

<sup>4</sup> The immediate relativistic version of the Schrödinger equation was derived by Gordon (1926) and Klein (1927), known as the Klein-Gordon equation. The Klein-Gordon equation is valid for the Bose-Einstein

particles, while the Dirac equation is valid for the Fermi-Dirac particles.

 *T* Θ<sup>D</sup> <sup>3</sup>

*CV* � <sup>12</sup>*π*<sup>4</sup>

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

<sup>2</sup> *eβh*¯ *<sup>ω</sup>* 

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

, (86)

, (88)

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

*<sup>E</sup>* <sup>=</sup> *Nn*0*h*¯ *<sup>ω</sup>* <sup>=</sup> *Nh*¯ *<sup>ω</sup>*

phonon system is

of the heat capacity as

and so the heat capacity of the system is

*CV* =

where Θ<sup>D</sup> is known as the Debye temperature.

**4. Matters as interacting liquids**

of dealing with the full details.

(Dirac, 1928a;b) for the relativistic ones.<sup>4</sup>

 *∂E ∂T V*

The procedure for the photonic subsystem is quite useful in computing thermodynamics for many kinds of elementary exciations, which are usually massless (Goldstone) bosons in a condensed matter system. It is predicted that any crystal must be completely ordered, and the atoms of each kind must occupy entirely definite positions, in a state of complete thermodynamic equilibrium (Andreev & Lifshitz, 1969; Leggett, 1970). It is well known that the ions vibrate even in the zero temperature with several vibration modes (Aschcroft & Mermin, 1976; Callaway, 1974; Jones & March, 1973a; Kittel, 2005; Landau & Lifshitz, 1980; Madelung, 1978; Pines, 1999). The energy spectrum of a phonon in a mode *j* and a wavevector **q** contains the zero vibration term

$$
\epsilon\_{i\mathbf{q}} = \left( n\_{j\mathbf{q}} + \frac{1}{2} \right) \hbar \omega\_j(\mathbf{q}) \,, \tag{80}
$$

where *nj***<sup>q</sup>** is the occupation number of the single-particle modes of *j* and **q**. The corresponding partition function is then written as

$$Z = \sum\_{\boldsymbol{\eta}\_{/} \mathbf{q}} \exp\left[-\beta \boldsymbol{\epsilon}\_{\boldsymbol{\jmath}} \left(\mathbf{q}\right)\right]$$

$$= \prod\_{\boldsymbol{\eta}, \mathbf{q}} \frac{\exp\left[-\beta \hbar \boldsymbol{\omega}\_{\boldsymbol{\jmath}} \left(\mathbf{q}\right)\right]}{1 - \exp\left[-\beta \hbar \boldsymbol{\omega}\_{\boldsymbol{\jmath}} \left(\mathbf{q}\right)\right]}.\tag{81}$$

The Helmholtz free energy of the phonon subsystem is

$$F = -\beta \ln Z = \beta \sum\_{j, \mathbf{q}} \ln \left\{ 2 \sinh \left[ \frac{\beta \hbar \omega\_j(\mathbf{q})}{2} \right] \right\}. \tag{82}$$

In the noninteracting phonon gas condition, the entropy *S*, internal energy *E*, and the volume constant specific heat *C* becomes

$$\begin{split} S &= -\left(\frac{\partial F}{\partial T}\right)\_V \\ &= k\_{\mathsf{B}} \sum\_{j,\mathbf{q}} \left[\frac{\beta}{2} \hbar \omega\_j(\mathbf{q}) \coth\left(\frac{\beta}{2} \hbar \omega\_j(\mathbf{q})\right) - \ln\left\{2 \sinh\left(\frac{\beta}{2} \hbar \omega\_j(\mathbf{q})\right)\right\}\right]. \end{split} \tag{83}$$

$$\begin{split} E &= F - TS - \left(\frac{\partial F}{\partial T}\right)\_V\\ &= \sum\_{j, \mathbf{q}} \left\{ \frac{\hbar \omega\_j \left(\mathbf{q}\right)}{2} + \frac{\hbar \omega\_j \left(\mathbf{q}\right)}{e^{\hbar \hbar \omega\_j \left(\mathbf{q}\right)} - 1} \right\}, \end{split} \tag{84}$$

$$\mathbf{C}\_{V} = \left(\frac{\partial E}{\partial T}\right)\_{V} = k\_{\mathbf{B}} \frac{\sum\_{j,\mathbf{q}} \left[\frac{\hbar \omega\_{l}(\mathbf{q})}{2}\right]^{2}}{\sinh\left[\frac{\hbar \omega\_{l}(\mathbf{q})}{2}\right]}.\tag{85}$$

Although the realistic phonon dispersion relations are complicated (Aschcroft & Mermin, 1976; Kittel, 2005), there are two useful model dispersion relations of phonons. Einstein (1907; 1911) modelled the density of states as constant frequences in each vibrational mode as D (*ω*) = *Nδ* (*ω* − *ω*0), where delta function is centered at *ω*0. This model is, in turn, useful 14 Will-be-set-by-IN-TECH

The procedure for the photonic subsystem is quite useful in computing thermodynamics for many kinds of elementary exciations, which are usually massless (Goldstone) bosons in a condensed matter system. It is predicted that any crystal must be completely ordered, and the atoms of each kind must occupy entirely definite positions, in a state of complete thermodynamic equilibrium (Andreev & Lifshitz, 1969; Leggett, 1970). It is well known that the ions vibrate even in the zero temperature with several vibration modes (Aschcroft & Mermin, 1976; Callaway, 1974; Jones & March, 1973a; Kittel, 2005; Landau & Lifshitz, 1980; Madelung, 1978; Pines, 1999). The energy spectrum of a phonon

where *nj***<sup>q</sup>** is the occupation number of the single-particle modes of *j* and **q**. The corresponding

−*β�<sup>j</sup>* (**q**)

−*βh*¯ *ω<sup>j</sup>* (**q**)

*βh*¯ *ω<sup>j</sup>* (**q**) 2

*β*

<sup>2</sup> *<sup>h</sup>*¯ *<sup>ω</sup><sup>j</sup>* (**q**)

, (84)

. (85)

−*βh*¯ *ω<sup>j</sup>* (**q**)

exp 

> exp

1 − exp

*j*,**q** ln 2 sinh 

*β*

*h*¯ *ω<sup>j</sup>* (**q**) *<sup>e</sup>βh*¯ *<sup>ω</sup>j*(**q**) <sup>−</sup> <sup>1</sup>

> *h*¯ *ωj*(**q**) 2 2

 *h*¯ *ωj*(**q**) 2

In the noninteracting phonon gas condition, the entropy *S*, internal energy *E*, and the volume

<sup>2</sup> *<sup>h</sup>*¯ *<sup>ω</sup><sup>j</sup>* (**q**)

Although the realistic phonon dispersion relations are complicated (Aschcroft & Mermin, 1976; Kittel, 2005), there are two useful model dispersion relations of phonons. Einstein (1907; 1911) modelled the density of states as constant frequences in each vibrational mode as D (*ω*) = *Nδ* (*ω* − *ω*0), where delta function is centered at *ω*0. This model is, in turn, useful

 − ln 2 sinh

*h*¯ *ω<sup>j</sup>* (**q**), (80)

. (81)

. (82)

. (83)

in a mode *j* and a wavevector **q** contains the zero vibration term

partition function is then written as

constant specific heat *C* becomes

*S* = −

= *k*<sup>B</sup> ∑ *j*,**q**

*E* = *F* − *TS* −

 *∂E ∂T V* = *k*<sup>B</sup>

= ∑ *j*,**q**

*CV* =

 *∂F ∂T V*

> *β*

*<sup>h</sup>*¯ *<sup>ω</sup><sup>j</sup>* (**q**)

*�i***<sup>q</sup>** = *nj***<sup>q</sup>** + 1 2 

*<sup>Z</sup>* <sup>=</sup> <sup>∑</sup>*nj*,**<sup>q</sup>**

= ∏ *j*,**q**

*<sup>F</sup>* = −*<sup>β</sup>* ln *<sup>Z</sup>* = *<sup>β</sup>*∑

<sup>2</sup> *<sup>h</sup>*¯ *<sup>ω</sup><sup>j</sup>* (**q**) coth

 *∂F ∂T V*

<sup>2</sup> <sup>+</sup>

∑*j*,**<sup>q</sup>**

sinh

The Helmholtz free energy of the phonon subsystem is

to treat high temperature phonon thermodynamics. The thermal energy of the noninteracting phonon system is

$$E = N n^{0} \hbar \omega = \frac{N \hbar \omega}{e^{\beta \hbar \omega} - 1},\tag{86}$$

and so the heat capacity of the system is

$$\mathbf{C}\_{V} = \left(\frac{\partial E}{\partial T}\right)\_{V} = Nk\_{\rm B} \left(\beta \hbar \omega\right)^{2} \frac{e^{\beta \hbar \omega}}{\left(e^{\beta \hbar \omega} - 1\right)^{2}}.\tag{87}$$

On the other hand, Debye (1912) approximated that the velocity of sound is taken as constant for each phonon mode with the dispersion relation *ω* ∼ *q*. In this case, the method used for the photon gas is directly applied. At low temperature limit, the Debye extracted the *T*<sup>3</sup> law of the heat capacity as

$$\mathbf{C}\_{V} \simeq \frac{12\pi^{4}}{5} \boldsymbol{n} k\_{\rm B} \left(\frac{T}{\Theta\_{\rm D}}\right)^{3},\tag{88}$$

where Θ<sup>D</sup> is known as the Debye temperature.

#### **4. Matters as interacting liquids**

The discussions in Sec. 3 has succeed to describe the thermodynamic properties of materials in many aspects, and hence such descriptions were treated in many textbooks. However, the oversimplified model fails the many important features on the material properties. One of the important origin of such failures is due to the ignorance of the *electromagnetic* interaction among the constituent particles; electrons and ions, which carry electric charges. However, the inclusion of interactions among the particles is enormously difficult to treat. To date the quantum field theory (QFT) is known as the standard method in dealing with the interacting particles. There are many good textbooks on the QFT (Berestetskii et al., 1994; Bjorken & Drell, 1965; Doniach & Sondheimer, 1982; Fetter & Walecka, 2003; Fradkin, 1991; Gross, 1999; Itzykson & Zuber, 1980; Mahan, 2000; Negele & Orland, 1988; Parisi, 1988; Peskin & Schroeder, 1995; Zinn-Justin, 1997) in treating the interacting particles systematically in various aspects. In this article, the idea of the treatments will be reviewed briefly, instead of dealing with the full details.

The idea of noninteracting particles inspires an idea to deal with the electronic subsystem as a sum of independent particles under a given potential field (Hartree, 1928), with the consideration of the effect of Pauli exclusion principle (Fock, 1930), which it is known as the exchange effect. This idea, known as the Hartree-Fock method, was mathematically formulated by introducing the Slater determinant (Slater, 1951) for the many-body electronic wave function. The individual wave function of an electron can be obtained by solving either Schrödinger equation (Schrödinger, 1926a;b;c;d) for the nonrelativistic cases or Dirac equation (Dirac, 1928a;b) for the relativistic ones.<sup>4</sup>

Since an electron carries a fundamental electric charge *e* in its motion, it is necessary to deal with electromagnetic waves or their quanta photons. Immediate necessity was arosen in order to deal with both electrons and photons in a single quantum theoretical framework in consideration of the Einstein's special theory of relativity (Einstein, 1905). Jordan & Pauli (1928) and Heisenberg & Pauli (1929) suggested that a new formalism to treat both the

<sup>4</sup> The immediate relativistic version of the Schrödinger equation was derived by Gordon (1926) and Klein (1927), known as the Klein-Gordon equation. The Klein-Gordon equation is valid for the Bose-Einstein particles, while the Dirac equation is valid for the Fermi-Dirac particles.

A quantum mechanical system is described equally well by specifying the function *K*, or by

Towards the Authentic *Ab Intio* Thermodynamics 559

Let us consider a situation that a particle propagates from 1 to 2 through 3 in a weak potential operator *U*ˆ (3), which differs from zero only for *t* between *t*<sup>1</sup> and *t*2. The kernel is expanded

To zeroth order in *U*ˆ , *K* is that for a free particle, *K*<sup>0</sup> (2, 1). Let us consider the situation if *U* differs from zero only for the infinitesimal time interval Δ*t*<sup>3</sup> between some time *t*<sup>3</sup> and *t*<sup>3</sup> + Δ*t*<sup>3</sup>

after solving the Schrödinger equation in Eq. (89). The particle at 2 then propagates *freely* from

We can decompose the Hamiltonian operator *H*ˆ by *H*ˆ <sup>0</sup> + *U*ˆ , where *H*ˆ <sup>0</sup> is the Hamiltonian

where 3� abbreviates **x**3, *t*<sup>3</sup> + Δ*t*3, by a free propagation. The difference of the wave function

where d3 = d3**x**3d*t*3. We can imagine that a particle travels as a free particle from point to point, but is scattered by the potential operator *U*ˆ at 3. The higher order terms are also

The analysis for the charged free Dirac particle gives a new interpretation of the antiparticle, which has the reversed charge of the particle; for example, a positron is the antiparticle of an electron. The Dirac equation (Dirac, 1928a;b) has negative energy states of an electron. Dirac interpreted himself that the negative energy states are fully occupied in vacuum, and an

*<sup>K</sup>* (2, 1) <sup>=</sup> *<sup>K</sup>*<sup>0</sup> (2, 1) <sup>+</sup> *<sup>K</sup>*(1) (2, 1) <sup>+</sup> *<sup>K</sup>*(2) (2, 1) <sup>+</sup> ··· . (94)

*K*<sup>0</sup> (3, 1) *ψ* (1) d3**x**1. (95)

*K*<sup>0</sup> (**x**2, *t*2; **x**3, *t*<sup>3</sup> + Δ*t*3) *ψ* (**x**3, *t*<sup>3</sup> + Δ*t*3). (97)

*U*ˆ (3) *ψ* (3) Δ*t*3. (98)

*K*<sup>0</sup> (2, 3) *U*ˆ (3) *K*<sup>0</sup> (3, 1) *ψ* (1) d3**x**3d<sup>3</sup>**x**1. (100)

*K*<sup>0</sup> (2, 3) *U*ˆ (3) *K*<sup>0</sup> (3, 1) d3, (101)

d3**x**3, (99)

*<sup>h</sup>*¯ <sup>Δ</sup>*t*<sup>3</sup>*ψ* (**x**, *t*3), (96)

specifying the Hamiltonian operator *H*ˆ from which it results.

for *t*<sup>1</sup> < *t*<sup>3</sup> < *t*2. The particle will propagate from 1 to 3 as a *free* particle,

*ψ* (3) =

For the short time interval Δ*t*3, the wave function will change to

operator of the free particle. The change in wave function by *U*ˆ will be

*ψ* (2) =

*i h*¯ 

<sup>Δ</sup>*<sup>ψ</sup>* <sup>=</sup> <sup>−</sup> *<sup>i</sup>*

The wave function at 2 is that of the propagated particle from *t*<sup>3</sup> + Δ*t*<sup>3</sup> to be

 *K*0 2, 3� *ψ* 3� 

*h*¯  *h*¯

*ψ* (**x**2, *t*2) =

Δ*ψ* (2) = −Δ*t*<sup>3</sup>

Therefore, the first order expansion of the kernel *K* is then

*<sup>K</sup>*(1) (2, 1) <sup>=</sup> <sup>−</sup> *<sup>i</sup>*

*<sup>ψ</sup>* (**x**, *<sup>t</sup>*<sup>3</sup> <sup>+</sup> <sup>Δ</sup>*t*3) <sup>=</sup> *<sup>e</sup>*−*<sup>i</sup> <sup>H</sup>*<sup>ˆ</sup>

in powers of *U*ˆ that

**x**3, *t*<sup>3</sup> + Δ*t*<sup>3</sup> as

at 2 is obtained by

analyzed in a similar way.

electrons and the radiations as quantized objects in such a way of a canonical transformation to the normal modes of their fields; this method is called the *second quantization* or the *field quantization*. The canonical transformation technique to the normal mode is a well established classical method for continuous media (Goldstein, 1980). The idea treats both the electrons and radiations as continuous fields and quantized them for their own normal modes (Heisenberg & Pauli, 1929; Jordan & Pauli, 1928). The practically available solutions was suggested for the nonrelativistic case by Bethe (1947) followed by the fully relativistic case by Dyson (1949a;b); Feynman (1949a;b); Schwinger (1948; 1949a;b); Tomonaga (1946), so called the *renormalization* for cancelling the unavoidable divergencies appeared in the quantum field theory. This kind of theory on the electrons and radiations is called as the quantum electrodynamics (QED), which is known to be the most precise theory ever achieved (Peskin & Schroeder, 1995) with the error between the theory and experiment to be less than a part per billion (ppb) (Gabrielse et al., 2006; 2007; Odom et al., 2006).

#### **4.1 The concepts of quantum field theory**

Feynman (1949b) visulized the underlying concept of the quantum field theory by reinterpreting the nonrelativistic Schrödinger equation with the Green's function concept. As in classical mechanics, a Hamiltonian operator *H*ˆ contains all the mechanical interactions of the system. The necessary physical information of the system is contained in the wave function *ψ*. The Schrödinger equation

$$i\hbar\frac{\partial\psi}{\partial t} = \hat{H}\psi,\tag{89}$$

describes the change in the wave function *ψ* in an infinitesimally time interval Δ*t* as due to the operation if an operator is *<sup>e</sup>*−*<sup>i</sup> <sup>H</sup>*<sup>ˆ</sup> *<sup>h</sup>*¯ Δ*t* . This description is equivalent to the description that the wave function *ψ* (**x**2, *t*2) at **x**<sup>2</sup> and *t*<sup>2</sup> is evolved one from the wave function *ψ* (**x**1, *t*1) at **x**<sup>1</sup> and *t*<sup>1</sup> through the equation

$$
\psi\left(\mathbf{x}\_{2},t\_{2}\right) = \int \mathbf{K}\left(\mathbf{x}\_{2},t\_{2};\mathbf{x}\_{1},t\_{1}\right)\psi\left(\mathbf{x}\_{1},t\_{1}\right)\mathbf{d}^{3}\mathbf{x}\_{1} \tag{90}
$$

where *K* is the kernel of the evolution and *t*<sup>2</sup> > *t*1. If *ψ* (**x**1, *t*1)is expanded in terms of the eigen function *φ<sup>n</sup>* with the eigenvalue *En* of the constant Hamiltonian operator *H*ˆ as ∑*<sup>n</sup> cnφ<sup>n</sup>* (**x**), one can find for *t*<sup>2</sup> > *t*<sup>1</sup>

$$K(2,1) = \sum\_{n} \phi\_{n}\left(\mathbf{x}\_{2}\right)\phi\_{n}^{\*}\left(\mathbf{x}\_{1}\right)e^{-i\frac{\mathbb{E}\_{\mu}}{\hbar}(t\_{2}-t\_{1})},\tag{91}$$

where we abbreviated 1 for **x**1, *t*<sup>1</sup> and 2 for **x**2, *t*<sup>2</sup> and define *K* (2, 1) = 0 for *t*<sup>2</sup> < *t*1. It is straightforward to show that *K* can be defined by that solution of

$$\left(i\hbar\frac{\partial}{\partial t\_2} - \hat{H}\_2\right) K\left(2, 1\right) = i\hbar\delta\left(2, 1\right),\tag{92}$$

where *<sup>δ</sup>* (2, 1) <sup>=</sup> *<sup>δ</sup>* (*t*<sup>2</sup> <sup>−</sup> *<sup>t</sup>*1) *<sup>δ</sup>* (*x*<sup>2</sup> <sup>−</sup> *<sup>x</sup>*1) *<sup>δ</sup>* (*y*<sup>2</sup> <sup>−</sup> *<sup>y</sup>*1) *<sup>δ</sup>* (*z*<sup>2</sup> <sup>−</sup> *<sup>z</sup>*1) and the subscript 2 on *<sup>H</sup>*<sup>ˆ</sup> <sup>2</sup> means that the operator acts on the variables of 2 of *K* (2, 1). The kernel *K* is now called as the *Green's function* and it is the total amplitude for arrival at **x**2, *t*<sup>2</sup> starting from **x**1, *t*1. The transition amplitude for finding a particle in state *χ* (2), if it was in *ψ* (1), is

$$
\int \chi^\* \left( 2 \right) K \left( 2, 1 \right) \psi \left( 1 \right) \mathbf{d}^3 \mathbf{x}\_1 \mathbf{d}^3 \mathbf{x}\_2. \tag{93}
$$

16 Will-be-set-by-IN-TECH

electrons and the radiations as quantized objects in such a way of a canonical transformation to the normal modes of their fields; this method is called the *second quantization* or the *field quantization*. The canonical transformation technique to the normal mode is a well established classical method for continuous media (Goldstein, 1980). The idea treats both the electrons and radiations as continuous fields and quantized them for their own normal modes (Heisenberg & Pauli, 1929; Jordan & Pauli, 1928). The practically available solutions was suggested for the nonrelativistic case by Bethe (1947) followed by the fully relativistic case by Dyson (1949a;b); Feynman (1949a;b); Schwinger (1948; 1949a;b); Tomonaga (1946), so called the *renormalization* for cancelling the unavoidable divergencies appeared in the quantum field theory. This kind of theory on the electrons and radiations is called as the quantum electrodynamics (QED), which is known to be the most precise theory ever achieved (Peskin & Schroeder, 1995) with the error between the theory and experiment to be less than

Feynman (1949b) visulized the underlying concept of the quantum field theory by reinterpreting the nonrelativistic Schrödinger equation with the Green's function concept. As in classical mechanics, a Hamiltonian operator *H*ˆ contains all the mechanical interactions of the system. The necessary physical information of the system is contained in the wave

describes the change in the wave function *ψ* in an infinitesimally time interval Δ*t* as due to

wave function *ψ* (**x**2, *t*2) at **x**<sup>2</sup> and *t*<sup>2</sup> is evolved one from the wave function *ψ* (**x**1, *t*1) at **x**<sup>1</sup> and

where *K* is the kernel of the evolution and *t*<sup>2</sup> > *t*1. If *ψ* (**x**1, *t*1)is expanded in terms of the eigen function *φ<sup>n</sup>* with the eigenvalue *En* of the constant Hamiltonian operator *H*ˆ as ∑*<sup>n</sup> cnφ<sup>n</sup>* (**x**), one

*φ<sup>n</sup>* (**x**2) *φ*<sup>∗</sup>

<sup>−</sup> *<sup>H</sup>*<sup>ˆ</sup> <sup>2</sup> 

where we abbreviated 1 for **x**1, *t*<sup>1</sup> and 2 for **x**2, *t*<sup>2</sup> and define *K* (2, 1) = 0 for *t*<sup>2</sup> < *t*1. It is

where *<sup>δ</sup>* (2, 1) <sup>=</sup> *<sup>δ</sup>* (*t*<sup>2</sup> <sup>−</sup> *<sup>t</sup>*1) *<sup>δ</sup>* (*x*<sup>2</sup> <sup>−</sup> *<sup>x</sup>*1) *<sup>δ</sup>* (*y*<sup>2</sup> <sup>−</sup> *<sup>y</sup>*1) *<sup>δ</sup>* (*z*<sup>2</sup> <sup>−</sup> *<sup>z</sup>*1) and the subscript 2 on *<sup>H</sup>*<sup>ˆ</sup> <sup>2</sup> means that the operator acts on the variables of 2 of *K* (2, 1). The kernel *K* is now called as the *Green's function* and it is the total amplitude for arrival at **x**2, *t*<sup>2</sup> starting from **x**1, *t*1. The transition

*<sup>n</sup>* (**x**1)*e*

<sup>−</sup>*<sup>i</sup> En <sup>h</sup>*¯ (*t*2−*t*1)

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> *<sup>H</sup>*<sup>ˆ</sup> *<sup>ψ</sup>*, (89)

. This description is equivalent to the description that the

*K* (**x**2, *t*2; **x**1, *t*1) *ψ* (**x**1, *t*1) d3**x**1, (90)

*K* (2, 1) = *ih*¯ *δ* (2, 1), (92)

*χ*<sup>∗</sup> (2) *K* (2, 1) *ψ* (1) d3**x**1d<sup>3</sup>**x**2. (93)

, (91)

*ih*¯ *∂ψ*

a part per billion (ppb) (Gabrielse et al., 2006; 2007; Odom et al., 2006).

*<sup>h</sup>*¯ Δ*t*

*ψ* (**x**2, *t*2) =

*<sup>K</sup>* (2, 1) <sup>=</sup> <sup>∑</sup>*<sup>n</sup>*

straightforward to show that *K* can be defined by that solution of

amplitude for finding a particle in state *χ* (2), if it was in *ψ* (1), is 

 *ih*¯ *<sup>∂</sup> ∂t*<sup>2</sup>

**4.1 The concepts of quantum field theory**

function *ψ*. The Schrödinger equation

the operation if an operator is *<sup>e</sup>*−*<sup>i</sup> <sup>H</sup>*<sup>ˆ</sup>

*t*<sup>1</sup> through the equation

can find for *t*<sup>2</sup> > *t*<sup>1</sup>

A quantum mechanical system is described equally well by specifying the function *K*, or by specifying the Hamiltonian operator *H*ˆ from which it results.

Let us consider a situation that a particle propagates from 1 to 2 through 3 in a weak potential operator *U*ˆ (3), which differs from zero only for *t* between *t*<sup>1</sup> and *t*2. The kernel is expanded in powers of *U*ˆ that

$$K(2,1) = K\_0(2,1) + K^{(1)}\left(2,1\right) + K^{(2)}\left(2,1\right) + \cdots \ . \tag{94}$$

To zeroth order in *U*ˆ , *K* is that for a free particle, *K*<sup>0</sup> (2, 1). Let us consider the situation if *U* differs from zero only for the infinitesimal time interval Δ*t*<sup>3</sup> between some time *t*<sup>3</sup> and *t*<sup>3</sup> + Δ*t*<sup>3</sup> for *t*<sup>1</sup> < *t*<sup>3</sup> < *t*2. The particle will propagate from 1 to 3 as a *free* particle,

$$
\psi\left(\mathbf{3}\right) = \int \mathbf{K}\_0\left(\mathbf{3}, \mathbf{1}\right) \psi\left(\mathbf{1}\right) \mathbf{d}^3 \mathbf{x}\_1. \tag{95}
$$

For the short time interval Δ*t*3, the wave function will change to

$$
\psi\left(\mathbf{x}, t\_3 + \Delta t\_3\right) = e^{-i\frac{H}{\hbar}\Delta t\_3} \psi\left(\mathbf{x}, t\_3\right),
\tag{96}
$$

after solving the Schrödinger equation in Eq. (89). The particle at 2 then propagates *freely* from **x**3, *t*<sup>3</sup> + Δ*t*<sup>3</sup> as

$$
\psi\left(\mathbf{x}\_{2},t\_{2}\right) = \int \mathbf{K}\_{0}\left(\mathbf{x}\_{2},t\_{2};\mathbf{x}\_{3},t\_{3} + \Delta t\_{3}\right) \psi\left(\mathbf{x}\_{3},t\_{3} + \Delta t\_{3}\right) \,. \tag{97}
$$

We can decompose the Hamiltonian operator *H*ˆ by *H*ˆ <sup>0</sup> + *U*ˆ , where *H*ˆ <sup>0</sup> is the Hamiltonian operator of the free particle. The change in wave function by *U*ˆ will be

$$
\Delta\psi = -\frac{i}{\hbar}\hat{\mathcal{U}}\left(\mathfrak{J}\right)\psi\left(\mathfrak{J}\right)\Delta t\_{\mathfrak{J}}.\tag{98}
$$

The wave function at 2 is that of the propagated particle from *t*<sup>3</sup> + Δ*t*<sup>3</sup> to be

$$
\psi\left(2\right) = \int K\_0\left(2\left,3'\right)\psi\left(3'\right)d^3\mathbf{x}\_{3'}\tag{99}
$$

where 3� abbreviates **x**3, *t*<sup>3</sup> + Δ*t*3, by a free propagation. The difference of the wave function at 2 is obtained by

$$
\Delta\Psi\left(\mathfrak{I}\right) = -\Delta t\_3 \frac{i}{\hbar} \int K\_0\left(\mathfrak{I}, \mathfrak{I}\right) \hat{\mathcal{U}}\left(\mathfrak{I}\right) K\_0\left(\mathfrak{I}, \mathfrak{I}\right) \psi\left(\mathfrak{I}\right) \mathrm{d}^3 \mathbf{x}\_3 \mathrm{d}^3 \mathbf{x}\_1. \tag{100}
$$

Therefore, the first order expansion of the kernel *K* is then

$$K^{(1)}\left(2,1\right) = -\frac{i}{\hbar} \int K\_0\left(2,3\right) \hat{\mathcal{U}}\left(3\right) K\_0\left(3,1\right) d\mathcal{S},\tag{101}$$

where d3 = d3**x**3d*t*3. We can imagine that a particle travels as a free particle from point to point, but is scattered by the potential operator *U*ˆ at 3. The higher order terms are also analyzed in a similar way.

The analysis for the charged free Dirac particle gives a new interpretation of the antiparticle, which has the reversed charge of the particle; for example, a positron is the antiparticle of an electron. The Dirac equation (Dirac, 1928a;b) has negative energy states of an electron. Dirac interpreted himself that the negative energy states are fully occupied in vacuum, and an

For a noninteracting particle, the eigen function will be a plane wave

**<sup>k</sup>** <sup>=</sup> *<sup>h</sup>*¯ *<sup>ω</sup>*<sup>0</sup>

*K*<sup>0</sup> (**k**, *ω*) =

the transition amplitude will be maximum around *ω* = *ω*<sup>0</sup>

the Fock (1930) consideration of the Pauli's exclusion principle.

the interaction from the noninteracting gas.

and the eigenvalue will be *�*<sup>0</sup>

which immediately yields

*K*0 **x**, *t*; **x**� , *t* � <sup>=</sup> <sup>1</sup> 2*π*<sup>4</sup> d3**k** ∞ −∞

*ω* = *ω*<sup>0</sup>

**4.2 Self-energy**

given in Eqs. (104) and (105) yields

*<sup>φ</sup>***<sup>k</sup>** (**x**) <sup>=</sup> <sup>1</sup>

d*ωei***k**·**<sup>x</sup>** *e* −*iω*(*t*−*t*� ) 

*θ* (*k* − *k*F) *<sup>ω</sup>* <sup>−</sup> *<sup>ω</sup>*<sup>0</sup>

**<sup>k</sup>** + *iη*

where *k*<sup>F</sup> is the Fermi wave vector. Since the ±*iη* were introduced to render the time integral in Eq. (104) convergent, it is convenient to take the limit *<sup>η</sup>* <sup>→</sup> <sup>0</sup>+. Eq. (106) diverges at

*free* particle remains in the state **k**: the particle will keep its momentum and hence its (kinetic) energy. This is nothing more than the celebrating statement of *inertial motion* by Galileo.

This idea that a particle propagates freely until it faced with the a scattering center, where the particles emit or absorb the interacting quanta, is nothing more than an extension of the model introduced by Drude (1900) for electrons in metals. We already obtained the thermodynamic information of noninteracting gases in Sec. 3. Hence, the remaining task is to see the effects of

Let us come back to the case of a particle propagating from 1 to 2 in the way given in Eq. (101) by considering the interaction process demonstrated in Eq. (103). The perturbation procedure for the interacting fermions includes, in its first-order expansion, two fundamental processes (Fetter & Walecka, 2003), which are the prototypes of the interactions of all order perturbation expansion. For the first case, the particle *a* propagates from 1 to 3, emits (absorbs) a boson propagating to 4, where the other particle *b* absorbs (emits) the boson. The particle *b* propagates after the absoption (emission) at 4 to the position 4 again, just before it absorbing (emitting) the boson. This process is known as the *vacuum polarization* and it is equivalent, for an electronic system, to the method of Hartree (1928). In terms of the classical electrodynamics, the process is nothing more than that an electron moves in a Coulomb potential generated by the neighboring charge density. Secondly, the particle propagates firstly from 1 to 3 and it emits (absorbs) a boson at 3 to change its state. The particle in the new state then propagates from 3 to 4, where the new state particle absorbs (emits) the boson propagated from 3, and change its states to the original one in propagating to 2. This process is known as the *exchange* and it is equivalent, for the electronic system, to

The energy spectrum of the particle itself will change during both the processes. Dyson (1949a;b) discussed that the changing the particle energy itself by the perturbative treatment of

+

**<sup>k</sup>** ± *iη*. If we compute the transition amplitude in Eq. (93) from the state *ψ* (1) to the state *ψ* (2) in the Fourier transformed space, the kernel function form of Eq. (106) implies that

<sup>√</sup>*<sup>V</sup>*

*n*, to be over **k**, in Eq. (91) becomes an integration and then the consideration of the identities

Towards the Authentic *Ab Intio* Thermodynamics 561

*ei***k**·**x**, (105)

+

*θ* (*k*<sup>F</sup> − *k*) *<sup>ω</sup>* <sup>−</sup> *<sup>ω</sup>*<sup>0</sup>

**<sup>k</sup>** − *iη*

, (106)

**<sup>k</sup>** <sup>±</sup> *<sup>i</sup>η*. In the limit *<sup>η</sup>* <sup>→</sup> <sup>0</sup>+, the

 ,

**<sup>k</sup>**. In the limit of an infinite volume, the summation over

*θ* (*k* − *k*F) *<sup>ω</sup>* <sup>−</sup> *<sup>ω</sup>*<sup>0</sup>

**<sup>k</sup>** − *iη*

*θ* (*k*<sup>F</sup> − *k*) *<sup>ω</sup>* <sup>−</sup> *<sup>ω</sup>*<sup>0</sup>

**<sup>k</sup>** + *iη*

elimination of one electron from the vacuum will carry a positive charge; the unoccupied state was interpreted as a *hole*. Feynman (1949b) reinterpreted that the *hole* is a positron, which is an electron propagting backward in time. The interpretation has the corresponding classical electrodynamic picture. If we record the trajectory of an electron moving in a magnetic field, the trajectory of the electron will be bent by the Lorentz force exerting on the electron. When we reverse the record in time of the electron in the magnetic field, the bending direction of the trajectory is that of the positively charged particle with the same mass to the electron. Therefore, we understand that a particle is propagting forward in time, while the corresponding antiparticle or the hole is propagating backward in time. Due to the negative energy nature of the hole or antiparticle, a particle-hole pair will be annihilated when the particle meet the hole at a position during their propagations in space-time coincidently. Reversely, vacuum can create the particle-hole pair from the vacuum fluctuations.

Now consider a system of two particles *a* and *b* propagate from 1 to 3 interacting at 5 for the particle *a* and from 2 to 4 interacting at 6 for the particle *b*. In the case of free particles, the kernel *K*<sup>0</sup> is a simple multiple of two free particle kernels *K*0*<sup>a</sup>* and *K*0*<sup>b</sup>* as

$$K\_0\left(3,4;1,2\right) = K\_{0a}\left(3,1\right)K\_{0b}\left(4,2\right). \tag{102}$$

When two particles are interacting through a two particle potential *U*ˆ , the first-order expansion term of *K* may be written (Feynman, 1949a) as

$$K^{(1)}\left(3,4;1,2\right) = -\frac{i}{\hbar} \iint K\_{0a}\left(3,5\right) K\_{0b}\left(4,6\right) \hat{\mathcal{U}}\left(5,6\right) K\_{0a}\left(5,1\right) K\_{0b}\left(6,2\right) d\mathbf{5d}\,\text{d}\,\text{f}\,\,\tag{103}$$

One important difference from Eq. (101) is that the interaction *U*ˆ at a specific space-time position 3 is replaced by the interaction field *U*ˆ (5, 6). The interaction field is also able to be quantized and the interaction field *U*ˆ (5, 6) is interpreted as an interaction field quantum propagating *freely* between 5 and 6. For the case of two-electron interaction, the particles are electrons interacting through the electromagnetic interaction. One electron *a* propagates from 1 to 3 and the other *b* propagates from 2 to 4. During their propagations, the electron *a* emits (aborbs) a photon at point 5, while the electron *b* absorbs (emits) the photon at 6. The wave function of each electron differs by the emission of the photon at 5 or 6 from its wave function at the origin of the propagation. The process of emission and absorption of the photons during the electrons propagations change the energy of the electronic subsystem. The process to compute the energy of the Fermi liquid in the perturbative treatment of the interaction requires the consideration of the essential many-body treatment available by the procedures suggested by Dyson (1949a;b); Feynman (1949a;b); Schwinger (1948; 1949a;b); Tomonaga (1946).

For the future reasons, it is useful to see the consequence of the step function behavior of the kernel *K*. As described above, *K* (2, 1) has its meaning as the solution of the Green's function Eq. (92) only if *t*<sup>2</sup> > *t*1. It is convenient to use multiply the step function *θ* (*t*<sup>2</sup> − *t*1) to the kernel *K* for implying the physical meaning. The step function has an integral representation

$$\theta\left(t - t'\right) = -\int\_{-\infty}^{\infty} \frac{\mathbf{d}\omega}{2\pi i} \frac{e^{-i\omega(t - t')}}{\omega + i\eta},\tag{104}$$

where *η* is an positive real number. If *t* > *t*� , then the contour must be closed in the lower-half *ω* plane, including the simple pole at *ω* = −*iη* with residue −1. If *t* < *t*� , then the contour must be closed in the upper-half *ω* plane and gives zero, because the integrand has no singularity for �*ω* > 0.

For a noninteracting particle, the eigen function will be a plane wave

$$\phi\_{\mathbf{k}}\left(\mathbf{x}\right) = \frac{1}{\sqrt{V}} e^{i\mathbf{k}\cdot\mathbf{x}},\tag{105}$$

and the eigenvalue will be *�*<sup>0</sup> **<sup>k</sup>** <sup>=</sup> *<sup>h</sup>*¯ *<sup>ω</sup>*<sup>0</sup> **<sup>k</sup>**. In the limit of an infinite volume, the summation over *n*, to be over **k**, in Eq. (91) becomes an integration and then the consideration of the identities given in Eqs. (104) and (105) yields

$$\mathcal{K}\_0\left(\mathbf{x},t;\mathbf{x}',t'\right) = \frac{1}{2\pi^4} \int \mathbf{d}^3 \mathbf{k} \int\_{-\infty}^{\infty} \mathbf{d}\omega e^{i\mathbf{k}\cdot\mathbf{x}} e^{-i\omega(t-t')} \left[ \frac{\theta\left(k-k\_\mathrm{F}\right)}{\omega - \omega\_\mathbf{k}^0 + i\eta} + \frac{\theta\left(k\_\mathrm{F}-k\right)}{\omega - \omega\_\mathbf{k}^0 - i\eta} \right] d\omega$$

which immediately yields

18 Will-be-set-by-IN-TECH

elimination of one electron from the vacuum will carry a positive charge; the unoccupied state was interpreted as a *hole*. Feynman (1949b) reinterpreted that the *hole* is a positron, which is an electron propagting backward in time. The interpretation has the corresponding classical electrodynamic picture. If we record the trajectory of an electron moving in a magnetic field, the trajectory of the electron will be bent by the Lorentz force exerting on the electron. When we reverse the record in time of the electron in the magnetic field, the bending direction of the trajectory is that of the positively charged particle with the same mass to the electron. Therefore, we understand that a particle is propagting forward in time, while the corresponding antiparticle or the hole is propagating backward in time. Due to the negative energy nature of the hole or antiparticle, a particle-hole pair will be annihilated when the particle meet the hole at a position during their propagations in space-time coincidently.

Reversely, vacuum can create the particle-hole pair from the vacuum fluctuations.

kernel *K*<sup>0</sup> is a simple multiple of two free particle kernels *K*0*<sup>a</sup>* and *K*0*<sup>b</sup>* as

expansion term of *K* may be written (Feynman, 1949a) as

*θ t* − *t* � = −

*ω* plane, including the simple pole at *ω* = −*iη* with residue −1. If *t* < *t*�

where *η* is an positive real number. If *t* > *t*�

*h*¯ 

*<sup>K</sup>*(1) (3, 4; 1, 2) <sup>=</sup> <sup>−</sup> *<sup>i</sup>*

Tomonaga (1946).

singularity for �*ω* > 0.

Now consider a system of two particles *a* and *b* propagate from 1 to 3 interacting at 5 for the particle *a* and from 2 to 4 interacting at 6 for the particle *b*. In the case of free particles, the

When two particles are interacting through a two particle potential *U*ˆ , the first-order

One important difference from Eq. (101) is that the interaction *U*ˆ at a specific space-time position 3 is replaced by the interaction field *U*ˆ (5, 6). The interaction field is also able to be quantized and the interaction field *U*ˆ (5, 6) is interpreted as an interaction field quantum propagating *freely* between 5 and 6. For the case of two-electron interaction, the particles are electrons interacting through the electromagnetic interaction. One electron *a* propagates from 1 to 3 and the other *b* propagates from 2 to 4. During their propagations, the electron *a* emits (aborbs) a photon at point 5, while the electron *b* absorbs (emits) the photon at 6. The wave function of each electron differs by the emission of the photon at 5 or 6 from its wave function at the origin of the propagation. The process of emission and absorption of the photons during the electrons propagations change the energy of the electronic subsystem. The process to compute the energy of the Fermi liquid in the perturbative treatment of the interaction requires the consideration of the essential many-body treatment available by the procedures suggested by Dyson (1949a;b); Feynman (1949a;b); Schwinger (1948; 1949a;b);

For the future reasons, it is useful to see the consequence of the step function behavior of the kernel *K*. As described above, *K* (2, 1) has its meaning as the solution of the Green's function Eq. (92) only if *t*<sup>2</sup> > *t*1. It is convenient to use multiply the step function *θ* (*t*<sup>2</sup> − *t*1) to the kernel *K* for implying the physical meaning. The step function has an integral representation

> ∞ −∞

must be closed in the upper-half *ω* plane and gives zero, because the integrand has no

d*ω* 2*πi* *e*−*iω*(*t*−*<sup>t</sup>* � )

*<sup>ω</sup>* <sup>+</sup> *<sup>i</sup><sup>η</sup>* , (104)

, then the contour

, then the contour must be closed in the lower-half

*K*<sup>0</sup> (3, 4; 1, 2) = *K*0*<sup>a</sup>* (3, 1) *K*0*<sup>b</sup>* (4, 2). (102)

*K*0*<sup>a</sup>* (3, 5) *K*0*<sup>b</sup>* (4, 6) *U*ˆ (5, 6) *K*0*<sup>a</sup>* (5, 1) *K*0*<sup>b</sup>* (6, 2) d5d6. (103)

$$K\_0\left(\mathbf{k},\omega\right) = \left[\frac{\theta\left(k - k\_{\rm F}\right)}{\omega - \omega\_{\rm k}^0 + i\eta} + \frac{\theta\left(k\_{\rm F} - k\right)}{\omega - \omega\_{\rm k}^0 - i\eta}\right],\tag{106}$$

where *k*<sup>F</sup> is the Fermi wave vector. Since the ±*iη* were introduced to render the time integral in Eq. (104) convergent, it is convenient to take the limit *<sup>η</sup>* <sup>→</sup> <sup>0</sup>+. Eq. (106) diverges at *ω* = *ω*<sup>0</sup> **<sup>k</sup>** ± *iη*. If we compute the transition amplitude in Eq. (93) from the state *ψ* (1) to the state *ψ* (2) in the Fourier transformed space, the kernel function form of Eq. (106) implies that the transition amplitude will be maximum around *ω* = *ω*<sup>0</sup> **<sup>k</sup>** <sup>±</sup> *<sup>i</sup>η*. In the limit *<sup>η</sup>* <sup>→</sup> <sup>0</sup>+, the *free* particle remains in the state **k**: the particle will keep its momentum and hence its (kinetic) energy. This is nothing more than the celebrating statement of *inertial motion* by Galileo.

#### **4.2 Self-energy**

This idea that a particle propagates freely until it faced with the a scattering center, where the particles emit or absorb the interacting quanta, is nothing more than an extension of the model introduced by Drude (1900) for electrons in metals. We already obtained the thermodynamic information of noninteracting gases in Sec. 3. Hence, the remaining task is to see the effects of the interaction from the noninteracting gas.

Let us come back to the case of a particle propagating from 1 to 2 in the way given in Eq. (101) by considering the interaction process demonstrated in Eq. (103). The perturbation procedure for the interacting fermions includes, in its first-order expansion, two fundamental processes (Fetter & Walecka, 2003), which are the prototypes of the interactions of all order perturbation expansion. For the first case, the particle *a* propagates from 1 to 3, emits (absorbs) a boson propagating to 4, where the other particle *b* absorbs (emits) the boson. The particle *b* propagates after the absoption (emission) at 4 to the position 4 again, just before it absorbing (emitting) the boson. This process is known as the *vacuum polarization* and it is equivalent, for an electronic system, to the method of Hartree (1928). In terms of the classical electrodynamics, the process is nothing more than that an electron moves in a Coulomb potential generated by the neighboring charge density. Secondly, the particle propagates firstly from 1 to 3 and it emits (absorbs) a boson at 3 to change its state. The particle in the new state then propagates from 3 to 4, where the new state particle absorbs (emits) the boson propagated from 3, and change its states to the original one in propagating to 2. This process is known as the *exchange* and it is equivalent, for the electronic system, to the Fock (1930) consideration of the Pauli's exclusion principle.

The energy spectrum of the particle itself will change during both the processes. Dyson (1949a;b) discussed that the changing the particle energy itself by the perturbative treatment of

A similar analysis can be carried out for the interaction between two particles, which always consists of the lowest-order interaction plus a series of proper expansion. The four

Towards the Authentic *Ab Intio* Thermodynamics 563

1 − Π (*q*) *U*<sup>0</sup> (*q*)

. (114)


*κ* (*q*) = 1 − *U*<sup>0</sup> (*q*) Π (*q*), (115)

*<sup>κ</sup>* (*q*) . (116)

*<sup>U</sup>* (*q*) <sup>=</sup> *<sup>U</sup>*<sup>0</sup> (*q*)

the *screening* of the lowest-order interaction by the polarization of the medium is obtained as

*<sup>U</sup>* (*q*) <sup>=</sup> *<sup>U</sup>*<sup>0</sup> (*q*)

Goldstone (1957) provided a new picture of the many-body systems with the quantum field theoretic point of view, presented above. Let us the free particle Hamiltonian *H*ˆ has a many-body eigenstate Φ, which is a determinant formed from *N* particles of the *ψn*, and which is able to be described by enumerating these *N* one-particle states. Suppose that *H*ˆ <sup>0</sup> has a non-degenerate ground state Φ<sup>0</sup> formed from the lowest *N* of the *ψn*. The states *ψ<sup>n</sup>* occupied in Φ<sup>0</sup> will be called unexcited states, and all the higher states *ψ<sup>n</sup>* will be called excited states. An eigenstate Φ of *H*ˆ <sup>0</sup> can be described by enumerating all the excited states which are occupied, and all the unexcited states which are not occupied. An unoccupied unexcited state is regarded as a *hole*, and the theory will deal with particles in excited states and holes in the unexcited states. In this treatment, the ground state Φ<sup>0</sup> is considered as a new *vacuum*, a particle is considered as an occupied states in the excited states, and the hole is essentially

different from the positrons, in which are the symmetric counterpart of the electrons.

*<sup>E</sup>* <sup>−</sup> *<sup>E</sup>*<sup>0</sup> <sup>=</sup> �Φ0<sup>|</sup> *<sup>H</sup>*<sup>ˆ</sup> <sup>1</sup> <sup>∑</sup>*<sup>n</sup>*

−<sup>1</sup>

*<sup>E</sup>*<sup>0</sup> <sup>−</sup> *<sup>H</sup>*<sup>ˆ</sup> <sup>0</sup>

Goldstone (1957) derived the energy difference between the system with and without

where the summation should do on the linked<sup>7</sup> terms of the perturbation. The noninteracting Hamiltonian *H*ˆ <sup>0</sup> in the denominator can be replaced by the corresponding eigenvalues, because Eq. (117) is interpreted by inserting a complete set of eigenstates of *H*ˆ <sup>0</sup> between each interaction *H*ˆ 1. The physical situation can be visualized as follows: (1) The interaction Hamiltonian *<sup>H</sup>*<sup>ˆ</sup> <sup>1</sup> operate on <sup>|</sup>Φ0� creates a state with two particles and two holes. This

or scatter the existing particle-hole pair and so on. (3) The final *H*ˆ <sup>1</sup> must then return the system to the ground state |Φ0�. This process gives the difference in energy of the interacting many-body system from the noninteracting one. By choosing the first-order perturbation in

<sup>6</sup> In the Dirac notation, a quantum state *<sup>n</sup>* is written in the Hilbert space of form <sup>|</sup>*n*� and the corresponding conjugate state is written as �*n*|. The wave function is the projection to the position space, such that *<sup>ψ</sup><sup>n</sup>* (**x**) <sup>=</sup> �**x**<sup>|</sup> *<sup>n</sup>*�. For an operator *<sup>A</sup>*ˆ, �*m*<sup>|</sup> *<sup>A</sup>*<sup>ˆ</sup> <sup>|</sup>*n*� is called as a *matrix element* to represent the probability for transtion from the state <sup>|</sup>*n*� by the operation *<sup>A</sup>*<sup>ˆ</sup> to the state �*m*|. Readers can see the details in Dirac

<sup>7</sup> During the expansion, there are terms describing pair creation and annihilations corresponding to the

 1 *<sup>E</sup>*<sup>0</sup> <sup>−</sup> *<sup>H</sup>*<sup>ˆ</sup> <sup>0</sup> *H*ˆ 1 *<sup>n</sup>*

. (2) The next *H*ˆ <sup>1</sup> can create more particle-hole pairs

dimensional Fourier transformation to the *q* coordinates yields

Introducing a generalized dielectric function

**4.3 Goldstone's theorem: the many-body formalism**

interactions, *H*ˆ 1, in the Dirac notation,<sup>6</sup> as

state propagates with

(1998); Sakurai (1994).

free particle Green's function.

the fermion interaction when we consider single particle propagation from 1 to 2. The energy of the particle differs from the noninteracting particle propagation, and it can be systematically included to the particle propagation kernel as

$$K(2,1) = K\_0(2,1) + \iint K(2,4) \, \Sigma \, (4,3) \, K(3,1) \, \text{d}3 \, \text{d}4,\tag{107}$$

where Σ is known as the *self-energy*. When we consider the all order perturbation, the exact single-particle propagation can be obtained by using the successive self-energy inclusion as

$$\begin{aligned} K(2,1) &= K\_0(2,1) + \iint K(2,4) \, \Sigma\_{\ }(4,3) \, K(3,1) \, \text{d}3 \, \text{d}4 \\ &+ \iiint K(2,6) \, \Sigma\_{\ }(6,5) \, K(5,4) \, \Sigma\_{\ }(4,3) \, K(3,1) \, \text{d}3 \, \text{d}4 \, \text{d}5 \, \text{d}6 + \cdots \, \end{aligned} \tag{108}$$

with the special care of the self-energy to be proper.<sup>5</sup> This equation is known as the *Dyson equation*.

The consequence of the Dyson equation can be seen easily if we perform a four dimensional Fourier transform on Eq. (108) with respect to the difference **y** − **x**, *t*<sup>2</sup> − *t*<sup>1</sup> into the momentum space to an algebraic form

$$K\left(k\right) = \mathcal{K}\_0\left(k\right) + \mathcal{K}\_0\left(k\right)\Sigma\left(k\right)\mathcal{K}\left(k\right),\tag{109}$$

where *k* abbreviates (**k**, *ω*). We can solve Eq. (109) as

$$K(k) = \frac{1}{K\_0^{-1}\left(k\right) - \Sigma\left(k\right)}.\tag{110}$$

Considering the self-energy Σ is complex, the *iη* in the free particle kernel Eq. (106) is no more relevant. Therefore, we obtain the solution of the Dyson equation as

$$K(k) = K\left(\mathbf{k}, \omega\right) = \frac{1}{\omega - \hbar^{-1}\epsilon\_\mathbf{k}^0 - \Sigma\left(\mathbf{k}, \omega\right)}.\tag{111}$$

The physical meaning of Eq. (111) is straightforward: an interacting particle propagates as the free particle does, but its excitation energy differs by a *dressing* term Σ (**k**, *ω*).

Lehmann (1954) and Galitskii & Migdal (1958) discussed the usefulness of Eq. (111) in the applications for many-body systems. In the Lehmann representation, the frequency *ω* is a complex number to be

$$
\hbar\hbar\omega = \varepsilon\_\mathbf{k} - i\gamma\_\mathbf{k'}\tag{112}
$$

where *γ***<sup>k</sup>** is the damping of the particle. The singularity of the exact Green's function *K* (**k**, *ω*), considered as a function of *ω*, determine both the excitation energy *�***<sup>k</sup>** of the system and its damping *γ***k**. Furthermore, the chemical potential can be determined as the point where �Σ (**k**, *ω*) changes the sign, because

$$\begin{array}{l}\text{G2}\left(\mathbf{k},\omega\right)\geq 0,\ \omega<\mu/\hbar\\\text{G2}\left(\mathbf{k},\omega\right)\leq 0,\ \omega>\mu/\hbar.\end{array}\tag{113}$$

<sup>5</sup> The proper implies the terms that cannot be disintegrated into the lower order expansion terms during the perturbation expansion.

20 Will-be-set-by-IN-TECH

the fermion interaction when we consider single particle propagation from 1 to 2. The energy of the particle differs from the noninteracting particle propagation, and it can be systematically

where Σ is known as the *self-energy*. When we consider the all order perturbation, the exact single-particle propagation can be obtained by using the successive self-energy inclusion as

with the special care of the self-energy to be proper.<sup>5</sup> This equation is known as the *Dyson*

The consequence of the Dyson equation can be seen easily if we perform a four dimensional Fourier transform on Eq. (108) with respect to the difference **y** − **x**, *t*<sup>2</sup> − *t*<sup>1</sup> into the momentum

> *<sup>K</sup>* (*k*) <sup>=</sup> <sup>1</sup> *K*−<sup>1</sup>

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

free particle does, but its excitation energy differs by a *dressing* term Σ (**k**, *ω*).

<sup>0</sup> (*k*) − Σ (*k*)

*<sup>ω</sup>* <sup>−</sup> *<sup>h</sup>*¯ <sup>−</sup>1*�*<sup>0</sup>

Considering the self-energy Σ is complex, the *iη* in the free particle kernel Eq. (106) is no more

The physical meaning of Eq. (111) is straightforward: an interacting particle propagates as the

Lehmann (1954) and Galitskii & Migdal (1958) discussed the usefulness of Eq. (111) in the applications for many-body systems. In the Lehmann representation, the frequency *ω* is a

where *γ***<sup>k</sup>** is the damping of the particle. The singularity of the exact Green's function *K* (**k**, *ω*), considered as a function of *ω*, determine both the excitation energy *�***<sup>k</sup>** of the system and its damping *γ***k**. Furthermore, the chemical potential can be determined as the point where

�Σ (**k**, *ω*) ≥ 0, *ω* < *μ*/¯*h*

<sup>5</sup> The proper implies the terms that cannot be disintegrated into the lower order expansion terms during

*K* (2, 4) Σ (4, 3) *K* (3, 1) d3d4

*K* (2, 4) Σ (4, 3) *K* (3, 1) d3d4, (107)

*K* (2, 6) Σ (6, 5) *K* (5, 4) Σ (4, 3) *K* (3, 1) d3d4d5d6 + ··· , (108)

*K* (*k*) = *K*<sup>0</sup> (*k*) + *K*<sup>0</sup> (*k*) Σ (*k*) *K* (*k*), (109)

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

*h*¯ *ω* = *�***<sup>k</sup>** − *iγ***k**, (112)

�<sup>Σ</sup> (**k**, *<sup>ω</sup>*) <sup>≤</sup> 0, *<sup>ω</sup>* <sup>&</sup>gt; *<sup>μ</sup>*/¯*h*. (113)

. (110)

. (111)

included to the particle propagation kernel as

*K* (2, 1) = *K*<sup>0</sup> (2, 1) +

+ 

where *k* abbreviates (**k**, *ω*). We can solve Eq. (109) as

relevant. Therefore, we obtain the solution of the Dyson equation as

*equation*.

space to an algebraic form

complex number to be

�Σ (**k**, *ω*) changes the sign, because

the perturbation expansion.

*K* (2, 1) = *K*<sup>0</sup> (2, 1) +

A similar analysis can be carried out for the interaction between two particles, which always consists of the lowest-order interaction plus a series of proper expansion. The four dimensional Fourier transformation to the *q* coordinates yields

$$\mathcal{U}\left(q\right) = \frac{\mathcal{U}\_0\left(q\right)}{1 - \Pi\left(q\right)\mathcal{U}\_0\left(q\right)}.\tag{114}$$

Introducing a generalized dielectric function

$$\kappa \left( \boldsymbol{q} \right) = 1 - \mathcal{U}\_0 \left( \boldsymbol{q} \right) \boldsymbol{\Pi} \left( \boldsymbol{q} \right) \, \tag{115}$$

the *screening* of the lowest-order interaction by the polarization of the medium is obtained as

$$\mathcal{U}\left(q\right) = \frac{\mathcal{U}\_0\left(q\right)}{\kappa\left(q\right)}.\tag{116}$$

#### **4.3 Goldstone's theorem: the many-body formalism**

Goldstone (1957) provided a new picture of the many-body systems with the quantum field theoretic point of view, presented above. Let us the free particle Hamiltonian *H*ˆ has a many-body eigenstate Φ, which is a determinant formed from *N* particles of the *ψn*, and which is able to be described by enumerating these *N* one-particle states. Suppose that *H*ˆ <sup>0</sup> has a non-degenerate ground state Φ<sup>0</sup> formed from the lowest *N* of the *ψn*. The states *ψ<sup>n</sup>* occupied in Φ<sup>0</sup> will be called unexcited states, and all the higher states *ψ<sup>n</sup>* will be called excited states. An eigenstate Φ of *H*ˆ <sup>0</sup> can be described by enumerating all the excited states which are occupied, and all the unexcited states which are not occupied. An unoccupied unexcited state is regarded as a *hole*, and the theory will deal with particles in excited states and holes in the unexcited states. In this treatment, the ground state Φ<sup>0</sup> is considered as a new *vacuum*, a particle is considered as an occupied states in the excited states, and the hole is essentially different from the positrons, in which are the symmetric counterpart of the electrons.

Goldstone (1957) derived the energy difference between the system with and without interactions, *H*ˆ 1, in the Dirac notation,<sup>6</sup> as

$$E - E\_0 = \left\langle \Phi\_0 \right| \hat{H}\_1 \sum\_n \left( \frac{1}{E\_0 - \hat{H}\_0} \hat{H}\_1 \right)^n \left| \Phi\_0 \right\rangle \tag{117}$$

where the summation should do on the linked<sup>7</sup> terms of the perturbation. The noninteracting Hamiltonian *H*ˆ <sup>0</sup> in the denominator can be replaced by the corresponding eigenvalues, because Eq. (117) is interpreted by inserting a complete set of eigenstates of *H*ˆ <sup>0</sup> between each interaction *H*ˆ 1. The physical situation can be visualized as follows: (1) The interaction Hamiltonian *<sup>H</sup>*<sup>ˆ</sup> <sup>1</sup> operate on <sup>|</sup>Φ0� creates a state with two particles and two holes. This state propagates with *<sup>E</sup>*<sup>0</sup> <sup>−</sup> *<sup>H</sup>*<sup>ˆ</sup> <sup>0</sup> −<sup>1</sup> . (2) The next *H*ˆ <sup>1</sup> can create more particle-hole pairs or scatter the existing particle-hole pair and so on. (3) The final *H*ˆ <sup>1</sup> must then return the system to the ground state |Φ0�. This process gives the difference in energy of the interacting many-body system from the noninteracting one. By choosing the first-order perturbation in

<sup>6</sup> In the Dirac notation, a quantum state *<sup>n</sup>* is written in the Hilbert space of form <sup>|</sup>*n*� and the corresponding conjugate state is written as �*n*|. The wave function is the projection to the position space, such that *<sup>ψ</sup><sup>n</sup>* (**x**) <sup>=</sup> �**x**<sup>|</sup> *<sup>n</sup>*�. For an operator *<sup>A</sup>*ˆ, �*m*<sup>|</sup> *<sup>A</sup>*<sup>ˆ</sup> <sup>|</sup>*n*� is called as a *matrix element* to represent the probability for transtion from the state <sup>|</sup>*n*� by the operation *<sup>A</sup>*<sup>ˆ</sup> to the state �*m*|. Readers can see the details in Dirac (1998); Sakurai (1994).

<sup>7</sup> During the expansion, there are terms describing pair creation and annihilations corresponding to the free particle Green's function.

For both fermionic and bosonic *statistics*, G is *periodic* over the range 2*βh*¯ and may expanded

Towards the Authentic *Ab Intio* Thermodynamics 565

*<sup>i</sup>ωnτ*<sup>G</sup> (**x**1, **<sup>x</sup>**2, *<sup>ω</sup>n*), (125)

*<sup>β</sup>h*¯ . (126)

<sup>d</sup>*τeiωnτ*<sup>G</sup> (**x**1, **<sup>x</sup>**2, *<sup>τ</sup>*), (127)

*<sup>β</sup>h*¯ , fermion . (128)

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

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

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

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

, (129)

. (130)

, (131)

, (132)

reduces to

2 

*<sup>n</sup>*, and the factor <sup>1</sup>

*<sup>β</sup>h*¯ <sup>∑</sup>*<sup>n</sup> e*

> *βh*¯ 0

(2*n*+1)*π*

solution of the Dyson equation in Eq. (110) for the temperature Green's function as

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

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

 1 0

d*λ λ*

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

where *s* is the spin of the particles, *η* is the infinitesimally small positive number, and *λ* is the

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

> *�*0 **<sup>k</sup>** <sup>+</sup> <sup>1</sup>

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

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

<sup>8</sup> When an interacting Hamiltonian *H*ˆ *<sup>I</sup>* is applied to the noninteracting system with the Hamiltonian *H*ˆ 0, the interaction strength is controlled by the parameter *λ* to yield Hamiltonian *H*ˆ = *H*ˆ <sup>0</sup> + *<sup>λ</sup>H*<sup>ˆ</sup> *<sup>I</sup>* with its *<sup>λ</sup>* parameterized eigenstate <sup>|</sup>Ψ<sup>0</sup> (*λ*)�. The energy difference *<sup>E</sup>* <sup>−</sup> *<sup>E</sup>*<sup>0</sup> is calculated by

The thermodynamic variables *N*, *E*, and Ω can by obtained by the Eq. (129):

∓*V* (2*s* + 1)

× 1 2 *h*¯ +

× 1 2 *h*¯ +

interaction strength control parameter.<sup>8</sup> The internal energy is

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

*<sup>β</sup>h*¯ , boson

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

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

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

1

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

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

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

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

The coordinates in **x** can be transformed to the momentum space **k**. There is a corresponding

2*nπ*

*<sup>ω</sup><sup>n</sup>* <sup>=</sup> *<sup>n</sup><sup>π</sup>*

<sup>G</sup> (**x**1, **<sup>x</sup>**2, *<sup>τ</sup>*) <sup>=</sup> <sup>1</sup>

. The Fourier coefficient in Eq. (125) is then

G (**x**1, **x**2, *ωn*) =

*ω<sup>n</sup>* =

with the noninteracting temperature Green's function of the form

Ω (*T*, *V*, *μ*) = Ω<sup>0</sup> (*T*, *V*, *μ*)

This shows that *<sup>e</sup>*−*iωnβh*¯ is equal to *<sup>e</sup>*−*in<sup>π</sup>* <sup>=</sup> (−1)

This is known as the *Matsubara frequency*.

in a Fourier seriese

1 ± (−1)

*n*

The grand potential is

 1 0 d*λ*

*<sup>λ</sup>* �Ψ<sup>0</sup> (*λ*)<sup>|</sup> *<sup>λ</sup>H*<sup>ˆ</sup> *<sup>I</sup>* <sup>|</sup>Ψ<sup>0</sup> (*λ*)�.

where

where

1 2 

the Goldstone scheme, it is successfully reproduced the Hartree-Fock equation. The method for the perturbation expansion in the Goldstone scheme is usually called as coupled cluster expansion.

#### **4.4 Matsubara method for finite temperatures**

For simplicity, consider a system with one kind of particles in dealing with the grand partition function defined in Eq. (2), *i.e.*,

$$\Xi = \sum\_{i} e^{-\beta(\epsilon\_{i} - \mu N)}.\tag{118}$$

It is convenient to introduce a *grand canonical Hamiltonian* operator

$$
\hat{\mathcal{K}} = \hat{H} - \mu \hat{\mathcal{N}}.\tag{119}
$$

The operator form forces the equivalent form of Eq. (118) to be a trace

$$
\Xi = e^{-\beta \Omega} = \text{Tr}e^{-\beta \hat{\mathcal{K}}},\tag{120}
$$

and its corresponding statistical density matrix in the form

$$
\mathfrak{H} = \frac{e^{-\beta \mathcal{K}}}{\Xi} = e^{\beta \left(\Omega - \mathcal{K}\right)}.\tag{121}
$$

Armed with Eqs. (120) and (121), all the thermodynamic information are calculable. However, difficulty is in its thermally evolving nature of those equations, instead of their time evolving nature of the Schrödinger equation.

In order to compute thermodynamic properties of ferromagnetic materials, Bloch (1932) suggested a statistical equation, resembling the ordinary Schrödinger equation, by treating the evolution of the statistical density with respect to temperature as

$$-\frac{\partial \rho}{\partial \beta} = \hat{H}\rho.\tag{122}$$

Matsubara (1955) and Gaudin (1960) proved that the Green's function formalism in quantum field theory also satisfies the Wick's theorem (Wick, 1950), hence the perturbation expansion in all order satisfies, when the "single-particle" field operator follows the Bloch equation by switching the real-time *t* to the *imaginary-time*, *τ* = *it*. Note that the term "single-particle" does not mean a real single-particle.

The procedure in Sec. 4.1 is applied also. The single-particle temperature Green's function G may follow the Green's function equation

$$
\left[\hbar\frac{\partial}{\partial\tau\_2} - \hat{\mathcal{K}}\right] \mathcal{G}\left(\mathbf{x}\_1\tau\_1; \mathbf{x}\_2\tau\_2\right) = \pm\hbar\delta(\mathbf{x}\_2 - \mathbf{x}\_1)\delta\left(\tau\_2 - \tau\_1\right),
\tag{123}
$$

where ± indicates the temperature Green's function for fermions and bosons, respectively. However, the temperature Green's function does not depend on *τ*<sup>1</sup> and *τ*<sup>2</sup> separately, but it depends on the difference *τ* ≡ *τ*<sup>1</sup> − *τ*2. Each imaginary time *τ* has *βh*¯ periodicity. The temperature Green's function is then in the form

$$\mathcal{G}\left(\mathbf{x}\_{1},\tau\_{1};\mathbf{x}\_{2},\tau\_{2}\right) = \mathcal{G}\left(\mathbf{x}\_{1},\mathbf{x}\_{2};\tau\right). \tag{124}$$

For both fermionic and bosonic *statistics*, G is *periodic* over the range 2*βh*¯ and may expanded in a Fourier seriese

$$\mathcal{G}\left(\mathbf{x}\_{1},\mathbf{x}\_{2},\tau\right) = \frac{1}{\beta\hbar} \sum\_{\mathbf{n}} \epsilon^{i\omega\_{\mathbf{n}}\tau} \mathcal{G}\left(\mathbf{x}\_{1},\mathbf{x}\_{2},\omega\_{\mathbf{n}}\right) \,,\tag{125}$$

where

22 Will-be-set-by-IN-TECH

the Goldstone scheme, it is successfully reproduced the Hartree-Fock equation. The method for the perturbation expansion in the Goldstone scheme is usually called as coupled cluster

For simplicity, consider a system with one kind of particles in dealing with the grand partition

<sup>−</sup>*β*<sup>Ω</sup> = Tr*e*

<sup>Ξ</sup> <sup>=</sup> *<sup>e</sup>*

Armed with Eqs. (120) and (121), all the thermodynamic information are calculable. However, difficulty is in its thermally evolving nature of those equations, instead of their time evolving

In order to compute thermodynamic properties of ferromagnetic materials, Bloch (1932) suggested a statistical equation, resembling the ordinary Schrödinger equation, by treating

Matsubara (1955) and Gaudin (1960) proved that the Green's function formalism in quantum field theory also satisfies the Wick's theorem (Wick, 1950), hence the perturbation expansion in all order satisfies, when the "single-particle" field operator follows the Bloch equation by switching the real-time *t* to the *imaginary-time*, *τ* = *it*. Note that the term "single-particle"

The procedure in Sec. 4.1 is applied also. The single-particle temperature Green's function G

where ± indicates the temperature Green's function for fermions and bosons, respectively. However, the temperature Green's function does not depend on *τ*<sup>1</sup> and *τ*<sup>2</sup> separately, but it depends on the difference *τ* ≡ *τ*<sup>1</sup> − *τ*2. Each imaginary time *τ* has *βh*¯ periodicity. The

<sup>−</sup> *∂ρ*

<sup>−</sup>*β*K<sup>ˆ</sup>

−*β*(*�i*−*μN*)

. (118)

, (120)

*<sup>β</sup>*(Ω−Kˆ). (121)

*∂β* <sup>=</sup> *<sup>H</sup>*<sup>ˆ</sup> *<sup>ρ</sup>*. (122)

G (**x**1*τ*1; **x**2*τ*2) = ±*h*¯ *δ*(**x**<sup>2</sup> − **x**1)*δ* (*τ*<sup>2</sup> − *τ*1), (123)

G (**x**1, *τ*1; **x**2, *τ*2) = G (**x**1, **x**2; *τ*). (124)

<sup>K</sup><sup>ˆ</sup> <sup>=</sup> *<sup>H</sup>*<sup>ˆ</sup> <sup>−</sup> *<sup>μ</sup>N*<sup>ˆ</sup> . (119)

Ξ = ∑ *i e*

It is convenient to introduce a *grand canonical Hamiltonian* operator

The operator form forces the equivalent form of Eq. (118) to be a trace

the evolution of the statistical density with respect to temperature as

and its corresponding statistical density matrix in the form

Ξ = *e*

*<sup>ρ</sup>*<sup>ˆ</sup> <sup>=</sup> *<sup>e</sup>*−*β*K<sup>ˆ</sup>

expansion.

**4.4 Matsubara method for finite temperatures**

function defined in Eq. (2), *i.e.*,

nature of the Schrödinger equation.

does not mean a real single-particle.

may follow the Green's function equation *h*¯ *∂ ∂τ*2

<sup>−</sup> <sup>K</sup><sup>ˆ</sup> 

temperature Green's function is then in the form

$$
\omega\_{\mathfrak{n}} = \frac{n\pi}{\beta\hbar}.\tag{126}
$$

This shows that *<sup>e</sup>*−*iωnβh*¯ is equal to *<sup>e</sup>*−*in<sup>π</sup>* <sup>=</sup> (−1) *<sup>n</sup>*, and the factor <sup>1</sup> 2 <sup>1</sup> <sup>±</sup> *<sup>e</sup>*−*iωnβh*¯ reduces to 1 2 1 ± (−1) *n* . The Fourier coefficient in Eq. (125) is then

$$\mathcal{G}\left(\mathbf{x}\_{1},\mathbf{x}\_{2},\omega\_{\hbar}\right) = \int\_{0}^{\beta\hbar} \mathbf{d}\tau e^{i\omega\_{\hbar}\tau} \mathcal{G}\left(\mathbf{x}\_{1},\mathbf{x}\_{2},\tau\right),\tag{127}$$

where

$$\omega\_{\boldsymbol{\eta}} = \begin{cases} \frac{2n\pi}{\beta\hbar}, & \text{boson} \\ \frac{(2n+1)\pi}{\beta\hbar}, & \text{fermion} \end{cases} \tag{128}$$

This is known as the *Matsubara frequency*.

The coordinates in **x** can be transformed to the momentum space **k**. There is a corresponding solution of the Dyson equation in Eq. (110) for the temperature Green's function as

$$\mathcal{G}\left(\mathbf{k},\omega\_{\rm n}\right) = \frac{1}{\mathcal{G}\_0^{-1} - \Sigma\left(\mathbf{k},\omega\_{\rm n}\right)},\tag{129}$$

with the noninteracting temperature Green's function of the form

$$\mathcal{G}\_0\left(\mathbf{k}, \omega\_{\hbar}\right) = \frac{1}{i\omega\_{\hbar} - \hbar^{-1}\left(\epsilon\_{\mathbf{k}}^0 - \mu\right)}.\tag{130}$$

The thermodynamic variables *N*, *E*, and Ω can by obtained by the Eq. (129): The grand potential is

$$\begin{split} \Omega\left(T, V, \mu\right) &= \Omega\_{0}\left(T, V, \mu\right) \\ &\quad \mp V\left(2s+1\right) \int\_{0}^{1} \frac{\mathbf{d}\lambda}{\lambda} \int \frac{\mathbf{d}^{3}k}{\left(2\pi\right)^{2}} \frac{1}{\beta\hbar} \sum\_{n} e^{i\omega\_{n}\eta} \\ &\quad \times \left[\frac{1}{2}\hbar + \frac{\frac{1}{2}\Sigma^{\lambda}\left(\mathbf{k}, \omega\_{n}\right)}{i\omega\_{n} - \hbar^{-1}\left(\epsilon\_{\mathbf{k}}^{0} - \mu\right) - \Sigma^{\lambda}\left(\mathbf{k}, \omega\_{n}\right)}\right], \end{split} \tag{131}$$

where *s* is the spin of the particles, *η* is the infinitesimally small positive number, and *λ* is the interaction strength control parameter.<sup>8</sup> The internal energy is

$$\begin{split} E\left(T, V, \mu\right) &= \mp V \left(2s + 1\right) \int \frac{\mathbf{d}^3 k}{\left(2\pi\right)^2} \frac{1}{\beta\hbar} \sum\_{n} e^{i\omega\_n \eta} \\ &\times \left[\frac{1}{2}\hbar + \frac{\epsilon\_\mathbf{k}^0 + \frac{1}{2}\Sigma^\lambda \left(\mathbf{k}, \omega\_n\right)}{i\omega\_n - \hbar^{-1}\left(\epsilon\_\mathbf{k}^0 - \mu\right) - \Sigma^\lambda \left(\mathbf{k}, \omega\_n\right)}\right], \end{split} \tag{132}$$

<sup>8</sup> When an interacting Hamiltonian *H*ˆ *<sup>I</sup>* is applied to the noninteracting system with the Hamiltonian *H*ˆ 0, the interaction strength is controlled by the parameter *λ* to yield Hamiltonian *H*ˆ = *H*ˆ <sup>0</sup> + *<sup>λ</sup>H*<sup>ˆ</sup> *<sup>I</sup>* with its *<sup>λ</sup>* parameterized eigenstate <sup>|</sup>Ψ<sup>0</sup> (*λ*)�. The energy difference *<sup>E</sup>* <sup>−</sup> *<sup>E</sup>*<sup>0</sup> is calculated by 1 0 d*λ <sup>λ</sup>* �Ψ<sup>0</sup> (*λ*)<sup>|</sup> *<sup>λ</sup>H*<sup>ˆ</sup> *<sup>I</sup>* <sup>|</sup>Ψ<sup>0</sup> (*λ*)�.

discrete set of points {*iωn*}, the thermodynamic variables Ω, *E*, and *N* are obtained. The spectral weight function is able to be derived from the temperature Green's function G. Then we can have the real-time Green's functions *G*¯ *<sup>R</sup>* and *G*¯ *<sup>A</sup>*, which give us the quasiparticle

Towards the Authentic *Ab Intio* Thermodynamics 567

Hence, we complete the equilibrium thermodynamics; information on macroscopic

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

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

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

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

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'

of magnetism, which we have not understood its underlying physics yet.

energy spectrum.

diagram.

**5. Irreversible processes**

thermodynamic variables as well as microstates.

society, see Kaufman (2001) for example.

understand the phase transitions.

theorem in Eq. (17) and the Bloch equation in Eq. (122).

**5.1 Linear response theory**

and the number of particles is

$$N\left(T, V, \mu\right) = \mp V \left(2s + 1\right) \int \frac{\mathbf{d}^3 k}{\left(2\pi\right)^2} \frac{1}{\beta\hbar} \sum\_{\mathbf{n}} \frac{e^{i\omega\_{\mathbf{n}}\eta}}{\left(\omega\_{\mathbf{n}} - \hbar^{-1}\left(\varepsilon\_{\mathbf{k}}^0 - \mu\right) - \Sigma^\lambda\left(\mathbf{k}, \omega\_{\mathbf{n}}\right)\right)}.\tag{133}$$

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 of the energy spectra is available by using the temperature Green's function.

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

$$\tilde{G}\left(\mathbf{k},\omega\right) = \frac{1}{1 \mp e^{-\beta\hbar\omega}} \tilde{G}^{R}\left(\mathbf{k},\omega\right) + \frac{1}{1 \mp e^{\beta\hbar\omega}} \tilde{G}^{A}\left(\mathbf{k},\omega\right). \tag{134}$$

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

$$\bar{\Gamma}\left(\mathbf{k},z\right) = \int\_{-\infty}^{\infty} \frac{\mathbf{d}\omega'}{2\pi} \frac{A\left(\mathbf{k},\omega'\right)}{z-\omega'},\tag{135}$$

where *A* (**k**, *ω*� ) is the spectral weight function (Galitskii & Migdal, 1958), which yields the integral representations of the *retarded* Green's function

$$\bar{\mathbf{G}}^{R}\left(\mathbf{k},\omega\right) = \bar{\Gamma}\left(\mathbf{k},\omega + i\eta\right) = \int\_{-\infty}^{\infty} \frac{\mathbf{d}\omega'}{2\pi} \frac{A\left(\mathbf{k},\omega'\right)}{\omega - \omega' + i\eta'}\tag{136}$$

and the *advanced* Green's function

$$\bar{\mathbf{G}}^{A}(\mathbf{k},\omega) = \bar{\Gamma}\left(\mathbf{k},\omega - i\eta\right) = \int\_{-\infty}^{\infty} \frac{\mathbf{d}\omega'}{2\pi} \frac{A\left(\mathbf{k},\omega'\right)}{\omega - \omega' - i\eta}.\tag{137}$$

The spectral weight function should satisfy the sum rule

$$\int\_{-\infty}^{\infty} \frac{\mathbf{d}\omega'}{2\pi} A\left(\mathbf{k}, \omega'\right) = 1\tag{138}$$

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

$$\bar{\Gamma}\left(\mathbf{k},\omega\right) \sim \frac{1}{\omega} \int\_{-\infty}^{\infty} \frac{\mathbf{d}\omega'}{2\pi} A\left(\mathbf{k},\omega'\right) \sim \frac{1}{\omega} \quad \text{for} \quad |\omega| \to \infty. \tag{139}$$

The spectral function constructs also the temperature Green's function

$$\mathcal{G}\left(\mathbf{k},\omega\_{\rm ll}\right) = \int\_{-\infty}^{\infty} \frac{\mathbf{d}\omega^{\prime}}{2\pi} \frac{A\left(\mathbf{k},\omega^{\prime}\right)}{i\omega\_{\rm ll}-\omega^{\prime}}.\tag{140}$$

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 discrete set of points {*iωn*}, the thermodynamic variables Ω, *E*, and *N* are obtained. The spectral weight function is able to be derived from the temperature Green's function G. Then we can have the real-time Green's functions *G*¯ *<sup>R</sup>* and *G*¯ *<sup>A</sup>*, which give us the quasiparticle energy spectrum.

Hence, we complete the equilibrium thermodynamics; information on macroscopic thermodynamic variables as well as microstates.
