A. Supercell method

In principle, exact enthalpies of formation of random phases can only be calculated with infinitely large SCs, where one can then properly invoke self-averaging without the need of configurational averaging over multiple instantiations. In practice, the enthalpies can only be calculated for finite SCs where the total number of atoms is not more than a few hundreds. To overcome this problem, Zunger et al.^{Reference Zunger, Wei, Ferreira and Bernard18} developed the so-called special quasirandom structure (SQS) method. The idea of this method is to mimic the random atomic configurations using much smaller SCs.

In the SQS method, a configuration (denoted by σ) is subdivided into a series of “figures”, or clusters, f = (k, m), where k is the number of vertices and m is the maximum distance between two vertices of f. In the binary case, each vertex (atomic site) i has a spin-like variable S _i, S _i = 1 if A (or B) occupies site i and S _i = −1 otherwise. For each figure, a product can be defined as ${\pi _f}\left( \sigma \right) = {D_f}\mathop \prod \limits_f {S_i}$, where D _f is the degeneracy of f. For a N-atom system, a correlation function for a figure can be defined as the average of π_f, i.e., ${\bar{\pi }_f}\left( \sigma \right) = {1 \over {N{D_f}}}{\pi _f}\left( \sigma \right)$. In the cluster expansion (CE) method, the correlation functions are used to expand physical properties. Using p _f to denote the effective cluster interaction, a configuration-dependent physical property can be expressed as $p\left( \sigma \right) = \mathop \sum \limits_f N{D_f}{\bar{\pi }_f}\left( \sigma \right){p_f}$, and the ensemble averaged physical property $p = \left\langle {p\left( \sigma \right)} \right\rangle = \sum\limits_f {N{D_f}} \left\langle {{{\bar{\pi }}_f}} \right\rangle {p_f}$.

Finite SCs can never have true random configuration. However, an error function can be defined to measure the difference between the randomness of a SQS SC and the true random configuration. Such an error function is usually defined as $\varepsilon \left( \sigma \right) = \sum\limits_f {N{D_f}} \left| {{{\bar{\pi }}_f}\left( \sigma \right) - {{\bar{\pi }}_{\rm{R}}}} \right|$. For binary cases, the ideal random configuration has ${\bar{\pi }_{\rm{R}}} = {\left( {2x - 1} \right)^k}$. Minimization of the error function yields SQS SCs. There are several codes that can generate SQS SCs, such as ATAT,^{Reference Van de Walle, Asta and Ceder19} UNCLE,^{Reference Lerch, Wieckhorst, Hart, Forcade and Muller20} CLUPAN,^{Reference Seko, Koyama and Tanaka21} etc.

In the case of multicomponent alloys, the correlation functions are defined differently. As described above, in the binary case, the correlation is computed by taking the product of spin variables, i.e., ${\pi _f}\left( \sigma \right) = {D_f}\prod\limits_f {\theta \left( {{S_i}} \right)}$ with cluster functions θ(S _i) = S _i. In the case of n-component alloys, it requires n − 1 cluster functions for the complete description of the correlations. Correspondingly, a figure with k vertices is connected with (n − 1)^k correlation functions. Sanchez et al.^{Reference Sanchez, Ducastelle and Gratias22} laid out the general multicomponent formalism, which was applied later by Wolverton and de Fontaine.^{Reference Wolverton and De Fontaine23} Readers are referred to these publications for more details.

A fully random solid solution retains the symmetry of the parent lattice. For example, an X-ray diffraction study of an ideal random fcc A-B solid solution alloy will only show the Bragg maxima of the underlying “average” fcc lattice. However, the usage of the SC approach breaks both the point-group symmetry and translational invariance. In addition, the electronic band structure of the SCs naturally contains multiple electron bands within a reduced Brillouin zone (BZ), making it difficult to compare the spectra signatures with the spectroscopic measurements, for example, angle resolved photoemission spectroscopy, which is typically defined within the BZ of the primitive cell. To overcome this difficulty, a state-of-the-art electron band unfolding method has recently been developed, to obtain the proper primitive cell spectral weight from the SC calculation. The idea was first proposed by Popescu et al.^{Reference Popescu and Zunger24} and the unfolding principle linking the SC band structure and the bands in the primitive cell reads as

$$A\left( {{\bi{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} }},E} \right) = \sum\limits_i {\left| {\left\langle {\Phi _{iK} |{\bi{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} }}} \right\rangle } \right|^2 } \delta \left( {\varepsilon _i - E} \right)\quad ,$$

where $A\left( {{\bi{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} }},E} \right)$ is the spectral function of a primitive cell, and ${\bi{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} }}$ and K are the wave vectors defined in a primitive cell and a SC, respectively. Φ_iK is the SC eigenstates for band index i. δ(ε_i − E) is the SC band structure. This method was first used to show the rapid disintegration of the valence band Bloch character in (In, Ga)N semiconductors and localized impurity band in Ga(N, P) semiconductors.^{Reference Popescu and Zunger24} Combined with the Wannier functions (treated as the basis function), this unfolding method is used to produce the proper spectral weight at very low computational cost.^{Reference Ku, Berlijn and Lee25} Employing this method, Ku and coworkers further demonstrated the gap opening due to orbital ordering in LaMnO₃^{Reference Ku, Berlijn and Lee25} and the important spectral change under Zn/Co doping in the Fe-based superconductor BaFe₂As₂.^{Reference Berlijn, Lin, Garber and Ku26} While the use of SCs to simulate the CSAs and HEAs are standard, this unfolding method has not been used so far to obtain the effects of disorder on the spectral function in HEAs.

The primary advantage of the SC method is that it takes the advantage of the standard electronic structure technology to obtain accurate total energies and to allow internal relaxation of the atomic positions that result from each atom having a different local configuration of surrounding atoms. Additionally, other effects beyond single site such as SRO, the spin disorder effect can be included in the SCs naturally, and the SC method provides a straightforward probe for various effects on the physical properties.

B. CPA and its derivatives

The simplest mean field theory for alloys is virtual crystal approximation (VCA).^{Reference Nordheim27,Reference Muto28} Within VCA, the inhomogeneous atomic potential is replaced by an effective potential, averaged over the local potentials for different species according to the compositions. Even though VCA sustains the translational invariance and seems to capture the electronic structure in dilute limit in disordered semiconductors, ferromagnets, and ferroelectric oxides,^{Reference de Gironcoli, Giannozzi and Baroni29–Reference Bellaiche and Vanderbilt32} it is not a realistic description of the scattering properties of disordered alloys. The reason is that VCA fails to capture the disorder-induced electronic band broadening and the proper inter-species hybridization effect that gives rise to the shifting of the bands.^{Reference Stocks, Temmerman, Szotek and Sterne33} Distinguished from the VCA, the CPA^{Reference Soven17} is derived from the configurational average of the scattering properties of different species embedded in an effective medium, with the self-consistency condition that this concentration-weighted average does not cause additional scattering beyond that of the effective CPA medium itself.^{Reference Gyorffy34} The basic idea is depicted in Fig. 2. When combined with multiple scattering Korringa–Kohn–Rostoker (KKR) Green’s function methods,^{Reference Korringa35,Reference Kohn and Rostoker36} the CPA condition can be represented as follows:

$${x_{\rm{A}}}{\tau ^{{\rm{A}},nn}} + {x_{\rm{B}}}{\tau ^{{\rm{B,}}nn}} = {\tau ^{{\rm{CPA,}}nn}}\quad ,$$

where x _A and x _B are concentrations of alloying elements A and B in a binary alloy. τ^A,nn is the scattering path operator which is defined from the Dyson equation, based on the single-site scattering matrix (t), and the free-space/host propagator (G ₀): τⁿⁿ′ = [t − 1 − G ₀]⁻¹; The scattering path operator for species A, transformed from a multiple scattering theory (MST) form of the Dyson equation, is τ^A,nn′ = [(t ^A)⁻¹ − (t ^CPA)⁻¹ + (τ^CPA)⁻¹]⁻¹. The component projected Green’s function $\left( {{G^\alpha }\left( {\vec{r},\vec{r}',E} \right)} \right)$ can be obtained from the corresponding scattering path operator (τ^A,nn′), and the total Green’s function $G\left( {\vec{r},\vec{r},E} \right)$ is then simply the concentration-weighted average of ${G^\alpha }\left( {\vec{r},\vec{r}',E} \right)$: $G\left( {\vec{r},\vec{r}',E} \right) = \mathop \sum \limits_\alpha {x_\alpha }{G^\alpha }\left( {\vec{r},\vec{r}',E} \right)$. It should be stressed that the CPA is a proper analytic theory, yielding a Green’s function and self-energy preserving both the translational invariance and the point-group symmetry of the underlying crystal structure. Furthermore, the knowledge of Green’s function allows us to calculate the physical observables of interest, such as density of state (n(E)), charge density (ρ(r)), and Bloch spectral function (BSF), A(E, k):

$$n\left( E \right) = - {1 \over \pi }\Im Tr\int {{\rm{d}}\vec{r}G\left( {\vec{r},\vec{r},E} \right)} \quad ,$$

$$\rho \left( r \right) = - {1 \over \pi }\Im Tr\int_{ - \infty }^{{E_{\rm{F}}}} {{\rm{d}}EG\left( {\vec{r},\vec{r},E} \right)} \quad ,$$

$$A\left( {E,k} \right) = - {1 \over \pi }\Im G\left( {E,k} \right)\quad .$$

FIG. 2.

Sketch of CPA: the embedding of either alloying components A (brown circles) or B (blue circles) within A_cB_1−c into the CPA effective medium should not cause additional scattering. 〈〉 means the statistical and configurationally averaged.

For more details of KKR-CPA, readers are referred to Ref. Reference Faulkner and Stocks37.

An important feature of the KKR-CPA formalism is that it involves a complex and energy-dependent self-energy, which describes both broadening of the electron bands and the band shift. Notably, this concept is absent in the standard SC technique. The KKR-CPA equation was first solved numerically by Stocks et al.,^{Reference Stocks, Temmerman and Gyorffy38} with applications to Ni–Cu alloys to determine the DOS, and further applied to study the electronic structure based on the concept of BSF—A(E, k) in disordered systems.^{Reference Faulkner and Stocks37} Figure 3(a) depicts the BSF at five different k points along the Δ direction in nonmagnetic Ni–Cu alloys. The peaks of the BSF indicate the positions of electronic bands which would be delta functions in a perfect crystal. The energy- and k-dependent broadening of the BSF is induced by the electron scattering due to disorder. The linewidth and variation of the BSF at the Fermi level is closely related to the k-dependent lifetime (τ_k), group velocity (v_k), and mean free path (λ_k) of electrons.^{Reference Stocks and Butler39} Figure 3(b) shows the BSF of Ni–Cu alloys, reproduced using SPR-KKR-CPA code,^{Reference Ebert, Ködderitzsch and Minár40,Reference Ebert41} along the high symmetry directions in the fcc lattice and gives the band structure in disordered alloys. The smearing of the bands indicates the electron scattering due to disorder.

FIG. 3.

(a) The calculated BSF A(E, k) for Cu_0.77Ni_0.23 alloys as a function of energy at five k points along Δ direction, using KKR-CPA. Adapted/Reproduced with permission from Ref. Reference Faulkner and Stocks37 (American Physical Society, 1980) (b) The BSF of Cu_0.77Ni_0.23 with Fermi level set to 0 eV. BSF is calculated using SPR-KKR-CPA.^{Reference Ebert, Ködderitzsch and Minár40,Reference Ebert41} Color encodes the magnitude of the BSF, given in arbitrary units.

There are two important energy scales for describing the strength of the disorder scattering: the band center mismatch among different species (δ) and the band width (W). For transition metals, the most relevant band center is simply the d-scattering resonance of the single-site potential, while the band width denotes how delocalized is the related orbital due to hybridization. The relative ratio between δ and W characterizes the strength of the disorder scattering. If δ ≪ W, the disorder scattering is weak and the electronic structure of the disordered alloy is in the VCA limit. This is the case for the disorder in delocalized sp bands where the corresponding band width is large. On the other hand, if δ ∼ W, then a strong disorder scattering of electrons is expected, and as a result, the electron bands are smeared out. If δ ≫ W, then the electron bands are no longer continuous and a split-band limit is reached. It can be shown that the CPA properly extrapolates between the weak and strong scattering limits. We note that both δ and W parameters are spin-dependent in general. For magnetic materials, the same argument is still valid within each spin channel. Take magnetic NiCo alloy as an illustration, as sketched in Fig. 4(a). Here, δ is much smaller than W in the majority spin channel, whilst δ is comparable with W in the minority spin channel. The reason for this behavior is that Ni and Co have different valence electron fillings and different local moments. The latter results in a different shift of the minority spin channels for two species. While the band centers in the majority spin channel are aligned, the spin-minority band center of Co is shifted, relative to that of Ni, by the difference of the onsite exchange splitting between two elements. Because of the Stoner-like proportionality between the moment and exchange splitting, band center mismatch in the spin minority channel is increased by approximately 1 eV (See Fig. 4(a), Adapted from Ref. Reference Mu, Samolyuk, Wimmer, Troparevsky, Khan, Mankovsky, Ebert and Stocks42).

FIG. 4.

(a) Energy scales for disorder scattering: The species-resolved density of state in NiCo is shown here. ε^σ_d denotes the 3d band centers for spin σ, where σ can be either ↑ or ↓. The band center mismatch δ_σ is defined as max{ε^σ_d (i)} − min{ε^σ_d (i)} with i denoting different species. The exchange splitting is Δ_Exch:= ε^↓_d − ε^↑_d. With elevated Δ_Exch from Ni to Co, the local moment is enhanced due to M = I × Δ_Exch, where the Stoner parameter I ≈ 1.0 μ_B/eV is a constant. (b) Fermi surface of NiCo (ΓXW slice) in the spin majority channel and spin minority channel: Sharper bands in the spin majority channel at the Fermi level while bands in the spin minority channel are more smeared out. The color indicates the intensity of the BSF, given in arbitrary units. Fermi surface is calculated using SPR-KKR-CPA.^{Reference Ebert, Ködderitzsch and Minár40,Reference Ebert41} Adapted/Reproduced with permission from Ref. Reference Mu, Samolyuk, Wimmer, Troparevsky, Khan, Mankovsky, Ebert and Stocks42.

Furthermore, the CPA not only provides a powerful and convenient probe of electronic structures of disordered alloys but also offers an opportunity to explore other forms of disorder that can also be mapped on to a random solid solution or lattice gas. One example is the disorder local moment (DLM), used to model the paramagnetic state of magnetic metals and alloys at infinite temperature.^{Reference Pindor, Staunton, Stocks and Winter43,Reference Gyorffy, Pindor, Staunton, Stocks and Winter44} The underlying idea of the DLM is to model the paramagnetic state of magnetic materials having very defined local moments by equally populated random distribution of spin-up and spin-down moments. The configurational average over random spin orientations is done within CPA. Along with the instantaneous spin disorder is the fluctuation of the size of the spin moment, usually giving the reduced local moment according to the Stoner theory of ferromagnetism. The DLM state, solved self-consistently, correctly describes the electronic structure and size variation of the local moments of Fe and Co in the paramagnetic state.^{Reference Pindor, Staunton, Stocks and Winter43,Reference Staunton, Gyorffy, Pindor, Stocks and Winter45} At finite temperatures lower than the transition temperature, the transverse spin fluctuations can still be treated in a mean field way using the effective Weiss field, and the spin distributions at arbitrary temperatures are described by a discrete set of orientations, each of which are treated as CPA components whose population is determined by the Weiss field,^{Reference Deák, Simon, Balogh, Szunyogh, Dos Santos Dias and Staunton46} as illustrated in Fig. 5(b). This spin disorder can be further combined with the chemical disorder to pursue the finite temperature electronic structure,^{Reference Belashchenko, Weerasinghe, Mu and Pujari47,Reference Mankovsky, Polesya, Chadova, Ebert, Staunton, Gruenbaum, Schoen, Back, Chen and Song48} electrical transport,^{Reference Mankovsky, Polesya, Chadova, Ebert, Staunton, Gruenbaum, Schoen, Back, Chen and Song48,Reference Ebert, Mankovsky, Chadova, Polesya, Minar and Koedderitzsch49} magnetic phase diagram,^{Reference Pujari, Larson, Antropov and Belashchenko50} and magnetic anisotropy^{Reference Staunton, Ostanin, Razee, Gyorffy, Szunyogh, Ginatempo and Bruno51–Reference Zhuravlev, Antropov and Belashchenko53} in magnetic metals and disordered alloys.

FIG. 5.

Configurational averaging for thermal atomic displacements (a) and thermal spin fluctuations (b). A discrete set of atomic displacements or spin orientation vectors at a given temperature are configurationally averaged using the CPA, leading to a periodic effective medium. Adapted/Reproduced with permission from Ref. Reference Ebert, Mankovsky, Chadova, Polesya, Minar and Koedderitzsch49 (American Physical Society, 2015).

Another application of the CPA concept is to model the thermal atomic displacements u_i(T). The amplitude of the thermal atomic displacement can be determined by the root mean square displacement at a certain temperature, obtained from the Debye model with the Debye temperature averaged from the Debye temperature of each alloying element according to

$$\left\langle {{u^2}} \right\rangle = {{3{\hbar ^2}} \over {4{\pi ^2}M{k_{\rm{B}}}{\Theta _{\rm{D}}}}}\left[ {{{\Phi \left( {{{{\Theta _{\rm{D}}}} \over T}} \right)} \over {{{{\Theta _{\rm{D}}}} \over T}}} + {1 \over 4}} \right]\quad ,$$

where k _B is the Boltzmann constant, Θ_D is the averaged Debye temperature, Φ(x) is the Debye function, and M is the averaged mass. Additionally, the thermal atomic displacement is randomly distributed over the space. Assuming uncorrelated site-to-site thermal displacements, one can replace the atomic displacement by a discrete set of {u_i} uniformly distributed over a solid angle 4π, with each discretized angular displacement treated as a CPA component [sketched in Fig. 5(a)]. This idea has been used in the study of the electrical resistivity and Gilbert damping parameter at finite temperatures.^{Reference Ebert, Mankovsky, Chadova, Polesya, Minar and Koedderitzsch49} The benefit of the CPA treatment of the thermal displacement is that the configurational average is obtained in the thermodynamic limit, thereby obviating the necessary usage of large SCs to include the full distribution of local displacements.

Another important application of the CPA method is to study the structural stability of alloys by calculating the total energies. In particular, the high computational efficiency of CPA renders it advantageous to study HEAs, where the large number of atomic species (≥5) and large chemical disorder present great challenges to the SC method. In practice, the CPA method is widely used in conjunction with the exact muffin-tin orbital (EMTO)^{Reference Vitos54,Reference Vitos, Abrikosov and Johansson55} method to study HEAs. The EMTO method adopts optimized overlapping muffin-tin (MT) potentials and uses full charge density technique to give a better description of metallic systems than the traditional MT KKR method. For example, using EMTO-CPA method, Ma et al.^{Reference Ma, Grabowski, Körmann, Neugebauer and Raabe56} calculated the fcc, bcc, and hcp phases of CoCrFeMnNi HEAs. By explicitly incorporating the free energy contributions from vibrational, electronic, and magnetic excitations, they found that the stable phase of CoCrFeMnNi transforms from the hcp to fcc structure as the temperature increases. Another advantage of the CPA method is that it only requires a small cell regardless of the concentration of the elements and therefore provides a convenient tool to study doping effects. For example, Tian et al.^{Reference Tian, Delczeg, Chen, Varga, Shen and Vitos57} investigated the structural stability of Al_xCoCrFeNi for different concentrations of Al. They predicted that Al_xCoCrFeNi would transform from the fcc to bcc structure as x increases from 0.6 to 1.2, which is in good agreement with experiment results. Other than structure stability, the total energies calculated can also be used to extract elastic constants from the energy–strain relations. For cubic systems, only 3 elastic constants, C ₁₁, C ₁₂, and C ₄₄ are independent, from which the bulk modulus, Young’s modulus, and Poisson’s ratio can then be determined. Using the EMTO + CPA method, the elastic constants of different HEAs have been calculated in a series of papers by Tian et al.^{Reference Tian, Delczeg, Chen, Varga, Shen and Vitos57–Reference Tian, Varga, Shen and Vitos60}

Despite the many successes of CPA, it is limited to systems without SRO due to its single-site nature. A path forward is to go beyond the single-site approximation as in the nonlocal CPA (NL-CPA). The basic idea of NL-CPA is similar to CPA except that NL-CPA performs configurational averaging over embedding atomic clusters rather than a single site within an effective medium. The SRO can be included by embedding a cluster with SRO explicitly. NL-CPA was first implemented within KKR formalism and was applied to study the effect of SRO parameters on the electronic structure of the classical alloy bcc Cu_0.5Zn_0.5.^{Reference Rowlands, Staunton, Györffy, Bruno and Ginatempo61,Reference Rowlands, Ernst, Györffy and Staunton62} Another interesting application is in Fe_0.5Pt_0.5 where the spin moment of Pt can be induced by increasing the number of nearby Fe atoms.^{Reference Ködderitzsch, Ebert, Rowlands and Ernst63} More details about the NL-CPA and other cluster-based CPA can be found in this review and references therein.^{Reference Rowlands64}

A. Boltzmann transport

Perhaps, the most intuitive treatment of electron transport in a dc field is the semiclassical Boltzmann transport equation (BTE).^{Reference Ashcroft and Mermin65} The conduction electrons are treated semiclassically, meaning that the motion of electrons with band index n can be described by a wave packet centered at position r with moment k, and the dynamics of electrons between successive scattering events under a slow-varying (as compared with the dimension of the wave packet) time-dependent field E(r, t) or H(r, t) are given by the semiclassical equations:

$$\dot{r} = {v_n}\left( k \right) = {{\partial {\varepsilon _n}\left( k \right)} \over {\hbar \partial k}}\quad ,$$

$$\hbar \dot{k} = - {\rm{e}}\left[ {E\left( {r,t} \right) + {1 \over c}{v_n}\left( k \right) \times H\left( {r,t} \right)} \right]\quad .$$

In the phase space, the electron distribution at a particular position and momentum can be represented by the electron distribution function $\left( {f\left( {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} ,t} \right)} \right)$, whose meaning is the number of electrons at time t around $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r} \to \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r} + {\rm{d}}\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r}$, with momentum $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} \to \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} + {\rm{d}}\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k}$, represented by $f\left( {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} ,t} \right){{{\rm{d}}\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r} {\rm{d}}\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} } / {4{\pi ^3}}}$. The electron distribution function can be further divided into the equilibrium contribution f _n(r, k) (Fermi-Dirac distribution) plus the nonequilibrium contribution which controls the transport properties. The evolution of the nonequilibrium electron distribution function, f(r, k, t) − f _n(r, k), is determined by diffusion effects, field-driven effects, and the scattering events (both scattering-in and scattering-out terms): ${{{\rm{d}}f\left( {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} ,t} \right)} \over {{\rm{d}}t}} = \left. {{{\partial f\left( {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} ,t} \right)} \over {\partial t}}} \right|_{{\rm{scatter}}} - \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over v} _k \left. {{{\partial f\left( {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} ,t} \right)} \over {\partial \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r} }}} \right|_{{\rm{diff}}} - \dot k\left. {{{\partial f\left( {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} ,t} \right)} \over {\partial k}}} \right|_{{\rm{field}}}$. We note that the band index n is assumed to be independent of the electron motion or scattering events. Therefore, the semiclassical theory ignores the inter-bands transition and is only suitable for transport under weak perturbations.

In the steady state, the time variation of f(r, k, t) under the above effects vanishes, leading to a simplified equation: ${\left. {{{\partial f\left( {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} ,t} \right)} \over {\partial t}}} \right|_{{\rm{scatter}}}} = {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over v} _k}{\left. {{{\partial f\left( {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} ,t} \right)} \over {\partial \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r} }}} \right|_{{\rm{diff}}}} + \dot{k}{\left. {{{\partial f\left( {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} ,t} \right)} \over {\partial k}}} \right|_{{\rm{field}}}}$. In a homogeneous dc electrical field, all the spatial gradient terms of f _n(r, k, t) also vanish. In relaxation-time approximation, a scattering rate $\left( {{1 \over {{\tau _k}}}} \right)$ is introduced, ${\left. {{{\partial f\left( {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over r} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over k} ,t} \right)} \over {\partial t}}} \right|_{{\rm{scatter}}}} = {1 \over {{\tau _k}}}\left( {f - {f_0}} \right)$, regardless of the details of scattering events. The resulting dc conductivity tensor (σ_αβ), derived from the BTE, reads as

$${\sigma _{{\rm{\alpha \beta }}}} = {{{{\rm{e}}^2}} \over {4{\pi ^3}\hbar }}\int {{\rm{d}}{S_{\rm{F}}}} {{{\tau _k}\left( {{\varepsilon _{\rm{F}}}} \right){{\left( {{v_k}} \right)}_{\rm{\alpha }}}{{\left( {{v_k}} \right)}_{\rm{\beta }}}} \over {\left| {{v_k}} \right|}}\quad ,$$

where τ_k(ε_F) is the relaxation time at the Fermi level for a particular k point, and the surface integral is performed only at the Fermi surface S _F.

The KKR-CPA method calculates the electronic structure [BSF A _B(k, ε)] of disordered alloys in both dilute limit and concentrated limit and provides a feasible way to calculate the above residual conductivity formula by estimating the relaxation time (τ_k(ε_F)) and Fermi velocity (v _k(ε_F)) from KKR-CPA electronic structure calculations. This approach was first applied to Ag_xPd_1−x to evaluate the residual resistivity at 0.2 ≤ x ≤ 0.8 with the scattering-in term neglected.^{Reference Stocks and Butler39} Good agreement with experiments was obtained, and an abnormal concentration-dependent lifetime and velocity was found that is due to the energy dependence of disorder self-energy at the top of d bands. However, as pointed out by Butler,^{Reference Butler66} in spite of excellent agreement with experiment, this method suffers from two problems: first, this method can only be applied to systems with well-defined bands at the Fermi level; second, the scattering-in term is neglected, which amounts to neglect the vertex correction, as shown by Butler.^{Reference Butler66} For the first problem, we should note that even though the conventional BTE is limited to the dilute limit, it can be extended to concentrated alloys by using the BSF for the electronic structure, provided that well-defined dispersion relationship, having small disorder broadening, exists at the Fermi level thereby allowing a relative unambiguous definition of group velocity. However, in the case of strong disorder scattering, the Fermi surface can be so smeared out that this theory is no longer applicable. For the second point, the neglect of the scattering-in term may lead to large errors when the electron scattering is dominated by s–p or p–d scattering.

The modified Boltzmann conductivity, including the scattering-in term, was first derived by Butler^{Reference Butler66} and was further implemented by Mertig and coworkers to study the residual resistivity in the 3d transition metal element doped nonmagnetic Cu or ferromagnetic Ni, Co (see Review Ref. Reference Mertig, Zeller and Dederichs67 and references therein). This approach gave very good agreement with experiment for 3d transition metal defects in ferromagnetic nickel.^{Reference Mertig, Zeller and Dederichs67} To get the residual resistivity of a ferromagnet, Mertig et al.^{Reference Mertig, Zeller and Dederichs67} used Mott’s two-channel model,^{Reference Mott68} assuming that the electrons travel independently in separate spin channels. This two-channel model neglects the spin-flipping scattering which might be induced by spin disorder or spin–orbit coupling. The spin-flip scattering time τ_sf was further included and applied to doped Cu systems.^{Reference Fedorov, Zahn, Gradhand and Mertig69} In addition, Fert et al.,^{Reference Fert and Campbell70,Reference Fert and Campbell71} in the extensive transport studies of additional 3d element defects in pure nickel, found that Cr defects, which antiferromagnetically couple to Ni, are very effective in raising the residual resistivity. This study further stimulated the discovery of giant magnetoresistance (GMR) in FM-AFM multilayers,^{Reference Baibich, Broto, Fert, Van Dau, Petroff, Etienne, Creuzet, Friederich and Chazelas72} leading to the Nobel prize.

B. Kubo formula and its derivatives

As mentioned, the Boltzmann transport theory cannot be used in systems with no clear definition of the group velocity and relaxation time at the Fermi surface. To overcome this problem, Butler^{Reference Butler66} proposed a more general approach based on Kubo’s linear response formula^{Reference Kubo73} for the current–current correlation functions. The original Kubo formula at T = 0 K leads to the Kubo–Greenwood (KG) equation, in which the conductivity tensor is expressed as

$${\sigma _{\mu \nu }} = {\hbar \over {\pi N\Omega }}Tr\left\{ {\left\langle {j_\alpha ^\mu {\mathop{\rm Im}\nolimits} {G^ + }\left( {{E_{\rm{F}}}} \right)j_\beta ^\nu {\mathop{\rm Im}\nolimits} {G^ + }\left( {{E_{\rm{F}}}} \right)} \right\rangle } \right\}\quad ,$$

where G ⁺ is the retarded Green function and the average 〈〉 refers to the configurational average. In this pioneering work by Butler, the KG equation for disordered alloys is formulated within KKR-CPA. The configurational average results in the decomposition of the conductivity tensor to the site-diagonal and off-diagonal components:

$${\sigma ^0}_{\mu \nu } = {{ - 4{m^2}} \over {\pi {\hbar ^3}\Omega }}\sum\limits_\alpha {{x_\alpha }} \sum\limits_{{L_1},{L_2},{L_3},{L_4}} {J_{{L_4}{L_1}}^{\alpha \mu }\tau _{{L_1}{L_2}}^{{\rm{CPA}},00}J_{{L_2}{L_3}}^{\alpha \nu }\tau _{{L_3}{L_4}}^{{\rm{CPA}},00}} \quad ,$$

$${\sigma ^1}_{\mu \nu } = {{ - 4{m^2}} \over {\pi {\hbar ^3}\Omega }}\sum\limits_{\alpha \beta } {{x_\alpha }{x_\beta }} \sum\limits_{{L_1},{L_2},{L_3},{L_4}} {J_{{L_4}{L_1}}^{\alpha \mu }{{\left[ {{{\left\{ {1 - \chi \omega } \right\}}^{ - 1}}\chi } \right]}_{{L_1},{L_2},{L_3},{L_4}}}J_{{L_2}{L_3}}^{\beta \nu }} \quad ,$$

where the quantity $J_{LL'}^{\alpha \mu }$ denotes the matrix elements of the μ-component of the current density operator j for species α with concentration x _α. L stands for the angular moment quantum number {l, m}. {1 − χω}⁻¹ corresponds to the scattering-in term, or the vertex correction. The KG formula is exact in the linear region, regardless of the scattering strength. This formula was first implemented in KKR-CPA by Swihart and coworkers^{Reference Swihart, Butler, Stocks, Nicholson and Ward74} and since been applied to many alloys. Notably, this is the first implementation of the alloy transport theory that included the effects of vertex correction. The impact of the vertex correction on the residual resistivity of Ag_xPd_1−x is largest on the Ag-rich side where the Fermi energy lies in the sp bands. We also note that the KG equation was also implemented by Turek et al.,^{Reference Turek, Kudrnovský, Drchal, Szunyogh and Weinberger75} with vertex correction included, in the framework tight-binding linearized muffin-tin orbital method (TB-LMTO)-CPA.

Since the seminal work by Butler, various derivatives of the Kubo-type formula were derived. The KG equation was implemented in the basis of NL-CPA by Tulip et al.^{Reference Tulip, Staunton, Lowitzer, Ködderitzsch and Ebert76} to explore the SRO effect on the electron transport in Ag_xPd_1−x [see Fig. 6(a)]. It turns out that the single-site CPA already worked excellently in Ag_xPd_1−x. Intuitively, SRO should reduce the residual resistivity in most scenarios. This reduction is apparent in Cu_xZn_1−x. However, there are some exceptions where the residual resistivity increases with SRO. An example is a series of so-called K-state alloys, usually in Ni, Cr containing alloys, such as NiCoCr, Ni25Cr, etc. The reason for this bizarre SRO-enhanced resistivity is attributed to the particular feature of the electronic structure under SRO: reduced amplitude of d electron hopping in the spin majority channel between the nearest atomic pair in the presence of SRO.^{Reference Lowitzer, Ködderitzsch, Ebert, Tulip, Marmodoro and Staunton77}

FIG. 6.

(a) Residual resistivity of Ag_cPd_1−c as a function of concentration c. ρ₀ was calculated using CPA (green) and nonlocal CPA (red) with the vertex correction included. Experimental data are shown (blue line with filled circles). Adapted/Reproduced with permission from Ref. Reference Tulip, Staunton, Lowitzer, Ködderitzsch and Ebert76 (American Physical Society, 2008) (b) Temperature-dependent longitudinal resistivity in Fe: The thermal atomic displacements are denoted as “vib”, and the thermal spin fluctuations are denoted as “fluct”. Spin fluctuations at a given temperature are obtained from the Monte Carlo simulation or experiments. Adapted/Reproduced with permission from Ref. Reference Ebert, Mankovsky, Chadova, Polesya, Minar and Koedderitzsch49 (American Physical Society, 2015).

In the above discussion, relativistic effects are neglected, which is consistent with the two-channel model. A corresponding full relativistic theory was explored by solving the KKR-CPA equation based on the Dirac Hamiltonian, which gives rise to the coupling between the two spin channels. The KG transport equation in the corresponding full relativistic frame was first implemented by Banhart and coworkers.^{Reference Banhart, Bernstein, Voitländer and Weinberger78} Applying it to magnetic NiFe alloys, Banhart et al.^{Reference Banhart, Vernes and Ebert79} found that residual resistivity is not sensitive to the scalar relativistic effect, yet is sensitive to the strength of spin–orbit coupling—which can be manipulated by scaling the speed of light. The importance of the spin–orbit coupling indicates the non-negligible spin-flipping conductance in NiFe and the breakdown of the two-channel model.

As mentioned in Sec. III.A, the effect of finite temperature lattice vibrations or of spin fluctuations can be included within single-site KKR-CPA, by employing the CPA analogy to configurational average of the uncorrelated thermal displacements, or thermal spin orientations (see Fig. 5). The thermal effect on the conductivity tensor of Ni and Fe were investigated by Ebert et al.^{Reference Ebert, Mankovsky, Chadova, Polesya, Minar and Koedderitzsch49} Figure 6(b) illustrates the longitudinal resistivity of Fe as a function of temperature. When the finite temperature moments are taken from the experiment, the calculated residual resistivity is in excellent agreement with the measured one. However, it is shown that the Monte Carlo simulation overestimate the local moments and the residual resistivity. In addition, the thermal displacement effect and the thermal spin disorder effect on electrical transport can be further separated out, allowing Ebert et al.^{Reference Ebert, Mankovsky, Chadova, Polesya, Minar and Koedderitzsch49} to demonstrate the violation of Matthiessen’s rule. Employing this method, Mankovsky et al.^{Reference Mankovsky, Polesya, Chadova, Ebert, Staunton, Gruenbaum, Schoen, Back, Chen and Song48} calculated the temperature-dependent transport properties of the ordered FeRh alloys and explained the abrupt resistivity drop at the antiferromagnetic to ferromagnetic phase transition at T = 320 K.

According to Neumann’s principle,^{Reference Neumann and Meyer80} all physical properties must sustain the symmetry of the parent lattice. For example, σ = D ⁻¹(R)σD(R), where D(R) is the matrix representation of the point group {R}. Note that the presence of magnetic ordering will lower the point group and usually leads to a subgroup of the parent group. Therefore, the symmetry of the conductivity tensor changes accordingly. For a nonmagnetic or paramagnetic material with underlying cubic symmetry, the σ tensor is diagonal and isotropic. If the system undergoes a magnetic phase transition, say, from paramagnetic to ferromagnetic with the spin quantization axis along z, then the conductivity tensor takes the form:

$$\sigma = \left( {\matrix{ {{\sigma _{xx}}} & {{\sigma _{xy}}} & 0 \cr { - {\sigma _{xy}}} & {{\sigma _{xx}}} & 0 \cr 0 & 0 & {{\sigma _{zz}}} \cr} } \right)\quad .$$

The diagonal terms of the conductivity tensor are distinguished by components parallel to the spin (σ_zz) and components perpendicular to the spin (σ_xx, σ_yy), as does the resistivity tensor. The anisotropic magnetoresistance (AMR) ratio is defined ${{\Delta \rho } \over \rho } = {{\hbox( {{\rho _ \bot } - {\rho _{\parallel} }} \hbox)} / {\hbox( {{2 / {3{\rho _ \bot }}} + {1 / {3{\rho _{\parallel} }}}} \hbox)}}$. We note that the AMR arises from the spin–orbit interaction, so a full relativistic treatment of the KKR-CPA is needed for this term. In addition, the off-diagonal term σ_xy is the anomalous Hall conductivity, which also stems from the spin–orbit interaction. Distinguished by the origins of scattering, σ_xy can be decomposed into the intrinsic contribution, which is related to the intrinsic electronic band structure and is equivalent to the Berry curvature^{Reference Nagaosa, Sinova, Onoda, MacDonald and Ong81} and also the extrinsic contributions, from side-jump scattering^{Reference Berger82} and skew scattering.^{Reference Smit and Volger83} The off-diagonal terms were solved using the more general Kubo–Streda equation,^{Reference Streda84,Reference Crépieux and Bruno85} which was then implemented in KKR-CPA by Lowitzer et al.^{Reference Lowitzer, Ködderitzsch and Ebert86} with applications to ferromagnetic FePd and NiPd alloys.

As another interesting transport phenomenon, the spin Hall effect^{Reference Sinova, Valenzuela, Wunderlich, Back and Jungwirth87} is receiving more and more attention due to potential applications in spintronic devices (spin–orbit torque^{Reference Liu, Pai, Li, Tseng, Ralph and Buhrman88}) to make conversions between the spins and the current, and eventually to reduce the energy consumption of the logic devices. The evaluation of the spin Hall conductivity in disordered alloys can be carried out by the Kubo–Streda equation as well. The detail is beyond the scope of this review. Here, we only refer to some outstanding studies,^{Reference Vernes, Györffy and Weinberger89,Reference Lowitzer, Gradhand, Ködderitzsch, Fedorov, Mertig and Ebert90} where the Spin Hall conductivity was derived from the modified Kubo–Streda equation and further implemented in the basis of KKR-CPA for disordered alloys.

C. Nonequilibrium Green’s function method

Besides Kubo formula, another approach to calculate the residual resistivity of metals and alloys is to use the Landauer–Büttiker formulation^{Reference Landauer91,Reference Büttiker92} for the conductance. Using this approach, the conductance can be obtained using SC models of different lengths with fixed lateral cross section. The conductivity or the resistivity can then be extracted by fitting the results to Ohm’s Law. A SC, whose geometry consists of a layer with periodic boundary condition along two lateral directions, is illustrated in Fig. 7(a). The layer geometry consists of a central active region that is embedded between semi-infinite leads. In the active region, various forms of disorders, such as chemical disorder, spin disorder, and atomic fluctuations, can be accounted for. The layer Green’s function is solved self-consistently. Assuming interactions between layers are short-ranged and limited within the two nearest layers, the Hamiltonian is tri-diagonal, which greatly reduces the computational expense. The conductance between spin channels (μ, ν), can be written as in terms of the transmission probability T, self-energy of the left lead (Γ^L) and the right lead (Γ^R), and the retarded and advanced Green’s functions (G ⁻) as

$${G_{\mu \nu }} \sim {{2{e^2}} \over \hbar }M{T_{\mu \nu }} = {{2{e^2}} \over \hbar }MTr\left\{ {\Gamma _\mu ^{\rm{R}}G_{\mu \nu }^ + \Gamma _\nu ^{\rm{L}}G_{\mu \nu }^ - } \right\}\quad .$$

FIG. 7.

(a) Slab geometry to calculate the Landauer–Büttiker conductance. Adapted/Reproduced with permission from Ref. Reference Wysocki, Sabirianov, van Schilfgaarde and Belashchenko93 (American Physical Society, 2009). (b) Conductance of Fe (black) and Ni (red) at a function of thickness of monolayers at T _c and 0.25T _c. Adapted/Reproduced with permission from Ref. Reference Wysocki, Sabirianov, van Schilfgaarde and Belashchenko93 (American Physical Society, 2009). (c) Conductance of ferromagnetic Fe and Gd with uncorrelated thermal atomic displacements included. Adapted/Reproduced with permission Ref. Reference Glasbrenner, Pujari and Belashchenko96 (American Physical Society, 2012).

Based on the TB-LMTO, this SC method has been applied to explore the spin disorder resistivity in Ni and Fe^{Reference Wysocki, Sabirianov, van Schilfgaarde and Belashchenko93} [as can be seen in Fig. 7(b)], temperature-dependent resistivity in transition metals, heavy element rare-earth metals^{Reference Glasbrenner, Belashchenko, Kudrnovský, Drchal, Khmelevskyi and Turek94–Reference Glasbrenner, Pujari and Belashchenko96} [as can be seen in Fig. 7(c)], inverse spin Hall effect in Pt/Py bilayers,^{Reference Wang, Wesselink, Liu, Yuan, Xia and Kelly97} and spin injections,^{Reference Belashchenko, Glasbrenner and Wysocki98} to name a few. The configurational average required in alloys is performed manually by averaging many SCs with different random distributions of the alloying elements. A large computational cost notwithstanding, this direct method can be applied to arbitrary complex systems and is straightforward to include different scattering effects, such as chemical disorder, short-range disorder, species-related displacement fluctuations, and complex spin structures. For example, magnetic SRO has been included in studies of spin-disorder resistivity of iron and nickel by Wysocki et al.^{Reference Wysocki, Sabirianov, van Schilfgaarde and Belashchenko93} It is found that while near neighbor magnetic SRO reduces the spin-disorder resistivity in both metals, it is a small effect in Fe but a sizable effect in Ni. Although this approach is a powerful approach to describe finite temperature transport by including the thermal atomic displacements and/or noncollinear spin textures, it has not yet been applied to disordered alloys. This method is also extensively implemented in the KKR formalism to investigate the GMR^{Reference Butler, Zhang, Nicholson and MacLaren99} and the tunneling magnetoresistance^{Reference Butler, Zhang, Schulthess and Maclaren100} in the classical trilayer system.

D. Applications to Cantor-Wu alloys

Using ab initio method, Mu et al.^{Reference Mu, Samolyuk, Wimmer, Troparevsky, Khan, Mankovsky, Ebert and Stocks42} explored the intriguing transport properties in the set of Cantor-Wu alloys. The ρ₀ measurement^{Reference Jin, Sales, Stocks, Samolyuk, Daene, Weber, Zhang and Bei8,Reference Zhang, Stocks, Jin, Lu, Bei, Sales, Wang, Béland, Stoller, Samolyuk, Caro, Caro and Weber15} of a subset of Cantor-Wu alloys shows that values of ρ₀ break into two subgroups of low (ρ₀ < 10 μΩ cm) and high (ρ₀ > 75 μΩ cm) resistivity alloys, with the least and most resistive alloys differing by almost two orders of magnitude (ρ₀ (NiCo) = 1.3 μΩ cm; ρ₀ (CrPdFeCoNi) = 124.8 μΩ cm). In particular, the high-ρ₀ of CrPdFeCoNi is close to the MIR limit, which is characterized by a value of the electron mean free path (λ_e) comparable to lattice spacing.^{Reference Ioffe and Regel101–Reference Gunnarsson, Calandra and Han103} This observation suggests that although the carrier densities are high, as in normal metals, λ_e can be controlled, from 10 to 10³ Å, by the numbers/types of alloying elements. Therefore, given the unusual behavior of the transport properties, it is important to understand the underlying scattering mechanisms that give rise to the large variation of λ_e in Cantor-Wu alloys. This is the subject of the remainder of this section.

Ultimately, all electron scattering arises from disorder-induced site-to-site potential fluctuations. However, evaluating electrical resistivity is a formidable task because of the multiple electron scattering mechanisms that underlie these fluctuations in the crystal potential, in particular: intrinsic site disorder, local atomic relaxation, magnetic moment scattering, and SRO. These effects are sketched in Fig. 8(a). In the single-site approximation, electron scattering due to site-disorder of the local potential (δ) can be thought of as that resulting from the random distribution of elements with different nuclear charges across the sites of a periodic lattice—giving an inhomogeneous local potential. In addition, in the presence of different magnetic moments across the sites, conduction electrons experience an inhomogeneous onsite exchange field (Δ_Exch), which gives rise to additional spin-channel dependent fluctuations in the site potential. The KKR-CPA method naturally includes the effects on the electronic structure and Fermi surface from both forms of site-disorder [δ, Δ_Exch].

FIG. 8.

(a) Sketch of various disorders in alloys: chemical disorder (colored circles for different species), atomic displacement away from its ideal site (white arrows) and local moments (black arrows). (b) Species-resolved atomic displacement distribution in selected Cantor-Wu alloys: NiFeCo, NiFeCoCr, and NiFeCoCrMn. (c) Species-resolved local moment distribution in selected Cantor-Wu alloys: NiFeCo, NiFeCoCr, and NiFeCoCrMn. Adapted/Reproduced with permission from Ref. Reference Mu, Samolyuk, Wimmer, Troparevsky, Khan, Mankovsky, Ebert and Stocks42.

Going beyond the idealized single-site approximation, this picture is modified by including the effects of displacement scattering resulting from the atomic displacements away from their ideal lattice sites. For alloys involving the mid period 3d-transition metal elements that all have potential to form complicated magnetic moment distributions, magnetic scattering becomes important. In actuality, extra magnetic scattering beyond single-sites can arise from both moment formation and how they couple amongst themselves—ferromagnetically, mixed ferromagnetic/antiferromagnetic, noncollinear, or spin-glass. The effects of the above scattering mechanisms on ρ₀ can be further explored by the extensions of KKR-CPA as discussed previously.

1. Onsite disorder [δ_σ, Δ_Exch]

Using ab initio KKR-CPA method and KG formula to calculate the residual resistivity of the full set of Cantor-Wu alloys, Mu et al.^{Reference Mu, Samolyuk, Wimmer, Troparevsky, Khan, Mankovsky, Ebert and Stocks42} characterized the most important scattering events as the site-diagonal disorder. Consistent with the transport measurements,^{Reference Jin, Sales, Stocks, Samolyuk, Daene, Weber, Zhang and Bei8,Reference Zhang, Stocks, Jin, Lu, Bei, Sales, Wang, Béland, Stoller, Samolyuk, Caro, Caro and Weber15} the Cantor-Wu alloys are characterized into two groups according to the magnitude of residual resistivity: low-ρ₀ alloys (NiPd, NiCo, NiFe, and NiFeCo), with ρ₀< 15 μΩ cm, and high-ρ₀ alloys (the remainder of the set of Cantor-Wu alloys) The ρ₀ of the full set of Cantor-Wu alloys is summarized in Fig. 9(a). It is interesting to note that alloying of Mn and/or Cr with combinations of Ni, Fe, and Co, always results in a high residual resistivity and an electron mean free path comparable to the lattice spacing.

FIG. 9.

(a) Effect of site-diagonal disorder: ρ₀ of Ni-based multicomponent alloys, from both LDA (blue bars) and GGA-PBE (green bars) exchange correlation functionals, compared with experimental data. Larger value of ρ₀ from GGA compared with LDA is due to the larger predicted local moments. (b) Modified resistivity by the displacement fluctuation (DF) in all considered alloys. (c) Effect of spin disorder (SD) on the residual resistivity in selected alloys. Adapted/Reproduced with permission from Ref. Reference Mu, Samolyuk, Wimmer, Troparevsky, Khan, Mankovsky, Ebert and Stocks42.

As mentioned in Sec. II, the effect of disorder scattering on the residual resistivity can be explained by the disorder-induced smearing of the electron structure. On the one hand, for alloys containing only elements from the post half-filled 3d shell-Fe, Co, and Ni, the majority-spin 3d band centers are aligned due to minimization of the kinetic energy. This results in a small δ_↑. Because different local moments form on different species and they couple ferromagnetically, the additional exchange splitting leads to a large band center mismatch (δ_↓) in the minority-spin channel and consequent large disorder scattering in that channel [see Fig. 4(a)]. As a result, the spin majority channel possesses a well-defined clean Fermi surface while the spin minority channel is “blocked” by the strong disorder scattering induced Fermi surface smearing. Because the resistivity in the majority spin channel is low, it acts as a short circuit with the resulting total residual resistivity also low. On the other hand, when alloying the late 3d elements with early 3d elements with lower electron filling—Cr and/or Mn—both majority-spin bands and minority-spin bands are shifted to high energies to realize charge neutrality. This leads to large δ in both spin channels (large disorder scattering) and much higher residual resistivity. In addition, the moments on Cr or Mn, if they exist and can couple either ferromagnetically or antiferromagnetically with the other spins, further modify δ but still result in large δ.

While accounting for the onsite disorder within the KKR-CPA calculations successfully separates the high-ρ₀ alloys and low-ρ₀ alloys, the theory significantly underestimates the residual resistivity in some alloys, particularly for NiCoCr, NiFeCoCrMn. To investigate possible causes of this underestimation, Mu et al.^{Reference Mu, Samolyuk, Wimmer, Troparevsky, Khan, Mankovsky, Ebert and Stocks42} considered the following additional electron scattering mechanisms.

A. Mean field theory

Even though the development of single-site mean field CPA successfully solves the electronic problem for disordered alloys, in a sense that the configurational average in thermodynamic limit is treated properly, its development and application to phonons has proved much more challenging. As a result, CPA techniques for calculating the effects of disorder on phonons (dispersion and lifetimes) are still much less developed and tested than for electrons—particularly in the availability of community codes that can treat general CSAs and HEAs. The reason lies in the fact that the phonon problem in disordered alloys involves both so-called diagonal disorder and off-diagonal disorder, and they are inseparable. The diagonal disorder is associated with the randomness of the onsite electronic energy (electrons) and the randomness of atomic mass (phonons), while the off-diagonal disorder, neglected by the original formulation of CPA, is associated with the differences in site-to-site hopping integrals (electrons) and force constants (phonons). It is clear that the off-diagonal disorder also plays a fundamental role in many nonlocal problems, and the neglect of it is problematic.

However, for electrons, the problem of neglecting the off-diagonal disorder is obviated by the reformulation of CPA with the help of MST as embodied by KKR-CPA and related techniques. In the KKR-CPA theory for electrons, onsite electronic energy is replaced by the single-site scattering t-matrix, which contains all of the information regarding the scattering properties of a single-site potential, while the off-diagonal site-to-site hopping becomes crystal free particle Green’s function which only depends on the crystal structure and thus is irrelevant to disorder. Therefore, single-site CPA describes the electron state well.

For phonons, no such simple reformulation has proved possible, with the result that generalizing CPA to include off-diagonal randomness is a prerequisite to making progress. Any perturbations, such as a single defect, exert large impacts on not only the mass matrix (onsite disorder) but also the force constant (off-site disorder). Both forms of disorders couple with each other, and the phonon—as the solution of the dynamical matrix—is a result of both onsite and offsite disorders. The original CPA theory for phonons in disordered alloys was developed by Taylor^{Reference Taylor106} and has been applied to systems such as Ge–Si alloy,^{Reference Wakabayashi107} Ni–Pt alloy,^{Reference Tsunoda, Kunitomi, Wakabayashi, Nicklow and Smith108} Al–Cu alloy,^{Reference Kaplan and Mostoller109} and Ni–Pd alloy.^{Reference Kamitakahara and Taylor110} As expected, no satisfactory results are obtained due to the neglect of off-diagonal randomness in the force constants. To include the off-diagonal disorder, successive extensions of CPA have been formulated,^{Reference Kaplan and Mostoller111–Reference Gonis and Garland114} but most of them proved to be problematic due to the lack of self-consistency (in the single-site effective medium sense), loss of the zero energy and wave-vector Goldstone mode, or incorrect analytic properties. It was not until 2002 when Ghosh et al.^{Reference Ghosh, Leath and Cohen115} developed the successful cluster generalization of CPA, called itinerant CPA (ICPA), includes both diagonal disorder and off-diagonal disorder correctly within an analytic and self-consistent theory.

ICPA has been applied to the Ni–Pd alloy with both mass disorder and force constant disorder.^{Reference Ghosh, Leath and Cohen115} Consistent with the experimental observations, the split bands for both LA and TA mode along (00ζ) direction and (0ζζ) direction are obtained from ICPA, providing a positive test for the theory. The phonon properties in other alloys, such as NiPt,^{Reference Ghosh, Leath and Cohen115,Reference Dutta and Ghosh116} FePd,^{Reference Alam, Ghosh and Mookerjee117,Reference Dutta and Ghosh118} and Cu₃Au,^{Reference Dutta, Bisht and Ghosh119} were also been studied by using the ICPA method. However, ICPA is only applicable for binary alloys and it only distinguishes the averaged force constant between different pairs (A–A, A–B, B–B), say, Φ_AA, Φ_BB, and Φ_AB, while neglecting the fluctuations within each type of pairwise interaction. This approximation needs to be further validated from the benchmark with both phonon dispersion and linewidth, measured from the inelastic neutron/X-ray scattering.

Besides ICPA, Alam and Mookerjee^{Reference Alam and Mookerjee120} reported a related but less used augmented space recursion (ASR) approach and further applied the ASR approach to study materials with force-constant dominant disorder, fcc Ni₈₈Cr₁₂ and BCC Fe₅₃Cr₄₇.^{Reference Alam, Ghosh and Mookerjee117,Reference Alam and Mookerjee120,Reference Alam, Chouhan and Mookerjee121} In particular, they explored the short-range order effect on the lattice dynamics and vibrational entropy of the BCC alloy Fe_1−xCr_x and found that the effect of SRO is not significant for the phonon entropy of mixing.^{Reference Alam, Chouhan and Mookerjee121}

B. Supercell phonon unfolding method

Similar to the electronic structure problem in disordered alloys, the phonon bands can also be solved by the unfolding method based on SQSs. Recently, a nonanalytic SC phonon unfolding (SPU) method^{Reference Ku, Berlijn and Lee25,Reference Boykin and Klimeck122–Reference Ikeda, Carreras, Seko, Togo and Tanaka125} was developed to map phonon bands of the SC BZ to the BZ of the primitive cell to directly compare with the measured phonon spectra. It has proven to be successful in semiconductor alloys,^{Reference Boykin, Kharche, Klimeck and Korkusinski123} dilute doped oxide,^{Reference Delaire, Al-Qasir, May, Li, Sales, Niedziela, Ma, Matsuda, Abernathy and Berlijn126} and concentrated metallic alloys, such as Au–Cu^{Reference Ikeda, Carreras, Seko, Togo and Tanaka125} alloys, Fe–Pd alloys,^{Reference Ikeda, Körmann, Dutta, Carreras, Seko, Neugebauer and Tanaka127} and bcc HEAs.^{Reference Körmann, Ikeda, Grabowski and Sluiter128} According to the formalism developed by Ikeda and coworkers,^{Reference Ikeda, Carreras, Seko, Togo and Tanaka125} the phonon spectral function A(k, ω) at the wave vector k and frequency ω is given as: $A\left( {k,\omega } \right) = \sum\limits_J {{{\left| {\left[ {\hat{p}\left( k \right)v\left( {K,J} \right)} \right]} \right|}^2}} \delta \left[ {\omega - \omega \left( {K,J} \right)} \right]$, where $\hat{p}\left( k \right)$ is the projection operator for wave vector k in the primitive BZ, v(K, J) is the eigenvector of the dynamical matrix in the SC cell at SC’s BZ for phonon mode J, and ω(K, J) is the corresponding eigenvalue. This SPU method, utilizing large SCs, naturally captures the effects of local structural relaxation, and couples straightforwardly with the first-principles calculation. On the negative side, it relies on the ansatz that configurational averaging is approximated by manually averaging the observables (spectral functions) over many special quasi-random supercells (SQS).^{Reference Zunger, Wei, Ferreira and Bernard18}

Combining the inelastic neutron scattering measurement and ab initio based SPU method, Delaire et al.^{Reference Delaire, Al-Qasir, May, Li, Sales, Niedziela, Ma, Matsuda, Abernathy and Berlijn126} demonstrated that the introduction of dilute heavy substitutional cation defects in FeSi leads to resonant impurity phonon bands in low-frequency regions. The impurity band hybridizes with the acoustic bands of the host and effectively enhances the phonon scattering while suppresses the lattice thermal conductivity of the thermoelectric.

Applying the SPU method to concentrated nonmagnetic Cu₃Au alloys with massive mass disorder and force constant disorder, the phonon spectra function can be obtained, which captures the experimental data quite well—as illustrated in Fig. 10. The SPU theory finds split bands for one transverse acoustic modes along (0ζζ) and further, species decomposition of the effective spectral function using the SPU method, reveals that this is due to nonoverlapping spectral weights between different types of pairwise interactions, such as Cu–Cu pair, Cu–Au pair, or Au–Au pair. The reason can be boiled down to the disparate of bond stiffness between different pairwise interactions and the mass disorder.

FIG. 10.

Unfolded phonon spectra of CuAu in the fcc lattice, based on 32-atom and 108-atom SC cells. Adapted/Reproduced with permission from Ref. Reference Ikeda, Carreras, Seko, Togo and Tanaka125 (American Physical Society, 2017).

More recently, SPU has been applied to a series of multicomponent bcc CSAs, from binary alloys to 5-component alloys, with special focus on the phonon broadening as a function of chemical complexity.^{Reference Körmann, Ikeda, Grabowski and Sluiter128} As expected, mass disorder and force constant disorder are both important for the phonon properties (phonon dispersion relationship and phonon linewidth). Neglecting both of them—which is equivalent to VCA, or merely including one type of disorder while averaging the other—gives erroneous phonon properties. By resolving the species-dependent contributions to the phonon spectra over the series of bcc alloys, Körmann and coworkers^{Reference Körmann, Ikeda, Grabowski and Sluiter128} found that the top of the phonon brands are dominated by the lighter elements and all the species resolved contributions are independent of actual alloys, with different elemental concentrations. They also found^{Reference Körmann, Ikeda, Grabowski and Sluiter128} that the change of vibration entropy associated with increased chemical complexity is material-specific and could be as large as the variation of the configurational entropy.

C. CLDM

The lattice vibrations or phonons in disordered alloys can also be obtained by parametrizing both the forces and the force constants. The atomic relaxation plays an important role in this effective CE of the forces/force constants and the thermodynamic properties of alloys. Distinguished by the range of interaction, the atomic relaxation can be separated into the long-ranged strain effect which arises from the acoustic sum rule (translational invariance) and the short-ranged local relaxation with respect to the homogeneous deformation of the ideal lattice. Therefore, the displacement of atom i can be written as ${{\bi{u}}_{\bi{i}}} = {u_{{\rm{\alpha \beta }}}}\left( {R_i^{\rm{\alpha }} - R_0^{\rm{\beta }}} \right) + {w_i}$, where u _αβ is the strain tensor and $R_i^{\rm{\alpha }}$ is the coordination of atom i, as compared with the origin coordinate $R_0^{\rm{\beta }}$. w _i is the internal part of atomic displacement, as opposed to the homogeneous deformation of the ideal lattice.

An approach to construct the effective CE for the forces in alloys is the so-called Kanzaki–Krivoglaz–Khachaturyan model (KKKM), which assumes configuration-independent force constants and Kanzaki forces, and the effective interaction therein is pairwise.^{Reference Kanzaki129,Reference Matsubara130} However, both assumptions are oversimplified and this model is generalized to the configuration-dependent lattice deformation model (CLDM) which captures the relaxation effect as well as many-body interactions.^{Reference Zhuravlev, An and Belashchenko131} The formation enthalpy is explicitly separated to the chemical part (H _chem) with all the atoms clamped on the ideal lattice site and the relaxation part (H _rel). Both parts depend on the local configurations (Σ). The latter part can be expressed as following, in the harmonic approximation:

$${H_{{\rm{rel}}}} = - \sum\limits_i {{u_i}{F_i}\left( {\rm{\Sigma }} \right)} + {1 \over 2}\sum\limits_{i,j} {{u_i}{A_{ij}}\left( {\rm{\Sigma }} \right){u_j}} \quad ,$$

where F _i(Σ) is the Kanzaki force and A _ij(Σ) is the configuration-dependent force constant. Minimizing the relaxation energy gives the expression of the forces and stress tensor as a function of configuration-dependent Kanzaki force, configuration-dependent force constant, atomic displacement and strain. In equilibrium, the atomic force and stress tensor vanish. The Kanzaki force is then parametrized according to the configurations and many-body interactions are introduced explicitly. In addition, the force constants are obtained from the linear response of the force to the atomic displacement within the harmonic approximation. All the parameters can be fitted based on the data set of calculated atomic forces within many different cells from first-principles. The relaxation energies, phonon spectra, concentration waves, and thermal dynamic phase transitions can then be obtained. The accuracy of the model can be systematically increased by including fitting parameters for interactions in longer range and more many-body interactions.

This approach was first applied to Cu–Au alloys and Fe–Pt alloys at three concentrations of 25, 50, and 75%. This approach shows that in all cases, the configuration-dependent force is comparable to the conventional Kanzaki force which is configuration independent. Figure 11 illustrates the fitting of the configuration-dependent force constants (blue dots), compared with conventional Kanzaki force fitting (open square). It is seen that the neglect of the configuration-dependent force is detrimental to the fitting for all data sets.

FIG. 11.

Fitting of the force constants referenced from the ideal fcc lattice for (a) Cu–Au and (b) Fe–Pt alloys. The forces predicted by the fit are plotted against the calculated forces. Each Cartesian component for each site is represented by one data point. Adapted/Reproduced with permission from Ref. Reference Zhuravlev, An and Belashchenko131 (American Physical Society, 2014).

With all the configuration-dependent force constants, it is feasible to evaluate the phonon spectra of disordered alloys. Zhuravlev et al.^{Reference Zhuravlev, An and Belashchenko131} calculated the phonon DOS by constructing the force constant for a 2048-atom SC based on the fitted force constants. The calculated DOS was compared with experiments and good overall agreement was reached. Ideally, the parametrized force constant model can be combined with the phonon unfolding method to get the phonon spectra in an underlying primitive cell.

D. Lattice dynamics at finite temperatures

The mean field phonon theory (ICPA), the SPU, and also the CLDM are theories at zero temperature, in which only harmonic constants are considered. In reality, even at zero temperature, anharmonic effects can be important particularly in alloys with large atomic-size-mismatches.^{Reference Zhuravlev, An and Belashchenko131} It is also interesting to note that phonon is not a good long-lived quasi-particle even at zero temperature, due to the local symmetry breaking by disorder. Therefore, the scattering of phonons in disorder environments cause the finite phonon lifetime, along with intrinsic phonon broadening.

For practical applications of the phonon theory, it is necessary to describe the lattice dynamics at finite temperatures when different interactions between various quasi-particles become important. For example, (i) phonon–phonon interaction leads to thermal expansion and renormalized phonon frequencies; (ii) electron–phonon interaction leads to anomalous phonon behavior such as Kohn anomalies or phonon drag in the electronic structure; (iii) magnon–phonon interaction is responsible for structure phase transition such as δ–γ phase transition in the paramagnetic Fe. How interactions affect the phonon properties in the underlying disorder environment is an interesting but difficult problem. All the attempts reviewed here to explore those interactions in disordered alloys are limited to the SC modeling of disorder.

Restricting the discussion to nonmagnetic metals and alloys, we first consider the temperature-dependent phonon spectra due to the anharmonic effect alone. Two major outcomes arise from the anharmonic effect: (i) the thermal expansion of the lattice; (ii) renormalization of the phonons. A straightforward way to include the anharmonic effect is ab initio molecular dynamics since MD is designed to include the temperature effect implicitly. However, ab initio MD suffers from the fact that obtaining an accurate free energy is extremely computational expensive, particularly for large cells. Here, we only review two effective first-principle approaches as to include the anharmonic effect in nonmagnetic metals and alloys.

One ab initio method to account for the thermal expansion effect is the quasi-harmonic approximation (QHA). The QHA requires both phonon calculations and total energy calculations at various constant volumes (V) and pressures (p). By minimizing the Gibbs free energy (see below) at a given temperature with respect to the volume, one can find the equilibrium volume at this particular temperature. Therefore, the evolution of the phonon dispersion relationship as a function of temperature is obtained by linking the temperature to the volume.

$$G\left( {T,p} \right) = {\rm{min}}_V \left[ {U\left( V \right) + F_{{\rm{ph}}} \left( {T;V} \right) + pV} \right]\quad ,$$

where F _ph(T; V) is the Helmholtz free energy defined as

$$F\left( T \right) = \int_0^\infty {{\rm{d}}\omega g\left( \omega \right)} \left[ {{{\hbar \omega } \over 2} + {k_{\rm{B}}}T\ln \left( {1 - {{\rm{e}}^{{{ - \hbar \omega } / {{k_{\rm{B}}}T}}}}} \right)} \right]\quad ,$$

where g(ω) is the phonon density of state at frequency ω. QHA has been applied to Ni–Al alloys to explore the structural phase transition and dynamic properties^{Reference Rubini and Ballone132} and is now widely used.

The renormalization of phonon frequencies by the anharmonic interaction is very common in metals. For example, nonmagnetic Ti, Zr, Hf in the bcc lattice are dynamically unstable within harmonic approximation, but the inclusion of anharmonic effects at high temperatures removes the soft modes and stabilizes the lattice. The self-consistent ab initio lattice dynamics (SCAILD) approach^{Reference Souvatzis, Eriksson, Katsnelson and Rudin133} is designed to efficiently evaluate the high-temperature phonons by using the SC method. The key idea is to iteratively solve the thermal atomic displacements and Hellman–Feynman forces from DFT, and then renormalize the phonon frequency accordingly until self-consistency is reached. Applying this method to study high-temperature phonons in Ti, Zr, and Hf, Souvatzis et al.^{Reference Souvatzis, Eriksson, Katsnelson and Rudin133} found stabilized phonon dispersions, consistent with the experimental data at the corresponding temperatures. However, this method provides no access to the phonon broadening that is the key feature of the phonon–phonon scattering. Moreover, apart from QHA, the SCAILD freezes the volume of the cells, and thereby no thermal expansion effect is included.

For magnetic metals and alloys, finite temperature induces the transverse spin fluctuations that perturb the spin orientation and longitudinal spin fluctuations that change the magnitude of the local moments. These effects can have huge impact on the phonon properties. Perhaps the most well-known example is the δ–γ phase transition of Fe in the paramagnetic state, which has been extensively studied in the past, using the dynamical mean field theory,^{Reference Leonov, Poteryaev, Anisimov and Vollhardt134} spin-space averaging (SSA) procedure,^{Reference Körmann, Dick, Grabowski, Hickel and Neugebauer135} and ab initio MD combined with DLM,^{Reference Steneteg, Alling and Abrikosov136} etc. Here we mainly focus on the SSA procedure and ab initio MD.

The SSA procedure, developed by Förmann et al.^{Reference Körmann, Dick, Grabowski, Hickel and Neugebauer135,Reference Körmann, Grabowski, Dutta, Hickel, Mauger, Fultz and Neugebauer137} for paramagnetic alloys, relies on a series of magnetic SQSs, and within each of them, the spin orientations are randomly and equally populated—as in the DLM model. The force constants of different magnetic SQSs are then averaged as the representative of the force constant in the magnetically disordered paramagnetic state. In this procedure, the averaging of the interatomic forces is performed in the adiabatic limit, in a sense that the spin dynamics is much faster that the atomic dynamics. Therefore, the interatomic force for one magnetic configuration is not representable and the averaging over various magnetic configurations is required. Employing this procedure, not only the phonons of fcc/bcc Fe in the paramagnetic state can be reproduced but also the phonon dispersion at arbitrary temperatures can also be pursued. The magnetic configurations at finite temperatures can also be modeled using Monte Carlo simulation of the Heisenberg model, which is capable to capture the magnetic SRO beyond T _c. The thermal expansion effect was also investigated simply by taking the experimental lattice parameters at different temperatures. In addition, Förmann et al.^{Reference Körmann, Grabowski, Dutta, Hickel, Mauger, Fultz and Neugebauer137} found that the phonon–magnon interaction for some phonon branches is an order of magnitude larger than the phonon shift due to the thermal expansion. This provides convincing evidence of the importance of the phonon–magnon interaction in magnetic alloys.

This SSA combined with magnetic SQS and SPU has been applied to the magnetic alloys Fe_0.75Pd_0.25, Fe_0.75 Pt_0.25 at finite temperatures.^{Reference Ikeda, Körmann, Dutta, Carreras, Seko, Neugebauer and Tanaka127} Since both alloys are Invar alloys, whose lattice parameters are not sensitive to the temperature, the thermal expansion effect is safely excluded from the consideration. Experiments found an unusual phonon hardening of a transverse acoustic mode in Fe–Pd and Fe–Pt at elevated temperatures. The calculated finite temperature phonon spectra in Fe_0.75Pd_0.25 and Fe_0.75 Pt_0.25 are shown in Fig. 12. The reason for the phonon hardening is traced down to the effect of spin disorder on the force constants. Magnetic SRO can also be included in the SSA technique, and it was found that the presence of magnetic SRO even above T _c significantly improves the agreement between theory and experiment in Fe_0.75Pd_0.25, Fe_0.75 Pt_0.25, and Fe.

FIG. 12.

Temperature-dependent phonon spectra for the lowest transverse phonon mode in (110) direction in (a) Fe_0.75Pd_0.25 and (b) Fe_0.75Pt_0.25. Blue and red corresponds to the phonon spectra in a ferromagnetic phase and paramagnetic phase, respectively. Green gives the phonon spectra at finite temperatures. Experimental data are also labeled in open circle. Adapted/Reproduced with permission from Ref. Reference Ikeda, Körmann, Dutta, Carreras, Seko, Neugebauer and Tanaka127 (Springer Nature, 2018, http://creativecommons.org/licenses/by/4.0/).

An ab initio MD combined with the DLM picture has been used to study the phonons in the paramagnetic phase.^{Reference Steneteg, Alling and Abrikosov136} This method started with collinear spin polarized MD based on a DLM state in a SC, i.e., with randomly oriented spins and vanishing total net moment. The time step for MD is denoted t _MD. After many ionic MD steps, the directions of all local moments are then completely rearranged at a time step t _SF. Thus, it is required that t _SF ≫ t _MD in those DLM-MD steps. After several DLM-MD steps, the value of observables of interest such as potential energy can be well-converged. This method has been applied to CrN to investigate the equation of state and bulk modulus in both antiferromagnetic phase and paramagnetic phase. The effect of electron correlation on the mixing thermodynamics in Cr_1−xAl_xN solid solutions is investigated using this DLM-MD technique.^{Reference Alling, Hultberg, Hultman and Abrikosov138} In addition, combined with the temperature-dependent effective potential method developed by Hellman et al.,^{Reference Hellman and Abrikosov139} DLM-MD technique could simulate the phonon band^{Reference Hellman, Abrikosov and Simak140} and vibration free energy of magnetic materials accurately,^{Reference Shulumba, Alling, Hellman, Mozafari, Steneteg, Odén and Abrikosov141} providing a reliable possibility to model the phase transition in magnetic materials.

A. Transport in bcc HEAs with larger atomic displacements

In Cantor-Wu alloys, the electron scattering due to displacement fluctuations is minute due to typical small atomic displacements since all alloying elements are of similar sizes, embedded in a close-packed fcc structure. It is interesting to study the electronic transport in the series of bcc HEAs and to explore the effect of atomic displacements on the electrical transport. The alloying elements from the conventional bcc HEAs (e.g., VNbMoTaW and VNbTiTaW) are 3d-, 4d-, 5d-, transition metal elements, from the left side of the periodic table. This means (i) the similar number of valence electrons of the alloying elements weakens the disorder scattering of the electrons; (ii) the displacement fluctuations are expected to be large due to the atomic size-mismatch effect. In addition, all of the bcc HEAs are nonmagnetic and thereby, the complex magnetic scattering is absent. Figure 13 compares the BSF of bcc VNbMoTaW and that of fcc NiFeCoCrMn HEAs (in the spin minority channel). As shown, BSF of VNbMoTaW exhibits weaker band smearing at the Fermi level and the residual resistivity of VNbMoTaW should be lower. Moreover, because atomic displacements are large in bcc alloys, their effects on transport properties of low-resistivity alloys are pronounced. Therefore, a series of bcc HEAs should display a completely different feature of transport properties compared with the fcc HEAs: the maximally chemical disorder-induced electron scattering is much weaker than that in fcc HEAs due to the nature of alloying elements while displacement fluctuations make an important contribution to residual resistivity.

FIG. 13.

(a) BSF for bcc VNbMoTaW and (b) majority-spin BSF fcc NiFeCoCrMn, given in arbitrary unit. BSF is calculated using SPR-KKR-CPA.^{Reference Ebert, Ködderitzsch and Minár40,Reference Ebert41} Adapted/Reproduced with permission from Ref. Reference Zhang, Stocks, Jin, Lu, Bei, Sales, Wang, Béland, Stoller, Samolyuk, Caro, Caro and Weber15 (Springer Nature, 2015, http://creativecommons.org/licenses/by/4.0/).

B. Effect of disorder on thermal transport

The current treatment of lattice vibrations in disordered alloys assumes well-defined phonon quasi-particles and also assumes that the presence of disorder only amounts to a perturbation to the phonon Hamiltonian. Within this approximation, the lattice-mediated thermal transport can be solved by the phonon gas model or BTE. However, this approximation is questionable since in many alloys that present force disorder and mass disorder, there is no well-defined dispersion relationship ω(k) at T = 0 K and phonons as well due to the disorder—particularly in CSAs. Allen and coworkers^{Reference Allen, Feldman, Fabian and Wooten142} first distinguished the lattice vibration in random alloys by the diffusons, propagons, and localons. Although both diffusons and propagons are delocalized, the difference between them lies in that propagons are wave-like quasi-particles with definite wave vectors while diffusons do not exhibit periodicity of waves. Moreover, localons are local vibrations around some atoms and decay quickly away from the centered atoms. Recently, Seyf et al.^{Reference Seyf, Yates, Bougher, Graham, Cola, Detchprohm, Ji, Kim, Dupuis, Lv and Henry143} studied thermal transport using the KG formula in disordered In_xGe_1−x As, and they found excellent agreement of the temperature-dependent thermal conductivity with the experiment. This method is essentially different from the phonon-gas model and can be used to sort out different contributions to the thermal conductivity from the diffusons, localons, and propagons. This study reveals that doping, even in dilute limit, effectively reduces the propagons character. In low-temperature regions, the thermal transport is however still governed by the propagons. With modest amounts of defects, diffusons become more important and dominate the thermal transport in high-temperature regions. This study invokes the rethinking of phonons in disordered systems and will encourage more fundamental studies on validation of phonons and the lattice-mediated thermal transport in disordered alloys.