text:
# Phonons in random alloys: the itinerant coherent-potential approximation

## I Introduction

Much research has been carried out over the past few decades on the nature of elementary excitations in disordered alloys. Many aspects of the lattice-vibrational, magnetic and electronic excitations in such systems have been intensively studied both theoretically and experimentally. Of them, the electronic problem has been covered in most detail in recent times with the emergence of first-principles techniques which have made it possible for the theories to attain a much higher degree of accuracy and reliability. Surprisingly, this is not true for phonons despite their being not only conceptually the simplest type of elementary excitation but also the most readily accessible to detailed experiment. From the early 60’s till the early 80’s there were many experimental investigations of phonons in random binary alloys expt ; expt1 ; expt2 ; expt3 ; expt4  by neutron scattering techniques. More recent experiments have been lacking, probably due to the absence of a reliable theory. The feature which makes the theory of phonon excitations difficult is the inseparability of diagonal and off-diagonal disorder. The reason for this is that the force-constant sum rule, i.e. the force constants between a site $`i`$ and its neighbors $`j`$ obey the relation $`𝚽_{ii}=_{ji}𝚽_{ij}`$, must be rigorously satisfied even if the system is disordered. In other words, a single defect at one site in the system perturbs even the diagonal hamiltonian on its neighbors as well, thereby imposing environmental disorder on the force-constants. Hence, any theory must include diagonal, off-diagonal, and environmental disorder as well in order to produce reliable results for phonon excitations in random alloys.

From the late 60’s there were many attempts to provide an adequate theory of phonons in random alloys. The first successful, self-consistent approximation was the coherent potential approximation(CPA) cpa . The CPA is a single-site, mean-field approximation generally capable of dealing only with diagonal disorder (mass disorder in the context of phonons). In the early 70’s there were several studies using the CPA cpa1 ; cpa2  which failed to establish it as a complete answer to the phonon problem in random alloys. The discrepancies with experiment confirmed this need for a theory which could include force-constant changes in addition to mass disorder. Several extensions of the CPA to include off-diagonal and environmental disorder were proposed over the next several years cpa3 ; cpa4 ; cpa5 ; cpa6 ; cpa7 ; cpa8  but only in certain very special cases, such as the separable cpa3  or the additive cpa4 ; cpa5  limits of off-diagonal and environmental disorder, were there successes. The more general approximations cpa6 ; cpa7 ; cpa8  produced Green’s functions which either failed to retain the necessary analytic properties, the translational invariance of the averaged system, or were not fully self-consistent. Moreover, all of these extensions failed to capture the effects of multisite or cluster scatterings which give rise to additional structures in quantities such as the spectral density functions. Later attempts which met with some success for real alloy systems included the recursion method hhk  which can handle large clusters and treats all kinds of disorder on an equal footing. However, the recursion method is neither self-consistent nor translationally invariant when used alone. Yussouf and Mookerjee mook  were able to provide a self-consistent generalization of the CPA to include 2-site scattering using a recursion method in conjunction with the augmented space formalism(ASF)asf . The ASF has the proper translational invariance, yields analytic Green’s functions, and can handle diagonal and off-diagonal disorder.

An alternative approach was provided by Kaplan, Leath, Gray and Diehl klgd (KLGD) which is also based on the ASF. This approach generalized the travelling cluster approximation of Mills and Ratanavararaksha tca  for diagonal disorder to include the other kinds of disorder and multisite effects. Using the diagram symmetry rule of Mills and Ratnavararaksha and the translational symmetry of the augmented-space operators, they presented a self-consistent multiple-scattering theory which allows one to work with a small number of atoms instead of treating large clusters as is done in recursion. It provides analytic, translationally-invariant approximations at all concentrations for diagonal, off-diagonal, and environmental disorder. It can be applied even to problems of charge transfer, lattice relaxation, and short-range order in the context of electronic excitations. However, they illustrated their method only with one-dimensional models and presented it in a very general and complex mathematical language.

In this paper, we present a simple, straightforward formulation of the KLGD method for single-site scattering of phonons in three-dimensional lattices and provide the first application of it to phonons in random alloys . We term it the itinerant coherent-potential approximation or ICPA ; it maintains translational invariance, unitarity, and analyticity of physical properties while including off-diagonal and environmental disorder. In addition to demonstrating its superiority over the single-site CPA and its previous extensions, we provide insight into the physics of force-constant disorder. Our results reveal the complex interplay of forces between various atomic species in a random environment, an important phenomenon which has never been addressed properly.

In Section II we describe the theory, introducing the augmented-space representation and its use in constructing the self-consistent scattering theory and the single-site itinerant coherent-potential approximation. In Section III we derive expressions for important physical quantities such as densities of states, spectral functions, inelastic scattering cross-sections and their sum rules in terms of the configuration-averaged Green’s function of the system. In Sections IV and V we present our results on Ni55Pd45 and Ni50Pt50 alloys as test cases and compare them with experiment. Concluding remarks are presented in Section VI.

## II Formalism

In this section, we briefly sketch the rationale behind augmented space, introduce its representations, and define the notation to be used throughout the paper. We present our discussions here only in the context of phonons. The formulation of the ICPA for other kinds of excitations is closely analogous.

### II.1 Augmented space and its representations

The description of disordered systems conventionally proceeds as follows: the dynamical behavior of a system is described by a Hamiltonian, whereas the statistical behavior of the disorder is imposed from outside. The Hamiltonian itself does not describe the full behavior of the random system, but has to be augmented with the distribution of the set of random potentials which are associated with the various configurations of the system. The physical properties are then obtained by ensemble averages over configurations. The CPA and its extensions employ this procedure.

An alternative procedure is that instead of looking at the excitations of the system as moving in a random array of disordered potentials, the excitations are considered to be moving in periodic potentials in the presence of a ‘field’ which specifies the disorder. The Hamiltonian, expanded to include the disorder field, then by itself completely describes the disordered system. Since the information on random configurations is already incorporated into the Hamiltonian, the configuration averaging is not a further process as in the mean-field approaches, but simply an evaluation of matrix elements. The idea of introducing a ‘disorder field’ to describe the random fluctuations in the system by extending the Hilbert space to include the disorder field and by representing the Hamiltonian in this new space constitutes the core of the augmented-space formalism. The extended Hilbert space which captures the random fluctuations is called the ‘augmented space’.

Here, we work only with a binary alloy A$`_{c_A}`$B$`_{c_B}`$. We assume that each lattice site is randomly occupied by an A atom or by a B atom. We wish to calculate the configuration-averaged values of the experimentally measurable physical quantities, for which we need a configuration-averaged Green’s function. In particular, we shall concentrate here on the configuration-averaged displacement-displacement (one-phonon) Green’s function rmp

$$G_{nm}^{\alpha \beta }\left(t\right)=\frac{1}{ih}u_n^\alpha \left(t\right);u_m^\beta \left(0\right),$$
(1)
or, after Fourier transformation to the frequency domain

$$𝐆\left(\omega ^2\right)=[𝐦\omega ^2𝚽]^1.$$
(2)
In Eqs.(1) and (2) $``$ stands for both configuration and thermodynamic averaging. In Eq.(1), $`m,n`$ specify lattice sites and $`\alpha \beta `$ the cartesian directions. $`u_n^\alpha (t)`$ is the displacement operator of an atom at the lattice site $`n`$ in the direction $`\alpha `$ at the time $`t`$. In Eq.(2) a bold symbol represents a matrix for which all indices are to be understood. The semi-colon ; denotes Bose time ordering. $`𝐦`$ is the mass operator, $`𝚽`$ is the force-constant operator, and $`\omega `$ is the frequency which contains a vanishingly small negative imaginary part. The masses are random,

$$m_{ij}^{\alpha \beta }=m_i\delta _{\alpha \beta }\delta _{ij},$$
(3)
with $`m_i`$ randomly taking on the value $`m^\mathrm{\Gamma }`$ if species $`\mathrm{\Gamma }`$=A,B is on site $`i`$. The force-constants take on the values $`\left(\varphi _{ij}^{\alpha \beta }\right)^{\mathrm{\Gamma }\mathrm{\Delta }}`$ if species $`\mathrm{\Gamma }`$ is on site $`i`$ and species $`\mathrm{\Delta }`$ is on site $`j`$.

It is $`𝐆`$ which carries all the dynamical informations of interest, and the essential difficulty of the theory of phonons in random systems arises from taking the configuration average of the inverse of the matrix $`𝐦\omega ^2𝚽`$. The augmented-space technique asf ; km  greatly facilitates this averaging. The displacements $`𝐮`$, masses $`𝐦`$, force-constants $`𝚽`$, and Greens function $`𝐆`$ are defined in the dynamical Hilbert space $`\mathrm{\Psi }`$ in which the Hamiltonian of the system operates. For a binary alloy, $`\mathrm{\Psi }`$ is augmented by the space $`\mathrm{\Theta }`$ of all possible atomic configurations of the system. The resulting augmented space $`\mathrm{\Omega }`$ is

$`\mathrm{\Omega }=\mathrm{\Psi }\mathrm{\Theta },`$
In $`\mathrm{\Omega }`$ or $`\mathrm{\Theta }`$ operators are represented by symbols with superposed carets. In the configuration representation within $`\mathrm{\Theta }`$, the state of site $`i`$ is specified by the single-site state $`|A_i`$ if A is on $`i`$ and by $`|B_i`$ if B is on $`i`$. With respect to these states, the occupation operators $`\widehat{\eta }_i^\mathrm{\Gamma }`$, $`\mathrm{\Gamma }`$=A,B,

$`\widehat{\eta }_i^A|A_i`$ $`=`$ $`|A_i,\widehat{\eta }_i^A|B_i=0,`$
$`\widehat{\eta }_i^B|B_i`$ $`=`$ $`|B_i,\widehat{\eta }_i^B|A_i=0`$ (4)
are represented by the matrices

$$\widehat{\eta }_i^A=\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right),\widehat{\eta }_i^B=\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right)=\widehat{𝐈}_i\widehat{\eta }_i^A.$$
(5)
The configuration of the entire system is specified by the direct product of all single-site states $`_i|\mathrm{\Gamma }_i`$, $`\mathrm{\Gamma }`$=A,B. The mass operator for site $`i`$ is given by,

$$\widehat{𝐦}_i^{}=m^A\widehat{\eta }_i^A+m^B\widehat{\eta }_i^B.$$
(6)
Similarly the force-constants for sites $`i`$ and $`j`$ are given by

$`\widehat{𝚽}_{ij}^{}`$ $`=`$ $`\varphi _{ij}^{AA}\widehat{\eta }_i^A\widehat{\eta }_j^A+\varphi _{ij}^{AB}\widehat{\eta }_i^A\widehat{\eta }_j^B`$ (7)
$`\text{ }+\varphi _{ij}^{BA}\widehat{\eta }_i^B\widehat{\eta }_j^A+\varphi _{ij}^{BB}\widehat{\eta }_i^B\widehat{\eta }_j^B,`$
with the cartesian indices understood.

Consider now a rotated representation for site $`i`$ in which the basis vectors for its configuration space are given by

$`|0_i`$ $`=`$ $`\sqrt{c_A}|A_i+\sqrt{c_B}|B_i,`$
$`|1_i`$ $`=`$ $`\sqrt{c_B}|A_i\sqrt{c_A}|B_i.`$ (8)
Constructing the configuration average of any operator $`\widehat{𝐀}`$ in $`\mathrm{\Theta }`$ can be carried out simply by taking the expectation value of $`\widehat{𝐀}`$ with the state

$$|f=\underset{i}{}|0_i,$$
(9)
Thus $`|0_i`$ is the site-average state (or the virtual-crystal state), $`|1_i`$ describes a fluctuation away from the average state on site $`i`$, and

$$|f_i=|1_i\underset{ji}{}|0_j.$$
(10)
is the state in which there is a fluctuation or a defect in the average state $`|f`$ only on site $`i`$. In this fluctuation representation the occupation operators $`\widehat{\eta }_i^A`$ and $`\widehat{\eta }_i^B`$ are transformed to

$`\widehat{\eta }_i^A`$ $`=`$ $`\left(\begin{array}{cc}c_A& \sqrt{c_Ac_B}\\ \sqrt{c_Ac_B}& c_B\end{array}\right)`$ (13)
$`\widehat{\eta }_i^B`$ $`=`$ $`\left(\begin{array}{cc}c_B& \sqrt{c_Ac_B}\\ \sqrt{c_Ac_B}& c_A\end{array}\right).`$ (16)
In transforming from the configuration representation to the fluctuation representation, $`\widehat{𝐦}^{}`$ goes to $`\widehat{𝐦}`$ and $`\widehat{𝚽}^{}`$ to $`\widehat{𝚽}`$, as given by Eqs.(6) and (7), respectively, with the $`\widehat{\eta }^\mathrm{\Gamma }`$ of Eq.(5) replaced by the $`\widehat{\eta }^\mathrm{\Gamma }`$ of Eq.(11). Thus the dynamical operators $`\widehat{𝐦}`$ and $`\widehat{𝚽}`$ are not diagonal with respect to the number of fluctuations or defects in the fluctuation representation and can create them, destroy them, or, in the case of $`\varphi _{ij}`$, cause them to travel, or ‘itinerate’. We refer to the movement of defects induced by the off-diagonal elements of the $`\widehat{\eta }_i^\mathrm{\Gamma }`$ as the itineration of fluctuations to distinguish it from the propagation of phonons. However, these operators are translationally invariant; the randomness in configuration is thus captured by translationally invariant operators in the configuration space $`\mathrm{\Theta }`$. The $`\widehat{\eta }^\mathrm{\Gamma }`$ operators constitute the disorder field referred to above.

Any operator $`\widehat{𝐀}`$ in this augmented space can be represented in block form,

$$\widehat{𝐀}=\left(\begin{array}{cc}\overline{𝐀}& 𝐀^{}\\ 𝐀_{}^{}{}_{}{}^{}& \stackrel{~}{𝐀}\end{array}\right),$$
(17)
where the bold notation A implies a matrix in the site and cartesian indices. The four elements of the block matrix are given by

$`\overline{𝐀}`$ $`=`$ $`𝐏\widehat{𝐀}𝐏,`$
$`𝐀^{}`$ $`=`$ $`𝐏\widehat{𝐀}\left(\mathrm{𝟏}𝐏\right),`$
$`𝐀_{}^{}{}_{}{}^{}`$ $`=`$ $`\left(\mathrm{𝟏}𝐏\right)\widehat{𝐀}𝐏,`$
$`\stackrel{~}{𝐀}`$ $`=`$ $`\left(\mathrm{𝟏}𝐏\right)\widehat{𝐀}\left(\mathrm{𝟏}𝐏\right),`$ (18)
where $`𝐏`$, the projection operator onto the virtual-crystal state, is given by $`𝐏=|ff|`$ . Thus,we see that $`\overline{𝐀}`$ is the configuration average of the quantity $`\widehat{𝐀}`$ while $`𝐀^{}`$,$`𝐀_{}^{}{}_{}{}^{}`$ generate the coupling between the average and the fluctuation states and $`\stackrel{~}{𝐀}`$ is that part of $`𝐀`$ entirely within the space of fluctuation states.

In the present paper we shall make the approximation of treating explicitly only single fluctuation states $`|f_l`$ in the fluctuation space $`\mathrm{\Theta }|ff|`$, although multiple-fluctuation states are treated implicitly via a self-consistency condition. States in $`\mathrm{\Omega }`$ can then be specified by $`|if`$ or $`|if_l`$ where $`i`$ is the site index of the dynamical variable in $`\mathrm{\Psi }`$, position or momentum, with the cartesian index understood. For the site indices of the corresponding matrix elements we shall often use the compact notation

$`if|\widehat{𝐀}|jf`$ $`=`$ $`\overline{A}_{ij},`$
$`if_l|\widehat{𝐀}|jf_l^{}`$ $`=`$ $`\stackrel{~}{A}_{ij}^{(l)(l^{})},`$
$`if|\widehat{𝐀}|jf_l`$ $`=`$ $`A_{ij}^{(l)},`$
$`if_l|\widehat{𝐀}|jf`$ $`=`$ $`A_{ij}^{(l)},`$ (19)
where $`l`$ and $`l^{}`$ denote the locations of the concentration fluctuation or defect. The parentheses around $`l`$ indicate that it is neither a site nor a cartesian direction index, but indicates instead the position of a fluctuation in the lattice.

### II.2 Multiple-scattering picture

A phonon propagating in a random alloy undergoes irreducible multiple scattering gonis  both repeatedly off a single fluctuation and successively off fluctuations on the different sites it encounters in the process. The CPA takes into account the former but not the latter. To illustrate how the treatment of this process of multiple scattering by fluctuations differs between the CPA and our formalism we employ a cartoon diagram (Fig.1). The top panel, a two-dimensional cross-section, illustrates the multiple-scattering process included in the CPA. There, the filled circle is a single ‘fluctuation site’ immersed in an average medium denoted by open circles. The arrow on the left is the direction of phonon propagation. When the phonon meets the fluctuation site, it undergoes irreducible multiple scattering at that site. In the CPA(diagonal disorder), the irreducible scattering by the defect site is confined to the defect site. The circle around the fluctuation site indicates the region of influence of the perturbation. None of the springs are affected by the presence of this defect since the force-constants are the same everywhere. One does an averaging over all the possible occupations of the single site. The phonon diagrams of the self-energy which describe this multiple scattering process completely are shown in Fig.2(a). There, the filled circles represent the fluctuation sites, the dotted lines represent successive scatterings from the fluctuation site, and the double solid line represents the self-consistent propagator.

The lower three panels in Fig.1 illustrate scattering sites in the ICPA. The difference from the CPA is that the region of influence is not only the site of fluctuation but also its neighboring environment around the fluctuation site. The figure shows an example (dotted contour) where the environment includes nearest neighbors only (The calculations could be extended to further neighbors as well). When the phonon interacts with the fluctuation site in the top panel of the three, it scatters also from all of its neighbors since their spring constants also undergo changes (denoted by the thick spring lines in contrast to the thin ones for the average medium). The whole cluster of atoms undergoes fluctuations in force-constants as the occupation of the fluctuation site changes. One has to keep in mind that the force-constant between the fluctuation site and its neighbor on the right, say, depends on the occupation of both sites, as is true for the next neighboring site on its right as well. So, one is led to include the irreducible scatterings by the fluctuation on all neighboring sites, which then requires inclusion of scattering by the fluctuations on its neighbors etc., until the irreducible scatterings extend throughout the entire sample. A simple example of this process is indicated in the middle and bottom panels. Indeed, Mills and Ratanavararakshatca  have shown that once there are non-diagonal terms in the scattering, the self-energy must include these migrations (itinerations) of the scatterer throughout the sample in order to attain unitarity and thereby guarantee that the average Green’s function will be properly analytic or Herglotz. The self-consistent scattering and the resulting coherent potential about a single defect thus itinerates from defect to defect throughout the sample, making it an itinerant coherent potential. The scattering could have started from any site in the sample so that the result is also fully translationally invariant, and the self-energy is $`q`$-dependent but diagonal in the $`q`$-space of the Brillouin zone of the underlying periodic lattice structure. Fig. 2(b) illustrates a typical self-energy diagram in the ICPA. The solid and dotted overlapping ellipses denote the multiple scattering by a single-site and its neighbors, i.e. by a cluster of atoms, and the subsequent itineration of this process. The thin dotted lines and the thick double lines are as in Fig.2(a).

In the multiple scattering framework, we calculate the self-energy $`𝚺\left(\omega ^2\right)`$, defined by

$$𝐆\left(\omega ^2\right)=[𝐆_{vca}^1\left(\omega ^2\right)𝚺\left(\omega ^2\right)]^1,$$
(20)
where $`𝐆_{vca}`$ is the unperturbed Green’s function,

$$𝐆_{vca}=\left(\overline{𝐦}\omega ^2\overline{𝚽}\right)^1,$$
(21)
and $`\overline{𝐦}`$ and $`\overline{𝚽}`$ are the configuration-averaged mass and force-constant operators respectively.

Our major task is to calculate the self-energy $`𝚺\left(\omega ^2\right)`$. Let us consider $`\widehat{𝐊}=\left(\widehat{𝐦}\omega ^2\widehat{𝚽}\right)=\widehat{𝐆}^1`$. Using the 2$`\times `$2 block representation of augmented-space operators of Eq.(12) we get

$$\widehat{𝐆}=\left(\begin{array}{cc}\overline{𝐆}& 𝐆^{}\\ 𝐆_{}^{}{}_{}{}^{}& \stackrel{~}{𝐆}\end{array}\right)=\left(\begin{array}{cc}\overline{𝐊}& 𝐊^{}\\ 𝐊_{}^{}{}_{}{}^{}& \stackrel{~}{𝐊}\end{array}\right)^1.$$
(22)
Using the relation for the inverse of an operator in 2$`\times `$2 block form mat , namely,

$$\widehat{𝐀}^1=\left(\begin{array}{cc}\left(\overline{𝐀}𝐀^{}\stackrel{~}{𝐀}^1𝐀^{}\right)^1& \left(\stackrel{~}{𝐀}𝐀_{}^{}{}_{}{}^{1}\overline{𝐀}𝐀^{}\right)^1\\ \left(\overline{𝐀}𝐀_{}^{}{}_{}{}^{1}\stackrel{~}{𝐀}𝐀^{}\right)^1& \left(\stackrel{~}{𝐀}𝐀^{}\overline{𝐀}^1𝐀^{}\right)^1\end{array}\right),$$
(23)
we get

$`\overline{𝐆}`$ $`=`$ $`[\left(\overline{𝐦}\omega ^2\overline{𝚽}\right)𝐊^{}\left(\stackrel{~}{𝐦}\omega ^2\stackrel{~}{𝚽}\right)^1𝐊^{}]^1`$ (24)
$`=`$ $`[𝐆_{vca}^1𝐊^{}\{𝐆_{vca}^1[(\overline{𝐦}\omega ^2\overline{𝚽})`$
$`(\stackrel{~}{𝐦}\omega ^2\stackrel{~}{𝚽})]\}^1𝐊^{}]^1`$
$`=`$ $`[𝐆_{vca}^1𝐊^{}\mathrm{𝐅𝐊}^{}]^1.`$
Therefore, the self-energy is given by

$$𝚺=𝐊^{}\mathrm{𝐅𝐊}^{},$$
(25)
where

$$𝐅=\stackrel{~}{𝐊}^1=\{𝐆_{vca}^1\stackrel{~}{𝐈}\stackrel{~}{𝐕}\}^1,$$
(26)
and where

$$\stackrel{~}{𝐕}=\left(\overline{m}\stackrel{~}{𝐈}\stackrel{~}{𝐦}\right)\omega ^2\left(\overline{\varphi }\stackrel{~}{𝐈}\stackrel{~}{𝚽}\right).$$
(27)
The quantity $`\stackrel{~}{𝐕}`$ denotes all perturbations to the average medium, and $`𝐅`$ contains the itineration of the fluctuation in the average medium.

Up to this point, the scattering formalism is exact. We now introduce the ICPA by restricting the states within the configuration space $`\mathrm{\Theta }|ff|`$ to the single-fluctuation states, the notation for which is given by Eq.(14). Making the site and cartesian indices explicit, we obtain for $`𝚺`$ in Eq.(20), under this restriction,

$$\mathrm{\Sigma }_{ij}^{\alpha \beta }=K_{\alpha i,\delta k}^{(m)}F_{\delta k,\gamma l}^{(m)(n)}K_{\gamma l,\beta j}^{(n)}.$$
(28)
The sumations are over the repeated indices, and the fluctuation itinerator $`𝐅`$ is given by a Dyson equation,

$$𝐅^{(i)(j)}=𝐆_{vca}[\delta _{(i)(j)}+\underset{l}{}\stackrel{~}{𝐕}^{(i)(l)}𝐅^{(l)(j)}],$$
(29)
where only the site index of the fluctuation is shown. The quantities in (23) are translationally invariant as follows:

$`K_{ik}^{(m)}`$ $`=`$ $`K_{im,km}^{(0)},`$
$`F_{kl}^{(m)(n)}`$ $`=`$ $`F_{km,lm}^{(0)(nm)}.`$ (30)
The single fluctuation in Eq.(23) can be considered to have been ‘created’ by $`𝐊^{(n)}`$ at site $`n`$, itinerated to site $`m`$ by $`𝐅^{(n)(m)}`$ and ‘destroyed’ by $`𝐊^{(m)}`$ at site $`m`$. The $`𝐊`$, $`𝐊^{}`$ and $`𝐅`$ matrices have elements which are non-zero only for site indices within the environment of the appropriate defects i.e. the indices $`i`$ and $`k`$ ($`l`$ and $`j`$) must be within the neighborhood perturbed by the defect at $`m`$($`n`$). The terms with more than one fluctuation(defect) present at a time correspond to coherent pair and ‘defect cluster’ scattering and are neglected in the single-site scattering considered here. All of these operators act in the augmented space. The Equations (20)-(24) define an itinerant single-site multiple scattering theory.

### II.3 Self-consistency

The restriction in Eq.(23) to states of $`\mathrm{\Theta }|ff|`$ containing only a single fluctuation is a very severe approximation. Multiple-fluctuation states are of course present in $`𝐅`$ and contribute to $`𝚺`$. In the spirit of the CPA, these are included approximately by introducing self-consistency. As in the CPA klgd ; rmp  we obtain self-consistency by replacing $`𝐆_{vca}`$ in $`𝐅`$ in Eq.(24) by a conditional propagator $`𝐆^{(i)}`$, identical to $`\overline{𝐆}=𝐆`$ except that all irreducible scatterings beginning or ending on site $`i`$ are omitted, so that $`𝐅`$ would then be given by

$$𝐅^{(i)(j)}=𝐆^{(i)}[\delta _{(i)(j)}+\underset{l}{}\stackrel{~}{𝐕}^{(i)(l)}𝐅^{(l)(j)}].$$
(31)
In parallel with Eq.(15), $`𝐆^{(i)}`$ contains a conditional self-energy $`𝚺^{(i)}`$ which, is like Eq.(23), except that it includes only those scatterings that neither start nor end on $`i`$,

$$𝐆^{(i)}=[\left(𝐆_{vca}\right)^1𝚺^{(i)}]^1,$$
(32)
$$𝚺^{(i)}=\underset{l,mi}{}𝐊^{(l)}𝐅^{(l)(m)}𝐊^{(m)}.$$
(33)
Referring to Fig.2a, the double line in the multiple-scattering graphs is the propagator $`𝐆^{(i)}`$ when the solid dot refers to site $`i`$. We obtain the itinerant CPA by allowing $`𝐊^{}`$,$`𝐊^{}`$, and $`\stackrel{~}{𝐕}`$ to include force-constant disorder as well and therefore defect itineration in Equations (26)-(28). This closed set of equations defines our single-site, self-consistent, multiple-scattering theory which, when solved, yields $`𝐅`$. Inserting $`𝐅`$ into Eq.(20) for $`𝚺`$ and the result into Eq.(15) then yields $`𝐆`$. It is already known that Eqs.(15),(20),(26)-(28) have a unique solution which yields a Herglotz average Green’s function klgd . A major difference between this and previous generalizations of the CPA is that for scattering from single-site fluctuations with off-diagonal and/or environmental disorder, as is considered here, the matrix representation of the operator $`\stackrel{~}{𝐕}`$ has elements which transfer or itinerate the fluctuation from site to site. This feature causes the self-energy to have nonzero off-diagonal elements in real-space extending across the sample and thus contributes importantly to such quantities as the two-particle vertex corrections in a way the CPA cannot chit .

It now remains to solve these equations, making use of the translational symmetry of the augmented-space operators. We accomplish this by Fourier transforms on the fluctuation-site labels,

$$A\left(\stackrel{}{q}\right)_{mn}=N^1\underset{l,l^{}}{}A_{l+m,l^{}+n}^{(l)(l^{})}e^{i\stackrel{}{q}\stackrel{}{R}_{ll^{}}},$$
(34)
and

$$A_{l+m,l^{}+n}^{(l)(l^{})}=N^1\underset{\stackrel{}{q}}{}A\left(\stackrel{}{q}\right)_{mn}e^{i\stackrel{}{q}\stackrel{}{R}_{ll^{}}},$$
(35)
where $`\stackrel{}{R}_{ll^{}}`$ is the lattice vector connecting the fluctuation sites $`l`$ and $`l^{}`$, $`m`$ and $`n`$ are neighbors of $`l`$ and $`l^{}`$ respectively, and the $`\stackrel{}{q}`$ sum is over the Brillouin zone.

We can also effect Fourier transforms on the site indices themselves. That of the self-energy is

$$\mathrm{\Sigma }\left(\stackrel{}{q}\right)=N^1\underset{ij}{}\mathrm{\Sigma }_{ij}e^{\stackrel{}{q}\stackrel{}{R}_{ij}}.$$
(36)
From Eqs.(20), (25) and (29),it follows that

$$\mathrm{\Sigma }\left(\stackrel{}{q}\right)=\underset{l,m,n,p}{}K_{lm}^{(0)}F\left(\stackrel{}{q}\right)_{mn}K_{np}^{(0)}e^{i\stackrel{}{q}\stackrel{}{R}_{lp}}.$$
(37)
In this notation, Eq (26) becomes

$$F\left(\stackrel{}{q}\right)_{mn}=G_{mn}^{(0)}+\underset{rp}{}G_{mr}^{(0)}\stackrel{~}{V}\left(\stackrel{}{q}\right)_{rp}F\left(\stackrel{}{q}\right)_{pn}.$$
(38)
The cartesian indices here are implicit so that each quantity is a 3$`\times `$3 matrix .

Since the range of interaction in real-space is finite, the perturbations $`𝐊^{(i)}`$ and $`\stackrel{~}{𝐕}\left(\stackrel{}{q}\right)`$ are finite matrices, nonzero only over a finite set of real sites. For example, if we consider nearest-neighbor perturbation only in a single-site approximation, the $`\stackrel{~}{𝐕}\left(\stackrel{}{q}\right)`$ and $`𝐊^{}`$ are 3$`(Z+1)\times `$ 3 $`\left(Z+1\right)`$ matrices where $`Z`$ is the number of nearest neighbors. This is the minimum matrix size necessary to exhibit all the impurity modes or states about each fluctuation site.

In full matrix notation, we obtain

$$𝐅\left(\stackrel{}{q}\right)=[𝐆_{}^{(0)}{}_{}{}^{1}\stackrel{~}{𝐕}\left(\stackrel{}{q}\right)]^1.$$
(39)
These matrices , for example, for an fcc lattice, are of dimension 39$`\times `$39.

In order to evaluate $`𝐆^{(0)}`$ we rewrite Eq. (27) as,

$$𝐆^{(0)}=[\left(𝐆_{vca}\right)^1𝚺^{(0)}]^1=[𝐆^1+\stackrel{~}{𝚺}^{(0)}]^1,$$
(40)
where, $`\stackrel{~}{𝚺}^{(0)}=\left(𝚺𝚺^{(0)}\right)`$. The conditional self-energy $`\stackrel{~}{𝚺}^{(0)}`$ contains only those scatterings which either start or end with a perturbation caused by a fluctuation at site $`0`$. Thus, to evaluate the self-consistent propagator $`𝐆^{(0)}`$, we need to know $`𝐆`$. But $`𝐆`$ is obtained from Eq.(15), which becomes

$`G\left(\stackrel{}{q}\right)`$ $`=`$ $`[G_{vca}\left(\stackrel{}{q}\right)^1\mathrm{\Sigma }\left(\stackrel{}{q}\right)]^1,`$
$`G_{ij}`$ $`=`$ $`N^1{\displaystyle \underset{\stackrel{}{q}}{}}G\left(\stackrel{}{q}\right)e^{i\stackrel{}{q}\stackrel{}{R}_{ij}}.`$ (41)
After reaching self-consistency by the procedure described below, we use these expressions to calculate densities of states and spectral functions.

The conditional self-energy $`\stackrel{~}{𝚺}^{(0)}`$ can be broken up into two contributions:

(i) Scattering that starts from a defect at site $`0`$ and ends at site $`j`$.

(ii) Scattering that starts at $`j`$ but ends at $`0`$.

This decomposition results in

$`\stackrel{~}{𝚺}^{(0)}`$ $`=`$ $`{\displaystyle \underset{j}{}}[𝐊^{(0)}𝐅^{(0)(j)}𝐊^{(j)}+𝐊^{(j)}𝐅^{(j)(0)}𝐊^{(0)}]`$ (42)
$`𝐊^{(0)}𝐅^{(0)(0)}𝐊^{(0)}.`$
The last term is subtracted to avoid overcounting when $`j`$=$`0`$.

In a block notation similar to that of Eq.(12), we have

$$\stackrel{~}{𝚺}^{(0)}=\left(\begin{array}{cc}𝚺_1& 𝚺_3\\ 𝚺_3^{}& 0\end{array}\right),$$
(43)
$$𝐆^{(0)}=\left(\begin{array}{cc}𝐆_1^{(0)}& 𝐆_3^{(0)}\\ 𝐆_3^{(0)}& 𝐆_2^{(0)}\end{array}\right),$$
(44)
$$𝐆=\left(\begin{array}{cc}𝐆_1& 𝐆_3\\ 𝐆_3^{}& 𝐆_2\end{array}\right).$$
(45)
where, for a genaral operator $`\widehat{𝐀}`$, $`𝐀_1`$ begins and ends with scattering about site $`0`$ , $`𝐀_2`$ neither begins nor ends with scattering about site $`0`$ and $`𝐀_3(𝐀_3^{})`$ begins(ends) with scattering at the site $`0`$ and ends(begins) with scattering about a site different from $`0`$. The term $`𝚺_2`$ is $`0`$ since $`\stackrel{~}{𝚺}^{(0)}`$ must begin or end at the site $`0`$. From Eq. (35), we have

$$𝐆^{(0)}=𝐆(𝐈+\stackrel{~}{𝚺}^{(0)}𝐆)^1,$$
(46)
which leads to

$`𝐆_1^{(0)}`$ $`=`$ $`\overline{𝐗}[𝐈+\left(𝚺_1𝚺_3𝐆_2𝚺_3^{}\right)\overline{𝐗}+𝚺_3𝐆_3^{}]^1,`$ (47)
where

$$\overline{𝐗}=\left(𝐈+𝐆_3𝚺_3^{}\right)^1𝐆_1.$$
(48)
after a lengthy algebraic analysis which was previously given in Ref.18.

In order to evaluate these expressions klgd , we need to calculate four terms: $`𝐆_1`$,$`𝚺_1`$, $`𝐆_3𝚺_3^{}`$ and $`𝚺_3𝐆_2𝚺_3^{}`$. The first term $`𝚺_1`$ is just a finite sum of finite matrices and can be evaluated directly, but the other two terms involve sums which range over all sites in the solid and must be evaluated by Fourier transforms. This is done in the following way,

$`\left(G_3\mathrm{\Sigma }_3^{}\right)_{t,t^{}}`$ $`=`$ $`{\displaystyle \underset{m}{}}{\displaystyle \underset{r,n,l}{}}G\left(\omega ^2\right)_{t,m+r}K_{m+r,m+n}^{(m)}`$
$`\times `$ $`F_{m+n,l}^{(m)(0)}K_{l,t^{}}^{(0)}{\displaystyle \underset{r}{}}G\left(\omega ^2\right)_{t,r}\stackrel{~}{\mathrm{\Sigma }}_{r,t^{}}^{(0)},`$
which becomes

$`\left(G_3\mathrm{\Sigma }_3^{}\right)_{t,t^{}}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{\stackrel{}{q}}{}}{\displaystyle \underset{r}{}}G\left(\stackrel{}{q}\right)e^{i\stackrel{}{q}\stackrel{}{R}_{tr}}M\left(\stackrel{}{q}\right)_{r,t^{}}`$ (49)
$`\text{ }\left(G_1\mathrm{\Sigma }_1\right)_{t,t^{}},`$
and,similarly,

$`\left(\mathrm{\Sigma }_3G_2\mathrm{\Sigma }_3^{}\right)_{t,t^{}}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{\stackrel{}{q},r,r^{}}{}}M\left(\stackrel{}{q}\right)_{tr}G\left(\stackrel{}{q}\right)e^{i\stackrel{}{q}\stackrel{}{R}_{rr^{}}}M\left(\stackrel{}{q}\right)_{r^{},t^{}}`$ (50)
$`(\mathrm{\Sigma }_1G_1\mathrm{\Sigma }_1+\mathrm{\Sigma }_1G_3\mathrm{\Sigma }_3^{}+\mathrm{\Sigma }_3G_3^{}\mathrm{\Sigma }_1)_{t,t^{}},`$
where

$$M\left(\stackrel{}{q}\right)_{r,t^{}}=\underset{nl}{}K_{rn}^{(0)}F\left(\stackrel{}{q}\right)_{nl}K_{lt^{}}^{(0)}.$$
(51)
In these equations, $`0`$ is the index of the single fluctuation-site in consideration; $`r,r^{},t,t^{},l,n`$ are the neighboring sites of $`0`$; and $`m,m^{}`$ are general sites in the sample. So, it is clear that one needs to work only on matrices of size 3$`(Z+1)\times `$3 $`\left(Z+1\right)`$ and use the Fourier transform of operators to handle the itineration of the fluctuation throughout the entire sample. An interesting point to note is that the quantities $`G_3\mathrm{\Sigma }_3^{}`$ and $`\mathrm{\Sigma }_3G_2\mathrm{\Sigma }_3^{}`$ represent the scattering and itineration of the disturbance including the effect of the off-diagonal and environmental disorder. In case of diagonal-disorder only, they vanish giving $`𝐆_1^{(0)}=𝐆_1\left(𝐈+𝚺_1\right)^1`$, which is the CPA self-consistent propagator, and the self-consistent set of equations reduces to the CPA equations.

The inputs to the self-consistency cycle are $`𝐆_{start}^{(0)}=𝐆_{vca}`$(or some better guess) , $`𝐊^{},𝐊^{}`$ and $`\stackrel{~}{𝐕}\left(\stackrel{}{q}\right)`$. The procedures for evaluating the latter three quantities are given in the Appendix. The cycle consists of the following steps :

1. calculation of $`𝐅\left(\stackrel{}{q}\right)`$ using Eqs.(33) and (34).
2. calculation of $`\mathrm{\Sigma }\left(\stackrel{}{q}\right)`$ using Eq.(32).
3. calculation of $`G\left(\stackrel{}{q}\right)`$ and $`G\left(\omega ^2\right)`$ using Eq.(36).
4. calculation of $`𝐆_1^{(0)}`$ using Eqs.(42),(43),(44) and (45)
5. If the results of steps 1.-4. are acceptably close to those of the previous cycle, stop. If not, use as input to step 1 and iterate.

The iterations are done till self-consistency is achieved for each $`\stackrel{}{q}`$-point in the Brillouin zone. In the process of achieving self-consistency, one calculates $`G`$ in both real-space and in $`\stackrel{}{q}`$-space; each is needed to obtain densities of states and spectral densities respectively. In the next section, we describe how these are used to calculate physical quantities of interest and discuss their significance.

## III Important quantities; sum rules

In this section we derive results for important physical quantities such as the densities of states (partial and total), spectral densities (partial and total) and inelastic scattering cross sections (coherent and incoherent). All of these are derivable from the real-space and the $`\stackrel{}{q}`$-space configuration-averaged Green’s function and enable us to make direct comparisons with experimental measurements.

### III.1 Densities of states

The total density of states for a 3-dimensional system is defined as,

$$\nu \left(\omega \right)=\frac{1}{3\pi N}m\{Tr\mathrm{𝐦𝐆}(\omega ^2)\},$$
(52)
where $`𝐦`$ is the mass matrix, and $`N`$ is the number of sites. In augmented space we have,

$`mG_{ii}`$ $`=`$ $`if|\widehat{𝐦}\widehat{𝐆}|if,`$ (53)
$`=`$ $`if|\widehat{𝐦}|ifif|\widehat{𝐆}|if`$
$`+if|\widehat{𝐦}|if_iif_i|\widehat{𝐆}|if,`$
$`=`$ $`\overline{m}\overline{G}_{00}+m^{}if_i|\widehat{𝐆}|if.`$
To evaluate the second term, we use the notation of Eq.(12) for the operators $`\widehat{𝐆}`$ and $`\widehat{𝐊}=\widehat{𝐆}^1`$. Then, using Eq.(18), we obtain

$$if_i|\widehat{𝐆}|if=𝐆^{}=\stackrel{~}{𝐊}^1𝐊^{}\overline{𝐆}=\mathrm{𝐅𝐊}^{}\overline{𝐆}.$$
(54)
We can, therefore, write

$`if_i|\widehat{𝐆}|if={\displaystyle \underset{l}{}}{\displaystyle \underset{j,n}{}}F_{ij}^{(i)(l)}K_{jn}^{(l)}G_{ni}.`$
Fourier transforming over the fluctuation site according to (29) gives

$`if_i|\widehat{G}|if`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{l,j,n}{}}{\displaystyle \underset{\stackrel{}{q}}{}}F\left(\stackrel{}{q}\right)_{0,jl}e^{i\stackrel{}{q}\stackrel{}{R}_{il}}`$
$`\text{ }\times K_{jl,nl}^{}G_{ni}.`$
The Fourier transform of $`𝐆`$ on the real-site index now gives,

$`if_i|\widehat{G}|if`$ $`=`$ $`{\displaystyle \frac{1}{N^2}}{\displaystyle \underset{l,j,n}{}}{\displaystyle \underset{\stackrel{}{q}\stackrel{}{q^{}}}{}}F\left(\stackrel{}{q}\right)_{0,jl}e^{i\stackrel{}{q}\stackrel{}{R}_{il}}`$
$`\text{ }\times K_{jl,nl}^{}G(\stackrel{}{q^{}})e^{i\stackrel{}{q}\stackrel{}{R}_{ni}}.`$
Finally we obtain,

$`if_i|\widehat{𝐆}|if`$ $`=`$ $`{\displaystyle \underset{mp}{}}{\displaystyle \underset{\stackrel{}{q}}{}}F\left(\stackrel{}{q}\right)_{0,m}K_{mp}^{}e^{i\stackrel{}{q}\stackrel{}{R}_p}G(\stackrel{}{q}),`$
where, $`m=jl`$, $`p=nl`$, the neighboring sites perturbed by the fluctuation. All the terms on the right hand side have been calculated already in the process of achieving self-consistency. The evaluation of the average density of states is thus straightforward.

The partial density of states for atoms of type $`s`$ is given by

$$\nu \left(\omega \right)_s=\frac{m_s}{3\pi N}m\{TrG(\omega ^2)^{ss}_{ii}\},$$
(55)
where

$$G^{ss}_{ii}=G^{ss}_{00}=0f|\widehat{\eta }_0^s\widehat{𝐆}|0f$$
(56)
because of translation invariance. We thus have

$`G^s_0=G^{ss}_{00}`$ $`=`$ $`0f|\widehat{\eta }_0^s|0f0f|\widehat{𝐆}|0f`$ (57)
$`+`$ $`if|\widehat{\eta }_0^s|0f_00f_0|\widehat{𝐆}|0f,`$
and, from Eq. (11), it follows that

$`\nu \left(\omega \right)_A`$ $`=`$ $`{\displaystyle \frac{m_A}{3\pi }}m[c_A\{0f|\widehat{𝐆}|0f\}`$
$`+\sqrt{c_Ac_B}\{0f_0|\widehat{𝐆}|0f\}],`$
$`\nu \left(\omega \right)_B`$ $`=`$ $`{\displaystyle \frac{m_B}{3\pi }}m[c_B\{0f|\widehat{𝐆}|0f\}`$
$`\sqrt{c_Ac_B}\{0f_0|\widehat{𝐆}|0f\}].`$
The elements of $`\widehat{𝐆}`$ in Eq.(53) were already evaluated while calculating the average density of states above.

The partial Green’s functions $`G^s_0`$ are used in calculating the incoherent scattering structure factor which is directly measured in the experiments,

$$S_{incoh}(\stackrel{}{Q},\omega )=\underset{s}{}b_s^2\stackrel{}{Q}mG^s\left(\omega \right)_0\stackrel{}{Q},$$
(59)
where $`b_s`$ is the incoherent scattering length for atoms of type $`s`$, and $`Q`$ is the phonon wavenumber.

### III.2 Spectral densities

The average spectral function is defined as,

$$𝒜_\lambda (\stackrel{}{q},\omega ^2)=\frac{1}{\pi }mG_\lambda (\stackrel{}{q},\omega ^2),$$
(60)
where $`\lambda `$ is a normal-mode branch index. More interesting quantities to calculate are the conditional or partial Green’s functions $`G^{ss^{}}(\stackrel{}{q},\omega ^2)`$ in $`\stackrel{}{q}`$-space because these enable one to calculate the coherent-scattering structure factors which are measured directly in the neutron-scattering experiments and are given by,

$$S_\lambda (\stackrel{}{q},\omega )_{coh}=\underset{ss^{}}{}d_sd_s^{}\frac{1}{\pi }mG_\lambda ^{ss^{}}(\stackrel{}{q},\omega ^2),$$
(61)
where, $`d_s`$ is the coherent scattering length for the species $`s`$.

The conditional Green’s functions are defined as

$`G^{ss^{}}(\stackrel{}{q},\omega ^2)={\displaystyle \frac{1}{N}}{\displaystyle \underset{ij}{}}G^{ss^{}}\left(\omega ^2\right)_{ij}e^{\stackrel{}{q}\stackrel{}{R}_{ij}},`$
$`G^{ss^{}}\left(\omega ^2\right)_{ij}=\widehat{\eta }_i^s\widehat{𝐆}\left(\omega ^2\right)\widehat{\eta }_j^s^{}=if|\widehat{\eta }_i^s\widehat{𝐆}\widehat{\eta }_j^s^{}|jf`$
$`=if|\widehat{\eta }_i^s|ifif|\widehat{𝐆}|jfjf|\widehat{\eta }_j^s^{}|jf`$
$`+if|\widehat{\eta }_i^s|if_iif_i|\widehat{𝐆}|jfjf|\widehat{\eta }_j^s^{}|jf`$
$`+if|\widehat{\eta }_i^s|ifif|\widehat{𝐆}|jf_jjf_j|\widehat{\eta }_j^s^{}|jf`$
$`+if|\widehat{\eta }_i^s|if_iif_i|\widehat{𝐆}|jf_jjf_j|\widehat{\eta }_j^s^{}|jf.`$ (62)
In Eq.(57) the index $`\lambda `$ is to be understood. These four terms include all the possible scattering processes when two different sites are occupied by two species. The four different terms involve calculations of the Green’s function under various circumstances of coupling between the average and the fluctuation states weighted by the appropriate concentrations .

We obtain from (57)

$`G^{ss^{}}_{ij}`$ $`=`$ $`c_sc_s^{}if|\widehat{𝐆}|jf+[c_s^{}\sqrt{c_s(1c_s)}`$ (63)
$`\times `$ $`(1)^{(1n^s)}if_i|\widehat{𝐆}|jf]+[c_s\sqrt{c_s^{}(1c_s^{})}`$
$`\times `$ $`(1)^{(1n^s^{})}if|\widehat{𝐆}|jf_j]+c_sc_s^{}`$
$`\text{ }\times \left(1\right)^{(n^s+n^s^{})}if_i|\widehat{𝐆}|jf_j.`$
The integer $`n^s`$ is equal to 1 if $`s`$=A and is equal to 0 if $`s`$=B.

These terms can be easily calculated using Fourier transforms as has been previously demonstrated for the density of states. The final forms of the conditional Green’s functions in $`\stackrel{}{q}`$-space are

$`G^{AA}(\stackrel{}{q},\omega ^2)`$ $`=`$ $`c_A^2G(\stackrel{}{q},\omega ^2)+c_A\sqrt{c_Ac_B}(T_1+`$
$`T_2)+c_Ac_BT_3,`$
$`G^{BB}(\stackrel{}{q},\omega ^2)`$ $`=`$ $`c_B^2G(\stackrel{}{q},\omega ^2)c_B\sqrt{c_Ac_B}(T_1+`$
$`T_2)+c_Ac_BT_3,`$
$`G^{AB}(\stackrel{}{q},\omega ^2)`$ $`=`$ $`c_Ac_BG(\stackrel{}{q},\omega ^2)+\sqrt{c_Ac_B}(c_BT_1`$
$`c_AT_2)c_Ac_BT_3,`$
$`G^{BA}(\stackrel{}{q},\omega ^2)`$ $`=`$ $`c_Ac_BG(\stackrel{}{q},\omega ^2)\sqrt{c_Ac_B}(c_AT_1`$ (64)
$`c_BT_2)c_Ac_BT_3,`$
where

$`T_1`$ $`=`$ $`{\displaystyle \underset{nm}{}}F_{0n}\left(\stackrel{}{q}\right)K_{nm}^{}e^{i\stackrel{}{q}\stackrel{}{R}_m}G(\stackrel{}{q},\omega ^2),`$
$`T_2`$ $`=`$ $`{\displaystyle \underset{nm}{}}G(\stackrel{}{q},\omega ^2)e^{i\stackrel{}{q}\stackrel{}{R}_n}K_{nm}^{}F_{m0}\left(\stackrel{}{q}\right),`$
$`T_3`$ $`=`$ $`F_{00}\left(\stackrel{}{q}\right)+{\displaystyle \underset{nm}{}}{\displaystyle \underset{lp}{}}F_{0n}\left(\stackrel{}{q}\right)K_{nm}^{}e^{i\stackrel{}{q}\stackrel{}{R}_m}`$ (65)
$`\text{ }\times G(\stackrel{}{q},\omega ^2)e^{i\stackrel{}{q}\stackrel{}{R}_l}K_{lp}^{}F_{p0}\left(\stackrel{}{q}\right),`$
and where $`n,m,l`$ and $`p`$ are the neighboring sites of the fluctuation site $`0`$ influenced by the perturbation. For the lattices with each site having inversion symmetry, $`T_1`$=$`T_2`$ holds because $`T_1`$ and $`T_2`$ are the contributions from two processes which are conjugate to one another. In that case, $`G^{AB}(\stackrel{}{q},\omega ^2)`$= $`G^{BA}(\stackrel{}{q},\omega ^2)`$ when $`c_A`$=$`c_B`$.

The sum rules for the conditional Green’s functions are derived the following way: Integrate $`\omega \widehat{𝐆}\left(\omega ^2\right)`$, Eq.(2), along the real axis, closing the contour above at infinity, obtaining

$`{\displaystyle 𝑑\omega \omega \widehat{𝐆}\left(\omega ^2\right)}`$ $`=`$ $`\widehat{𝐦}^1\pi i,`$
or

$`{\displaystyle _0^{\mathrm{}}}𝑑\omega 2\omega m\widehat{𝐆}\left(\omega ^2\right)`$ $`=`$ $`\widehat{𝐦}^1\pi .`$ (66)
Similarly, using Eq.(57), taking the Fourier transform of Eq.(58), carrying out the contour integral and inserting (61) yields the sum rule for the partial spectral functions

$$_0^{\mathrm{}}d\omega 2\omega mG^{ss^{}}(\stackrel{}{q},\omega ^2)=\pi \frac{c_s}{m_s}\delta _{ss^{}}.$$
(67)
For the total Green’s function we obtain

$$_0^{\mathrm{}}d\omega 2\omega mG(\stackrel{}{q},\omega ^2)=\pi (\frac{c_A}{m_A}+\frac{c_B}{m_B}).$$
(68)
The experimental dispersion curves are obtained from the wave-vector dependence of the peak frequencies of the structure factors as measured, after a deconvolution of the experimental resolution function. The question is whether the dispersion curves so obtained, which incorporate the effect of the coherent scattering lengths, differ significantly from those obtained from the peak frequencies of the Green’s function itself which gives, in principle, a proper description of the dynamics but does not contain the scattering length weighting. To answer that question one needs to recognize that the peak positions in $`m𝐆`$ are very closely related to the zeroes of $`e𝐆^1`$ at a given wave vector. If we diagonalize the Hermitian $`e𝐆^1`$ both with respect to mode and species index, each of the two species components of $`e𝐆^1`$ will have a zero. Correspondingly, each of the two components of $`m𝐆`$ will have a peak, if $`m𝚺`$ does not wipe it out. So, in the species representation, the different matrix elements $`mG^{ss^{}}`$ will thus all have these peaks at nearly the same frequencies. Thus, the weighting of the $`G^{ss^{}}`$ by the scattering lengths will not shift the peak positions significantly even when the scattering lengths differ appreciably. The intensities and the lineshapes of $`S_{coh}`$ and $`m𝐆`$ may differ significantly, but the peak positions will generally differ little. In summary, the structure factor fairly accurately reflects the phonon dynamics contained in $`𝐆`$ with regard to the dispersion curves, an important fact illustrated below in the next two sections where we present our calculations on Ni55Pd45 and Ni50Pt50 alloys.

## IV Application to Ni55Pd45; weak force-constant disorder

In the next two sections we explore the relative importance of force-constant disorder and mass disorder for the vibrational properties of random alloys in two specific alloys Ni55Pd45 and Ni50Pt50. In the former alloy, the mass disorder is much larger than the force-constant disorder. The mass ratio $`m_{Pd}/m_{Ni}`$ is 1.812, whereas the Pd force-constants are only about 15$`\%`$ larger than those of Ni dbm ; bvk . In the latter alloy, both the mass disorder and the force-constant disorder are large, providing an interesting contrast between the two materials. For both the cases we have done our calculations on 200 $`\omega `$-points and have used a small imaginary frequency part of -0.01 in the Green’s function. For the Brillouin zone integration 356 $`\stackrel{}{q}`$-points in the irreducible 1/48-th of the zone produced well converged results. The simplest linear-mixing scheme was used to accelerate the convergence. For both the cases the number of iterations ranged from 3 to 13 depending on the frequency $`\omega `$.

For Ni55Pd45, we compare the results of virtual crystal(VCA), CPA, and ICPA computations, using the ICPA force-constants to construct the averages used in the VCA and CPA and compare the results with experiment. We make a distinction between that use of the VCA and of ‘mean crystal’ models in which the average mass is employed and a set of ‘mean-crystal’ force-constants are fitted as parameters to the experimental data.

Kamitakahara and Brockhouse expt2  investigated Ni55Pd45 by inelastic neutron scattering and reported a strange observation. A theoretical calulation based on a mean crystal model having the average mass and fitted force-constants between those of Ni and Pd agreed closely with the experimental dispersion curves. This was quite a puzzle because it suggested that the large mass disorder had little effect. There were theoretical studies on this system using recursionmook1  and the average t-matrix approximation kt , but no theoretical results for the frequencies were available. In an attempt to solve this puzzle, we have carried out calculations with the CPA, the ICPA, and the VCA, as well as with the mean-crystal model used in Ref. 3 .

For the ICPA calculation, we assumed that the explicit scattering caused by the force-constant disorder was confined to the nearest-neighbors. This assumption is justified because the nearest-neighbor force-constants are an order of magnitude larger than those of the further neighbors so that the nearest-neighbors feel the effect of disorder most strongly. For the virtual-crystal or the average medium into which the scattering was embedded, we kept terms in the Hamiltonian up to the fourth neighbor, which turned out to be sufficient. The problem with force-constant disorder scattering calculations is the general absence of prior information about species-dependent force-constants. We note that Pd is the larger atom here. In the alloy, the Ni-Pd separation is larger than the Ni-Ni separation. As a result, the Ni-Pd force-constants should be less than the Ni-Ni ones. Using this intuitive argument and for simplicity in this illustration, we kept the $`\varphi _{NiNi}^{\alpha \beta }`$ and $`\varphi _{PdPd}^{\alpha \beta }`$ the same as those of the pure materials dbm  and reduced the $`\varphi _{NiPd}^{\alpha \beta }`$ below the $`\varphi _{NiNi}^{\alpha \beta }`$ by an $`\alpha \beta `$ independent factor. The dispersion curves were obtained from the ICPA calculations using $`\varphi _{NiPd}^{\alpha \beta }`$=0.7$`\varphi _{NiNi}^{\alpha \beta }`$ (solid lines) for the nearest neighbors and using the force-constants of Ref.3 for the higher neighbors. They are compared in the top panel of Fig.3(a) with the experimental results expt2  (open circles) and the VCA for the same force-constants(dashed lines). These ICPA dispersion curves were constructed by numerically determining the peaks in the coherent scattering structure factor $`S_\lambda (\stackrel{}{q},\omega )_{coh}`$ given by Eq.(56), which was calculated using the partial spectral functions of Eqs.(59) and (60) and weighting them with the coherent scattering lengths for Ni and Pd. We could thus make a direct comparison with the experimental results because the neutron data observed in the experiments inherently incorporates the effect of the scattering lengths of the species.

Excellent agreement of the ICPA with experiment was obtained for all three symmetry directions and for each branch by varying only one parameter in the force-constant matrix. This suggests that the force-constant disorder is weak and the system is dominated by the mass-disorder, as one would expect from the numerical values of the parameters. It is confirmed by the results of the CPA calculations shown in the bottom panel of the Fig. 3(a), using the same force-constants. As in the top panel, the solid lines are the CPA results, the circles are the experimental points, and the dashed lines are the VCA results. The agreement with the experiment again suggests the dominance of the mass-disorder, but there are more interesting points to note. In the long-wavelength (low $`\stackrel{}{q}`$) regime, the VCA, the CPA, and the ICPA curves are indistinguishable because the self-averaging of both mass and force-constants over a single wavelength reduces both the CPA and the ICPA to the VCA. But, as we move to high wave-vectors , the VCA deviates to frequencies below the experimentally observed ones. This fact is due to the use of an average mass in the Hamiltonian. In the high-wave-vector region, the lighter atoms, i.e. Ni in this case, dominate and push the frequencies up. That is why the CPA and the ICPA agree very well across the Brillouin zone while the VCA fails for the high wave-vectors. The reason that Kamitakahara and Brockhouse got a very good fit to the experimental points in Ref. 3 by using their mean-crystal model is that they obtained parametrized force-constants which were higher than those calculated in the VCA. Though they had used the average mass in their calculations the higher values of the force-constants (They used $`\overline{𝚽}=c_A𝚽_{AA}+c_B𝚽_{BB}`$ rather than $`\overline{𝚽}=c_{AA}^2𝚽_{AA}+c_{BB}^2𝚽_{BB}+2c_Ac_B𝚽_{AB}`$) compensated for their omission of the effect of the mass fluctuations. This is illustrated in Fig.3(b). There, the dashed lines are their mean-crystal model calculations, the circles are the experimental points, and the solid lines represent a CPA calculation with the force-constants used in Ref. 3. Here we see that the CPA yields frequencies that are too high in the large wave-vector region. The CPA captures the effect of mass-fluctuation and the domination of the Ni atoms for higher wave-vectors, but the higher values of the assumed mean force-constants pushes the frequencies up further, thereby worsening the agreement with the experiment. But, there was yet another puzzle to solve. Another striking feature is that in spite of incorporating the scattering lengths in our calculations in determining the frequencies there was little change in the ICPA results with respect to the CPA results even though the coherent scattering lengths of Ni and Pd differ significantly( The coherent scattering length for Ni is 1.03 while that of Pd is 0.6). This lack of change can be understood from a comparison between the partial and total spectral functions (Fig. 4(a)) and the partial and the total coherent structure factors (Fig 4(b)). In these figures, we have shown examples of ICPA spectral functions and structure factors along the \[$`\zeta `$,0,0\] direction, $`\zeta =\frac{|\stackrel{}{q}|}{|\stackrel{}{q}_{max}|}`$, for a low $`\zeta `$ , one in the middle of the Brillouin zone and one at the edge of the Brilluoin zone. In Fig. 5(a), the solid lines are the total spectral function while the dotted lines are the Ni-Ni spectra, the long-dashed lines are the Pd-Pd spectra and the dot-dashed lines are the Ni-Pd contribution. In each case, the peaks corresponding to the dominating species and that in the Ni-Pd curves occur at the same general positions. For example, in the \[.3,0,0\]-L curves, the peak in the spectral function is mostly that of Pd atoms while for the \[.5,0,0\]-L and -L curves, the contributions are from Ni atoms, the Pd-Pd contribution here is much less and that too is almost completely neutralised by the Ni-Pd contribution in the low frequency region. The occurence of the peaks of the Ni-Pd spectral functions and that of the Ni-Ni or Pd-Pd spectral functions almost at the same position across the Brillouin zone suggests that the inclusion of scattering lengths would primarily alter the relative weights of various contributions and thereby altering only the line shapes, while the dispersion curves wouldn’t undergo significant change. This is demonstrated in Fig. 4(b). There, the solid lines are the coherent structure factors $`S_\lambda (\stackrel{}{q},\omega )`$, the dotted lines are the Ni-Ni spectral functions weighted using the scattering length of Ni, the long-dashed lines are the Pd-Pd spectral functions weighted using the scattering length of Pd and the dot-dashed lines are the Ni-Pd spectral functions weighted by the scattering lengths of Ni and Pd according to Eq.(56). One can see that the weighting affects primarily the peak heights. These explicit numerical results confirm the qualitative argument given at the end of section III.

The disorder-induced widths are important because the effect of disorder is often manifested in them more directly than in the frequencies. Kamitakahara and Brockhouse extracted full widths at half maxima(FWHM) from their neutron groups by assuming that the observed line shape could be adequately approximated by the convolution of a Gaussian resolution function(representing the experimental resolution) with a Lorentzian natural line shape. Thus, for a comparison of our results with theirs, we have fitted our structure factors to a Lorentzian to extract the widths. The results are shown in Fig.5. The solid lines are the widths (FWHM) obtained in the ICPA, and the dotted lines are those obtained in the CPA. The circles are the experimental points, the filled circles having better resolution expt2  were those in which Kamitakahara and Brockhouse had more confidence. Generally, there is little difference between the widths obtained in the CPA and the ones obtained in the ICPA. In all three symmetry directions and for all branches, the ICPA performs slightly better than the CPA for high wave-vectors. The worst agreement with the experiment is for high wave-vectors in the $`[\zeta ,0,0]`$ and $`[\zeta ,\zeta ,\zeta ]`$ longitudinal branches and the $`[\zeta ,\zeta ,\zeta ]`$ transverse branch. In these cases, the low values of the widths in the theoretical calculations can be understood from the shape of the structure factors. From the examples in Fig.4 one can see that the agreement with experiment is good when we have a symmetric line shape, for example, for the \[.5,0,0\]L mode. On the other hand, the worst agreements with the experimental widths are for cases where we obtain a highly asymmetric line shape, for example, for the L mode. Fitting Lorentzians to such asymmetric line shapes is not condusive to meaningful values of the FWHMs. Because they obtain higher widths than the theories in those particular cases where the observations have been made with worse resolution(open circles in Fig.5), it is also not clear how trustworthy their treatment of the resolution function is.

The discussion above clearly tells us that for Ni55Pd45, the dominant effect is mass disorder. That alloy therefore does not provide a proper test of the ICPA. Nevertheless, our discussions show how a mean-crystal model can compensate for the neglect of mass fluctuations in alloys with little force-constant disorder through the introduction of erroneous mean force-constants, a classic case of cancellation of errors.

## V Ni50Pt50; strong mass and force-constant disorder

The mass ratio $`m_{Pt}/m_{Ni}`$ is 3, quite large compared to that in the NiPd system. The force-constants of Pt are, on an average, 55$`\%`$ larger dbm  than those of Ni. This makes it a potential example of strong force-constant disorder. Tsunoda et al.expt3  investigated NixPt1-x by inelastic neutron scattering and compared their observations with the CPA. Here, for illustration, we have considered $`x`$=0.5 only because that makes it a concentrated alloy and the failure of CPA was, qualitatively, very prominent at this concentration. They compared their incoherent scattering data with that of the CPA which predicted a split band separating out Ni and Pt contributions with a gap between them. The experiments did not reveal a split-band, and it was very clear that the inter-species forces play a significant role. We performed calculations with the CPA and the ICPA. As before, we used the ICPA force-constants in the CPA. The choice of ICPA force-constants was more difficult than for the NiPd because of the larger size difference between Ni and Pt. In this alloy, the Ni-Pt separation is also larger than the Ni-Ni separation. As a result, the Ni-Pt force-constants should also be less than those of Ni-Ni. Moreover, a pair of Ni atoms would find themselves in a cage partly made of larger Pt atoms which would therefore reduce the Ni-Ni force-constants relative to their values in the pure material. Similarly, the bigger Pt atoms find themselves compressed between much smaller Ni atoms, which would increase the Pt-Pt force-constants with respect to their values in pure Pt. Using this intuitive argument, we found that the following guesses for the force-constants worked well: $`\varphi _{NiNi}^{xy}`$, $`\varphi _{PtPt}^{xy}`$, $`\varphi _{NiNi}^{zz}`$ and $`\varphi _{PtPt}^{zz}`$ are kept the same as those of the pure materials dbm  and

$`\varphi _{NiNi}^{yy}=\varphi _{NiNi}^{xx}=0.9\varphi _{Ni(pure)}^{xx},`$
$`\varphi _{PtPt}^{yy}=\varphi _{PtPt}^{xx}=1.1\varphi _{Pt(pure)}^{xx},`$
$`\varphi _{NiPt}^{\alpha \beta }=0.8\varphi _{NiNi}^{\alpha \beta },\text{ for all }\alpha \text{,}\beta \text{.}`$
In Fig.6, we compare the ICPA results for the incoherent neutron structure factor (Eq.(54)) with those of the CPA and the experiment expt3 . In Fig. 6(a), the solid line stands for the ICPA results while the dotted line stands for the CPA. In Fig.6(b), the solid line is the ICPA results and the dotted line is the experimental curve. The CPA results suggest a split-band behaviour in the middle of the band clearly separating the Pt contribution in the low frequency region from the Ni contribution in the high frequency region. The overall contribution from the low frequency region is much less than that of the high frequency region in this system because the low frequency region is dominated by the heavier atom Pt which has a much lower incoherent scattering length expt3  than Ni, 0.1 in comparison to 4.5 for Ni. Including only mass fluctuations and ignoring the Ni-Pt correlated motion induced by the environmental disorder gives rise to this spurious gap in the CPA results. This point is further discussed below in connection with our analysis of the coherent scattering results. On the other hand, by incorporating the force-constant disorder, as is done in the ICPA, we get rid of this spurious gap and obtain good agreement with the experimental results, including the position of the right band-edge. The influence of the force-constant disorder is demonstrated more prominently in the dispersion curves and the line shapes.

In Fig.7, we compare the dispersion curves and widths obtained in the ICPA from the coherenet scattering structure factors, using the force constants as given above, with those in the CPA, using the averages of the same force-constants, and with the experimental results expt3 . The top three panels are the results obtained in the ICPA for the three symmetry directions. The procedure has already been discussed in the previous section. The bottom panels are the CPA results. The circles and error bars in the left panels are the experimental frequencies and the widths(FWHMs). The filled circles are experimental results with more accuracy and greater resolution. The shaded regions in the top panels and in the leftmost panel in the bottom span the calculated FWHM’s. The FWHM’s for the middle and the right panels in the bottom are indicated by the thin dotted lines. For the $`[\zeta ,0,0]`$ and the $`[\zeta ,\zeta ,\zeta ]`$ directions, the solid lines represent the longitudinal modes and the dotted lines the transverse modes. For the $`[\zeta ,\zeta ,0]`$ direction, the solid lines show the longitudinal mode while the long-dashed and the dot-dashed curves stand for the T1 and the T2 transverse modes, respectively. The experimental results are available only for the $`[\zeta ,0,0]`$ directions in this system. The ICPA agrees much better with the experiments than the CPA for both the longitudinal and the transverse branches. The CPA frequencies are generally below the experimental ones at low frequencies and above the experimental ones at high frequencies. The discrepancy gets worse as we move from the middle of the zone towards the zone edge. This can be explained the following way: the high wave-vector region is dominated by the lighter atoms. The use of the average force-constants by the CPA coupled with the domination of the lighter mass pushes the frequencies further up, thus producing a significant deviation from the experimental ones. The severity of this effect can be understood from the widths as well. In the CPA, the experimental points stay well outside the disorder-induced widths centered at the peak frequencies. The discrepancy is substantially reduced by the inclusion of force-constant disorder, as is seen from the ICPA results. Its inclusion changes the dispersion curves qualitatively as well. In the CPA, the bands extend fully across the Brilluoin zone for all symmetry directions while in the ICPA, the Pt-dominated peaks disappear at high-$`\zeta `$ and the Ni-dominated peaks wash out at low-$`\zeta `$ for all modes and symmetry directions, an effect observed in the experiments. This is a clear consequence of the force-constant disorder which can be understood by inspecting the spectral line shapes. In Fig.8(a) we present the partial and the total spectral densities for three different wave-vectors, one on the lower side of the zone, one in the middle and one at the boundary. The solid lines are the total spectral functions, the dotted lines are the Ni-Ni partial spectral functions, the long-dashed lines are the Pt-Pt contributions while the dot-dashed lines are the Ni-Pt contributions. In Fig.8(b) we present the partial and the total coherent scattering structure factors i.e. the spectral functions weighted by the coherent scattering lengths of the species. The coherent scattering lengths of Ni and Pt differ by only 7$`\%`$(the scattering length of Ni is 1.03 while that of Pt is 0.95), much less than do those of Ni and Pd. However, even this small difference produces significant changes in the line shapes and in the peak frequencies. In Fig.8(b), the solid lines give the total structure factor, the dotted lines are for weighted Ni-Ni contributions, the long-dashed lines give the weighted Pt-Pt contributions and the dot-dashed lines stand for weighted Ni-Pt contributions. A close inspection of the various contributions reveals the fact that unlike in NiPd, the Ni-Pt contribution plays the key role in determining the weight in the middle of the band (and in obtaining the merged bands in Fig.6) as well as adding or subtracting weights to the Ni-Ni or Pt-Pt contributions, thus elevating or suppressing one of the peaks. For example, in Fig.8(a), in the \[.5,0,0\]T curves, the Ni-Pt contribution adds weight to the total spectral function on top of the Pt-Pt peak at the low frequenies while it subtracts weight from the Ni-Ni contribution at higher frequencies thereby causing a weakly defined peak at high frequencies. In the \[.5,0,0\]L and in the T curves, the Ni-Pt contribution adds weight between the Pt-Pt and Ni-Ni peaks, thereby removing the gap in the spectrum. The Ni-Pt contribution is totally due to inclusion of force-constant disorder, since diagonal disorder produces no such contribution. Thus, the CPA produces spectral functions having two well-defined peaks corresponding to the Ni-Ni and the Pt-Pt contributions with a gap in between resulting in extended dispersion curves and split bands. The effect of incorporation of the difference in scattering lengths can be seen from these two figures as well. For example, in the T curves, there are two well-defined peaks in the total spectral functions, whereas the low-frequency peak is transformed into a shoulder in the total structure factor. This is because Ni has the larger scattering length which enhances the weight associated with the Ni-Pt contribution thereby cancelling more effectively the contribution from the Pt-Pt part. Similar effects are seen in the \[.5,0,0\]L and L curves. Moreover, this weighting sometimes produces a weakly defined peak whose FWHM cannot be well determined, which explains the observed washing out of the dispersion curves noted above. The effect of the small difference in scattering lengths is amplified by the force-constant disorder through the Ni-Pt structure factor.

In sum, the force-constant disorder plays a significant role in Ni50Pt50, and a theory with mass disorder only, like the CPA, fails both qualitatively and quantitatively in such cases. On the other hand, the ICPA successfully explains the effects of force-constant disorder through its effect on the partial structure factors. It also demonstrates the relative importance of the contributions of various atomic species to the coherent and incoherent structure factors which the CPA cannot. The ICPA and the Ni50Pt50 system therefore provide a proper test case for force-constant disorder and show that the ICPA can form a basis for understanding the lattice dynamics of other binary alloys.

## VI Conclusions

We have presented a straight-forward and tractable formulation of the KLGD klgd  method for single-site scattering of phonons in three dimensional lattices. We have demonstrated how this multiple-scattering based formalism captures the effects of off-diagonal and environmental disorder. The use of augmented-space to keep track of the configurations of the system has made the formalism simple yet powerful. The resulting translational invariance makes it numerically tractable as well. A significant contribution beyond Ref.18 is the derivation of the partial Green’s functions in real-space and the derivation of the partial spectral functions as well as their sum rules. This enables one to make direct comparison with neutron scattering experiments because of the incorporation of the scattering lengths of the different species. We have applied the formalism to real random alloys for the first time. In Ni55Pd45 we have demonstrated that mass disorder plays the prominent role, and the CPA consequently does a rather good job whereas the mean-crystal model requires erroneous fitted force-constants. Our partial structure factors enable us to understand the insensitivity of the normal modes towards the difference in the coherent scattering lengths of the two species despite the significant difference of 43$`\%`$ in this system. The Ni50Pt50 results demonstrate the prominence of force-constant disorder even in a case where the mass ratio is 3. The ICPA agrees well with both the coherent and the incoherent scattering experiments, whereas the CPA fails, both quantitatively and qualitatively. We are able to establish the role of the force-constant differences between the species in great detail with the help of the partial spectral functions and the partial structure factors. We have consequently clearly demonstrated that for systems like NiPt, where the force-constants are strongly species dependent, the determination of their values is crucial. However, we had no prior information about the species dependence of the force-constants. Intuitive arguments led to a set of force-constants which turned out to be quite good. A better understanding of the role of disorder in the lattice dynamics of random alloys could be achieved with prior information about the force-constants. These could be obtained, e.g., from first principles calculations on a set of ordered alloys.

###### Acknowledgements.

We thank Professors David Vanderbilt, Gabriel Kotliar and Karin Rabe for useful discussions and for the use of their computer systems. We thank W. Kamitakahara and R. Nicklow for useful communications regarding their experimental results.

## Appendix A Matrix elements

For the calculations in the nearest-neighbor approximation, one needs to evaluate the $`\{3\left(Z+1\right)\}^2`$ matrix elements of the operators $`𝐊^{}`$, $`𝐊^{}`$, $`\stackrel{~}{𝐕}\left(\stackrel{}{q}\right)`$ and $`𝐆_{vca}^1`$ and use them as inputs. These evaluations are done in augmented space using Eq.(14). The symmetry of the lattice structure is used to reduce the number of matrix elements evaluated. Here, we give results only for an fcc lattice. All the matrices are, therefore, of dimension $`39\times 39`$.

In an fcc system, each atom has 12 nearest neighbors with coordinates $`(\pm \frac{1}{2},\pm \frac{1}{2},0)`$, $`(\pm \frac{1}{2},0,\pm \frac{1}{2})`$, and $`(0,\pm \frac{1}{2},\pm \frac{1}{2})`$ with respect to the coordinates of the reference atom at $`(0,0,0)`$. The force-constants, between the atom $`0`$ and its neighbors satisfy the following cubic symmetry relation.

$$\varphi _{0j}^{\alpha \beta }=\varphi _{0j}^{\beta \alpha }=\varphi _{j0}^{\alpha \beta }=\varphi _{0k}^{\alpha \beta },$$
(69)
where $`\stackrel{}{R}_{0j}=\stackrel{}{R}_{0k}`$ and $`k`$ and $`j`$ are two neighbors on opposite sides of site $`0`$. For example, the force-constant matrix between the atoms $`(0,0,0)`$ and $`(\frac{1}{2},\frac{1}{2},0)`$ is of the form,

$$\varphi _{(000,\frac{1}{2}\frac{1}{2}0)}=\left(\begin{array}{ccc}a& b& 0\\ b& a& 0\\ 0& 0& g\end{array}\right).$$
(70)
The force-constant matrices between the atom $`0`$ and its other neighbors can easily be calculated from (A2) via the cubic symmetry operations. The results are:

$`\left(G_{vca}^1\right)_{ij}^{\alpha \beta }`$ $`=`$ $`\overline{m}\omega ^28D_1^{xx}4D_1^{xx}\text{ }\text{if }i=j\text{}\alpha =\beta ,`$ (71)
$`=`$ $`0\text{ }\text{if }i=j\text{}\alpha \beta ,`$
$`=`$ $`D_1^{xx}\text{ }\text{ if }i=0\text{}j=n=112\text{}\alpha =\beta \text{}R_j^\alpha =0,`$
$`=`$ $`D_1^{xx}\text{ }\text{if }i=0\text{}j=n=112\text{}\alpha =\beta \text{}R_j^\alpha 0\text{ },`$
$`=`$ $`2R_j^\alpha \times 2R_j^\beta \times D_1^{xy}\text{ }\text{if }i=0\text{}j=n=112\text{}\alpha \beta ,`$
$`=`$ $`0\text{ }\text{otherwise},`$
and

$`\left(G_{vca}^1\right)_{ij}^{\alpha \beta }=\left(G_{vca}^1\right)_{ji}^{\alpha \beta }=\left(G_{vca}^1\right)_{ij}^{\beta \alpha },\text{ }\text{for all }i\text{,}j\text{,}\alpha \text{, and }\beta .`$
Also, we find

$`\left(K^{}\right)_{ij}^{\alpha \beta }`$ $`=`$ $`m^{}\omega ^28D_2^{xx}4D_2^{xx}\text{ }\text{if }i=j=0\text{}\alpha =\beta ,`$ (72)
$`=`$ $`0\text{ }\text{if }i=j=0\text{}\alpha \beta ,`$
$`=`$ $`D_2^{xx}\text{ }\text{ if }i=0\text{}j=n=112\text{}\alpha =\beta \text{}R_j^\alpha =0,`$
$`=`$ $`D_2^{xx}\text{ }\text{if }i=0\text{}j=n=112\text{}\alpha =\beta \text{}R_j^\alpha 0\text{ },`$
$`=`$ $`2R_j^\alpha \times 2R_j^\beta \times D_2^{xy}\text{ }\text{if }i=0\text{}j=n=112\text{}\alpha \beta ,`$
$`=`$ $`\left(K^{}\right)_{0j}^{\alpha \beta }\text{ }\text{if }i0\text{}j=i=n=112,`$
$`=`$ $`0\text{ }\text{otherwise},`$
and

$`\left(K^{}\right)_{ij}^{\alpha \beta }=\left(K^{}\right)_{ji}^{\alpha \beta }=\left(K^{}\right)_{ij}^{\beta \alpha },\text{ }\text{for all }i\text{,}j\text{,}\alpha \text{, and }\beta .`$
Similarly, we obtain

$`\stackrel{~}{V}\left(\stackrel{}{q}\right)_{ij}^{\alpha \beta }`$ $`=`$ $`L^\alpha \text{ }\text{if }i=j=0\text{}\alpha =\beta ,`$ (73)
$`=`$ $`4D_3^{xy}Sin(q_\alpha )Sin(q_\beta )\text{ }\text{if }i=j=0\text{}\alpha \beta ,`$
$`=`$ $`D_5\text{ }\text{ if }i=0\text{}j=n=112\text{}\alpha =\beta \text{}R_j^\alpha =0,`$
$`=`$ $`D_6\text{ }\text{if }i=0\text{}j=n=112\text{}\alpha =\beta \text{}R_j^\alpha 0\text{ },`$
$`=`$ $`2R_j^\alpha \times 2R_j^\beta \times D_7\text{ }\text{if }i=0\text{}j=n=112\text{}\alpha \beta ,`$
$`=`$ $`D_4^{xx}D_1^{xx}\text{ }\text{ if }i0\text{}j=i=n=112\text{}\alpha =\beta \text{}R_i^\alpha =0,`$
$`=`$ $`D_4^{xx}D_1^{xx}\text{ }\text{if }i0\text{}j=i=n=112\text{}\alpha =\beta \text{}R_i^\alpha 0\text{ },`$
$`=`$ $`2R_j^\alpha \times 2R_j^\beta \times \left(D_4^{xy}D_1^{xy}\right)\text{ }\text{if }i0\text{}j=i=n=112\text{}\alpha \beta ,`$
$`=`$ $`D_3^{xx}e^{i\stackrel{}{q}\stackrel{}{R}_i}\text{ }\text{if }ji=n=112\text{}\stackrel{}{R}_j=\stackrel{}{R}_i\text{}\alpha =\beta \text{}R_j^\alpha =0,`$
$`=`$ $`D_3^{xx}e^{i\stackrel{}{q}\stackrel{}{R}_i}\text{ }\text{if }ji=n=112\text{}\stackrel{}{R}_j=\stackrel{}{R}_i\text{}\alpha =\beta \text{}R_j^\alpha 0,`$
$`=`$ $`D_3^{xy}e^{i\stackrel{}{q}\stackrel{}{R}_i}\times 2R_i^\alpha \times 2R_i^\beta \text{ }\text{if }ji=n=112\text{}\stackrel{}{R}_j=\stackrel{}{R}_i\text{}\alpha \beta ,`$
$`=`$ $`0\text{ }\text{otherwise},`$
and

$`\stackrel{~}{V}\left(\stackrel{}{q}\right)_{ij}^{\alpha \beta }`$ $`=`$ $`\left(\stackrel{~}{V}\left(\stackrel{}{q}\right)^{}\right)_{ji}^{\alpha \beta }=\stackrel{~}{V}\left(\stackrel{}{q}\right)_{ij}^{\beta \alpha },\text{ }\text{for all }i\text{,}j\text{,}\alpha \text{, and }\beta .`$
In these evaluations, we have used the notation $`R_{0j}=R_j`$, $`n=112`$ to represent the 12 nearest neighbors, and the notations

$`\overline{m}`$ $`=`$ $`c_Am^A+c_Bm^B,`$
$`m^{}`$ $`=`$ $`\sqrt{c_Ac_B}\left(m^Am^B\right),`$
$`\stackrel{~}{m}`$ $`=`$ $`c_Bm^A+c_Am^B,`$
$`D_1^{xx}`$ $`=`$ $`c_A^2a_{AA}+c_B^2a_{BB}+2c_Ac_Ba_{AB},`$
$`D_1^{xx}`$ $`=`$ $`c_A^2g_{AA}+c_B^2g_{BB}+2c_Ac_Bg_{AB},`$
$`D_1^{xy}`$ $`=`$ $`c_A^2b_{AA}+c_B^2b_{BB}+2c_Ac_Bb_{AB},`$
$`D_2^{xx}`$ $`=`$ $`\sqrt{c_Ac_B}\{c_Aa_{AA}c_Ba_{BB}+(c_Bc_A)a_{AB}\},`$
$`D_2^{xx}`$ $`=`$ $`\sqrt{c_Ac_B}\{c_Ag_{AA}c_Bg_{BB}+(c_Bc_A)g_{AB}\},`$
$`D_2^{xy}`$ $`=`$ $`\sqrt{c_Ac_B}\{c_Ab_{AA}c_Bb_{BB}+(c_Bc_A)b_{AB}\},`$
$`D_3^{xx}`$ $`=`$ $`c_Ac_B\left(a_{AA}+a_{BB}2a_{AB}\right),`$
$`D_3^{xx}`$ $`=`$ $`c_Ac_B\left(g_{AA}+g_{BB}2g_{AB}\right),`$
$`D_3^{xy}`$ $`=`$ $`c_Ac_B\left(b_{AA}+b_{BB}2b_{AB}\right),`$
$`D_4^{xx}`$ $`=`$ $`c_Ac_B\left(a_{AA}+a_{BB}\right)+\left(c_A^2+c_B^2\right)a_{AB},`$
$`D_4^{xx}`$ $`=`$ $`c_Ac_B\left(g_{AA}+g_{BB}\right)+\left(c_A^2+c_B^2\right)g_{AB},`$
$`D_4^{xy}`$ $`=`$ $`c_Ac_B\left(b_{AA}+b_{BB}\right)+\left(c_A^2+c_B^2\right)b_{AB},`$
$`D_5`$ $`=`$ $`\left(D_1^{xx}D_4^{xx}\right)+D_3^{xx}e^{i\stackrel{}{q}\stackrel{}{R}_j},`$
$`D_6`$ $`=`$ $`\left(D_1^{xy}D_4^{xy}\right)+D_3^{xy}e^{i\stackrel{}{q}\stackrel{}{R}_j},`$
$`D_7`$ $`=`$ $`\left(D_1^{xx}D_4^{xx}\right)+D_3^{xx}e^{i\stackrel{}{q}\stackrel{}{R}_j},`$
$`L^\alpha `$ $`=`$ $`(\overline{m}\stackrel{~}{m})\omega ^28(D_1^{xx}D_4^{xx})4(D_1^{xx}D_4^{xx})4D_3^{xx}\{Cosq_\alpha (Cosq_\gamma +Cosq_\delta )`$ (74)
$`\text{ }4D_3^{xx}Cosq_\gamma Cosq_\delta ;\text{ }\gamma ,\delta \alpha .`$
The symmetries of the force-constant matrices are reflected in the operators as well. The effect of itineration is captured in $`\stackrel{~}{V}\left(\stackrel{}{q}\right)_{ij}^{\alpha \beta }`$ through the quantities $`D_3^{\alpha \beta }`$. When there is no force-constant disorder the $`D_3`$ terms vanish, and $`\stackrel{~}{𝐕}`$ becomes independent of $`\stackrel{}{q}`$. A $`\stackrel{}{q}`$-independent self-energy results, and we arrive at the CPA equations. To illustrate how to obtain the various matrix elements of the operators, we present the calculation of $`K_{01}^{(0)}`$ where $`R_1=(\frac{1}{2},\frac{1}{2},0)`$.

$`K_{01}^{(0)}`$ $`=`$ $`0f|\widehat{𝐊}|1f_0,`$
$`=`$ $`0f|\left(\widehat{𝐦}\omega ^2\widehat{𝚽}\right)|1f_0.`$
Using Eqs.(6)and (11),

$`0f|\widehat{𝐦}|1f_0`$ $`=`$ $`0,`$
and using Eqs.(7) and (11),

$`0f|\widehat{𝚽}|1f_0`$ $`=`$ $`f|\widehat{𝚽}_{01}|f_0,`$
$`=`$ $`(\sqrt{c_A}A_0|+\sqrt{c_B}B_0|),(\sqrt{c_A}A_1|+\sqrt{c_B}B_1|)|\{\varphi _{01}^{AA}\widehat{\eta }_0^A\widehat{\eta }_1^A+\varphi _{01}^{BB}\widehat{\eta }_0^B\widehat{\eta }_1^B+\varphi _{01}^{AB}\widehat{\eta }_0^A\widehat{\eta }_1^B+\varphi _{01}^{BA}\widehat{\eta }_0^B\widehat{\eta }_1^A\}|`$
$`\text{ }(\sqrt{c_B}|A_0\sqrt{c_A}|B_0),(\sqrt{c_A}|A_1+\sqrt{c_B}|B_1).`$
If we use the cartesian coordinates explicitly, then the xx component is, for example,

$`f|\mathrm{\Phi }_{01}^{xx}|f_0`$ $`=`$ $`\sqrt{c_Ac_B}\{c_Aa_{AA}c_Ba_{BB}+(c_Bc_A)a_{AB}\}`$
$`=`$ $`D_2^{xx}.`$
The other components can be calculated similarly.