Theory of Electron Correlation in Disordered Crystals

Abstract

This paper presents a new method of describing the electronic spectrum and electrical conductivity of disordered crystals based on the Hamiltonian of electrons and phonons. Electronic states of a system are described by the tight-binding model. Expressions for Green’s functions and electrical conductivity are derived using the diagram method. Equations are obtained for the vertex parts of the mass operators of the electron–electron and electron–phonon interactions. A system of exact equations is obtained for the spectrum of elementary excitations in a crystal. This makes it possible to perform numerical calculations of the energy spectrum and to predict the properties of the system with a predetermined accuracy. In contrast to other approaches, in which electron correlations are taken into account only in the limiting cases of an infinitely large and infinitesimal electron density, in this method, electron correlations are described in the general case of an arbitrary density. The cluster expansion is obtained for the density of states and electrical conductivity of disordered systems. We show that the contribution of the electron scattering processes to clusters is decreasing, along with increasing the number of sites in the cluster, which depends on a small parameter.

1. Introduction

Advances in the description of disordered systems are mainly due to the development of the pseudopotential method [1]. However, due to the nonlocal nature of the pseudopotential, there is a problem of portability of the pseudopotential. It is impossible to use nuclear potentials, determined by the properties of some systems, in order to describe other systems. The use of the theory of Vanderbilt ultra-soft potentials [2,3] and the method of projector-augmented waves proposed by Blochl [4,5], allowed for achieving fundamental progress in investigating the electronic structure and the properties of the system. In the augmented-wave projector method, the wave function of the valence states of an electron (all-electron orbital) is expressed in terms of a pseudo-wave function. The pseudo wave function is expanded in a series of pseudo partial wave functions. The wave function is expanded in a series of partial wave functions with the same coefficients as in the expression for the pseudo wave function. Partial wave functions are described by the Schrödinger equation for non-interacting atoms. The expression for the pseudo Hamiltonian, as an equation for the pseudo wave function, is derived by minimizing the full energy functional. This approach was further developed through the use of the generalized gradient approximation proposed in [6,7,8,9,10]. The paper [10] describes the application of this method for calculating the electronic structure of crystals and molecules using the VASP and GAUSSIAN software packages, respectively.

It should be noted that in [11,12,13,14,15,16,17,18,19], the crystals electronic structure was carried out, including the Coulomb long-range interaction between electrons of different sites on the crystal lattice, thanks to a method based on the tight-binding model [20,21] and the functional density theory. However, such methods are suitable only for describing crystals characterized by ideal ordering. In disordered crystals, effects associated with localized electronic states occur. These effects cannot be described in a model where the crystal is treated as an ideal one.

Calculations of the electronic structure of an alloy are based on using the self-consistent method of the Korringa–Kohn–Rostoker-coherent potential approximations are made in work [22,23,24]. In [25], a virtual crystal approximation was proposed to study the properties of alloys by the density functional method. This approach is applied in the Vanderbilt ultra-soft pseudopotential scheme to predict the properties of Pb(Zr_0.5Ti_0.5)O₃ alloys in its paraelectric and ferroelectric phases.

In our work [26], we present a new method of describing the electronic spectrum and electrical conductivity of disordered crystals based on the Hamiltonian of electrons and phonons. Electronic states of a system are described by the tight-binding model. Calculations of two-time Green’s functions are based on temperature Green’s functions [26]. This uses a known relation between spectral representation for two-time and temperature Green’s function [27]. The calculation of the temperature Green’s functions of disordered crystal based on diagram technics are analogous to the diagram technique for a homogeneous system [27]. A system of exact equations is obtained for the spectrum of elementary excitations in a crystal. This makes it possible to perform numerical calculations of the energy spectrum and to predict the properties of the system with a predetermined accuracy.

2. Hamiltonian of an Electron–Phonon System for a Disordered Crystal

The Hamiltonian of the disordered system (alloy, disordered semiconductor, and disordered dielectric) consists of the sum of the Hamiltonian of electrons in the nucleus field, the Hamiltonian of electron–electron interaction, and the Hamilton of nucleus. The motion of the ion subsystem is reduced to nucleus oscillations near the equilibrium position under the influence of the nuclei interaction force, and their indirect interaction through electrons. In the Wannier representation, the system Hamiltonian is as follows [26]:

$$H=H_0+H_{\mathrm{int}}$$

(1)

where zero-order Hamiltonian

$$H_0=H_e^{(0)}+H_{ph}^{(0)}$$

(2)

consists of the Hamiltonian of the electrons in the field of the cores of the atom’s ideal ordered crystal

$$H_e^{(0)}={\displaystyle \sum _{\begin{array}{c}ni\mathsf{\gamma }\\ n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }\end{array}}{h}_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{\left(0\right)}{a}_{ni\mathsf{\gamma }}^{+}{a}_{n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}}$$

(3)

and the harmonic phonon Hamiltonian for the motion of the cores of the atom’s ideal ordered crystal

$$H_{ph}^{(0)}={\displaystyle \sum _{ni\mathsf{\alpha }}\frac{P_{ni\mathsf{\alpha }}^2}{2M_i}}+\frac{1}{2}{\displaystyle \sum _{\begin{array}{c}ni\mathsf{\alpha }\\ n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }\end{array}}{\Phi }_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}^{\left(0\right)}{u}_{ni\mathsf{\alpha }}}{u}_{n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}$$

(4)

Symbol n denotes the number of a unit cell, i denotes the number of a node in a unit cell, and γ denotes all of the other quantum numbers for the orbital, including spin. The symbol h⁽⁰⁾ denotes the “hopping integral” that connects the respective orbitals. For the phonon Hamiltonian, α is a spatial direction (x, y, or z), $P_{ni\mathsf{\alpha }}$ is the core momentum, $M_i$ is the mass of the core, $u_{ni\mathsf{\alpha }}$ is the deviation of the core from the equilibrium position of the lattice site, and ${\Phi }_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}^{\left(0\right)}$ is the corresponding spring-constant matrix.

The interaction Hamiltonian in Equation (1) is the perturbation of the system due to all of the effects we will be including. It is composed of six pieces:

$$H_{\mathrm{int}}=\mathsf{\delta }\Phi +H_{ec}+H_{eph}+H_{ee}+H_{phc}+H_{phph}$$

(5)

$\mathsf{\delta }\Phi$ is the modification of the core–core Coulomb interaction due to the disordered atoms added to the system; it is the difference between the original core–core repulsion Hamiltonian and the new one. The electronic Hamiltonian is modified by the term

$$H_{ec}={\displaystyle \sum _{\begin{array}{c}ni\mathsf{\gamma }\\ n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }\end{array}}{w}_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}{a}_{ni\mathsf{\gamma }}^{+}{a}_{n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}}$$

(6)

which is the difference between the new hopping Hamiltonian and the original periodic one. The electron–phonon interaction is given by

$$H_{eph}={\displaystyle \sum _{\begin{array}{c}ni\mathsf{\gamma }\\ n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }\end{array}}{v}_{ ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{\prime }{a}_{ni\mathsf{\gamma }}^{+}{a}_{n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}}$$

(7)

It is described in more detail below. The Hamiltonian of the Coulomb interaction between electrons is given by the term

$$\begin{array}{l}{H}_{ee}=\frac{1}{2}{\displaystyle \sum _{\begin{array}{l}{n}_1,n_2,\\ n_3,n_4\end{array}}{v}_{n_3,n_4}^{\left(2\right)n_1,n_2}{a}_{n_1}^{+}{a}_{n_2}^{+}{a}_{n_3}{a}_{n_4}},\\ n=\left(ni\mathsf{\gamma }\right).\end{array}$$

(8)

The modification of the interaction of the phonons with the cores caused by the disordering of the atoms is given by

$$\begin{array}{l}{H}_{phc}=\frac{1}{2}{\displaystyle \sum _{\begin{array}{c}ni\mathsf{\alpha }\\ n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }\end{array}}\Delta M_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}^{-1}{P}_{ni\mathsf{\alpha }}}{P}_{n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}+\\ +\frac{1}{2}{\displaystyle \sum _{\begin{array}{c}ni\mathsf{\alpha }\\ n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }\end{array}}\Delta {\Phi }_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}{u}_{ni\mathsf{\alpha }}}{u}_{n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }},\end{array}$$

(9)

where

$$\Delta M_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}^{-1}=\left(\frac{1}{M_{n{i}^{\prime }}}-\frac{1}{M_i}\right){\mathsf{\delta }}_{n{n}^{\prime }}{\mathsf{\delta }}_{i{i}^{\prime }}{\mathsf{\delta }}_{\mathsf{\alpha }{\mathsf{\alpha }}^{\prime }}$$

(10)

$\Delta {\Phi }_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}={\Phi }_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}-{\Phi }_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}^{\left(0\right)}$, and $M_{ni}$ and $M_i$ are the masses of the atoms at site $(ni)$ for the disordered and ordered alloy, respectively.

We also include the cubic anharmonic potential terms for the phonons (under the assumption that they remain small and can be treated perturbatively via

$$\begin{array}{l}{H}_{phph}=\frac{1}{3!}{\displaystyle \sum _{\begin{array}{l}ni\mathsf{\alpha }\\ n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }\\ n^{″}{i}^{″}{\mathsf{\alpha }}^{″}\end{array}}{\Phi }_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime },n^{″}{i}^{″}{\mathsf{\alpha }}^{″}}^{(0)}}{u}_{ni\mathsf{\alpha }}\times \\ \times u_{n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}{u}_{n^{″}{i}^{″}{\mathsf{\alpha }}^{″}}.\end{array}$$

(11)

The operators $a_{ni\mathsf{\gamma }}^{+}$ and $a_{ni\mathsf{\gamma }}^{}$ create and destroy electrons in the state described by Vane’s function ${\mathsf{\varphi }}_{ni\mathsf{\gamma }}\left(\mathsf{\xi }\right)=\langle \mathsf{\xi }|ni\mathsf{\gamma }\rangle$, where $\mathsf{\xi }=\left(\mathbf{r},{\mathsf{\sigma }}^{\prime }\right)$ are the spatial and z-components of the spin coordinates of the wave function.

To construct the Wannier functions, we use analytical expressions for the wave functions of an electron in the field of atomic nuclei of type λ localized at the lattice sites $(ni)$ of an ideally ordered crystal:

$$\begin{array}{l}{\mathsf{\psi }}_{ni\mathsf{\delta }}\left(\mathbf{r}-{\mathbf{r}}_{ni}\right)=R_{\tilde{\mathsf{\epsilon }}l}\left(\left|\mathbf{r}-{\mathbf{r}}_{ni}\right|\right)Y_{lm}\left(\widehat{\mathbf{r}-{\mathbf{r}}_{ni}}\right),\\ Y_{lm}\left(\widehat{\mathbf{r}-{\mathbf{r}}_{ni}}\right)=Y_{lm}\left(\mathsf{\theta },\mathsf{\phi }\right),\end{array}$$

(12)

where $\mathsf{\theta },\mathsf{\phi }$ are the angular spherical coordinates of the vector $\mathbf{r}-{\mathbf{r}}_{ni}.$

Here, $\mathsf{\delta }=\tilde{\mathsf{\epsilon }}lm$ is an index that incorporates the quantum numbers for the energy value $\tilde{\mathsf{\epsilon }}$, the angular momentum quantum numbers are l and m, r is the electron position vector, and ${\mathbf{r}}_{ni}$ is the position vector for the atom at site (ni) in equilibrium

$$\begin{array}{l}{\mathbf{r}}_{ni}={\mathbf{r}}_n+{\mathsf{\rho }}_i,\\ {\mathbf{r}}_n={\displaystyle \sum _{\mathsf{\nu }}{l}_{\mathsf{\nu }}{\mathbf{a}}_{\mathsf{\nu }}},\end{array}$$

(13)

${\mathbf{r}}_n$ is the position vector of the unit cell $n$ of the crystal lattice, and ${\mathsf{\rho }}_i$ is the vector of the relative position of the node of the sublattice i in the unit cell $n$. The coordinates $l_{\mathsf{\nu }}$ of the radius vector ${\mathbf{r}}_n$ of the unit cell $n$ of the crystal lattice are integers. The number $\mathsf{\nu }$ takes on values $\mathsf{\nu }=1,2,3$ for three-dimensional crystals, $\mathsf{\nu }=1,2$ for two-dimensional crystals, and $\mathsf{\nu }=1$ for one-dimensional crystals.

Basis orthogonalization is performed with the Lowdin method [28]

$$|{\tilde{\mathsf{\psi }}}_{ni\mathsf{\delta }}\rangle =S^{-1/2}|{\mathsf{\psi }}_{ni\mathsf{\delta }}\rangle ,S_{ni\mathsf{\delta },n^{\prime }{i}^{\prime }{\mathsf{\delta }}^{\prime }}=\langle {\mathsf{\psi }}_{ni\mathsf{\delta }}|{\mathsf{\psi }}_{n^{\prime }{i}^{\prime }{\mathsf{\delta }}^{\prime }}\rangle ,$$

(14)

where $S_{ni\mathsf{\delta },n^{\prime }{i}^{\prime }{\mathsf{\delta }}^{\prime }}$ are the overlapping matrix.

Vane’s functions ${\mathsf{\varphi }}_{ni\mathsf{\gamma }}\left(\mathbf{r},{\mathsf{\sigma }}^{\prime }\right)$, on which the Hamiltonian of the system are represented as in Equation (1), are defined from equation:

$${\mathsf{\varphi }}_{ni\mathsf{\gamma }}\left(\mathbf{r},{\mathsf{\sigma }}^{\prime }\right)={\tilde{\mathsf{\psi }}}_{ni\mathsf{\delta }}\left(\mathbf{r}-{\mathbf{r}}_{ni}\right){\chi }_{\mathsf{\sigma }}\left({\mathsf{\sigma }}^{\prime }\right)$$

(15)

where ${\chi }_{\mathsf{\sigma }}({\mathsf{\sigma }}^{\prime })$ is the spin part of wave function, $\mathsf{\gamma }=\mathsf{\delta }\mathsf{\sigma }$.

The orthogonalized wave function can be represented as:

$${\tilde{\mathsf{\psi }}}_{n_1{i}_1{\mathsf{\delta }}_1}\left(r_1,{\mathsf{\theta }}_1,{\mathsf{\phi }}_1\right)={\displaystyle \sum _{n_2,i_2{\mathsf{\delta }}_2}{S}_{n_2{i}_2{\mathsf{\delta }}_2,n_1{i}_1{\mathsf{\delta }}_1}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{R}_{{\tilde{\mathsf{\epsilon }}}_2{l}_2}\left(r_2\right)Y_{l_2{m}_2}\left({\mathsf{\theta }}_2,{\mathsf{\phi }}_2\right).$$

(16)

In expression (16):

$$\begin{array}{l}{\mathbf{r}}_1=\mathbf{r}-{\mathbf{r}}_{n_1{i}_1},{\mathbf{r}}_2=\mathbf{r}-{\mathbf{r}}_{n_2{i}_2}={\mathbf{r}}_1-{\mathbf{r}}_{n_2{i}_2{n}_1{i}_1},\\ r_2={\left({\left(x^1-x_{n_2{i}_2{n}_1{i}_1}^1\right)}^2+{\left(x^2-x_{n_2{i}_2{n}_1{i}_1}^2\right)}^2+{\left(x^3-x_{n_2{i}_2{n}_1{i}_1}^3\right)}^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.},\\ x^1=r_1\sin {\mathsf{\theta }}_1\cos {\mathsf{\phi }}_1,x^2=r_1\sin {\mathsf{\theta }}_1\sin {\mathsf{\phi }}_1,x^3=r_1\cos {\mathsf{\theta }}_1,\\ x_{n_2{i}_2{n}_1{i}_1}^{\mathsf{\alpha }}={\displaystyle \sum _{\mathsf{\nu }}\left(l_{\mathsf{\nu }}^{\left(2\right)}-l_{\mathsf{\nu }}^{\left(1\right)}\right)}{a}_{\mathsf{\nu }}^{\mathsf{\alpha }}+{\mathsf{\rho }}_{i_2}^{\mathsf{\alpha }}-{\mathsf{\rho }}_{i_1}^{\mathsf{\alpha }},\end{array}$$

(17)

$$\cos {\mathsf{\theta }}_2=\frac{r_1\cos {\mathsf{\theta }}_1-x_{n_2{i}_2{n}_1{i}_1}^3}{r_2},$$

(18)

$${\mathsf{\phi }}_2=\arccos \frac{r_1\sin {\mathsf{\theta }}_1\cos {\mathsf{\phi }}_1-x_{n_2{i}_2{n}_1{i}_1}^1}{r_2{\left(1-{\cos }^2{\mathsf{\theta }}_2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}.$$

(19)

Summation over $n_2{i}_2$ in expression (16) means summation over ${\mathbf{r}}_{n_2{i}_2}$, in accordance with Formula (13).

The overlap matrix $S_{ni\mathsf{\delta },n^{\prime }{i}^{\prime }{\mathsf{\delta }}^{\prime }}$ is found from the equation:

$$\begin{array}{l}{S}_{n_1{i}_1{\mathsf{\delta }}_1,n_2{i}_2{\mathsf{\delta }}_2}=\\ {\displaystyle \iiint R_{{\tilde{\mathsf{\epsilon }}}_1{l}_1}\left(r_1\right)Y_{l_1{m}_1}^{\ast }\left({\mathsf{\theta }}_1,{\mathsf{\phi }}_1\right)}{R}_{{\tilde{\mathsf{\epsilon }}}_2{l}_2}\left(r_2\right)Y_{l_2{m}_2}\left({\mathsf{\theta }}_2,{\mathsf{\phi }}_2\right)r_1^2\sin {\mathsf{\theta }}_{1}d{r}_{1}d{\mathsf{\theta }}_{1}d{\mathsf{\phi }}_1\end{array}$$

(20)

where $r_2,{\mathsf{\theta }}_2,{\mathsf{\phi }}_2$ are expressed through $r_1,{\mathsf{\theta }}_1,{\mathsf{\phi }}_1$ in accordance with Formulas (17)–(19).

To find matrix $S_{n_2{i}_2{\mathsf{\delta }}_2,n_1{i}_1{\mathsf{\delta }}_1}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$ in expression (16), we find the Fourier transform of the matrix (20):

$$S_{i_1{\mathsf{\delta }}_1,i_2{\mathsf{\delta }}_2}\left(\mathbf{k}\right)={\displaystyle \sum _{n_2}{S}_{n_1{i}_1{\mathsf{\delta }}_1,n_2{i}_2{\mathsf{\delta }}_2}}{e}^{i\mathbf{k}({\mathbf{r}}_{n_2{i}_2}-{\mathbf{r}}_{n_1{i}_1})}.$$

(21)

The vector $\mathbf{k}$ is defined by the expression

$$\begin{array}{l}\mathbf{k}={\displaystyle \sum _{\mathsf{\nu }}{k}_{\mathsf{\nu }}{\mathbf{b}}_{\mathsf{\nu }}},\\ \left({\mathbf{a}}_{\mathsf{\nu }}{\mathbf{b}}_{{\mathsf{\nu }}^{\prime }}\right)=2\pi {\mathsf{\delta }}_{\mathsf{\nu }{\mathsf{\nu }}^{\prime }},\end{array}$$

(22)

${\mathbf{b}}_{\mathsf{\nu }}$ is the basis vector of the translations of the reciprocal lattice.

Summing over $n_2$ on the right-hand side of Formula (21) is easy to do if we replace it according to (13) and use

$$\mathbf{k}\left({\mathbf{r}}_{n_2{i}_2}-{\mathbf{r}}_{n_1{i}_1}\right)={\displaystyle \sum _{\mathsf{\alpha }=1}^3{\displaystyle \sum _{{\mathsf{\nu }}^{\prime }}{k}_{{\mathsf{\nu }}^{\prime }}}}{b}_{{\mathsf{\nu }}^{\prime }}^{\mathsf{\alpha }}\left({\displaystyle \sum _{\mathsf{\nu }}\left(l_{\mathsf{\nu }}^{(2)}-l_{\mathsf{\nu }}^{(1)}\right)a_{\mathsf{\nu }}^{\mathsf{\alpha }}}+{\mathsf{\rho }}_{i_2}^{\mathsf{\alpha }}-{\mathsf{\rho }}_{i_1}^{\mathsf{\alpha }}\right).$$

(23)

As the matrix element $S_{n_1{i}_1{\mathsf{\delta }}_1,n_2{i}_2{\mathsf{\delta }}_2}$ decreases with the distance between the nodes $n_1{i}_1$, $n_2{i}_2$, in numerical calculations, when summing over $n_2$ in expression (21), it is sufficient to restrict ourselves to a few coordination spheres. In this case, summation over $n_2$ is reduced to summation over $l_{\mathsf{\nu }}^{(2)}$.

The matrix $S_{n_1{i}_1{\mathsf{\delta }}_1,n_2{i}_2{\mathsf{\delta }}_2}$ has an infinite rank. The rank of the matrix $S_{i_1{\mathsf{\delta }}_1,i_2{\mathsf{\delta }}_2}\left(\mathbf{k}\right)$ is finite, which allows for finding the matrix $S_{i_1{\mathsf{\delta }}_1,i_2{\mathsf{\delta }}_2}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\left(\mathrm{k}\right)$. The matrix $S_{n_2{i}_2{\mathsf{\delta }}_2,n_1{i}_1{\mathsf{\delta }}_1}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$ in expression (16) is found from equation:

$$S_{n_1{i}_1{\mathsf{\delta }}_1,n_2{i}_2{\mathsf{\delta }}_2}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}=\frac{1}{N}{\displaystyle \sum _{\mathbf{k}}{S}_{i_1{\mathsf{\delta }}_1,i_2{\mathsf{\delta }}_2}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\left(\mathbf{k}\right)}{e}^{-\mathrm{i}\mathbf{k}({\mathbf{r}}_{n_2{i}_2}-{\mathbf{r}}_{n_1{i}_1})}.$$

(24)

The values $h_{n_1{i}_1{\mathsf{\gamma }}_1,n_2{i}_2{\mathsf{\gamma }}_2}^{\left(0\right)}$ in Equation (3) are the matrix elements of the kinetic and potential energy ${\displaystyle \sum _{ni}{v}^{\mathsf{\lambda }i}\left(\mathbf{r}-{\mathbf{r}}_{ni}\right)}$ of the electron in the field of the cores of the atom’s ideal ordered crystal. The values $h_{n_1{i}_1{\mathsf{\gamma }}_1,n_2{i}_2{\mathsf{\gamma }}_2}^{\left(0\right)}$ are defined by the expression:

$$h_{n_1{i}_1{\mathsf{\gamma }}_1,n_2{i}_2{\mathsf{\gamma }}_2}^{\left(0\right)}={\displaystyle \sum _{n_3{i}_3{\mathsf{\delta }}_3}{E}_{i_1{\tilde{\mathsf{\epsilon }}}_1}{S}_{n_3{i}_3{\mathsf{\delta }}_3,n_1{i}_1{\mathsf{\delta }}_1}^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\ast }{S}_{n_3{i}_3{\mathsf{\delta }}_3,n_2{i}_2{\mathsf{\delta }}_2}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{\mathsf{\delta }}_{{\mathsf{\sigma }}_1,{\mathsf{\sigma }}_2}+{\displaystyle \sum _{n_3{i}_3\ne n_1{i}_1}{v}_{n_1{i}_1{\mathsf{\gamma }}_1,n_2{i}_2{\mathsf{\gamma }}_2}^{n_3{i}_3}},\mathsf{\gamma }=\mathsf{\delta }\mathsf{\sigma }.$$

(25)

In the Formula (25)

$$\begin{gathered} \begin{array}{l}{v}_{n_1{i}_1{\mathsf{\gamma }}_1,n_2{i}_2{\mathsf{\gamma }}_2}^{n_3{i}_3}=\\ {\displaystyle \iiint {\tilde{\mathsf{\psi }}}_{n_1{i}_1{\mathsf{\delta }}_1}^{\ast }\left(r_1,{\mathsf{\theta }}_1,{\mathsf{\phi }}_1\right)}{v}^{i_3}\left(r_3\right){\tilde{\mathsf{\psi }}}_{n_2{i}_2{\mathsf{\delta }}_2}\left(r_2,{\mathsf{\theta }}_2,{\mathsf{\phi }}_2\right)r_1^2\sin {\mathsf{\theta }}_{1}d{r}_{1}d{\mathsf{\theta }}_{1}d{\mathsf{\phi }}_1\times {\mathsf{\delta }}_{{\mathsf{\sigma }}_1,{\mathsf{\sigma }}_2}\end{array}, \\ {E}_{i_1{\tilde{\mathsf{\epsilon }}}_1}=-\frac{m{e}^4{\left(Z_{i_1}\right)}^2}{2{\hslash }^2{\tilde{\mathsf{\epsilon }}}_1^2},{\tilde{\mathsf{\epsilon }}}_1=1,2,3,\dots , \end{gathered}$$

(26)

Here, $r_2,{\mathsf{\theta }}_2,{\mathsf{\phi }}_2$ is expressed through $r_1,{\mathsf{\theta }}_1,{\mathsf{\phi }}_1$ in accordance with Formulas (17)–(19). The expression for $r_3$ is obtained from expression (17) for $r_2$ replacement $x_{n_2{i}_2{n}_1{i}_1}^{\mathsf{\alpha }}$ by $x_{n_3{i}_3{n}_1{i}_1}^{\mathsf{\alpha }}$. Summation over $n_3{i}_3$ in expression (25) means summation over ${\mathbf{r}}_{n_3{i}_3}$, in accordance with Formula (13).

In Formula (26) and $e$ are the mass and charge of the electron, respectively, and $Z_i$ is the ordinal number of an atom of the sort $\mathsf{\lambda }$ located in the site $ni$ of an ideally ordered crystal. $\hslash$ denotes the Planck’s constant.

The matrix element of the electron–ion interaction Hamiltonian in Equation (6) is given by

$$w_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}={\displaystyle \sum _{n^{″}{i}^{″}}{w}_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{n^{″}{i}^{″}}}$$

(27)

where

$$w_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{n^{″}{i}^{″}}={\displaystyle \sum _{\mathsf{\lambda }}{c}_{n^{″}{i}^{″}}^{\mathsf{\lambda }} w_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{\mathsf{\lambda }{n}^{″}{i}^{″}}},$$

(28)

$$w_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{\mathsf{\lambda }{n}^{″}{i}^{″}}=v_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{\mathsf{\lambda }{n}^{″}{i}^{″}}+\Delta v_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{\mathsf{\lambda }{n}^{″}{i}^{″}}-v_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{i^{″}}$$

(29)

$v_{n_1{i}_1{\mathsf{\gamma }}_1,n_2{i}_2{\mathsf{\gamma }}_2}^{\mathsf{\lambda }{n}^{″}{i}^{″}}$ is a matrix element of the potential of the core of the atom $v^{\mathsf{\lambda }}\left(\mathbf{r}-{\mathbf{r}}_{n^{″}{i}^{″}}\right)$.

The expression for $v_{n_1{i}_1{\mathsf{\gamma }}_1,n_2{i}_2{\mathsf{\gamma }}_2}^{\mathsf{\lambda }{n}_3{i}_3}$ is obtained from Formula (26) by replacing $v^{i_3}\left(r_3\right)$ with $v^{\mathsf{\lambda }}\left(r_3\right)$.

In Equation (28), $c_{ni}^{\mathsf{\lambda }}$ is a discrete binary random number taking the values of 1 or 0, depending on whether an atom of type $\mathsf{\lambda }$ is at site (ni) or not, respectively. The symbol $\Delta v_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{\mathsf{\lambda }{n}^{″}{i}^{″}}$ will be defined next.

The expression for the electron–phonon interaction in Equation (7) is found through derivatives of the potential energy of the electrons in the ion core field due to a displacement of the atom by vector ${\mathrm{u}}_{ni}$. In Equation (7), the value of $v_{ ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{\prime }$ is given by

$$v_{ ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{\prime }={\displaystyle \sum _{n^{″}{i}^{″}\mathsf{\alpha }}{v}^{\prime }{}_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{n^{″}{i}^{″}\mathsf{\alpha }}}{u}_{n^{″}{i}^{″}\mathsf{\alpha }},$$

(30)

where

$$v^{\prime }{}_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{n^{″}{i}^{″}\mathsf{\alpha }}={\displaystyle \sum _{\mathsf{\lambda }}{c}_{n^{″}{i}^{″}}^{\mathsf{\lambda }} v^{\prime }{}_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{\mathsf{\lambda }{n}^{″}{i}^{″}\mathsf{\alpha }}}$$

(31)

with $v^{\prime }{}_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{\mathsf{\lambda }{n}^{″}{i}^{″}\mathsf{\alpha }}$ the matrix elements of the following operator:

$$\begin{gathered} -{\mathrm{e}}_{n^{″}{i}^{″}\mathsf{\alpha }}\frac{d}{d\left|\mathbf{r}-{\mathbf{r}}_{n^{″}{i}^{″}}\right|}{v}^{\mathsf{\lambda }{i}^{″}}\left(\left|\mathbf{r}-{\mathbf{r}}_{n^{″}{i}^{″}}\right|\right), \\ {\mathrm{e}}_{n^{″}{i}^{″}}=\frac{\mathbf{r}-{\mathbf{r}}_{n^{″}{i}^{″}}}{\left|\mathbf{r}-{\mathbf{r}}_{n^{″}{i}^{″}}\right|}. \end{gathered}$$

(32)

The expression for $v^{\prime }{}_{n_1{i}_1{\mathsf{\gamma }}_1,n_2{i}_2{\mathsf{\gamma }}_2}^{\mathsf{\lambda }{n}_3{i}_3\mathsf{\alpha }}$ is obtained from Formula (26) by replacing in it $v^{i_3}\left(r_3\right)$ with

$$-\frac{\left(x_{\mathsf{\alpha }}-x_{n_3{i}_3{n}_1{i}_1}^{\mathsf{\alpha }}\right)}{r_3}\frac{d}{d{r}_3}{v}^{\mathsf{\lambda }{n}_3{i}_3}\left(r_3\right).$$

(33)

$\Delta v_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{\mathsf{\lambda }{n}^{″}{i}^{″}}$ in Equation (29) describes electron scattering on the static displacement of the atoms, and is defined by the equation

$$\Delta v_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{\mathsf{\lambda }{n}^{″}{i}^{″}}={\displaystyle \sum _{\mathsf{\alpha }}{v}^{\prime }{}_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{\mathsf{\lambda }{n}^{″}{i}^{″}\mathsf{\alpha }}}{u}_{n^{″}{i}^{″}\mathsf{\alpha }}^{s,\mathsf{\lambda }}$$

(34)

where $u_{n^{″}{i}^{″}\mathsf{\alpha }}^{s,\mathsf{\lambda }}$ is the $\mathsf{\alpha }$ projection of the static displacement of the atom of type $\mathsf{\lambda }$ in the site, and $n^{″}{i}^{″}$ I caused by the difference in the atomic radii of the components of the disordered crystal.

Upon receipt of expressions (27)–(34), it was taken into account that the potential energy operator of the electron in the field of the atoms core can be expressed as:$v^{ni}\left(\mathrm{r}-{{\mathrm{r}}^{\prime }}_{ni}\right)$, $\mathbf{r}{\prime }_{ni}={\mathbf{r}}_{ni}+{\mathbf{u}}_{ni}^s+{\mathbf{u}}_{ni}$, with r being the electron’s radius vector, ${\mathbf{r}}_{ni}$ the radius-vector of atom’s equilibrium position in the site of the crystal lattice (ni), ${\mathbf{u}}_{ni}^s$ the vector of atom’s static displacement from equilibrium position in site (ni), and ${\mathbf{u}}_{ni}$ the atom’s displacement operator in site (ni). Expanding $v^{ni}\left(\mathrm{r}-{{\mathrm{r}}^{\prime }}_{ni}\right)$ in the series in powers $u_{ni\mathsf{\alpha }}$ and restricting ourselves to linear terms, we arrive at expressions (27)–(34).

The matrix of the force constants arising from the direct Coulomb interaction of the ionic cores has the form:

$$\begin{array}{l}{\Phi }_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}=-\frac{Z_{ni}{Z}_{n^{\prime }{i}^{\prime }}{e}^2}{4\pi {\mathsf{\epsilon }}_0{\left|{\mathbf{r}}_n+{\mathsf{\rho }}_i-{\mathbf{r}}_{n^{\prime }}-{\mathsf{\rho }}_{i^{\prime }}\right|}^5}\times \\ \times [3\left(r_{n\mathsf{\alpha }}+{\mathsf{\rho }}_{i\mathsf{\alpha }}-r_{n^{\prime }\mathsf{\alpha }}-{\mathsf{\rho }}_{i^{\prime }\mathsf{\alpha }}\right)\left(r_{n{\mathsf{\alpha }}^{\prime }}+{\mathsf{\rho }}_{i{\mathsf{\alpha }}^{\prime }}-r_{n^{\prime }{\mathsf{\alpha }}^{\prime }}-{\mathsf{\rho }}_{i^{\prime }{\mathsf{\alpha }}^{\prime }}\right)-\\ -{\left|{\mathbf{r}}_n+{\mathsf{\rho }}_i-{\mathbf{r}}_{n^{\prime }}-{\mathsf{\rho }}_{i^{\prime }}\right|}^2{\mathsf{\delta }}_{\mathsf{\alpha }{\mathsf{\alpha }}^{\prime }}],ni\ne n^{\prime }{i}^{\prime }.\end{array}$$

(35)

where $Z_{ni}$ is the serial number of the atom located in the lattice site $ni$ of the disordered crystal, which is given by the expression

$$Z_{ni}={\displaystyle \sum _{\mathsf{\lambda }}{c}_{ni}^{\mathsf{\lambda }} }{Z}_i.$$

(36)

This matrix ${\Phi }_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}$ satisfies the following constraint:

$${\displaystyle \sum _{n^{\prime }{i}^{\prime }}{\Phi }_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}}=0.$$

(37)

Multicenter integrals $v_{n_3,n_4}^{\left(2\right)n_1,n_2},n=\left(ni\mathsf{\gamma }\right)$ in Formula (8) can be represented as

$$\begin{array}{l}{v}_{n_3{i}_3{\mathsf{\gamma }}_3,n_4{i}_4{\mathsf{\gamma }}_4}^{n_1{i}_1{\mathsf{\gamma }}_1,n_2{i}_2{\mathsf{\gamma }}_2}=e^2{\mathsf{\delta }}_{{\mathsf{\sigma }}_1{\mathsf{\sigma }}_4}{\mathsf{\delta }}_{{\mathsf{\sigma }}_2{\mathsf{\sigma }}_3}{\displaystyle \iint \frac{1}{\left|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{″}\right|}}\\ \times {\tilde{\mathsf{\psi }}}_{n_1{i}_1{\mathsf{\delta }}_1}^{\ast }\left(r_1^{\prime },{\mathsf{\theta }}_1^{\prime },{\mathsf{\phi }}_1^{\prime }\right){\tilde{\mathsf{\psi }}}_{n_2{i}_2{\mathsf{\delta }}_2}^{\ast }\left(r_1^{″},{\mathsf{\theta }}_1^{″},{\mathsf{\phi }}_1^{″}\right){\tilde{\mathsf{\psi }}}_{n_3{i}_3{\mathsf{\delta }}_3}\left(r_2^{″},{\mathsf{\theta }}_2^{″},{\mathsf{\phi }}_2^{″}\right){\tilde{\mathsf{\psi }}}_{n_4{i}_4{\mathsf{\delta }}_4}\left(r_2^{\prime },{\mathsf{\theta }}_2^{\prime },{\mathsf{\phi }}_2^{\prime }\right)\\ \times d^3{r}_1^{\prime }{d}^3{r}_1^{″}.\end{array}$$

(38)

In Formula (38)

$$\left|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{″}\right|={\left({\displaystyle \sum _{\mathsf{\alpha }}{\left(x^{\prime }{}^{\mathsf{\alpha }}-x^{″}{}^{\mathsf{\alpha }}-x_{n_2{i}_2{n}_1{i}_1}^{\mathsf{\alpha }}\right)}^2}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.},$$

(39)

$$d^3{r}_1^{\prime }=r_1^{\prime }{}^2\sin {\mathsf{\theta }}_1^{\prime }d{r}_1^{\prime }d{\mathsf{\theta }}_1^{\prime }d{\mathsf{\phi }}_1^{\prime },$$

(40)

When integrating over $r_1^{\prime },{\mathsf{\theta }}_1^{\prime },{\mathsf{\phi }}_1^{\prime }$ in Formula (38), $r_2^{\prime },{\mathsf{\theta }}_2^{\prime },{\mathsf{\phi }}_2^{\prime }$ should be expressed through $r_1^{\prime },{\mathsf{\theta }}_1^{\prime },{\mathsf{\phi }}_1^{\prime }$, in accordance with Formulas (17)–(19), in which I is necessary to replace $x_{n_2{i}_2{n}_1{i}_1}^{\mathsf{\alpha }}$ with $x_{n_4{i}_4{n}_1{i}_1}^{\mathsf{\alpha }}$. When integrating over $r_1^{″},{\mathsf{\theta }}_1^{″},{\mathsf{\phi }}_1^{″}$ in Formula (38), $r_2^{″},{\mathsf{\theta }}_2^{″},{\mathsf{\phi }}_2^{″}$ should be expressed through $r_1^{″},{\mathsf{\theta }}_1^{″},{\mathsf{\phi }}_1^{″}$, in accordance with Formulas (17)–(19), in which it is necessary to replace $x_{n_2{i}_2{n}_1{i}_1}^{\mathsf{\alpha }}$ with $x_{n_3{i}_3{n}_2{i}_2}^{\mathsf{\alpha }}$.

So, Formulas (17)–(19) describe the procedure for calculating the matrix elements $h_{n_1{i}_1{\mathsf{\gamma }}_1,n_2{i}_2{\mathsf{\gamma }}_2}^{\left(0\right)}$, $v_{n_3{i}_3{\mathsf{\gamma }}_3,n_4{i}_4{\mathsf{\gamma }}_4}^{\left(2\right)n_1{i}_1{\mathsf{\gamma }}_1,n_2{i}_2{\mathsf{\gamma }}_2}$ Hamiltonian (1), containing one-electron and two-electron integrals.

3. Green’s Functions of Electrons and Phonons

We employ a Green’s function-based formalism to perform the calculations. Ultimately, we need the real-time retarded $G_r^{AB}(t,t^{\prime })$ and advanced $G_a^{AB}(t,t^{\prime })$ Green’s functions are each defined as follows [26]:

$$\begin{gathered} G_r^{AB}(t,t^{\prime })=-\frac{i}{\hslash }\mathsf{\theta }(t-t^{\prime })<[\tilde{A}(t),\tilde{B}(t^{\prime })]>, \\ {G}_a^{AB}(t,t^{\prime })=\frac{i}{\hslash }\mathsf{\theta }(t^{\prime }-t)<[\tilde{A}(t),\tilde{B}(t^{\prime })]>. \end{gathered}$$

(41)

Here, the operators are expressed in the Heisenberg representation

$$\tilde{A}(t)=e^{i\mathsf{H} t/\hslash }A e^{-i\mathsf{H} t/\hslash },$$

(42)

where $\hslash$ is Planck’s constant, $\mathsf{H}=H-{\mathsf{\mu }}_{\mathrm{e}}{N}_e$, ${\mathsf{\mu }}_{\mathrm{e}}$ is the chemical potential of the electronic subsystem, and $N_e$ is the electron number operator given by

$$N_e={\displaystyle \sum _{ni\mathsf{\gamma }}{a}_{ni\mathsf{\gamma }}^{+}{a}_{ni\mathsf{\gamma }}^{}}.$$

(43)

In addition, the commutator or anticommutator is defined via

$$[A,B]=AB\mp BA,$$

(44)

where the commutator is used for Bose operators (−) and the anticommutator is used for Fermi operators (+). The symbol $\mathsf{\theta }(t)$ is Heaviside’s unit step function. The angle brackets $\langle \dots \rangle$ denote the thermal averaging with respect to the density matrix ρ

$$<A>=Tr (\mathsf{\rho }A), \mathsf{\rho }=e^{(\Omega -\mathsf{H})/\Theta },$$

(45)

where $\Omega$ is the thermodynamic potential of the system given by exp(Ω/Θ) = Trexp(−H/Θ) and $\Theta =k_{b}T$, with k_b Boltzmann’s constant and $T$ the temperature.

The thermal Green’s functions are defined by

$$G_{}^{AB}(\mathsf{\tau },{\mathsf{\tau }}^{\prime })=-<T_{\mathsf{\tau }}\tilde{A}(\mathsf{\tau })\tilde{B}({\mathsf{\tau }}^{\prime })>,$$

(46)

where the imaginary-time operator $\tilde{A}(\mathsf{\tau })$ is derived from the real-time Heisenberg representation and the substitution $t=-i\hslash \mathsf{\tau }$. Hence,

$$\tilde{A}(\mathsf{\tau })=e^{\mathsf{H} \mathsf{\tau }}A e^{-\mathsf{H} \mathsf{\tau }}.$$

(47)

In addition, the time-ordering operator satisfies

$$\begin{array}{l}{T}_{\mathsf{\tau }}\tilde{A}(\mathsf{\tau })\tilde{B}({\mathsf{\tau }}^{\prime })=\mathsf{\theta }(\mathsf{\tau }-{\mathsf{\tau }}^{\prime })\tilde{A}(\mathsf{\tau })\tilde{B}({\mathsf{\tau }}^{\prime })+\\ \pm \mathsf{\theta }({\mathsf{\tau }}^{\prime }-\mathsf{\tau })\tilde{B}({\mathsf{\tau }}^{\prime })\tilde{A}(\mathsf{\tau })\end{array},$$

(48)

where the plus sign is used for Bose operators and the minus sign for Fermi operators.

We next go to the interaction representation by introducing the operator

$$\mathsf{\sigma }(\mathsf{\tau })=e^{{\mathsf{H}}_0 \mathsf{\tau }}{e}^{-\mathsf{H} \mathsf{\tau }},$$

(49)

with $\mathsf{H}={\mathsf{H}}_0+H_{\mathrm{int}}$ and ${\mathsf{H}}_0=H_0-{\mathsf{\mu }}_e{N}_e$.

Differentiating the expression for $\mathsf{\sigma }(\mathsf{\tau })$ in Equation (68) with respect to $\mathsf{\tau }$ and then integrating starting from 0, with the boundary condition $\mathsf{\sigma }(0)=1$, we obtain:

$$\mathsf{\sigma }(\mathsf{\tau })=T_{\mathsf{\tau }}\exp \left[-{\displaystyle \underset{0}{\overset{\mathsf{\tau }}{\int }}{H}_{\mathrm{int}}({\mathsf{\tau }}^{\prime })d{\mathsf{\tau }}^{\prime }}\right],$$

(50)

where $H_{\mathrm{int}}(\mathsf{\tau })=e^{{\mathsf{H}}_0 \mathsf{\tau }}{H}_{\mathrm{int}} e^{-{\mathsf{H}}_0 \mathsf{\tau }}$. Employing this result yields

$$\tilde{A}(\mathsf{\tau })={\mathsf{\sigma }}^{-1}(\mathsf{\tau })A(\mathsf{\tau })\mathsf{\sigma }(\mathsf{\tau }),$$

(51)

with A(τ) in the Heisenberg representation with respect to the noninteracting Hamiltonian. Substituting these results into the definition of the thermal Green’s function creates the alternate interaction-representation form for the Green’s function, given by

$$G^{AB}(\mathsf{\tau },{\mathsf{\tau }}^{\prime })=-\frac{<T_{\mathsf{\tau }}A(\mathsf{\tau })B({\mathsf{\tau }}^{\prime })\mathsf{\sigma }(1/\Theta ){>}_0}{<\mathsf{\sigma }(1/\Theta ){>}_0}$$

(52)

where all time dependence is with respect to the noninteracting Hamiltonian and the trace over all states is with respect to the noninteracting states

$$<A{>}_0=Tr({\mathsf{\rho }}_{0}A), {\mathsf{\rho }}_0=e^{({\Omega }_0^{}-H_0)/\Theta }.$$

(53)

Next, we expand expression (51) for $\mathsf{\sigma }(1/\mathsf{\theta })$ in a series in powers of the interaction Hamiltonian $H_{\mathrm{int}}^{}(\mathsf{\tau })$ and substitute this expression in Formula (53).

The diagrammatic method is generated by expanding σ(τ) in a power series in terms of $H_{\mathrm{int}}^{}(\mathsf{\tau })$, and then using Wick’s theorem to evaluate the resulting operator averages.

The numbers of quantum states for different operators in the interaction Hamiltonian $H_{\mathrm{int}}^{}(\mathsf{\tau })$ (5)–(11), (51) are different, and the values of the argument $\mathsf{\tau }$ are the same.

Each operator can be assigned a quantum state number and an argument number $\mathsf{\tau }$, if in expression (51) for $\mathsf{\sigma }(1/\mathsf{\theta })$ the operator $H_{\mathrm{int}}(\mathsf{\tau })$ is replaced by an operator $H_{\mathrm{int}}(\mathsf{\tau },{\mathsf{\tau }}_1,\dots ,{\mathsf{\tau }}_k)$ in which the values of the argument $\mathsf{\tau }$ for operators with different quantum states are different, the matrix elements differ from the matrix elements of the operator $H_{\mathrm{int}}(\mathsf{\tau })$ by a factor $\mathsf{\delta }(\mathsf{\tau }-{\mathsf{\tau }}_1)\dots \mathsf{\delta }(\mathsf{\tau }-{\mathsf{\tau }}_k)$, and the single integral over $\mathsf{\tau }$ is replaced by the integral over $\mathsf{\tau },{\mathsf{\tau }}_1,\dots ,{\mathsf{\tau }}_k$ multiplicity $k+1$. The multiplicity of the integral is different for different types of interaction. In expression (53) for $G^{AB}(\mathsf{\tau },{\mathsf{\tau }}^{\prime })$, the term of the series for $\mathsf{\sigma }(1/\mathsf{\theta })$ (51) forms a multiple sum over quantum states and an integral over $\mathsf{\tau }$ of the mean T-product of operators $H_{\mathrm{int}}(\mathsf{\tau },{\mathsf{\tau }}_1,\dots ,{\mathsf{\tau }}_k)$ multiplied by an operator $A\left(\mathsf{\tau }\right)B\left({\mathsf{\tau }}^{\prime }\right)$. The T-product of operators is averaged over the occupation numbers of the quantum states of the system of noninteracting electrons and phonons, in accordance with Formula (53). The numbers of the quantum states for the operators in the indicated T-product are ordered by the matrix elements of the interaction operators $H_{\mathrm{int}}(\mathsf{\tau },{\mathsf{\tau }}_1,\dots ,{\mathsf{\tau }}_k)$, in accordance with Formulas (5)–(11), in such a way that pairs of operators are formed. This is due to the fact that among the average T-products of operators, only those in which the number of operators is even for the electron subsystem and the phonon subsystem are nonzero. The quantum state for each operator of the pair, except for the operators $A\left(\mathsf{\tau }\right)$, $B\left({\mathsf{\tau }}^{\prime }\right)$, coincides with one of the quantum states for the corresponding matrix element of the interaction operators $H_{\mathrm{int}}(\mathsf{\tau },{\mathsf{\tau }}_1,\dots ,{\mathsf{\tau }}_k)$ in the given product.

Let us give the averaging technique in expression (71) a simpler form. For this, in the T-product of each pair of operators $a_{n_1}\left({\mathsf{\tau }}_1\right)a_{n_2}^{+}\left({\mathsf{\tau }}_2\right),n=\left(ni\mathsf{\gamma }\right)$ for the electron subsystem and $u_{n_1}\left({\mathsf{\tau }}_1\right)u_{n_2}\left({\mathsf{\tau }}_2\right)$, $P_{n_1}\left({\mathsf{\tau }}_1\right)P_{n_2}\left({\mathsf{\tau }}_2\right)$, ${\mathrm{u}}_{n_1}\left({\mathsf{\tau }}_1\right)P_{n_2}\left({\mathsf{\tau }}_2\right)$, $P_{n_1}\left({\mathsf{\tau }}_1\right)u_{n_2}\left({\mathsf{\tau }}_2\right)$, $n=\left(ni\mathsf{\alpha }\right)$ for the phonon subsystem, in the Hamiltonian of the system of noninteracting electrons and phonons ${\mathsf{H}}_0$ (2)–(4), (53), we compare the sum of the products of pairs of operators ${\mathsf{H}}_{0n_1{n}_2}$, the numbers of quantum states of which coincide with the numbers of quantum states depending on $\mathsf{\tau }$ the operators of the pair.

Provided that the numbers of quantum states for the operators in the T-product are ordered by the matrix elements of the interaction operators $H_{\mathrm{int}}(\mathsf{\tau },{\mathsf{\tau }}_1,\dots ,{\mathsf{\tau }}_k)$ standing in it, the operators $\exp \left(-{\mathrm{H}}_{0n_1{n}_2}/\mathrm{\theta }\right)$, $\exp \left(-{\mathrm{H}}_{0n_1{n}_2}\mathsf{\tau }\right)$ change places with the products depending on $\mathsf{\tau }$ other pairs of operators. It follows from this that the average of the T-product of several operators in expression (53) is equal to the product of the average T-products of pairs of operators that determine the matrices of the Green’s functions for the zero-order Hamiltonian ${\mathsf{H}}_0$. This statement also extends to the case when the quantum states for the operators of a pair coincide with the quantum states for the operators of other pairs. This follows from the fact that the distribution function of a system of an infinite number of particles over the occupation numbers of quantum states has a sharp maximum, and the most probable value of a physical quantity is equal to its average value. The quantum state $n_i$ and the argument ${\mathsf{\tau }}_i$ for each operator of the pair, except for the operators $A\left(\mathsf{\tau }\right)$, $B\left({\mathsf{\tau }}^{\prime }\right)$, coincides with one of the quantum states $n_i$ and arguments ${\mathsf{\tau }}_i$ for the corresponding matrix element of the interaction operators $H_{\mathrm{int}}(\mathsf{\tau },{\mathsf{\tau }}_1,\dots ,{\mathsf{\tau }}_k)$ in the given product. In expression (71), the Green’s function $G^{AB}(\mathsf{\tau },{\mathsf{\tau }}^{\prime })$ is summed over quantum states $n_i$ and integrated over arguments ${\mathsf{\tau }}_i$.

The averaging technique described above in expression (71) for the Green’s function $G^{AB}(\mathsf{\tau },{\mathsf{\tau }}^{\prime })$ is the essence of Wick’s theorem. This technique then generalizes the approach used for the homogeneous system [27].

The technique for calculating the Green’s function $G^{AB}(\mathsf{\tau },{\mathsf{\tau }}^{\prime })$ (53) becomes clearer if the terms on the right-hand side of Equation (53) are represented in the form of diagrams. If the Green’s function of the system is expressed as a series only over connected diagrams [27], then the denominator in Formula (53) will cancel out with the same factor in the numerator. So, the thermal Green’s function is expanded in terms of connected diagrams. The indicated diagrammatic series can be easily summed up, which makes it possible to go beyond the framework of the first approximations of the perturbation theory and obtain equations for the Green’s functions of the system.

Summing up the indicated series, using the standard relation between the spectral representations of the temperature and real-time Green’s functions and performing an analytical continuation on the real axis, we obtain the following equations for the retarded Green’s functions [26] (hereinafter the dependence on r is suppressed):

$$\begin{gathered} G^{a{a}^{+}}(\mathsf{\epsilon })=G_0^{a{a}^{+}}(\mathsf{\epsilon })+G_0^{a{a}^{+}}(\mathsf{\epsilon }) \left(w+{\Sigma }_{eph}(\mathsf{\epsilon })+{\Sigma }_{ee}(\mathsf{\epsilon })\right) G_{}^{a{a}^{+}}(\mathsf{\epsilon }) \\ \begin{array}{l}{G}^{uu}(\mathsf{\epsilon })=G_0^{uu}(\mathsf{\epsilon })+G_0^{uu}(\mathsf{\epsilon }) \left(\Delta \Phi +{\Sigma }_{phe}(\mathsf{\epsilon })+{\Sigma }_{phph}(\mathsf{\epsilon })\right)\times \\ \times G_{}^{uu}(\mathsf{\epsilon })+G_0^{uP}(\mathsf{\epsilon })\Delta M^{-1}{G}_{}^{Pu}(\mathsf{\epsilon }),\end{array} \\ \begin{array}{l}{G}^{PP}(\mathsf{\epsilon })=G_0^{PP}(\mathsf{\epsilon })+G_0^{PP}(\mathsf{\epsilon })\Delta M^{-1}{G}_{}^{PP}(\mathsf{\epsilon })+G_0^{Pu}(\mathsf{\epsilon })\times \\ \times \left(\Delta \Phi +{\Sigma }_{phe}(\mathsf{\epsilon })+{\Sigma }_{phph}(\mathsf{\epsilon })\right) G_{}^{uP}(\mathsf{\epsilon }),\end{array} \\ \begin{array}{l}{G}^{uP}(\mathsf{\epsilon })=G_0^{uP}(\mathsf{\epsilon })+G_0^{uP}(\mathsf{\epsilon })\Delta M^{-1}{G}_{}^{PP}(\mathsf{\epsilon })+G_0^{uu}(\mathsf{\epsilon })\times \\ \times \left(\Delta \Phi +{\Sigma }_{phe}(\mathsf{\epsilon })+{\Sigma }_{phph}(\mathsf{\epsilon })\right) G_{}^{uP}(\mathsf{\epsilon }),\end{array} \\ \begin{array}{l}{G}^{Pu}(\mathsf{\epsilon })=G_0^{Pu}(\mathsf{\epsilon })+G_0^{Pu}(\mathsf{\epsilon }) \left(\Delta \Phi +{\Sigma }_{phe}(\mathsf{\epsilon })+{\Sigma }_{phph}(\mathsf{\epsilon })\right)\times \\ \times G_{}^{uu}(\mathsf{\epsilon })+G_0^{PP}(\mathsf{\epsilon })\Delta M^{-1}{G}_{}^{Pu}(\mathsf{\epsilon }),\end{array} \end{gathered}$$

(54)

where $\mathsf{\epsilon }=\hslash \mathsf{\omega }$. Here, $G^{a{a}^{+}}(\mathsf{\epsilon })$, $G^{uu}(\mathsf{\epsilon })$, $G^{PP}(\mathsf{\epsilon })$, $G^{uP}(\mathsf{\epsilon })$, $G^{Pu}(\mathsf{\epsilon })$ are the real-frequency representation of the single-particle Green’s function of the electrons, the coordinate-coordinate, momentum–momentum, coordinate–momentum, and momentum–coordinate Green’s functions of the phonons, respectively; and ${\Sigma }_{eph}(\mathsf{\epsilon })$, ${\Sigma }_{phe}(\mathsf{\epsilon })$, ${\Sigma }_{ee}(\mathsf{\epsilon })$, ${\Sigma }_{phph}(\mathsf{\epsilon })$ are the corresponding self-energies (mass operators) for the electron–phonon, phonon–electron, electron–electron, and phonon–phonon interactions.

The real-time and real-frequency Green’s functions are related by standard Fourier transform relations given by

$$G_{r,a}^{AB}(t)=\frac{1}{2\pi }{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{G}_{r,a}^{AB}(\mathsf{\omega })} e^{-i\mathsf{\omega }t}d\mathsf{\omega }$$

(55)

and

$$G_{r,a}^{AB}(\mathsf{\omega })={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{G}_{r,a}^{AB}(t)} e^{i\mathsf{\omega }t}dt.$$

(56)

The thermal Green’s functions are periodic (bosons) or antiperiodic (fermions) on the interval −1/Θ ≤ τ < 1/Θ, and hence have a Fourier series representation in terms of their Matsubara frequencies, as follows:

$$G^{AB}(\mathsf{\tau })=\Theta {\displaystyle \sum _{{\mathsf{\omega }}_n}{G}^{AB}({\mathsf{\omega }}_n) e^{-i{\mathsf{\omega }}_n\mathsf{\tau }}}$$

(57)

and

$$G^{AB}({\mathsf{\omega }}_n)=\frac{1}{2}{\displaystyle \underset{-1/\Theta }{\overset{1/\Theta }{\int }}{G}^{AB}(\mathsf{\tau })} e^{i{\mathsf{\omega }}_n\mathsf{\tau }}d\mathsf{\tau },$$

(58)

where the Matsubara frequencies satisfy

$$\begin{array}{l}{\mathsf{\omega }}_n=\{\begin{array}{cc}\hfill 2n\pi \Theta \hfill & \mathrm{for} \mathrm{Bose} \mathrm{particles},\hfill \\ \hfill (2n+1)\pi \Theta \hfill & \mathrm{for} \mathrm{Fermi} \mathrm{particles},\hfill \end{array} \\ n=0, \pm 1, \pm 2,\dots \end{array}$$

(59)

The electronic Green’s functions are infinite matrices with indices given by the lattice site n, the basis site i, and the other quantum numbers γ. Similarly, the phonon Green’s functions also are infinite matrices with the same lattice and basis site dependence, plus a dependence on the spatial coordinate direction α. Using the equations of motion for Green’s functions, one can obtain simple expressions for the zero-order Green’s functions, namely [26]:

$$G_0^{a{a}^{+}}(\mathsf{\epsilon })={[\mathsf{\epsilon }-H_0^{(1)}]}^{-1},$$

(60)

with

$$H_0^{(1)}=\Vert h_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{(0)}\Vert ,$$

(61)

$$G_0^{uu}(\mathsf{\epsilon })={\left[{\mathsf{\omega }}^2{\mathrm{M}}^{(0)}-{\Phi }^{(0)}\right]}^{-1},$$

(62)

with

$${\Phi }^{(0)}=\Vert {\Phi }^{(0)}{}_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}\Vert$$

(63)

and

$$M^{(0)}=\Vert M_i{\mathsf{\delta }}_{n{n}^{\prime }}{\mathsf{\delta }}_{i{i}^{\prime }}{\mathsf{\delta }}_{\mathsf{\alpha }{\mathsf{\alpha }}^{\prime }}\Vert .$$

(64)

Here, the double lines denote a matrix.

When the perturbations are small, given by

$$\frac{{\left(\frac{{\mathsf{\epsilon }}^2}{{\hslash }^2}\Delta M+\Delta \Phi +{\Sigma }_{phe}(\mathsf{\epsilon })+{\Sigma }_{phph}(\mathsf{\epsilon })\right)}_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}}{{\Phi }_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}^{(0)}}<<1,$$

(65)

then the solution of the system of equations in Equation (55) becomes

$$G^{a{a}^{+}}(\mathsf{\epsilon })={\left[\mathsf{\epsilon }-H_0^{\left(1\right)}-\left(w+{\Sigma }_{eph}(\mathsf{\epsilon })+{\Sigma }_{ee}(\mathsf{\epsilon })\right)\right]}^{-1},$$

(66)

$$G^{uu}(\mathsf{\epsilon })={\left[{\mathsf{\omega }}^2{M}^{\left(0\right)}-{\Phi }^{\left(0\right)}-\left(\frac{{\mathsf{\epsilon }}^2}{{\hslash }^2}\Delta M+\Delta \Phi +{\Sigma }_{phe}(\mathsf{\epsilon })+{\Sigma }_{phph}(\mathsf{\epsilon })\right)\right]}^{-1},$$

(67)

$$G_{}^{PP}(\mathsf{\epsilon })=\frac{{\mathsf{\epsilon }}^2}{{\hslash }^2}{(M^{(0)})}^2 G_{}^{uu}(\mathsf{\epsilon }),$$

(68)

where

$$\Delta M=\Vert (M_i-M_{ni}){\mathsf{\delta }}_{n{n}^{\prime }}{\mathsf{\delta }}_{i{i}^{\prime }}{\mathsf{\delta }}_{\mathsf{\alpha }{\mathsf{\alpha }}^{\prime }}\Vert ,\mathsf{\epsilon }=\hslash \mathsf{\omega }.$$

(69)

The mass operator of the Green’s function of electrons for the electron–phonon interaction ${\Sigma }_{eph}(\mathsf{\tau },{\mathsf{\tau }}^{\prime })$ is described by the diagram in Figure 1. The mass operator of the Green’s function of electrons for the electron–phonon interaction ${\Sigma }_{eph}(\mathsf{\tau },{\mathsf{\tau }}^{\prime })$ is described by the diagram in [26,29,30].

Solid lines in Figure 1 correspond to the Green’s function of electrons $G_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{a{a}^{+}}(\mathsf{\tau },{\mathsf{\tau }}^{\prime })$ and dashed lines correspond to the Green’s function of phonons $G_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}^{uu}(\mathsf{\tau },{\mathsf{\tau }}^{\prime })$. The vertex part ${\Gamma }_{ni\mathsf{\gamma },n_1{i}_1{\mathsf{\gamma }}_1}^{n_2{i}_2{\mathsf{\alpha }}_2}({\mathsf{\tau }}_2,\mathsf{\tau },{\mathsf{\tau }}_1)$ of the mass operator of the Green’s function is described by the diagrams in Figure 2.

The not shaded triangle in Figure 2 corresponds to equation

$${\Gamma }_{0 ni\mathsf{\gamma },n_1{i}_1{\mathsf{\gamma }}_1}^{n_2{i}_2{\mathsf{\alpha }}_2}({\mathsf{\tau }}_2,\mathsf{\tau },{\mathsf{\tau }}_1)=v^{\prime }{}_{ni\mathsf{\gamma },n_1{i}_1{\mathsf{\gamma }}_1}^{n_2{i}_2{\mathsf{\alpha }}_2}\mathsf{\delta }(\mathsf{\tau }-{\mathsf{\tau }}_2) \mathsf{\delta }(\mathsf{\tau }-{\mathsf{\tau }}_1).$$

(70)

In Figure 1 and Figure 2, summation for internal points $\tilde{n}\mathsf{\gamma }$ is carried out. Summation of $\tilde{n}\mathsf{\gamma }$ provides summation of $ni\mathsf{\gamma }$ and integration over $\mathsf{\tau }$. Expressions that correspond to each diagram are attributed to multiplier ${(-1)}^{n+F}$, where $n$ is the diagram’s order (namely the number of vertices ${\Gamma }_0$ in the diagram), and $F$ is the number of lines for the Green’s function of electrons $G_{}^{a{a}^{+}}$. This function goes out and goes in in the same vertices.

Explicitly, the electron–phonon self-energy becomes

$$\begin{array}{l}{\Sigma }_{\mathrm{e}ph ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}(\mathsf{\epsilon })=-\frac{1}{4\pi i}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{\mathsf{\epsilon }}^{\prime }\cot \mathrm{h}\left(\frac{{\mathsf{\epsilon }}^{\prime }}{2\Theta }\right) {\Gamma }_{ ni\mathsf{\gamma }, n_3{i}_3{\mathsf{\gamma }}_3}^{\left(0\right)n_1{i}_1{\mathsf{\alpha }}_1}\times }\\ \times \left[G_{n_1{i}_1{\mathsf{\alpha }}_1,n_2{i}_2{\mathsf{\alpha }}_2}^{uu}({\mathsf{\epsilon }}^{\prime })-G_{n_1{i}_1{\mathsf{\alpha }}_1,n_2{i}_2{\mathsf{\alpha }}_2}^{uu *}({\mathsf{\epsilon }}^{\prime })\right] G_{n_3{i}_3{\mathsf{\gamma }}_3,n_4{i}_4{\mathsf{\gamma }}_4}^{a{a}^{+}}\times \\ \times (\mathsf{\epsilon }-{\mathsf{\epsilon }}^{\prime }){\Gamma }_{ n_4{i}_4{\mathsf{\gamma }}_4,n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{n_2{i}_2{\mathsf{\alpha }}_2},\end{array}$$

(71)

$${\Gamma }_{ ni\mathsf{\gamma }, n_3{i}_3{\mathsf{\gamma }}_3}^{\left(0\right)n_1{i}_1{\mathsf{\alpha }}_1}=v^{\prime }{}_{ni\mathsf{\gamma }, n_3{i}_3{\mathsf{\gamma }}_3}^{n_1{i}_1{\mathsf{\alpha }}_1}.$$

(72)

where repeated indices are summed over.

Phonon–electron interaction is described by the diagram in Figure 3. Phonon–electron interaction is described by the diagram in [29,30].

The designation in Figure 3 corresponds to designations in Figure 1 and Figure 2.

The self-energy of the phonon due to the phonon–electron interaction is given by

$$\begin{array}{l}{\Sigma }_{ph\mathrm{e} ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}(\mathsf{\epsilon })=\frac{1}{2\pi i}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{\mathsf{\epsilon }}^{\prime } f({\mathsf{\epsilon }}^{\prime }){\Gamma }_{n_2{i}_2{\mathsf{\gamma }}_2,n_1{i}_1{\mathsf{\gamma }}_1}^{ \left(0\right)ni\mathsf{\alpha }}}\times \\ \times \{\left[G_{n_1{i}_1{\mathsf{\gamma }}_1,n_3{i}_3{\mathsf{\gamma }}_3}^{a{a}^{+}}(\mathsf{\epsilon }+{\mathsf{\epsilon }}^{\prime })-G_{n_1{i}_1{\mathsf{\gamma }}_1,n_3{i}_3{\mathsf{\gamma }}_3}^{a{a}^{+}*}(\mathsf{\epsilon }+{\mathsf{\epsilon }}^{\prime })\right]\times \\ \times G_{n_4{i}_4{\mathsf{\gamma }}_4,n_2{i}_2{\mathsf{\gamma }}_2}^{a{a}^{+}*}({\mathsf{\epsilon }}^{\prime })+G_{n_1{i}_1{\mathsf{\gamma }}_1,n_3{i}_3{\mathsf{\gamma }}_3}^{a{a}^{+}}(\mathsf{\epsilon }+{\mathsf{\epsilon }}^{\prime })\times \\ \times \left[G_{n_4{i}_4{\mathsf{\gamma }}_4,n_2{i}_2{\mathsf{\gamma }}_2}^{a{a}^{+}}({\mathsf{\epsilon }}^{\prime })-G_{n_4{i}_4{\mathsf{\gamma }}_4,n_2{i}_2{\mathsf{\gamma }}_2}^{a{a}^{+}*}({\mathsf{\epsilon }}^{\prime })\right]\}\times \\ \times {\Gamma }_{ n_3{i}_3{\mathsf{\gamma }}_3,n_4{i}_4{\mathsf{\gamma }}_4}^{ n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}.\end{array}$$

(73)

where $f(\mathsf{\epsilon })$ is the so-called Fermi–Dirac distribution function.

Diagrams for the mass operator ${\Sigma }_{ee}(\mathsf{\tau },{\mathsf{\tau }}^{\prime })$ that describe the electron–electron interaction, are shown in Figure 4. Diagrams for the mass operator ${\Sigma }_{ee}(\mathsf{\tau },{\mathsf{\tau }}^{\prime })$ that describe electron–electron interaction, are shown in [26,29,30].

The vertex parts ${\Gamma }_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{n_2{i}_2{\mathsf{\gamma }}_2,n_1{i}_1{\mathsf{\gamma }}_1}({\mathsf{\tau }}_2,{\mathsf{\tau }}_1\mathsf{\tau },{\mathsf{\tau }}^{\prime })$ are shown in diagrams in Figure 5.

Not shaded triangle in Figure 5 corresponds to equation

$${\Gamma }_{0 ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{n_2{i}_2{\mathsf{\gamma }}_2,n_1{i}_1{\mathsf{\gamma }}_1}({\mathsf{\tau }}_2,{\mathsf{\tau }}_1\mathsf{\tau },{\mathsf{\tau }}^{\prime })={\tilde{v}}_{n_1{i}_1{\mathsf{\gamma }}_1,n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{(2)ni\mathsf{\gamma },n_2{i}_2{\mathsf{\gamma }}_2}\times \mathsf{\delta }(\mathsf{\tau }-{\mathsf{\tau }}_2)\mathsf{\delta }(\mathsf{\tau }-{\mathsf{\tau }}_1)\mathsf{\delta }(\mathsf{\tau }-{\mathsf{\tau }}^{\prime }),$$

(74)

$${\tilde{v}}_{n_1{i}_1{\mathsf{\gamma }}_1,n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{(2)ni\mathsf{\gamma },n_2{i}_2{\mathsf{\gamma }}_2}=v_{n_1{i}_1{\mathsf{\gamma }}_1,n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{(2)ni\mathsf{\gamma },n_2{i}_2{\mathsf{\gamma }}_2}-v_{n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime },n_1{i}_1{\mathsf{\gamma }}_1}^{(2)ni\mathsf{\gamma },n_2{i}_2{\mathsf{\gamma }}_2}$$

(75)

The mass operator that describes the electron–electron interaction is:

$${\Sigma }_{ee ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}(\mathsf{\epsilon })={\Sigma }_{ee ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{(1)}+{\Sigma }_{ee ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{(2)}(\mathsf{\epsilon }),$$

(76)

$${\Sigma }_{\mathrm{ee} n,n^{\prime }}^{(1)}=-\frac{1}{2\pi i}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{\mathsf{\epsilon }}^{\prime } f\left({\mathsf{\epsilon }}^{\prime }\right) {\Gamma }_{ n, n^{\prime }}^{ \left(0\right) n_2,n_1}}\left[G_{n_1,n_2}^{a{a}^{+}}({\mathsf{\epsilon }}^{\prime })-G_{n_1,n_2}^{a{a}^{+} *}({\mathsf{\epsilon }}^{\prime })\right],$$

(77)

$$\begin{gathered} {\Sigma }_{\mathrm{ee}\hspace{0.17em} n,n^{\prime }}^{(2)}\left(\mathsf{\epsilon }\right)={\left(\frac{1}{2\pi i}\right)}^2{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{\mathsf{\epsilon }}_1}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{\mathsf{\epsilon }}_2 f\left({\mathsf{\epsilon }}_1\right)f\left({\mathsf{\epsilon }}_2\right){\Gamma }_{ n_2,n_1}^{ \left(0\right) n,n_3} } \\ \begin{array}{l}\times [G_{n_2,n_5}^{a{a}^{+} }(\mathsf{\epsilon }-{\mathsf{\epsilon }}_1-{\mathsf{\epsilon }}_2)G_{n_1,n_4}^{a{a}^{+}\ast }({\mathsf{\epsilon }}_1)-\\ -G_{n_2,n_5}^{a{a}^{+}\ast }(\mathsf{\epsilon }-{\mathsf{\epsilon }}_1-{\mathsf{\epsilon }}_2)G_{n_1,n_4}^{a{a}^{+} }({\mathsf{\epsilon }}_1)]\times \end{array} \\ \begin{array}{l}\times \left[G_{n_6,n_3}^{a{a}^{+}}\left({\mathsf{\epsilon }}_2\right)-G_{n_6,n_3}^{a{a}^{+} *}\left({\mathsf{\epsilon }}_2\right)\right]\\ -\left[G_{n_2,n_5}^{a{a}^{+}}\left(\mathsf{\epsilon }-{\mathsf{\epsilon }}_1-{\mathsf{\epsilon }}_2\right)-G_{n_2,n_5}^{a{a}^{+} *}\left(\mathsf{\epsilon }-{\mathsf{\epsilon }}_1-{\mathsf{\epsilon }}_2\right)\right]\times \end{array} \\ \begin{array}{l}\times [G_{n_1,n_4}^{a{a}^{+}}\left({\mathsf{\epsilon }}_1\right)G_{n_6,n_3}^{a{a}^{+}}\left({\mathsf{\epsilon }}_2\right)\\ -G_{n_1,n_4}^{a{a}^{+} *}\left({\mathsf{\epsilon }}_1\right)G_{n_6,n_3}^{a{a}^{+} *}\left({\mathsf{\epsilon }}_2\right)]\}{\Gamma }_{ n_4, n^{\prime }}^{ n_5,n_6},\end{array} \end{gathered}$$

(78)

$${\Gamma }_{n_2,n_1}^{\left(0\right)n,n_3}={\tilde{v}}_{n_2,n_1}^{(2)n,n_3}=v_{n_2,n_1}^{(2)n,n_3}-v_{n_1,n_2}^{(2)n,n_3}n=ni\mathsf{\gamma }$$

(79)

A similar result for the contribution to the phonon self-energy ${\Sigma }_{phph}(\mathsf{\epsilon })$ from phonon–phonon coupling is given in [26].

In deriving the expressions in Equations (72), (74), (78), and (79), we employed the standard resummation techniques for any function $\mathsf{\varphi }(z)$ that is analytic in the region covered by contour C, which encloses all of the Matsubara frequencies. Namely, we have

$$\begin{array}{l}\Theta {\displaystyle \sum _{{\mathsf{\omega }}_n}\mathsf{\varphi }(i{\mathsf{\omega }}_n)}=\frac{1}{4\pi i}{\displaystyle \underset{C}{\oint }dz \coth \left(\frac{z}{2\Theta }\right)\mathsf{\varphi }\left(z\right)} \\ ({\mathsf{\omega }}_n=2\mathrm{n}\pi \Theta )\end{array},$$

(80)

for the Bosonic case, and

$$\begin{array}{l}\Theta {\displaystyle \sum _{{\mathsf{\omega }}_n}\mathsf{\varphi }(i{\mathsf{\omega }}_n)}=-\frac{1}{2\pi i}{\displaystyle \underset{C}{\oint }dz \tilde{f}\left(\frac{z}{\Theta }\right)\mathsf{\varphi }\left(z\right)} \\ ({\mathsf{\omega }}_n=\left(2\mathrm{n}+1\right)\pi \Theta )\end{array},$$

(81)

for the Fermionic case, with

$$\tilde{f}\left(\frac{z}{\Theta }\right)=\frac{1}{\exp \left(\frac{z}{\Theta }\right)+1}.$$

(82)

We comment that for the many-body Green’s functions described here, it is customary to have the chemical potential situated at zero frequency, as dine here.

In general, the renormalization of the vertex of the mass operator of the Green’s functions in expressions (72), (74), and (79) can be performed using Figure 2 and Figure 5. The diagrams in Figure 2 and Figure 5 correspond to the equation

$$\begin{array}{l}{\Gamma }_{ n_3{i}_3{\mathsf{\gamma }}_3,n_4{i}_4{\mathsf{\gamma }}_4}^{n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}={\Gamma }_{ n_3{i}_3{\mathsf{\gamma }}_3,n_4{i}_4{\mathsf{\gamma }}_4}^{\left(0\right)n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}-\frac{1}{2\pi i}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d\mathsf{\epsilon }f\left(\mathsf{\epsilon }\right)}{\Gamma }_{ n_5{i}_5{\mathsf{\gamma }}_5,n_6{i}_6{\mathsf{\gamma }}_6}^{\left(0\right)n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}\\ \times \left[G_{n_6{i}_6{\mathsf{\gamma }}_6,n_7{i}_7{\mathsf{\gamma }}_7}^{a{a}^{+}}\left(\mathsf{\epsilon }\right)G_{n_8{i}_8{\mathsf{\gamma }}_8,n_5{i}_5{\mathsf{\gamma }}_5}^{a{a}^{+}}\left(\mathsf{\epsilon }\right)-G_{n_6{i}_6{\mathsf{\gamma }}_6,n_7{i}_7{\mathsf{\gamma }}_7}^{a{a}^{+}\ast }\left(\mathsf{\epsilon }\right)G_{n_8{i}_8{\mathsf{\gamma }}_8,n_5{i}_5{\mathsf{\gamma }}_5}^{a{a}^{+}\ast }\left(\mathsf{\epsilon }\right)\right]\\ \times {\Gamma }_{ n_7{i}_7{\mathsf{\gamma }}_7,n_8{i}_8{\mathsf{\gamma }}_8}^{\left(0\right)n_9{i}_9{\mathsf{\alpha }}_9}{G}_{n_9{i}_9{\mathsf{\alpha }}_9,n_{10}{i}_{10}{\mathsf{\alpha }}_{10}}^{uu}\left(0\right){\Gamma }_{n_3{i}_3{\mathsf{\gamma }}_3,n_4{i}_4{\mathsf{\gamma }}_4}^{n_{10}{i}_{10}{\mathsf{\alpha }}_{10}}\end{array}$$

(83)

and

$$\begin{array}{l}{\Gamma }_{ n_4, n^{\prime }}^{ n_5,n_6}={\Gamma }_{ n_4, n^{\prime }}^{\left(0\right) n_5,n_6}-\frac{1}{2\pi i}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d\mathsf{\epsilon }f\left(\mathsf{\epsilon }\right)}{\Gamma }_{ n_4, n_8}^{\left(0\right) n_5,n_7}\\ \times \left[G_{n_7,n_9}^{a{a}^{+}}\left(\mathsf{\epsilon }\right)G_{n_8,n_{10}}^{a{a}^{+}\ast }\left(\mathsf{\epsilon }\right)-G_{n_7,n_9}^{a{a}^{+}\ast }\left(\mathsf{\epsilon }\right)G_{n_8,n_{10}}^{a{a}^{+}}\left(\mathsf{\epsilon }\right)\right]\\ \times {\Gamma }_{ n_{10}, n^{\prime }}^{ n_9,n_6},n=ni\mathsf{\gamma }.\end{array}$$

(84)

Summation is implied over repeated indices in expressions (84) and (85).

The Fermi level ${\mathsf{\epsilon }}_F\equiv {\mathsf{\mu }}_{\mathrm{e}}$ of the system is determined by equation:

$$<Z>={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}f(\mathsf{\epsilon }) g_e(\mathsf{\epsilon }) d\mathsf{\epsilon }},$$

(85)

$$f(\mathsf{\epsilon })=\frac{1}{\exp (\frac{\mathsf{\epsilon }-{\mathsf{\epsilon }}_F}{\Theta } )+1}$$

(86)

where $<Z>$ is the average number of electrons per atom and $g_e(\mathsf{\epsilon })$ is the many-body electronic density of states, which satisfies

$$g_e\left(\mathsf{\epsilon }\right)=-\frac{1}{\pi \mathsf{\nu }N} \mathrm{Im} Tr{\langle G^{a{a}^{+}}\left(\mathsf{\epsilon }\right)\rangle }_c$$

(87)

Here, ${\langle \dots \rangle }_c$ denotes configurational averaging over the disorder, N is the number of primitive lattice cells, and $\mathsf{\nu }$ is the number of atoms per primitive cell. We drop the letter c on the configurational averaging for simplicity. In Equation (86), $\langle Z\rangle$ is the average number of electrons per atom.

It should be noted that the first term in the electron self-energy due to electron–electron interactions, ${\Sigma }_{ee ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{(1)}$ in Equation (77), describes the Coulomb and exchange electron–electron interactions in the Hartree–Fock approximation. The second term, ${\Sigma }_{ee ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{(2)}(\mathsf{\epsilon })$, which is caused by corrections beyond Hartree–Fock, describes the effects of electron correlations. As opposed to the procedures used in [13,14,15,16,17,18,19], the long-range Coulomb interaction of electrons located at different lattice sites of the crystal is described by taking into account an arbitrary number of energy bands.

The expression for the Green’s function in Equations (67) and (68) differs from the corresponding expressions for the Green’s function of a single-particle Hamiltonian of a disordered system, only from the different self-energy contributions. Hence, we solve for the Green’s function using the well-known methods of disordered systems theory [31,32].

To find the density of states by Formula (88), it is necessary to find the average values of the Green’s functions defined by Formulas (67) and (68).

4. Localized Magnetic Moments

As we will be working with magnetic moments in the remainder of the paper, we now slightly modify our notation so that the symbol $\mathsf{\gamma }=\mathsf{\delta }\mathsf{\sigma }=\tilde{\mathsf{\epsilon }}lm\mathsf{\sigma }$ now refers to all other quantum numbers, except for spin, and we introduce the spin quantum number σ explicitly in all of the following equations. The electron–electron self-energy in Equation (67) requires the occupation number $Z_{ni\mathsf{\delta }\mathsf{\sigma }}^{ \mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}$ of the different electronic states $(ni\mathsf{\delta }\mathsf{\sigma })$, where we are explicitly including the dependence on σ. The explicit values for $Z_{ni\mathsf{\delta }\mathsf{\sigma }}^{ \mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}$ are calculated from Equation (86), where the total electronic density of states $g_{\mathrm{e}}(\mathsf{\epsilon })$ is replaced by the partial density of states $g_{ni\mathsf{\delta }\mathsf{\sigma }}^{\mathsf{\lambda } m_{\mathsf{\lambda }\mathrm{i}}}(\mathsf{\epsilon })$ for energy band $\mathsf{\delta }$ and spin projection $\mathsf{\sigma }$, to allow for magnetic solutions. The occupation numbers $Z_{ni\mathsf{\delta }\mathsf{\sigma }}^{ \mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}$ and the partial density of states $g_{ni\mathsf{\delta }\mathsf{\sigma }}^{\mathsf{\lambda } m_{\mathsf{\lambda }\mathrm{i}}}(\mathsf{\epsilon })$ then satisfy

$$Z_{ni\mathsf{\delta }\mathsf{\sigma }}^{\begin{array}{l} \\ \mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}\end{array}}={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}f(\mathsf{\epsilon }) g_{ni\mathsf{\delta }\mathsf{\sigma }}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}(\mathsf{\epsilon }) d\mathsf{\epsilon }},$$

(88)

$$g_{ni\mathsf{\delta }\mathsf{\sigma }}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}(\mathsf{\epsilon })=-\frac{1}{\pi }{\mathrm{Im}\langle G_{ni\mathsf{\delta }\mathsf{\sigma },ni\mathsf{\delta }\mathsf{\sigma }}^{a{a}^{+}}(\mathsf{\epsilon })\rangle |}_{(ni)\in \mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}$$

(89)

Note that disorder averaging is done under the assumption that an atom of type λ is located at the site $(ni)$, and its projection of the localized magnetic moment onto the z-axis is equal to $m_{\mathsf{\lambda }i}$. The probability of this configuration is $P_{ni}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}$, and we have the obvious constraint that

$${\displaystyle \sum _{\mathsf{\lambda },m_{\mathsf{\lambda }\mathrm{i}}}{P}_{ni}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}=1}.$$

(90)

The total charge and magnetization for each orbital on a site are given by

$$Z_{ni\mathsf{\delta }}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}=Z_{ni\mathsf{\delta }\mathsf{\sigma }}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}+Z_{ni\mathsf{\delta },-\mathsf{\sigma }}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}, m_{\mathsf{\lambda }i\mathsf{\delta }}=Z_{ni\mathsf{\delta }\mathsf{\sigma }}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}-Z_{ni\mathsf{\delta },-\mathsf{\sigma }}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}},$$

(91)

and by

$$Z_{ni\mathsf{\delta }\mathsf{\sigma }}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}=\frac{Z_{ni\mathsf{\delta }}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}+m_{\mathsf{\lambda }i\mathsf{\delta }}}{2}, Z_{ni\mathsf{\delta },-\mathsf{\sigma }}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}=\frac{Z_{ni\mathsf{\delta }}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}-m_{\mathsf{\lambda }i\mathsf{\delta }}}{2},$$

(92)

respectively. We need to sum over all other quantum numbers to get the totals:

$$Z_{ni}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}={\displaystyle \sum _{\mathsf{\delta }}{Z}_{ni\mathsf{\delta }}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}}, m_{\mathsf{\lambda }i}={\displaystyle \sum _{\mathsf{\delta }}{m}_{\mathsf{\lambda }i\mathsf{\delta }}}.$$

(93)

5. Density of Electronic and Phononic States

In Equations (67) and (68), by introducing the mass operator as the sum of one site operators and selecting it as a zero approximation, the effective medium Green’s function cluster expansion for Green’s functions $G^{a{a}^{+}}(\mathsf{\epsilon })$, $G^{uu}(\mathsf{\epsilon })$, is performed. The specified expansion is a generalization of the cluster expansion for the Green’s function $G^{a{a}^{+}}(\mathsf{\epsilon })$ of single-particle Hamiltonian. Green’s functions of the effective environment are defined by the expressions:

$${\tilde{G}}^{a{a}^{+}}(\mathsf{\epsilon })={\left[\mathsf{\epsilon }-h^{(0)}-{\tilde{\Sigma }}_{eph}(\mathsf{\epsilon })-{\tilde{\Sigma }}_{ee}(\mathsf{\epsilon })-{\mathsf{\sigma }}_e\left(\mathsf{\epsilon }\right)\right]}^{-1},$$

(94)

$${\tilde{G}}^{uu}(\mathsf{\epsilon })={\left[\frac{{\mathsf{\epsilon }}^2}{{\hslash }^2}{M}^{\left(0\right)}-{\Phi }^{\left(0\right)}-{\tilde{\Sigma }}_{phe}(\mathsf{\epsilon })-{\mathsf{\sigma }}_{ph}\left(\mathsf{\epsilon }\right)\right]}^{-1}.$$

(95)

Expressions for the operators ${\tilde{\Sigma }}_{eph}(\mathsf{\epsilon })$, ${\tilde{\Sigma }}_{phe}(\mathsf{\epsilon })$, and ${\tilde{\Sigma }}_{ee}(\mathsf{\epsilon })$ are obtained from the expressions for the mass operators ${\Sigma }_{eph}(\mathsf{\epsilon })$, ${\Sigma }_{phe}(\mathsf{\epsilon })$, and ${\Sigma }_{ee}(\mathsf{\epsilon })$ (72)–(80) by replacing the vertex parts ${\Gamma }_{ ni\mathsf{\gamma }, n_3{i}_3{\mathsf{\gamma }}_3}^{\left(0\right)n_1{i}_1{\mathsf{\alpha }}_1}$, ${\Gamma }_{n_2,n_1}^{\left(0\right)n,n_3}$, $n\equiv ni\mathsf{\gamma }$ by their expressions for ideally ordered crystals, and replacing the Green’s functions $G^{a{a}^{+}}(\mathsf{\epsilon })$, $G^{uu}(\mathsf{\epsilon })$ with the Green’s functions of the effective medium ${\tilde{G}}^{a{a}^{+}}(\mathsf{\epsilon })$, ${\tilde{G}}^{uu}(\mathsf{\epsilon })$. Expressions for operators ${\mathsf{\sigma }}_e\left(\mathsf{\epsilon }\right)$ and ${\mathsf{\sigma }}_{ph}\left(\mathsf{\epsilon }\right)$ in Formulas (95) and (96) will be given below.

The Green’s functions in Equations (67) and (68) satisfy a Dyson equation that can be expressed in terms of a T-matrix via:

$$G(\mathsf{\epsilon })=\tilde{G}(\mathsf{\epsilon })+\tilde{G}(\mathsf{\epsilon }) T(\mathsf{\epsilon }) \tilde{G}(\mathsf{\epsilon }),$$

(96)

where the T-matrix T is represented by a series, in which each term describes the scattering of clusters with different numbers of nodes expressed schematically by

$$T={\displaystyle \sum _{(n_1{i}_1)}{t}^{n_1{i}_1}}+{\displaystyle \sum _{(n_1{i}_1)\ne (n_2{i}_2)}^{}{T}^{(2) n_1{i}_1,n_2{i}_2}}+\dots$$

(97)

Here, we have

$$T^{(2) n_1{i}_1,n_2{i}_2}={\left[I-t^{n_1{i}_1}\tilde{G}{t}^{n_2{i}_2}\tilde{G}\right]}^{-1}{t}^{n_1{i}_1}\tilde{G}{t}^{n_2{i}_2}\left[I+\tilde{G}{t}^{n_1{i}_1}\right],$$

(98)

where $t^{n_1{i}_1}$ is the on-site scattering operator, which is given by

$$t^{n_1{i}_1}={\left[I-({\Sigma }^{n_1{i}_1}-{\mathsf{\sigma }}^{n_1{i}_1})\tilde{G}\right]}^{-1}({\Sigma }^{n_1{i}_1}-{\mathsf{\sigma }}^{n_1{i}_1}).$$

(99)

The self-energy employed in Equation (67), ${\Sigma }_e^{n_1{i}_1}(\mathsf{\epsilon })$, satisfies

$$w+{\Sigma }_{eph}(\mathsf{\epsilon })+{\Sigma }_{ee}(\mathsf{\epsilon })-{\tilde{\Sigma }}_{eph}(\mathsf{\epsilon })-{\tilde{\Sigma }}_{ee}(\mathsf{\epsilon })={\displaystyle \sum _{\left(n_1{i}_1\right)}{\Sigma }_e^{n_1{i}_1}(\mathsf{\epsilon })},$$

(100)

for the electrons. For the phonons, we have

$$\begin{array}{l}\frac{{\mathsf{\epsilon }}^2}{{\hslash }^2}\Delta M+\Delta \Phi +{\Sigma }_{phe}(\mathsf{\epsilon })+{\Sigma }_{phph}(\mathsf{\epsilon })-\\ -{\tilde{\Sigma }}_{phe}(\mathsf{\epsilon })-{\tilde{\Sigma }}_{phph}(\mathsf{\epsilon })={\displaystyle \sum _{\left(n_1{i}_1\right)}{\Sigma }_{ph}^{n_1{i}_1}(\mathsf{\epsilon })}\end{array}$$

(101)

Using Equations (72), (77)–(79), and (101), we obtain the expression for the intrinsic energy part ${\Sigma }_{eni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{\mathsf{\lambda }{n}_1{i}_1}(\mathsf{\epsilon })$, which describes the scattering of electrons:

$$\begin{array}{l}{\Sigma }_{e{n}_1{i}_1{\mathsf{\gamma }}_1,n_2{i}_2{\mathsf{\gamma }}_2}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}ni}(\mathsf{\epsilon })=w_{n_1{i}_1{\mathsf{\gamma }}_1,n_2{i}_2{\mathsf{\gamma }}_2}^{\mathsf{\lambda }ni}+\\ +{\displaystyle \sum _{\begin{array}{l}{n}_3{i}_3{\mathsf{\gamma }}_3\\ n_4{i}_4{\mathsf{\gamma }}_4\end{array}}{{\displaystyle \tilde{v}}}_{n_3{i}_3{\mathsf{\gamma }}_3,n_2{i}_2{\mathsf{\gamma }}_2}^{(2)n_1{i}_1{\mathsf{\gamma }}_1,n_4{i}_4{\mathsf{\gamma }}_4}\left(Z_{n_3{i}_3{\mathsf{\gamma }}_3,n_4{i}_4{\mathsf{\gamma }}_4}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}ni}-{\tilde{Z}}_{n_3{i}_3{\mathsf{\gamma }}_3,n_4{i}_4{\mathsf{\gamma }}_4}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}ni}\right),}\end{array}$$

(102)

where

$$Z_{n_3{i}_3{\mathsf{\gamma }}_3,n_4{i}_4{\mathsf{\gamma }}_4}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}ni}=-\frac{1}{\pi }{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}f(\mathsf{\epsilon },{\mathsf{\epsilon }}_F) {\mathrm{Im}\langle G_{n_3{i}_3{\mathsf{\gamma }}_3,n_4{i}_4{\mathsf{\gamma }}_4}^{a{a}^{+}}(\mathsf{\epsilon })\rangle |}_{(ni)\in \mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}} d\mathsf{\epsilon }}.$$

(103)

The value of ${\tilde{Z}}_{n_3{i}_3{\mathsf{\gamma }}_3,n_4{i}_4{\mathsf{\gamma }}_4}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}ni}$ in Equation (103) is derived from Equation (104) by replacing the full Green’s function by the effective medium Green’s function. The diagonal elements of the matrix $Z_{n_3{i}_3{\mathsf{\gamma }}_3,n_4{i}_4{\mathsf{\gamma }}_4}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}ni}$ in Equation (104) are equal to the occupation numbers of the electron states $Z_{ni\mathsf{\delta }\mathsf{\sigma }}^{ \mathsf{\lambda }{m}_{\mathsf{\lambda }\mathrm{i}}}$ in Equation (89).

Using Equations (10), (70), (74), and (102), we obtain the expression for the intrinsic energy part ${\Sigma }_{phni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}^{\mathsf{\lambda }{n}_1{i}_1}(\mathsf{\epsilon })$, which describes the scattering of phonons:

$${\Sigma }_{phni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}^{\mathsf{\lambda }{n}_1{i}_1}(\mathsf{\epsilon })=\frac{{\mathsf{\epsilon }}^2}{{\hslash }^2}(M_{i_1}-M_{\mathsf{\lambda }}){\mathsf{\delta }}_{n{n}^{\prime }}{\mathsf{\delta }}_{i{i}^{\prime }}{\mathsf{\delta }}_{\mathsf{\alpha }{\mathsf{\alpha }}^{\prime }}$$

(104)

It should be noted that, in the limit of an infinite crystal, on the right-hand side of Equations (103) and (105), the terms inversely proportional to the number of lattice sites are neglected.

We require the fulfillment of the condition

$$\langle t^{0i_1}\rangle =0,$$

(105)

from which follows the system of coupled equations for the operator, in Formulas (95) and (96):

$$\begin{array}{l}{\mathsf{\sigma }}_e^{0i_1}(\mathsf{\epsilon })={\langle {[1-({\Sigma }_e^{0i_1}(\mathsf{\epsilon })-{\mathsf{\sigma }}_e^{0i_1}(\mathsf{\epsilon })){\tilde{G}}_{}^{a{a}^{+}}(\mathsf{\epsilon })]}^{-1}\rangle }^{-1}\times \\ \times \langle {[1-({\Sigma }_e^{0i_1}(\mathsf{\epsilon })-{\mathsf{\sigma }}_e^{0i_1}(\mathsf{\epsilon })){\tilde{G}}_{}^{a{a}^{+}}(\mathsf{\epsilon })]}^{-1}{\Sigma }_e^{0i_1}(\mathsf{\epsilon })\rangle \end{array},$$

(106)

$$\begin{array}{l}{\mathsf{\sigma }}_{ph}^{0i_1}(\mathsf{\epsilon })={\langle {[1-({\Sigma }_{ph}^{0i_1}(\mathsf{\epsilon })-{\mathsf{\sigma }}_{ph}^{0i_1}(\mathsf{\epsilon })){\tilde{G}}_{}^{uu}(\mathsf{\epsilon })]}^{-1}\rangle }^{-1}\times \\ \times \langle {[1-({\Sigma }_{ph}^{0i_1}(\mathsf{\epsilon })-{\mathsf{\sigma }}_{ph}^{0i_1}(\mathsf{\epsilon })){\tilde{G}}_{}^{uu}(\mathsf{\epsilon })]}^{-1}{\Sigma }_{ph}^{0i_1}(\mathsf{\epsilon })\rangle \end{array}.$$

(107)

The matrix elements of the Green’s function of the electron subsystem of the effective medium can be calculated using Fourier transformation

$${\tilde{G}}_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{a{a}^{+}}(\mathsf{\epsilon })=\frac{1}{N}{\displaystyle \sum _{\mathbf{k}}{\tilde{G}}_{i\mathsf{\gamma },i^{\prime }{\mathsf{\gamma }}^{\prime }}^{a{a}^{+}}(\mathbf{k},\mathsf{\epsilon })}{e}^{i\mathbf{k}({\mathbf{r}}_{ni}-{\mathbf{r}}_{n^{\prime }{i}^{\prime }})},$$

(108)

$${\tilde{G}}^{a{a}^{+}}(\mathbf{k},\mathsf{\epsilon })={\left(\mathsf{\epsilon }-\tilde{H}\left(\mathbf{k},\mathsf{\epsilon }\right)\right)}^{-1},$$

(109)

where

$$\tilde{H}\left(\mathbf{k},\mathsf{\epsilon }\right)=h^{(0)}\left(\mathbf{k}\right)+{\tilde{\Sigma }}_{eph}\left(\mathbf{k},\mathsf{\epsilon }\right)+{\tilde{\Sigma }}_{ee}\left(\mathbf{k},\mathsf{\epsilon }\right)+{\mathsf{\sigma }}_e\left(\mathbf{k},\mathsf{\epsilon }\right),$$

(110)

$N$ is the number of primitive unit cells.

We do a similar procedure for the effective medium phonon Green’s function, which satisfies

$${\tilde{G}}_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}^{uu}(\mathsf{\epsilon })=\frac{1}{N}{\displaystyle \sum _{\mathbf{k}}{\tilde{G}}_{i\mathsf{\alpha },i^{\prime }{\mathsf{\alpha }}^{\prime }}^{uu}(\mathbf{k},\mathsf{\epsilon })e^{i\mathbf{k}({\mathbf{r}}_{ni}-{\mathbf{r}}_{n^{\prime }{i}^{\prime }})}},$$

(111)

$${\tilde{G}}^{uu}(\mathbf{k},\mathsf{\epsilon })={\left(\frac{{\mathsf{\epsilon }}^2}{{\hslash }^2}{M}^{(0)}-\tilde{\Phi }\left(\mathbf{k},\mathsf{\epsilon }\right)\right)}^{-1},$$

(112)

where we have

$$\tilde{\Phi }\left(\mathbf{k},\mathsf{\epsilon }\right)={\Phi }^{(0)}\left(\mathbf{k}\right)+{\tilde{\Sigma }}_{phe}\left(\mathbf{k},\mathsf{\epsilon }\right)+{\tilde{\Sigma }}_{phph}\left(\mathbf{k},\mathsf{\epsilon }\right)+{\mathsf{\sigma }}_{ph}\left(\mathsf{\epsilon }\right),$$

(113)

$$M_{i\mathsf{\alpha },i^{\prime }{\mathsf{\alpha }}^{\prime }}^{(0)}=M_i{\mathsf{\delta }}_{i{i}^{\prime }}{\mathsf{\delta }}_{\mathsf{\alpha }{\mathsf{\alpha }}^{\prime }}.$$

(114)

Note that wave vector k varies within the first Brillouin zone. Furthermore, we have that ${\tilde{\Sigma }}_{eph}\left(\mathbf{k},\mathsf{\epsilon }\right)$ is the Fourier transformation of the matrix ${\Sigma }_{\mathrm{e}ph ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}(\mathsf{\epsilon })$ given in Equation (72), for which the terms $v^{\prime }{}_{ni\mathsf{\gamma }, n_3{i}_3{\mathsf{\gamma }}_3}^{n_1{i}_1{\mathsf{\alpha }}_1}$ are replaced by the values for a pure crystal and the corresponding Green’s functions are those of the effective medium. The other self-energies given by ${\tilde{\Sigma }}_{ee}\left(\mathbf{k},\mathsf{\epsilon }\right)$, ${\tilde{\Sigma }}_{phe}\left(\mathbf{k},\mathsf{\epsilon }\right)$, and ${\tilde{\Sigma }}_{phph}\left(\mathbf{k},\mathsf{\epsilon }\right)$ are defined similarly. In Equation (114), ${\Phi }^{(0)}\left(\mathbf{k}\right)$ is the Fourier transform of the matrix ${\Phi }^{(0)}{}_{ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}$, which describes the atomic nucleus repulsion. The self-energy ${\tilde{\Sigma }}_{phe}\left(\mathbf{k},\mathsf{\epsilon }\right)$ describes the attractive interaction between the atomic nuclei and the electrons.

The matrix ${\Sigma }_{phni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}^{\mathsf{\lambda }{n}_1{i}_1}(\mathsf{\epsilon })$ in expression (105) is diagonal. From expression (108), it follows that the matrix ${\mathsf{\sigma }}_{ph}\left(\mathsf{\epsilon }\right)$ is also diagonal, and its Fourier transform ${\mathsf{\sigma }}_{phi\mathsf{\alpha },i^{\prime }{\mathsf{\alpha }}^{\prime }}\left(\mathsf{\epsilon }\right)={\mathsf{\sigma }}_{phi}\left(\mathsf{\epsilon }\right){\mathsf{\delta }}_{i{i}^{\prime }}{\mathsf{\delta }}_{\mathsf{\alpha }{\mathsf{\alpha }}^{\prime }}$ in expression (114) does not depend on the wave vector $\mathbf{k}$. In the diagonal disorder approximation, the matrix $w_{ni\mathsf{\gamma },n^{\prime }{i}^{\prime }{\mathsf{\gamma }}^{\prime }}^{\mathsf{\lambda }{n}_1{i}_1}$ in expression (103) is diagonal in indices $n,n^{\prime }$. Neglecting the second term on the right-hand side of Equation (103), we obtain from Equation (107) that in this approximation the matrix ${\mathsf{\sigma }}_e\left(\mathsf{\epsilon }\right)$ is also diagonal in indices $n,n^{\prime }$, and its Fourier transform ${\mathsf{\sigma }}_{ei\mathsf{\gamma },i^{\prime }{\mathsf{\gamma }}^{\prime }}\left(\mathsf{\epsilon }\right)$ in expression (111) does not depend on wave vector $\mathbf{k}$. The Fourier transform of the mass operator of electron–phonon interaction has the form:

$$\begin{array}{l}{\tilde{\Sigma }}_{\mathrm{e}ph i\mathsf{\gamma },i^{\prime }{\mathsf{\gamma }}^{\prime }}(\mathbf{k},\mathsf{\epsilon })=-\frac{1}{4\pi i}\frac{1}{N}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{\mathsf{\epsilon }}_1\cot \mathrm{h}\left(\frac{{\mathsf{\epsilon }}_1}{2\Theta }\right)}\\ \times {\displaystyle \sum _{{\mathbf{k}}_1} {\Gamma }_{ i\mathsf{\gamma }, i_3{\mathsf{\gamma }}_3}^{\left(0\right)i_1{\mathsf{\alpha }}_1}\left(-\mathbf{k},\mathbf{k}-{\mathbf{k}}_1\right)}\left[{\tilde{G}}_{i_1{\mathsf{\alpha }}_1,i_2{\mathsf{\alpha }}_2}^{uu}({\mathbf{k}}_1{\mathsf{\epsilon }}_1)-{\tilde{G}}_{i_1{\mathsf{\alpha }}_1,i_2{\mathsf{\alpha }}_2}^{uu *}({\mathbf{k}}_1{\mathsf{\epsilon }}_1)\right]\\ \times {\tilde{G}}_{i_3{\mathsf{\gamma }}_3,i_4{\mathsf{\gamma }}_4}^{a{a}^{+}}\left(\mathbf{k}-{\mathbf{k}}_1,\mathsf{\epsilon }-{\mathsf{\epsilon }}_1\right){\tilde{\Gamma }}_{ i_4{\mathsf{\gamma }}_4,i^{\prime }{\mathsf{\gamma }}^{\prime }}^{i_2{\mathsf{\alpha }}_2}\left(-\mathbf{k}+{\mathbf{k}}_1,\mathbf{k}\right).\end{array}$$

(115)

The Fourier transform of the phonon–electron interaction mass operator is:

$$\begin{array}{l}{\tilde{\Sigma }}_{ph\mathrm{e} ni\mathsf{\alpha },n^{\prime }{i}^{\prime }{\mathsf{\alpha }}^{\prime }}(\mathbf{k},\mathsf{\epsilon })=\frac{1}{2\pi i}\frac{1}{N}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{\mathsf{\epsilon }}_1 f({\mathsf{\epsilon }}_1)}{\displaystyle \sum _{}{\Gamma }_{i_2{\mathsf{\gamma }}_2,i_1{\mathsf{\gamma }}_1}^{ \left(0\right)i\mathsf{\alpha }}}\left(-{\mathbf{k}}_1,\mathbf{k}+{\mathbf{k}}_1\right)\\ \times \{\left[{\tilde{G}}_{i_1{\mathsf{\gamma }}_1,i_3{\mathsf{\gamma }}_3}^{a{a}^{+}}(\mathbf{k}+{\mathbf{k}}_1,\mathsf{\epsilon }+{\mathsf{\epsilon }}_1)-{\tilde{G}}_{i_1{\mathsf{\gamma }}_1,i_3{\mathsf{\gamma }}_3}^{a{a}^{+}*}(\mathbf{k}+{\mathbf{k}}_1,\mathsf{\epsilon }+{\mathsf{\epsilon }}_1)\right]\\ \times {\tilde{G}}_{i_4{\mathsf{\gamma }}_4,i_2{\mathsf{\gamma }}_2}^{a{a}^{+}*}({\mathbf{k}}_1,{\mathsf{\epsilon }}_1)+{\tilde{G}}_{i_1{\mathsf{\gamma }}_1,i_3{\mathsf{\gamma }}_3}^{a{a}^{+}}(\mathbf{k}+{\mathbf{k}}_1,\mathsf{\epsilon }+{\mathsf{\epsilon }}_1)\\ \times \left[{\tilde{G}}_{i_4{\mathsf{\gamma }}_4,i_2{\mathsf{\gamma }}_2}^{a{a}^{+}}({\mathbf{k}}_1,{\mathsf{\epsilon }}_1)-{\tilde{G}}_{i_4{\mathsf{\gamma }}_4,i_2{\mathsf{\gamma }}_2}^{a{a}^{+}*}({\mathbf{k}}_1,{\mathsf{\epsilon }}_1)\right]\}{\tilde{\Gamma }}_{ i_3{\mathsf{\gamma }}_3,i_4{\mathsf{\gamma }}_4}^{ i^{\prime }{\mathsf{\alpha }}^{\prime }}\left(-\mathbf{k}-{\mathbf{k}}_1,{\mathbf{k}}_1\right).\end{array}$$

(116)

The vertex parts of the mass operators of electron–phonon and phonon–electron interactions are determined by the equation:

$$\begin{array}{l}{\tilde{\Gamma }}_{ i_3{\mathsf{\gamma }}_3,i_4{\mathsf{\gamma }}_4}^{n{i}^{\prime }{\mathsf{\alpha }}^{\prime }}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)={\Gamma }_{ i_3{\mathsf{\gamma }}_3,i_4{\mathsf{\gamma }}_4}^{\left(0\right)i^{\prime }{\mathsf{\alpha }}^{\prime }}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)-\frac{1}{2\pi i}\frac{1}{N}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d\mathsf{\epsilon }f\left(\mathsf{\epsilon }\right)}\\ \times {\displaystyle \sum _{}{\Gamma }_{ i_5{\mathsf{\gamma }}_5,i_6{\mathsf{\gamma }}_6}^{\left(0\right)i^{\prime }{\mathsf{\alpha }}^{\prime }}\left({\mathbf{k}}_1+{\mathbf{k}}_2-{\mathbf{k}}_5,{\mathbf{k}}_5\right)}\\ \times [{\tilde{G}}_{i_6{\mathsf{\gamma }}_6,i_7{\mathsf{\gamma }}_7}^{a{a}^{+}}\left({\mathbf{k}}_5,\mathsf{\epsilon }\right){\tilde{G}}_{i_8{\mathsf{\gamma }}_8,i_5{\mathsf{\gamma }}_5}^{a{a}^{+}}\left(-{\mathbf{k}}_1-{\mathbf{k}}_2+{\mathbf{k}}_5\mathsf{\epsilon }\right)\\ -{\tilde{G}}_{i_6{\mathsf{\gamma }}_6,i_7{\mathsf{\gamma }}_7}^{a{a}^{+}\ast }\left({\mathbf{k}}_5,\mathsf{\epsilon }\right){\tilde{G}}_{i_8{\mathsf{\gamma }}_8,i_5{\mathsf{\gamma }}_5}^{a{a}^{+}\ast }\left(-{\mathbf{k}}_1-{\mathbf{k}}_2+{\mathbf{k}}_5\mathsf{\epsilon }\right)]\\ \times {\Gamma }_{ i_7{\mathsf{\gamma }}_7,i_8{\mathsf{\gamma }}_8}^{\left(0\right)i_9{\mathsf{\alpha }}_9}\left(-{\mathbf{k}}_5,-{\mathbf{k}}_1-{\mathbf{k}}_2+{\mathbf{k}}_5\right){\tilde{G}}_{i_9{\mathsf{\alpha }}_9,i_{10}{\mathsf{\alpha }}_{10}}^{uu}\left({\mathbf{k}}_1+{\mathbf{k}}_2,0\right)\\ \times {\tilde{\Gamma }}_{i_3{\mathsf{\gamma }}_3,i_4{\mathsf{\gamma }}_4}^{i_{10}{\mathsf{\alpha }}_{10}}\left({\mathbf{k}}_1,{\mathbf{k}}_1+{\mathbf{k}}_2\right).\end{array}$$

(117)

In expressions (116)–(118)

$$\begin{array}{l}{\Gamma }_{ i_1{\mathsf{\gamma }}_1,i_2{\mathsf{\gamma }}_2}^{\left(0\right)i\mathsf{\alpha }}\left({\mathbf{k}}_1,{\mathbf{k}}_2\right)=\\ {\displaystyle \sum _{n_1,n_2}{v}^{\prime }{}_{n_1{i}_1{\mathsf{\gamma }}_1,n_2{i}_2{\mathsf{\gamma }}_2}^{ni\mathsf{\alpha }}}\exp \left(i{\mathbf{k}}_1\left({\mathbf{r}}_{n_1{i}_1}-{\mathbf{r}}_{ni}\right)+i{\mathbf{k}}_2\left({\mathbf{r}}_{n_2{i}_2}-{\mathbf{r}}_{ni}\right)\right).\end{array}$$

(118)

The Fourier transform of the mass operator of the electron–electron interaction can be represented as:

$${\tilde{\Sigma }}_{ee i\mathsf{\gamma },i^{\prime }{\mathsf{\gamma }}^{\prime }}(\mathbf{k},\mathsf{\epsilon })={\tilde{\Sigma }}_{ee i\mathsf{\gamma },i^{\prime }{\mathsf{\gamma }}^{\prime }}^{(1)}\left(\mathbf{k}\right)+{\tilde{\Sigma }}_{ee i\mathsf{\gamma },i^{\prime }{\mathsf{\gamma }}^{\prime }}^{(2)}(\mathbf{k},\mathsf{\epsilon }),$$

(119)

$$\begin{array}{l}{\tilde{\Sigma }}_{\mathrm{ee} \mathrm{i}\mathsf{\gamma },i^{\prime }{\mathsf{\gamma }}^{\prime }}^{(1)}\left(\mathbf{k}\right)=-\frac{1}{2\pi i}\frac{1}{N}\\ \times {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{\mathsf{\epsilon }}_1 f\left({\mathsf{\epsilon }}_1\right){\displaystyle \sum _{{\mathbf{k}}_1}{\Gamma }_{i\mathsf{\gamma },i^{\prime }{\mathsf{\gamma }}^{\prime }}^{ \left(0\right) {\mathrm{i}}_2{\mathsf{\gamma }}_2,i_1{\mathsf{\gamma }}_1}\left(-\mathbf{k},-{\mathbf{k}}_1,{\mathbf{k}}_1\right)}}\\ \times \left[{\tilde{G}}_{i_1{\mathsf{\gamma }}_1,i_2{\mathsf{\gamma }}_2}^{a{a}^{+}}({\mathbf{k}}_1,{\mathsf{\epsilon }}_1)-{\tilde{G}}_{i_1{\mathsf{\gamma }}_1,i_2{\mathsf{\gamma }}_2}^{a{a}^{+} *}({\mathbf{k}}_1,{\mathsf{\epsilon }}_1)\right]\end{array},$$

(120)

$$\begin{gathered} \begin{array}{l}{\tilde{\Sigma }}_{\mathrm{ee}\hspace{0.17em}i\mathsf{\gamma },i^{\prime }{\mathsf{\gamma }}^{\prime }}^{(2)}\left(\mathbf{k},\mathsf{\epsilon }\right)={\left(\frac{1}{2\pi i}\right)}^2\frac{1}{N^2}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{\mathsf{\epsilon }}_1}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{\mathsf{\epsilon }}_2}\\ \times f\left({\mathsf{\epsilon }}_1\right)f\left({\mathsf{\epsilon }}_2\right){\displaystyle \sum _{{\mathbf{k}}_1,{\mathbf{k}}_2}{\Gamma }_{{ \mathrm{i}}_2{\mathsf{\gamma }}_2,i_1{\mathsf{\gamma }}_1}^{ \left(0\right) \mathrm{i}\mathsf{\gamma },i_3{\mathsf{\gamma }}_3}\left(-\mathbf{k},-{\mathbf{k}}_1-{\mathbf{k}}_2+\mathbf{k},{\mathbf{k}}_1\right)}\end{array} \\ \begin{array}{l}\times \{[G_{i_2{\mathsf{\gamma }}_2,i_5{\mathsf{\gamma }}_5}^{a{a}^{+} }(\mathbf{k}-{\mathbf{k}}_1-{\mathbf{k}}_2,\mathsf{\epsilon }-{\mathsf{\epsilon }}_1-{\mathsf{\epsilon }}_2)G_{i_1{\mathsf{\gamma }}_1,i_4{\mathsf{\gamma }}_4}^{a{a}^{+}\ast }({\mathbf{k}}_1,{\mathsf{\epsilon }}_1)-\\ -G_{i_2{\mathsf{\gamma }}_2,i_5{\mathsf{\gamma }}_5}^{a{a}^{+}\ast }(\mathbf{k}-{\mathbf{k}}_1-{\mathbf{k}}_2,\mathsf{\epsilon }-{\mathsf{\epsilon }}_1-{\mathsf{\epsilon }}_2)G_{i_1{\mathsf{\gamma }}_1,i_4{\mathsf{\gamma }}_4}^{a{a}^{+}}({\mathbf{k}}_1,{\mathsf{\epsilon }}_1)]\end{array} \\ \begin{array}{l}\times \left[{\tilde{G}}_{i_6{\mathsf{\gamma }}_6,i_3{\mathsf{\gamma }}_3}^{a{a}^{+}}\left({\mathbf{k}}_2,{\mathsf{\epsilon }}_2\right)-{\tilde{G}}_{i_6{\mathsf{\gamma }}_6,i_3{\mathsf{\gamma }}_3}^{a{a}^{+}\ast }\left({\mathbf{k}}_2,{\mathsf{\epsilon }}_2\right)\right]\\ -[{\tilde{G}}_{i_2{\mathsf{\gamma }}_2,i_5{\mathsf{\gamma }}_5}^{a{a}^{+}}\left(\mathbf{k}-{\mathbf{k}}_1-{\mathbf{k}}_2,\mathsf{\epsilon }-{\mathsf{\epsilon }}_1-{\mathsf{\epsilon }}_2\right)\\ -{\tilde{G}}_{i_2{\mathsf{\gamma }}_2,i_5{\mathsf{\gamma }}_5}^{a{a}^{+}\ast }\left(\mathbf{k}-{\mathbf{k}}_1-{\mathbf{k}}_2,\mathsf{\epsilon }-{\mathsf{\epsilon }}_1-{\mathsf{\epsilon }}_2\right)]\end{array} \\ \begin{array}{l}\times [{\tilde{G}}_{i_1{\mathsf{\gamma }}_1,i_4{\mathsf{\gamma }}_4}^{a{a}^{+}}\left({\mathbf{k}}_1,{\mathsf{\epsilon }}_1\right){\tilde{G}}_{i_6{\mathsf{\gamma }}_6,i_3{\mathsf{\gamma }}_3}^{a{a}^{+}}\left({\mathbf{k}}_2,{\mathsf{\epsilon }}_2\right)\\ -{\tilde{G}}_{i_1{\mathsf{\gamma }}_1,i_4{\mathsf{\gamma }}_4}^{a{a}^{+}\ast }\left({\mathbf{k}}_1,{\mathsf{\epsilon }}_1\right){\tilde{G}}_{i_6{\mathsf{\gamma }}_6,i_3{\mathsf{\gamma }}_3}^{a{a}^{+}\ast }\left({\mathbf{k}}_2,{\mathsf{\epsilon }}_2\right)]\}\\ \times {\tilde{\Gamma }}_{ {\mathrm{i}}_4{\mathsf{\gamma }}_4,i^{\prime }{\mathsf{\gamma }}^{\prime }}^{ {\mathrm{i}}_5{\mathsf{\gamma }}_5,i_6{\mathsf{\gamma }}_6}\left({\mathbf{k}}_1+{\mathbf{k}}_2-\mathbf{k},-{\mathbf{k}}_2,{\mathbf{k}}_1\right).\end{array} \end{gathered}$$

(121)

The vertex part of the mass operator of the electron–electron interaction is determined by the equation:

$$\begin{array}{l}{\tilde{\Gamma }}_{ {\mathrm{i}}_4{\mathsf{\gamma }}_4,i^{\prime }{\mathsf{\gamma }}^{\prime }}^{ {\mathrm{i}}_5{\mathsf{\gamma }}_5,i_6{\mathsf{\gamma }}_6}\left({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3\right)={\Gamma }_{ {\mathrm{i}}_4{\mathsf{\gamma }}_4,i^{\prime }{\mathsf{\gamma }}^{\prime }}^{ \left(0\right){\mathrm{i}}_5{\mathsf{\gamma }}_5,i_6{\mathsf{\gamma }}_6}\left({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3\right)\\ -\frac{1}{2\pi i}\frac{1}{N}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d\mathsf{\epsilon }f\left(\mathsf{\epsilon }\right){\displaystyle \sum _{{\mathbf{k}}_4}{\Gamma }_{{ \mathrm{i}}_4{\mathsf{\gamma }}_4,i_8{\mathsf{\gamma }}_8}^{\left(0\right) {\mathrm{i}}_5{\mathsf{\gamma }}_5,i_7{\mathsf{\gamma }}_7}}}\left({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_4\right)\\ \times [{\tilde{G}}_{i_7{\mathsf{\gamma }}_7,i_9{\mathsf{\gamma }}_9}^{a{a}^{+}}\left({\mathbf{k}}_4,\mathsf{\epsilon }\right){\tilde{G}}_{i_8{\mathsf{\gamma }}_8,i_{10}{\mathsf{\gamma }}_{10}}^{a{a}^{+}\ast }\left(-{\mathbf{k}}_1-{\mathbf{k}}_2-{\mathbf{k}}_4,\mathsf{\epsilon }\right)\\ -{\tilde{G}}_{i_7{\mathsf{\gamma }}_7,i_9{\mathsf{\gamma }}_9}^{a{a}^{+}\ast }\left({\mathbf{k}}_4,\mathsf{\epsilon }\right){\tilde{G}}_{i_8{\mathsf{\gamma }}_8,i_{10}{\mathsf{\gamma }}_{10}}^{a{a}^{+}}\left(-{\mathbf{k}}_1-{\mathbf{k}}_2-{\mathbf{k}}_4,\mathsf{\epsilon }\right)]\\ \times {\tilde{\Gamma }}_{ {\mathrm{i}}_{10}{\mathsf{\gamma }}_{10},i^{\prime }{\mathsf{\gamma }}^{\prime }}^{ {\mathrm{i}}_9{\mathsf{\gamma }}_9,i_6{\mathsf{\gamma }}_6}\left({\mathbf{k}}_1+{\mathbf{k}}_2+{\mathbf{k}}_4,-{\mathbf{k}}_4,{\mathbf{k}}_3\right).\end{array}$$

(122)

In expression (123)

$$\begin{array}{l}{\Gamma }_{ i_1{\mathsf{\gamma }}_1,i\mathsf{\gamma }}^{\left(0\right)i_2{\mathsf{\gamma }}_2,i_3{\mathsf{\gamma }}_3}\left({\mathbf{k}}_1,{\mathbf{k}}_2,{\mathbf{k}}_3\right)=\\ {\displaystyle \sum _{n_1,n_2,n_3}{\tilde{v}}_{n_3{i}_3{\mathsf{\gamma }}_3,ni\mathsf{\gamma }}^{\left(2\right) {\mathrm{n}}_1{i}_1{\mathsf{\gamma }}_1,{\mathrm{n}}_2{i}_2{\mathsf{\gamma }}_2}}\exp \left(i{\mathbf{k}}_1\left({\mathbf{r}}_{n_1{i}_1}-{\mathbf{r}}_{ni}\right)\right)\\ \times \exp \left(i{\mathbf{k}}_2\left({\mathbf{r}}_{n_2{i}_2}-{\mathbf{r}}_{ni}\right)+i{\mathbf{k}}_3\left({\mathbf{r}}_{n_3{i}_3}-{\mathbf{r}}_{ni}\right)\right).\end{array}$$

(123)

Cluster decomposition for the Green’s function of electrons and phonons of disordered crystal can be obtained from Equations (95)–(100). The density of the electron and phonon states are presented as infinite series. Here, processes of scattering on clusters with different numbers of atoms are described by each term. It is shown that the contribution of the scattering processes of electrons and phonons in clusters decreases with increasing the number of atoms in the cluster by a small parameter

$$p(\mathsf{\epsilon })=\frac{1}{r\mathsf{\nu }}\left|{\displaystyle \sum _{(n_2{i}_2)\ne (n_1{i}_1),i,\mathsf{\gamma }}{\langle t^{n_1{i}_1}(\mathsf{\epsilon })\tilde{G}(\mathsf{\epsilon })t^{n_2{i}_2}(\mathsf{\epsilon })\tilde{G}(\mathsf{\epsilon })\rangle }_{0i\mathsf{\gamma },0i\mathsf{\gamma }}}\right|,$$

(124)

where $r$ is the total number of energy bands included in the calculation.

We have shown previously [26,30,31,32] that this parameter remains small when many parameters of the system are changed, except possibly for narrow energy intervals near the band edges.

By neglecting the contribution of processes of electron scattering in clusters consisting of three or more atoms that are small by the above parameter in Equation (125) for the density of electronic states, we obtain:

$$g_e\left(\mathsf{\epsilon }\right)=\frac{1}{v}{\displaystyle \sum _{i,\mathsf{\delta },\mathsf{\sigma },\mathsf{\lambda },m_{\mathsf{\lambda }i}}{P}_{0i}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}}{g}_{0\hspace{0.17em}i\mathsf{\delta }\hspace{0.17em}\mathsf{\sigma }}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}}}\left(\mathsf{\epsilon }\right),$$

(125)

$$\begin{array}{l}{g}_{0\hspace{0.17em}i\mathsf{\delta }\hspace{0.17em}\mathsf{\sigma }}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}}(\mathsf{\epsilon })=-\frac{1}{\pi } \mathrm{Im} \{\tilde{G}+\tilde{G} t^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}\hspace{0.17em}\tilde{G}+{\displaystyle \sum _{\begin{array}{c}\left(lj\right)\ne \left(0i\right)\\ {\mathsf{\lambda }}^{\prime },m_{{\mathsf{\lambda }}^{\prime }j}\end{array}}{P}_{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}lj\hspace{0.17em}\hspace{0.17em}0\hspace{0.17em}i}^{{\mathsf{\lambda }}^{\prime }\hspace{0.17em}{m}_{{\mathsf{\lambda }}^{\prime }j}/\mathsf{\lambda }\hspace{0.17em}{m}_{\mathsf{\lambda }i}}}\\ \times \tilde{G} \left[t^{{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}\hspace{0.17em}+T^{(2)\mathsf{\lambda }{m}_{\mathsf{\lambda }i} 0i,{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}\right]{{\tilde{G}}^{\hspace{0.17em}}\}}^{0i\mathsf{\delta }\mathsf{\sigma },0i\mathsf{\delta }\mathsf{\sigma }},\end{array}$$

(126)

$$\begin{array}{l}{T}^{(2)\mathsf{\lambda }{m}_{\mathsf{\lambda }i} 0i,{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}={\left[I-t^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i} 0i}\tilde{G}{t}^{{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}\tilde{G}\right]}^{-1}\\ \times t^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i} 0i}\tilde{G}{t}^{{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}\left[I+\tilde{G}{t}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i} 0i}\right]\end{array}$$

(127)

where $\tilde{G}={\tilde{G}}^{a{a}^{+}}(\mathsf{\epsilon })$.

Similarly averaging of the phonon Green’s function $G^{uu}(\mathsf{\epsilon })$ yields the phononic density of states:

$$g_{ph}\left(\mathsf{\epsilon }\right)=\frac{1}{\mathsf{\nu }}{\displaystyle \sum _{i,\mathsf{\alpha },\mathsf{\lambda }}{P}_{0i}^{\mathsf{\lambda }}{g}_{0i\mathsf{\alpha }\hspace{0.17em}}^{\mathsf{\lambda }}}\left(\mathsf{\epsilon }\right),$$

(128)

$$\begin{array}{l}{g}_{\hspace{0.17em}0i\mathsf{\alpha }}^{\mathsf{\lambda }}(\mathsf{\epsilon })=-\frac{1}{\pi } 2\frac{\mathsf{\epsilon }}{{\hslash }^2}{M}_i\mathrm{Im} \{\tilde{G}+\tilde{G} t^{\mathsf{\lambda }\hspace{0.17em}0i}\hspace{0.17em}\tilde{G}+{\displaystyle \sum _{\begin{array}{c}\left(lj\right)\ne \left(0i\right)\\ {\mathsf{\lambda }}^{\prime }\end{array}}{P}_{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}lj\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0i}^{{\mathsf{\lambda }}^{\prime }/\mathsf{\lambda }} }\\ \times \tilde{G} \left[t^{{\mathsf{\lambda }}^{\prime }\hspace{0.17em}lj}\hspace{0.17em}+T^{(2)\mathsf{\lambda } 0i,{\mathsf{\lambda }}^{\prime }lj}\right]{{\tilde{G}}^{\hspace{0.17em}}\}}^{0i\mathsf{\alpha },0i\mathsf{\alpha }},\end{array}$$

(129)

where $\tilde{G}={\tilde{G}}^{uu}(\mathsf{\epsilon })$.

In Equation (127), $P_{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}lj\hspace{0.17em}\hspace{0.17em}0\hspace{0.17em}i}^{{\mathsf{\lambda }}^{\prime }\hspace{0.17em}{m}_{{\mathsf{\lambda }}^{\prime }j}/\mathsf{\lambda }\hspace{0.17em}{m}_{\mathsf{\lambda }i}}$ is the conditional probability to find an atom of type λ′ at site (lj) for the atom with magnetic moment m_λ′j, provided that the sites in the unit cell at the origin (0i) have an atom of type λ with a magnetic moment m_λi. Here, $t_{ni}^{\mathsf{\lambda }\hspace{0.17em}{m}_{\mathsf{\lambda }i}}$ is the value of the matrix element of a single-center operator for scattering in the case where an atom of type λ is located at site (ni) and has a magnetic moment m_λi.

When the system is disordered, we need to consider a random arrangement of the disordered atomic sites. Hence, in Equation (129), the probability of an atom of type λ to be at site (0i) is given by

$$P_{0i}^{\mathsf{\lambda }}=<c_{0i}^{\mathsf{\lambda }}>,$$

(130)

where $c_{ni}^{\mathsf{\lambda }}$ is a discrete binary random number taking the values of 1 or 0, depending on whether an atom of type $\mathsf{\lambda }$ is at site (ni) or not, respectively (36). The joint probabilities in Equations (126), (127), (129), and (130) are defined by the following:

$$\begin{gathered} P_{lj\hspace{0.17em}\hspace{0.17em}0i}^{{\mathsf{\lambda }}^{\prime }\hspace{0.17em}\mathsf{\lambda }}=P_{0i}^{\mathsf{\lambda }\hspace{0.17em}}{P}_{lj\hspace{0.17em}\hspace{0.17em}0i}^{{\mathsf{\lambda }}^{\prime }/\mathsf{\lambda }}=<c_{lj}^{{\mathsf{\lambda }}^{\prime }}{c}_{0i}^{\mathsf{\lambda }}>, \\ {P}_{0i}^{\mathsf{\lambda }\hspace{0.17em}{m}_{\mathsf{\lambda }i}}=P_{0i}^{\mathsf{\lambda }\hspace{0.17em}}\hspace{0.17em}{P}_{0i}^{\hspace{0.17em}{m}_{\mathsf{\lambda }i}},P_{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}lj\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0i}^{{\mathsf{\lambda }}^{\prime }\hspace{0.17em}{m}_{{\mathsf{\lambda }}^{\prime }j}/\hspace{0.17em}\mathsf{\lambda }{m}_{\mathsf{\lambda }i}}=P_{\hspace{0.17em}lj\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0i}^{{\mathsf{\lambda }}^{\prime }\hspace{0.17em}/\mathsf{\lambda }}\hspace{0.17em}\hspace{0.17em}{P}_{\hspace{0.17em}\hspace{0.17em}lj\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0i}^{\hspace{0.17em}{m}_{{\mathsf{\lambda }}^{\prime }j}/\hspace{0.17em}{m}_{\mathsf{\lambda }i}} \\ {P}_{lj\hspace{0.17em}\hspace{0.17em}0i}^{m_{{\mathsf{\lambda }}^{\prime }j}\hspace{0.17em}{m}_{\mathsf{\lambda }i}}=P_{0i}^{m_{\mathsf{\lambda }i}}{P}_{lj\hspace{0.17em}\hspace{0.17em}0i}^{m_{{\mathsf{\lambda }}^{\prime }j}/m_{\mathsf{\lambda }i}}=<c_{lj}^{m_{{\mathsf{\lambda }}^{\prime }j}}{c}_{0i}^{m_{\mathsf{\lambda }i}}>. \end{gathered}$$

(131)

The probabilities are determined by the interatomic pair correlations ${\mathsf{\epsilon }}_{lj\hspace{0.17em}\hspace{0.17em}0i}^{\hspace{0.17em}B\hspace{0.17em}B}$, ${\mathsf{\epsilon }}_{\hspace{0.17em}\hspace{0.17em}lj\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0i}^{\hspace{0.17em}{\mathsf{\mu }}_{{\mathsf{\lambda }}^{\prime }j}^{-}\hspace{0.17em}{\mathsf{\mu }}_{\mathsf{\lambda }i}^{-}}$ via [30]:

$$P_{lj\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0i}^{{\mathsf{\lambda }}^{\prime }/\mathsf{\lambda }}=P_{lj}^{{\mathsf{\lambda }}^{\prime }}+\frac{{\mathsf{\epsilon }}_{lj\hspace{0.17em}\hspace{0.17em}0i}^{B\hspace{0.17em}B}}{P_{0i}^{\mathsf{\lambda }}}({\mathsf{\delta }}_{{\mathsf{\lambda }}^{\prime }B}-{\mathsf{\delta }}_{{\mathsf{\lambda }}^{\prime }A})({\mathsf{\delta }}_{\mathsf{\lambda }B}-{\mathsf{\delta }}_{\mathsf{\lambda }A}),\begin{array}{l}{P}_{lj\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0i}^{m_{{\mathsf{\lambda }}^{\prime }j}/m_{\mathsf{\lambda }i}}=P_{lj}^{m_{{\mathsf{\lambda }}^{\prime }j}}+\frac{{\mathsf{\epsilon }}_{\hspace{0.17em}\hspace{0.17em}lj\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0i}^{\hspace{0.17em}{\mathsf{\mu }}_{{\mathsf{\lambda }}^{\prime }j}^{-}\hspace{0.17em}{\mathsf{\mu }}_{\mathsf{\lambda }i}^{-}}}{P_{0i}^{m_{\mathsf{\lambda }i}}}({\mathsf{\delta }}_{m_{{\mathsf{\lambda }}^{\prime }j},{\mathsf{\mu }}_j^{-}}-\\ -{\mathsf{\delta }}_{m_{{\mathsf{\lambda }}^{\prime }j},{\mathsf{\mu }}_j^{+}})\left({\mathsf{\delta }}_{m_{\mathsf{\lambda }i},{\mathsf{\mu }}_i^{-}}-{\mathsf{\delta }}_{m_{\mathsf{\lambda }i},{\mathsf{\mu }}_i^{+}}\right)\end{array},$$

(132)

where δ is the Kronecker delta function. Note that the interatomic pair correlations also satisfy

$${\mathsf{\epsilon }}_{lj\hspace{0.17em}\hspace{0.17em}0i}^{BB}=<(c_{lj}^B-c_j^B)(c_{0i}^B-c_i^B)>, {\mathsf{\epsilon }}_{\hspace{0.17em}\hspace{0.17em}lj\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0i}^{\hspace{0.17em}{\mathsf{\mu }}_{{\mathsf{\lambda }}^{\prime }j}^{-}\hspace{0.17em}{\mathsf{\mu }}_{\mathsf{\lambda }i}^{-}}=<(c_{lj}^{{\mathsf{\mu }}_{{\mathsf{\lambda }}^{\prime }j}^{-}}-c_j^{{\mathsf{\mu }}_{{\mathsf{\lambda }}^{\prime }j}^{-}})(c_{0i}^{{\mathsf{\mu }}_{\mathsf{\lambda }i}^{-}}-c_i^{{\mathsf{\mu }}_{\mathsf{\lambda }i}^{-}})>.$$

(133)

The notations $P_{0i}^{\hspace{0.17em}{m}_{\mathsf{\lambda }i}}$ and $P_{\hspace{0.17em}\hspace{0.17em}lj\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0i}^{\hspace{0.17em}{m}_{{\mathsf{\lambda }}^{\prime }j}/\hspace{0.17em}{m}_{\mathsf{\lambda }i}}$ indicate the probabilities of the static fluctuations of the magnetization.

As an example, when we have a binary alloy, consisting of two sublattices, and two types of atoms A and B, we obtain

$$P_{0i}^A=x_A-\frac{{\mathsf{\nu }}_2}{\mathsf{\nu }}{\mathsf{\eta }}_a$$

(134)

for the first sublattice and

$$P_{0i}^A=x_A+\frac{{\mathsf{\nu }}_1}{\mathsf{\nu }}{\mathsf{\eta }}_a$$

(135)

for the second sublattice, with

$$P_{0i}^B=1-P_{0i}^A.$$

(136)

Here, ν = ν₁ + ν₂ is the total number of sublattice sites, x_A, and x_B = 1 − x_A are the concentrations of the atomic components A and B in the alloy, and η_a is the parameter that measures the long-range atomic order.

The two values $m_{\mathsf{\lambda }i}^{}={\mathsf{\mu }}_{\mathsf{\lambda }i}^{+}$ and ${\mathsf{\mu }}_{\mathsf{\lambda }i}^{-}$ represent the projections of the localized magnetic moment onto the z axis. The probability $P_{0i}^{\hspace{0.17em}{m}_{\mathsf{\lambda }i}}$ is connected with the long-range magnetic parameter η_m via the expressions

$$P_{0i}^{{\mathsf{\mu }}_{\mathsf{\lambda }i}^{+}}=x_{{\mathsf{\mu }}_{\mathsf{\lambda }}^{+}}-\frac{{\mathsf{\nu }}_2}{\mathsf{\nu }}{\mathsf{\eta }}_m$$

(137)

for sublattice 1 and

$$P_{0i}^{{\mathsf{\mu }}_{\mathsf{\lambda }i}^{+}}=x_{{\mathsf{\mu }}_{\mathsf{\lambda }}^{+}}+\frac{{\mathsf{\nu }}_1}{\mathsf{\nu }}{\mathsf{\eta }}_m$$

(138)

for sublattice 2, with

$$P_{0i}^{{\mathsf{\mu }}_{\mathsf{\lambda }i}^{-}}=1-P_{0i}^{{\mathsf{\mu }}_{\mathsf{\lambda }i}^{+}}.$$

(139)

Here, $x_{{\mathsf{\mu }}_{\mathsf{\lambda }}^{+}}$ and $x_{{\mathsf{\mu }}_{\mathsf{\lambda }}^{-}}=1-x_{{\mathsf{\mu }}_{\mathsf{\lambda }}^{+}}$ are equal to the relative number of lattice sites with localized magnetic moment projections ${\mathsf{\mu }}_{\mathsf{\lambda }i}^{+}$ and ${\mathsf{\mu }}_{\mathsf{\lambda }i}^{-}$, respectively.

For an ideally ordered crystal, the Green’s function in Equation (97) is:

$$G\left(\mathsf{\epsilon }\right)=\lim \tilde{G}\left(\mathsf{\epsilon }\right),\mathsf{\sigma }\left(\mathsf{\epsilon }\right)\to 0,$$

(140)

where the Green’s function $\tilde{G}\left(\mathsf{\epsilon }\right)$ is given by Formulas (95) and (96). The energies of the electrons and phonons of the crystal are determined from the equations for the poles of the Green’s functions:

$$\det \Vert \mathsf{\epsilon }{\mathsf{\delta }}_{i{i}^{\prime }}{\mathsf{\delta }}_{\mathsf{\gamma }{\mathsf{\gamma }}^{\prime }}-{\tilde{H}}_{i\mathsf{\gamma },i^{\prime }{\mathsf{\gamma }}^{\prime }}\left(\mathbf{k},\mathsf{\epsilon }\right)\Vert =0,$$

(141)

$$\det \Vert \frac{{\mathsf{\epsilon }}^2}{{\hslash }^2}{M}_i{\mathsf{\delta }}_{i{i}^{\prime }}{\mathsf{\delta }}_{\mathsf{\alpha }{\mathsf{\alpha }}^{\prime }}-{\tilde{\Phi }}_{i\mathsf{\alpha },i^{\prime }{\mathsf{\alpha }}^{\prime }}\left(\mathbf{k},\mathsf{\epsilon }\right)\Vert =0,$$

(142)

where ${\tilde{H}}_{i\mathsf{\gamma },i^{\prime }{\mathsf{\gamma }}^{\prime }}\left(\mathbf{k},\mathsf{\epsilon }\right)$, ${\tilde{\Phi }}_{i\mathsf{\alpha },i^{\prime }{\mathsf{\alpha }}^{\prime }}\left(\mathbf{k},\mathsf{\epsilon }\right)$ are given by Formulas (111) and (114).

6. Free Energy

The Gibbs free energy or, in other words, the thermodynamic potential of the system, satisfies [27]:

$$\Omega =-\Theta \ln \mathrm{Tr}(e^{-H/\Theta }).$$

(143)

The Hamiltonian H is defined in Equation (1). To perform the trace, we need to sum over all of the band states, but we also need to take into account the disorder averaging. The latter is commonly handled via a configurational average [26]. Using Formulas (50) and (144), we represent the thermodynamic potential in the form:

$$\Omega =\langle \mathsf{\delta }\Phi \rangle -\Theta S_c+{\Omega }_e^{(0)}+{\Omega }_{ph}^{(0)}+{\Omega }^{\prime },$$

(144)

where ${\Omega }_e^{(0)},{\Omega }_{ph}^{(0)}$ are the thermodynamic potentials for the electrons and the phonons in the field of the ionic cores, respectively. ${\Omega }^{\prime }$ is the component of the thermodynamic potential that is caused by the mutual scattering of electrons and phonons; it is defined by

$${\Omega }^{\prime }=-\Theta \ln \langle <\mathsf{\sigma }(1/\Theta ){>}_0\rangle ,$$

(145)

with σ given in Equation (50) for the interaction picture.

In addition, $S_c=-<\ln P_c>$ is the configurational entropy, where $P_c$ denotes the distribution function for atoms with a specific z-component of the magnetic moment on a given lattice site. The angular brackets $\langle \dots \rangle$ denote the configurational averaging over different disorder configurations for a given density of disorder.

Next, we use the “integration over the coupling constant” method to simplify the results further. By replacing the interacting Hamiltonian $H_{\mathrm{int}}$ (defined in Equation (5)) by H_int(λ) = λH_int, differentiating the expression for the piece of the thermodynamic potential ${\Omega }^{\prime }(\mathsf{\lambda })$ in Equation (146) with respect to $\mathsf{\lambda }$ and then integrating (with the boundary conditions ${\Omega }^{\prime }(0)=0$, ${\Omega }^{\prime }(1)={\Omega }^{\prime }$), we obtain the following after a long derivation:

$$\begin{array}{l}{\Omega }^{\prime }=-\frac{1}{\pi \mathsf{\nu }N}\mathrm{Im}{\displaystyle \underset{0}{\overset{1}{\int }}\frac{d\mathsf{\lambda }}{\mathsf{\lambda }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d\mathsf{\epsilon }}}[f(\mathsf{\epsilon })\times \\ \times Tr\langle \left(w(\mathsf{\lambda })+{\Sigma }_{eph}(\mathsf{\epsilon },\mathsf{\lambda })+{\Sigma }_{ee}(\mathsf{\epsilon },\mathsf{\lambda })\right)G^{a{a}^{+}}(\mathsf{\epsilon },\mathsf{\lambda })\rangle +\\ +\frac{1}{2}\coth \left(\frac{\mathsf{\epsilon }}{2\Theta }\right)Tr\langle \Delta M^{-1}(\mathsf{\lambda })G^{PP}(\mathsf{\epsilon },\mathsf{\lambda })+\\ +\left(\Delta \Phi (\mathsf{\lambda })+{\Sigma }_{phph}(\mathsf{\epsilon },\mathsf{\lambda })\right)G^{uu}(\mathsf{\epsilon },\mathsf{\lambda })\rangle ]\end{array}.$$

(146)

The contribution to the thermodynamic potential from the electrons (in the field of the ionic cores) is also simple to find. It is given by

$${\Omega }_e^{(0)}=-\Theta {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\ln \left(1+e^{({\mathsf{\mu }}_e-\mathsf{\epsilon })/\Theta }\right)} g_e^{(0)}(\mathsf{\epsilon }) d\mathsf{\epsilon }.$$

(147)

Similarly, the contribution to the thermodynamic potential from the phonons (in the field of the ionic cores) is given by

$${\Omega }_{ph}^{(0)}=\Theta {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\ln \left(1-e^{-\mathsf{\epsilon }/\Theta }\right)} g_{ph}^{(0)}(\mathsf{\epsilon }) d\mathsf{\epsilon }.$$

(148)

The values $g_e^{(0)}(\mathsf{\epsilon })$ and $g_{ph}^{(0)}(\mathsf{\epsilon })$ in Equations (148) and (149) are given by Formulas (126)–(130), in which one should put: $t^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}=t^{\mathsf{\lambda }0i}=0,\tilde{G}(\mathsf{\epsilon })=G_0(\mathsf{\epsilon })$ (61)–(65).

Finally, the configurational entropy can be represented as [26]:

$$S_c=-\left[{\displaystyle \sum _{\mathsf{\lambda },m_{\mathsf{\lambda }i},ni}{P}_{ni}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}}\ln P_{ni}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}}}+\frac{1}{2}{\displaystyle \sum _{\begin{array}{l}\mathsf{\lambda },m_{\mathsf{\lambda }i},ni\\ {\mathsf{\lambda }}^{\prime },m_{{\mathsf{\lambda }}^{\prime }j},lj\\ (ni)\ne (lj)\end{array}}{P}_{nilj}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}}\ln \frac{P_{nilj}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}}}{P_{ni}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}}{P}_{lj}^{{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}}}}+\cdots \right]$$

(149)

Ultimately, we are interested in determining the Helmholz free energy, $F$, as a function of the volume $V$,the temperature $T$, the number of electrons N_e, and the parameters of interatomic and magnetic correlations $({\epsilon }_{n_1{i}_1{n}_2{i}_2},\eta )$. The Helmholz free energy can be found directly from the thermodynamic potential. Namely, it satisfies $F=\Omega +{\mu }_e<N_e>$. The free energy per atom, can be approximated by [26]:

$$F=\langle \delta \Phi \rangle -\Theta S_c+{\Omega }_e+{\Omega }_{ph}+{\mu }_e<Z>,$$

(150)

where ${\Omega }_e$ and ${\Omega }_{ph}$ are given by Equations (148) and (149), but with $g_e^{(0)}(\epsilon )$, $g_{ph}^{(0)}(\epsilon )$ replaced by $g_e(\epsilon )$, $g_{ph}(\epsilon )$ (see Equations (126)–(130)).The values of the parameters of the interatomic and magnetic correlations $({\epsilon }_{n_1{i}_1{n}_2{i}_2},\eta )$ are found from the condition for the minimum free energy $F$ (151).

7. Electrical Conductivity

Assuming the system to be driven not too far from equilibrium, we are allowed then to make use of the linear response formalism of Kubo for the electrical conductivity tensor [33],

$${\mathsf{\sigma }}_{\mathsf{\alpha }\mathsf{\beta }}\left(\mathsf{\omega }\right)={\displaystyle \underset{0}{\overset{1/\Theta }{\int }}{\displaystyle \underset{0}{\overset{\infty }{\int }}{e}^{i\mathsf{\omega }t-\mathsf{\delta }t}}\langle {\tilde{J}}_{\mathsf{\beta }}\left(0\right){\tilde{J}}_{\mathsf{\alpha }}\left(t+i\hslash \mathsf{\tau }\right)\rangle d\mathsf{\tau }dt}.$$

(151)

In this equation, $J_{\mathsf{\alpha }}$ is the current operator along the $\mathsf{\alpha }$ spatial direction. The real part of the conductivity, called the optical conductivity, can then be represented in terms of the imaginary part of the retarded response function, or equivalently as

$$\mathrm{Re}{\mathsf{\sigma }}_{\mathsf{\alpha }\mathsf{\beta }}\left(\mathsf{\omega }\right)=\frac{i}{2\mathsf{\omega }}\left[G_r^{J_{\mathsf{\alpha }}{J}_{\mathsf{\beta }}}\left(\mathsf{\omega }\right)-G_a^{J_{\mathsf{\alpha }}{J}_{\mathsf{\beta }}}\left(\mathsf{\omega }\right)\right],$$

(152)

in terms of the retarded and advanced response functions. The current operator is just the number operator for the electrons, multiplied by their velocity and the electric charge, and then summed over all states. It is compactly represented via

$$J_{\mathsf{\alpha }}\left(t\right)=e{\displaystyle \int {\Psi }^{+}\left(\mathsf{\xi },t\right)v_{\mathsf{\alpha }}\Psi \left(\mathsf{\xi },t\right)}d\mathsf{\xi },$$

(153)

where ${\Psi }^{+}\left(\mathsf{\xi },t\right)$ and $\Psi \left(\mathsf{\xi },t\right)$ are the field operators for the creation and annihilation of electrons, respectively, ${\mathsf{\nu }}_{\mathsf{\alpha }}$ is the operator of the $\mathsf{\alpha }$ component of the band velocity, and $e$ is the electron charge. The integration over $\mathsf{\xi }$ sums over all states.

To get the retarded response function on the real frequency axis, we must analytically continue the thermal response functions. The thermal current–current response function is defined to be

$$G^{J_{\mathsf{\alpha }}{J}_{\mathsf{\beta }}}(\mathsf{\tau },{\mathsf{\tau }}^{\prime })=\frac{e^2}{N{V}_1}{\displaystyle \sum _{n_1{n}_2{n}_3{n}_4}{v}_{\mathsf{\alpha }{n}_4{n}_2}{v}_{\mathsf{\beta }{n}_3{n}_1}{G}^{″}}\left(n_1{\mathsf{\tau }}^{\prime },n_2\mathsf{\tau },n_3{\mathsf{\tau }}^{\prime },n_4\mathsf{\tau }\right),$$

(154)

where V₁ is the volume of the primitive unit cell, and the two-particle thermal Green’s function is given by the following time-ordered expectation value:

$$G^{″}(n_1{\mathsf{\tau }}^{\prime },n_2\mathsf{\tau },n_3{\mathsf{\tau }}^{\prime },n_4\mathsf{\tau })=\langle T_{\mathsf{\tau }}{a}_{n_1}({\mathsf{\tau }}^{\prime })a_{n_2}(\mathsf{\tau })a_{n_2}^{+}(\mathsf{\tau })a_{n_4}^{+}(\mathsf{\tau })\mathsf{\sigma }(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathsf{\theta }$}\right.)\rangle {\langle \mathsf{\sigma }(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathsf{\theta }$}\right.)\rangle }^{-1},(n=ni\mathsf{\gamma }).$$

(155)

The two-particle Green’s function from Equation (156) is described by the diagram in Figure 6.

The numbers of Figure 6 correspond to point numbers, e.g., 1 corresponds to $(n_1{i}_1{\mathsf{\gamma }}_1{\mathsf{\tau }}_1).$

Using the diagram technique for two-particle temperature Green’s function and neglecting the contributions of scattering processes on clusters of three or more sites for the conductivity tensor, we can get:

$$\begin{array}{l}\mathrm{Re}{\mathsf{\sigma }}_{\mathsf{\alpha }\mathsf{\beta }}\left(\mathsf{\omega }\right)=\frac{e^2\hslash }{4\pi V_1\mathsf{\epsilon }}\{{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{\mathsf{\epsilon }}_1\left[f\left({\mathsf{\epsilon }}_1+\mathsf{\epsilon }\right)-f\left({\mathsf{\epsilon }}_1\right)\right]}{\displaystyle \sum _{s,s^{\prime }=+,-}(2{\mathsf{\delta }}_{s{s}^{\prime }}-1)}{\displaystyle \sum _{\mathsf{\gamma },i}\{[v_{\mathsf{\beta }}\tilde{K}\left({\mathsf{\epsilon }}_1^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1^{s^{\prime }}+\mathsf{\epsilon }\right)]}\\ +{\displaystyle \sum _{\mathsf{\lambda },m_{\mathsf{\lambda }i}}{P}_{0i}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}}}\tilde{K}({\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon },v_{\mathsf{\beta }},{\mathsf{\epsilon }}_1{}^s)(t^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^s)\tilde{K}({\mathsf{\epsilon }}_1{}^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon })t^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon })\\ +{\displaystyle \sum _{\mathsf{\lambda },m_{\mathsf{\lambda }i}}{P}_{0i}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}}}{\displaystyle \sum _{\begin{array}{l}lj\ne 0i,\\ {\mathsf{\lambda }}^{\prime },m_{{\mathsf{\lambda }}^{\prime }j}\end{array}}{P}_{lj\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0i}^{{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}/\mathsf{\lambda }{m}_{\mathsf{\lambda }i}}[[\tilde{K}({\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon },v_{\mathsf{\beta }},{\mathsf{\epsilon }}_1{}^s)v_{\mathsf{\alpha }}\tilde{G}({\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon })]}\\ \times T^{(2)\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i,{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}({\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon })\\ +[\tilde{K}({\mathsf{\epsilon }}_1{}^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon })v_{\mathsf{\beta }}\tilde{G}({\mathsf{\epsilon }}_1{}^s)]T^{(2)\mathsf{\lambda } m_{\mathsf{\lambda }i}0i,{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}({\mathsf{\epsilon }}_1{}^s)\\ +\tilde{K}({\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon },v_{\mathsf{\beta }},{\mathsf{\epsilon }}_1{}^s)[(t^{{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}({\mathsf{\epsilon }}_1{}^s)\tilde{K}({\mathsf{\epsilon }}_1{}^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon })t^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon })\\ +\left(t^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^s)+t^{{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}({\mathsf{\epsilon }}_1{}^s)\right)\tilde{K}({\mathsf{\epsilon }}_1{}^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon })T^{(2)\mathsf{\lambda } m_{\mathsf{\lambda }i}0i,{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}({\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon })\\ +T^{(2){\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj,\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^s)\tilde{K}({\mathsf{\epsilon }}_1{}^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon })t^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon })\\ +T^{(2){\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj,\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^s)\tilde{K}({\mathsf{\epsilon }}_1{}^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon })T^{(2)\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i,{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}({\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon })\\ +{T^{(2){\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj,\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^s)\tilde{K}({\mathsf{\epsilon }}_1{}^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon })T^{(2){\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj,\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^{s^{\prime }}+\mathsf{\epsilon })]]\}}^{{}^{0i\mathsf{\gamma },0i\mathsf{\gamma }}}\\ +\frac{1}{N}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{\mathsf{\epsilon }}_{1}d{\mathsf{\epsilon }}_2}f\left({\mathsf{\epsilon }}_1\right)f\left({\mathsf{\epsilon }}_2\right)\Delta G_{\mathsf{\alpha }\mathsf{\beta }}^{II}\left({\mathsf{\epsilon }}_1,{\mathsf{\epsilon }}_2;\mathsf{\epsilon }\right)}\},\\ \\ \end{array}$$

(156)

where

$$\begin{gathered} \tilde{K}\left({\mathsf{\epsilon }}_1^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1^{s^{\prime }}+\mathsf{\epsilon }\right)={\tilde{G}}^{a{a}^{+}}\left({\mathsf{\epsilon }}_1^s\right)v_{\mathsf{\alpha }}{\tilde{G}}^{a{a}^{+}}\left({\mathsf{\epsilon }}_1^{s^{\prime }}+\mathsf{\epsilon }\right), \\ {\tilde{G}}^{a{a}^{+}}\left({\mathsf{\epsilon }}_1^{+}\right)={\tilde{G}}_r^{a{a}^{+}}\left({\mathsf{\epsilon }}_1\right), \\ {\tilde{G}}^{a{a}^{+}}\left({\mathsf{\epsilon }}_1^{-}\right)={\tilde{G}}_a^{a{a}^{+}}\left({\mathsf{\epsilon }}_1\right)={\left({\tilde{G}}_r^{a{a}^{+}}\right)}^{*}\left({\mathsf{\epsilon }}_1\right), \end{gathered}$$

(157)

And the two-particle interaction term denoted by $\Delta G^{II}({\mathsf{\epsilon }}_1,{\mathsf{\epsilon }}_2),$ is given by the equation:

$$\begin{gathered} \Delta G_{\alpha \beta }^{II}\left({\epsilon }_1,{\epsilon }_2;\epsilon \right)=\frac{i}{2\pi }{v}_{\alpha n_4{n}_2}{v}_{\beta n_3{n}_1}\{[{\tilde{G}}_{r{n}_1{n}_6}^{a{a}^{+}}\left({\epsilon }_1\right)-{\tilde{G}}_{a{n}_1{n}_6}^{a{a}^{+}}\left({\epsilon }_1\right)]\times \\ \times [{\tilde{G}}_{r{n}_2{n}_5}^{a{a}^{+}}\left({\epsilon }_2\right)-{\tilde{G}}_{a{n}_2{n}_5}^{a{a}^{+}}\left({\epsilon }_2\right)][{\tilde{G}}_{a{n}_7{n}_4}^{a{a}^{+}}\left({\epsilon }_2-\epsilon \right){\tilde{G}}_{r{n}_8{n}_3}^{a{a}^{+}}\left({\epsilon }_1+\epsilon \right)- \\ -{\tilde{G}}_{r{n}_7{n}_4}^{a{a}^{+}}\left({\epsilon }_2-\epsilon \right){\tilde{G}}_{a{n}_8{n}_3}^{a{a}^{+}}\left({\epsilon }_1+\epsilon \right)]+{\tilde{G}}_{a{n}_1{n}_6}^{a{a}^{+}}\left({\epsilon }_1-\epsilon \right)[{\tilde{G}}_{r{n}_2{n}_5}^{a{a}^{+}}\left({\epsilon }_2\right)-{\tilde{G}}_{a{n}_2{n}_5}^{a{a}^{+}}\left({\epsilon }_2\right)]\times \\ \times {\tilde{G}}_{a{n}_7{n}_4}^{a{a}^{+}}\left({\epsilon }_2-\epsilon \right)[{\tilde{G}}_{r{n}_8{n}_3}^{a{a}^{+}}\left({\epsilon }_1\right)-{\tilde{G}}_{a{n}_8{n}_3}^{a{a}^{+}}\left({\epsilon }_1\right)]-{\tilde{G}}_{r{n}_1{n}_6}^{a{a}^{+}}\left({\epsilon }_1-\epsilon \right)\times \\ \times [{\tilde{G}}_{r{n}_2{n}_5}^{a{a}^{+}}\left({\epsilon }_2\right)-{\tilde{G}}_{a{n}_2{n}_5}^{a{a}^{+}}\left({\epsilon }_2\right)]{\tilde{G}}_{r{n}_7{n}_4}^{a{a}^{+}}\left({\epsilon }_2-\epsilon \right)[{\tilde{G}}_{r{n}_8{n}_3}^{a{a}^{+}}\left({\epsilon }_1\right)-{\tilde{G}}_{a{n}_8{n}_3}^{a{a}^{+}}\left({\epsilon }_1\right)]+ \\ +\left[{\tilde{G}}_{a{n}_1{n}_6}^{a{a}^{+}}\left({\epsilon }_1-\epsilon \right){\tilde{G}}_{r{n}_2{n}_5}^{a{a}^{+}}\left({\epsilon }_2+\epsilon \right)-{\tilde{G}}_{r{n}_1{n}_6}^{a{a}^{+}}\left({\epsilon }_1-\epsilon \right){\tilde{G}}_{a{n}_2{n}_5}^{a{a}^{+}}\left({\epsilon }_2+\epsilon \right)\right]\times \\ \times [{\tilde{G}}_{r{n}_7{n}_4}^{a{a}^{+}}\left({\epsilon }_2\right)-{\tilde{G}}_{a{n}_7{n}_4}^{a{a}^{+}}\left({\epsilon }_2\right)][{\tilde{G}}_{r{n}_8{n}_3}^{a{a}^{+}}\left({\epsilon }_1\right)-{\tilde{G}}_{a{n}_8{n}_3}^{a{a}^{+}}\left({\epsilon }_1\right)]+ \\ +[{\tilde{G}}_{r{n}_1{n}_6}^{a{a}^{+}}\left({\epsilon }_1\right)-{\tilde{G}}_{a{n}_1{n}_6}^{a{a}^{+}}\left({\epsilon }_1\right)]{\tilde{G}}_{r{n}_2{n}_5}^{a{a}^{+}}\left({\epsilon }_2+\epsilon \right)[{\tilde{G}}_{r{n}_7{n}_4}^{a{a}^{+}}\left({\epsilon }_2\right)-{\tilde{G}}_{a{n}_7{n}_4}^{a{a}^{+}}\left({\epsilon }_2\right)] \\ \times {\tilde{G}}_{r{n}_8{n}_3}^{a{a}^{+}}\left({\epsilon }_1+\epsilon \right)-[{\tilde{G}}_{r{n}_1{n}_6}^{a{a}^{+}}\left({\epsilon }_1\right)-{\tilde{G}}_{a{n}_1{n}_6}^{a{a}^{+}}\left({\epsilon }_1\right)]{\tilde{G}}_{a{n}_2{n}_5}^{a{a}^{+}}\left({\epsilon }_2+\epsilon \right)\times \\ \times [{\tilde{G}}_{r{n}_7{n}_4}^{a{a}^{+}}\left({\epsilon }_2\right)-{\tilde{G}}_{a{n}_7{n}_4}^{a{a}^{+}}\left({\epsilon }_2\right)]{\tilde{G}}_{a{n}_8{n}_3}^{a{a}^{+}}\left({\epsilon }_1+\epsilon \right)\}{\tilde{\Gamma }}_{n_5{n}_8}^{n_6{n}_7}, \left(n=ni\mathsf{\gamma }\right) \end{gathered}$$

(158)

Summation over repeated indices in expression (159) is implied.

For the static conductivity tensor, we can get (ω→0):

$$\begin{array}{l}{\mathsf{\sigma }}_{\mathsf{\alpha }\mathsf{\beta }}=\frac{e^2\hslash }{4\pi V_1}\{{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}d{\mathsf{\epsilon }}_1\frac{\partial f}{\partial {\mathsf{\epsilon }}_1}}{\displaystyle \sum _{s,s^{\prime }=+,-}(2{\mathsf{\delta }}_{s{s}^{\prime }}-1)}{\displaystyle \sum _{\mathsf{\gamma },i}\{[v_{\mathsf{\beta }}\tilde{K}\left({\mathsf{\epsilon }}_1^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1^{s^{\prime }}\right)]}\\ +{\displaystyle \sum _{\mathsf{\lambda },m_{\mathsf{\lambda }i}}{P}_{0i}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}}}\tilde{K}({\mathsf{\epsilon }}_1{}^{s^{\prime }},v_{\mathsf{\beta }},{\mathsf{\epsilon }}_1{}^s)(t^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^s)\tilde{K}({\mathsf{\epsilon }}_1{}^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1{}^{s^{\prime }})t^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^{s^{\prime }})\\ +{\displaystyle \sum _{\mathsf{\lambda },m_{\mathsf{\lambda }i}}{P}_{0i}^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}}}{\displaystyle \sum _{\begin{array}{l}lj\ne 0i,\\ {\mathsf{\lambda }}^{\prime },m_{{\mathsf{\lambda }}^{\prime }j}\end{array}}{P}_{lj\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0i}^{{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}/\mathsf{\lambda }{m}_{\mathsf{\lambda }i}}[[\tilde{K}({\mathsf{\epsilon }}_1{}^{s^{\prime }},v_{\mathsf{\beta }},{\mathsf{\epsilon }}_1{}^s)v_{\mathsf{\alpha }}\tilde{G}({\mathsf{\epsilon }}_1{}^{s^{\prime }})]}\\ \times T^{(2)\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i,{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}({\mathsf{\epsilon }}_1{}^{s^{\prime }})\\ +[\tilde{K}({\mathsf{\epsilon }}_1{}^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1{}^{s^{\prime }})v_{\mathsf{\beta }}\tilde{G}({\mathsf{\epsilon }}_1{}^s)]T^{(2)\mathsf{\lambda } m_{\mathsf{\lambda }i}0i,{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}({\mathsf{\epsilon }}_1{}^s)\\ +\tilde{K}({\mathsf{\epsilon }}_1{}^{s^{\prime }},v_{\mathsf{\beta }},{\mathsf{\epsilon }}_1{}^s)[(t^{{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}({\mathsf{\epsilon }}_1{}^s)\tilde{K}({\mathsf{\epsilon }}_1{}^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1{}^{s^{\prime }})t^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^{s^{\prime }})\\ +\left(t^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^s)+t^{{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}({\mathsf{\epsilon }}_1{}^s)\right)\tilde{K}({\mathsf{\epsilon }}_1{}^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1{}^{s^{\prime }})T^{(2)\mathsf{\lambda } m_{\mathsf{\lambda }i}0i,{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}({\mathsf{\epsilon }}_1{}^{s^{\prime }})\\ +T^{(2){\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj,\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^s)\tilde{K}({\mathsf{\epsilon }}_1{}^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1{}^{s^{\prime }})t^{\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^{s^{\prime }})\\ +T^{(2){\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj,\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^s)\tilde{K}({\mathsf{\epsilon }}_1{}^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1{}^{s^{\prime }})T^{(2)\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i,{\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj}({\mathsf{\epsilon }}_1{}^{s^{\prime }})\\ +{T^{(2){\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj,\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^s)\tilde{K}({\mathsf{\epsilon }}_1{}^s,v_{\mathsf{\alpha }},{\mathsf{\epsilon }}_1{}^{s^{\prime }})T^{(2){\mathsf{\lambda }}^{\prime }{m}_{{\mathsf{\lambda }}^{\prime }j}lj,\mathsf{\lambda }{m}_{\mathsf{\lambda }i}0i}({\mathsf{\epsilon }}_1{}^{s^{\prime }})]]\}}^{0i\mathsf{\gamma },0i\mathsf{\gamma }}\}\\ \end{array}$$

(159)

The electron velocity satisfies the conventional definition

$$v_{\mathsf{\alpha }}\left(\mathbf{k}\right)=\frac{1}{\hslash }\frac{\partial H_0^{(1)}\left(\mathbf{k}\right)}{\partial k_{\mathsf{\alpha }}}.$$

(160)

When deriving expression (160), the last small term in the expression for electrical conductivity (157) is neglected.

The method developed in this work was applied in [28] to study the effect of an impurity on the energy spectrum and electrical conductivity of carbon nanotubes.

8. Energy Spectrum of Graphene with Adsorbed Potassium Atoms

To calculate the electron spectrum of graphene with adsorbed potassium atoms, we chose the wave functions of the 2s and 2p states of neutral noninteracting carbon atoms as the basis. In the calculation of the matrix elements of the Hamiltonian, we took three first coordination spheres. The energy spectrum of graphene was calculated for the temperature T = 0 K. In calculations, we neglect the renormalization of vertices of the mass operator of the electron–electron interaction. The dependence of the energy of an electron on the wave vector for graphene is calculated from the equation for Green’s function poles for the electron subsystem, defined in Equation (142).

In Figure 7a, we show the dependence of the electron energy ε in graphene with adsorbed potassium atoms on the wave vector k. Vector k is directed from the Brillouin zone center (point Γ) to the Dirac point (point K).

In Figure 7, the structural periodic distance from a potassium atom to a carbon atom is 0.28 nm. It is seen from Figure 7 that, at the ordered arrangement of potassium atoms, a gap in the energy spectrum of graphene arises. Its value depends on the concentration of adsorbed potassium atoms, their location in the unit cell, and the distance to carbon atoms. We established that, at the potassium concentration such that the unit cell includes two carbon atoms and one potassium atom, the latter being placed on the graphene surface above a carbon atom at a distance of 0.286 nm, the energy gap is ~0.25 eV (see Figure 7b). A more complex dependence of the electron energy on the wave vector in the region of the energy gap in comparison with that previously investigated in [34,35,36] in a simple two-band model is due to the effect of band hybridization. The location of the Fermi level in the energy spectrum depends on the potassium concentration and is in the energy interval $-0.36 \mathrm{Ry}\le {\mathsf{\epsilon }}_F\le Ry0.36$. Such a situation is realized if graphene is placed on a potassium support.

9. Conclusions

A novel approach to the description of the electronic spectrum, the thermodynamic potential, and the electrical conductivity of disordered crystals, based on the Hamiltonian of electrons and phonons, constitutes the main issue of the present work. Expressions for Green’s functions, thermodynamic potential, and electrical conductivity are derived using the diagram method. Equations are obtained for the vertex parts of the mass operators of electron–electron and electron–phonon interactions. A system of exact equations is obtained for the spectrum of elementary excitations in a crystal. This makes it possible to perform numerical calculations of the energy spectrum and the properties of the system with a predetermined accuracy. In contrast to other approaches, in which electron correlations are taken into account only in the limiting cases of an infinitely large and infinitesimal electron density, in this method, electron correlations are described in the general case of an arbitrary density.

It was found that a gap appears in the energy spectrum of graphene with an ordered arrangement of potassium atoms. Its value depends on the concentration of adsorbed potassium atoms, their location in the unit cell, and the distance to carbon atoms. It was found that at such a concentration of potassium, the unit cell includes two carbon atoms and one potassium atom, the latter being located on the graphene surface above the carbon atom at a distance of 0.286 nm, and the band gap is ~0.25 eV. Such a situation is realized if graphene is placed on a potassium support. A more complex dependence of the electron energy on the wave vector in the region of the energy gap in comparison with that previously investigated in [34,35,36] in a simple two-band model is due to the effect of band hybridization.