[1]#1 [1]*~#1 *g_loc *T

# Dynamical mean-field theory of the Anderson–Hubbard model with local and
non-local disorder in tensor formulation

A. Weh *  M<EMAIL_ADDRESS>Theoretical Physics II,
Institute of Physics, University of Augsburg, 86135 Augsburg, Germany  Y.
Zhang *  M Kavli Institute of Theoretical Sciences, University of Chinese
Academy of Sciences, Beijing, 100190, China  A. Östlin *  M Theoretical
Physics III, Center for Electronic Correlations and Magnetism, Institute of
Physics, University of Augsburg, 86135 Augsburg, Germany  Augsburg Center for
Innovative Technologies, University of Augsburg, 86135 Augsburg, Germany  H.
Terletska  Department of Physics and Astronomy, Middle Tennessee State
University, Murfreesboro, TN 37132, USA  D. Bauernfeind  Center for
Computational Quantum Physics, Flatiron Institute, 162 5th Avenue, New York,
NY 10010, USA  K.-M. Tam  Department of Physics & Astronomy, Louisiana State
University, Baton Rouge, LA 70803, USA  Center for Computation & Technology,
Louisiana State University, Baton Rouge, LA 70803, USA  H. G. Evertz *  M
Institute of Theoretical and Computational Physics, Graz University of
Technology, 8010 Graz, Austria  K. Byczuk *  M Institute of Theoretical
Physics, Faculty of Physics, University of Warsaw, ul. Pasteura 5, 02-093
Warszawa, Poland  D. Vollhardt  Theoretical Physics III, Center for Electronic
Correlations and Magnetism, Institute of Physics, University of Augsburg,
86135 Augsburg, Germany  L. Chioncel *  M Theoretical Physics III, Center for
Electronic Correlations and Magnetism, Institute of Physics, University of
Augsburg, 86135 Augsburg, Germany  Augsburg Center for Innovative
Technologies, University of Augsburg, 86135 Augsburg, Germany

(August 28, 2024)

###### Abstract

To explore correlated electrons in the presence of local and non-local
disorder, the Blackman–Esterling–Berk method for averaging over off-diagonal
disorder is implemented into dynamical mean-field theory using tensor
notation. The impurity model combining disorder and correlations is solved
using the recently developed fork tensor-product state solver, which allows
one to calculate the single particle spectral functions on the real-frequency
axis. In the absence of off-diagonal hopping, we establish exact bounds of the
spectral function of the non-interacting Bethe lattice with coordination
number $Z$. In the presence of interaction, the Mott insulating paramagnetic
phase of the one-band Hubbard model is computed at zero temperature in alloys
with site- and off-diagonal disorder. When the Hubbard $U$ parameter is
increased transitions from an alloy band-insulator through a correlated metal
into a Mott insulating phase are found to take place.

##  I Introduction

The electronic structure and transport properties of real materials are
strongly influenced by the Coulomb interaction between the electrons and the
presence of disorder Lee and Ramakrishnan (1985); Mott (1990); Belitz and
Kirkpatrick (1994). In particular, both electronic correlations and randomness
are driving forces behind a transition from a metallic to an insulating state
due to the localization of electrons (“metal-insulator transition” (MIT)).
While the Mott–Hubbard MIT is caused by the repulsive interaction between the
electrons Mott (1949, 1990); Imada  _et al._ (1998), the Anderson MIT is a
result of coherent backscattering of non-interacting electrons from randomly
distributed impurities Anderson (1958); Lee and Ramakrishnan (1985); Evers and
Mirlin (2008); Abrahams (2010). The interplay between interactions and static
disorder gives rise to many unusual and often unexpected phenomena Lee and
Ramakrishnan (1985); Belitz and Kirkpatrick (1994); Abrahams  _et al._ (2001);
Abrahams (2010); Dobrosavljevic  _et al._ (2012); Nandkishore and Huse (2015);
Abanin  _et al._ (2019). The simplest model of disordered interacting
electrons is the Anderson–Hubbard model, obtained by supplementing a single-
band Hubbard model with local and/or non-local disorder. If the disorder acts
only locally, i.e., via random local potentials (“diagonal disorder”), this
model is able to describe substitutionally disordered binary alloys. However,
in general disorder also affects the amplitudes for hopping between two sites
— especially when the bandwidths of the host and dopant are very different —
leading to additional “off-diagonal disorder”. In analytic calculations local
disorder is easier to treat and was studied extensively Ziman (1979). In
particular, the coherent potential approximation (CPA) Soven (1967); Taylor
(1967); Velický  _et al._ (1968); Elliott  _et al._ (1974); Economou (2006);
Gonis (1992) provides the best single-site approximation for non-interacting
systems with local disorder. For that reason the simultaneous investigation of
diagonal disorder within the CPA and of interacting electrons with local
(Hubbard) interaction within the dynamical mean-field theory (DMFT) Metzner
and Vollhardt (1989); Georges and Krauth (1992); Georges  _et al._ (1996);
Kotliar and Vollhardt (2004) fit together particularly well, since both DMFT
Metzner and Vollhardt (1989); Georges  _et al._ (1996) and CPA Vlaming and
Vollhardt (1992); Janiš and Vollhardt (1992) become exact in the limit of
infinite spatial dimensions or lattice coordination number.

The treatment of the Anderson–Hubbard model with off-diagonal disorder
received somewhat less attention Milovanović  _et al._ (1989); Dobrosavljević
and Kotliar (1993, 1994); Ulmke and Scalettar (1997); Denteneer  _et al._
(2001); Miranda and Dobrosavljevic (2012). In particular, Dobrosavljević and
Kotliar Dobrosavljević and Kotliar (1993, 1994) investigated this model within
DMFT by employing a functional integral representation for quantum averages
and the replica method for disorder averaging. Thereby they were able to study
Hubbard models with arbitrary disorder on Bethe lattices, as well as models on
an arbitrary lattice with a special distribution of the off-diagonal disorder.
In this way they studied the formation of local moments and the Mott
transition in disordered systems.

In 1971 Blackman, Esterling, and Berk (BEB) Blackman  _et al._ (1971);
Esterling (1975) showed that off-diagonal disorder can, in principle, be
incorporated into the CPA framework, such that both diagonal and off-diagonal
disorder are tractable within a single–site approximation. In the absence of
electronic interactions the BEB formalism was incorporated into the dynamical
cluster approximation within the typical medium cluster theory Terletska  _et
al._ (2014) and applied to multi-band systems Zhang  _et al._ (2016).

In our paper we extend these investigations of disordered systems by including
a local (“Hubbard”) interaction between the electrons. To this end we
investigate the Anderson–Hubbard model with diagonal and off-diagonal disorder
within the CPA in the BEB formulation Blackman  _et al._ (1971); Esterling
(1975) in a tensor formulation Koepernik  _et al._ (1997, 1998), while
correlation effects are treated within the DMFT Georges  _et al._ (1996). We
compute the spectral function and discuss the occurrence of the MIT in the
presence of diagonal and off-diagonal disorder. For these purposes an accurate
zero-temperature many-body solver on the real energies is used, the so-called
fork tensor-product state solver Bauernfeind  _et al._ (2017). We observe
successive alloy-insulator to metal and metal to Mott-insulator transitions
with increasing values of the Hubbard $U$ parameter. A similar transition
scenario was previously discussed for models with diagonal disorder solved
within DMFT, but using finite temperature solvers such as the Hirsch-Fye
algorithm Ulmke  _et al._ (1995) or the perturbative non-crossing
approximation Lombardo  _et al._ (2006). Contrary to the diagonal disorder
model in which the CPA solution provides a common bandwidth for all alloy
components, the presence of off-diagonal elements causes the formation of
different effective bandwidths for alloy components.

The paper is organized as follows. In Section II the BEB theory of the multi-
component Anderson–Hubbard model is formulated, and the computational scheme
is discussed. In Section III numerical results for a single-band Bethe lattice
are presented. Comparison with earlier results on DMFT+CPA allow us to
identify effects specifically due to off-diagonal disorder. Finally,
conclusions and a summarizing discussion are presented in Section IV.

##  II Blackman–Esterling–Berk theory for the Anderson–Hubbard model

In the simplest case an alloy consists of two types of atoms, $\mathsf{A}$ and
$\mathsf{B}$, with diagonal substitutional disorder, such that only the site-
diagonal elements of the Hamiltonian vary stochastically according to the
atomic species occupying the given site Soven (1967); Taylor (1967); Velický
_et al._ (1968). The physical quantities of interest, for instance the
spectral function, are those averaged over the possible disorder realizations.
Therefore, the idea of the CPA is to replace the ensemble with random
configurations by a periodic system with “average” atoms, whose properties are
determined self-consistently. The CPA finds a natural description in the
language of scattering theory. Assuming the origin to be occupied by atoms of
type $\mathsf{A}$ or $\mathsf{B}$ and all other sites by average atoms, the
scattering by the atom at the origin is easily computed. The self-consistency
condition of the CPA, which expresses that the scattering at the origin
vanishes on average, then allows one to compute the coherent Green’s function
for the average atom.

The BEB formalism Blackman  _et al._ (1971); Esterling (1975) is a
generalization of the CPA such that it becomes applicable also to off-diagonal
disorder. Similar to CPA it was formulated for a tight-binding model in which
the hopping matrix elements depend on the species of atoms occupying the two
sites connected by the hopping. For the above binary alloy example, the BEB
hopping matrix elements are
$t_{\mathsf{A}\mathsf{A}},t_{\mathsf{A}\mathsf{B}}$ and
$t_{\mathsf{B}\mathsf{B}}$. When the hopping matrix elements are equal the BEB
formalism reduces to the CPA Soven (1967); Taylor (1967), and for the binary
alloy case the scalar CPA equation becomes a $2\times 2$ matrix equation. An
in-depth analysis of the BEB method was performed in a tight-binding formalism
by Gonis and Garland Gonis and Garland (1977) using locators, propagators, and
a variational technique proving the analyticity of the BEB-CPA Green’s
function. A realistic multiband formulation of the BEB-CPA was introduced more
than three decades ago by Papaconstantopoulos, Gonis, and Laufer
Papaconstantopoulos  _et al._ (1989). Then Koepernik _et al._ Koepernik  _et
al._ (1997) developed the BEB-CPA extension within a full potential local-
orbital approach, and more recently a similar implementation was made within a
pseudopotential approach Herbig  _et al._ (2017). Shvaika Shvaika (2003) found
a connection between the Falicov–Kimball model with correlated hopping and the
BEB-CPA by rewriting the Hamiltonian as a $2\times 2$ matrix. For Hamiltonians
with interactions the BEB-CPA was employed by Burdin and Fulde Burdin and
Fulde (2007) to study the interplay between the Kondo effect and disorder. In
an attempt to address localization in strongly disordered electronic systems,
the typical medium theory Dobrosavljević and Kotliar (1993) was combined with
the dynamical cluster approximation including effects induced by off-diagonal
disorder Terletska  _et al._ (2014). However, this approach did not include
electronic interactions.

In the present paper we extend the BEB formalism to _interacting_ electrons
such that it can be applied to the multi-component Anderson–Hubbard model;
localization effects will not be addressed. The model is defined in Section
II.2. The corresponding DMFT equations are solved using the recently developed
fork tensor-product state solver Bauernfeind  _et al._ (2017).

###  II.1 Configurational averages and notation

In the conventional approach to systems with random variables (diagonal and/or
off-diagonal) the Green’s function is first expanded and then an average over
an appropriate set of terms is performed. By contrast, the BEB method treats
both diagonal- and off-diagonal randomness on equal footing by employing an
extended representation, which will be discussed below. The Green’s functions
are then evaluated using conventional expansion techniques. The formalism
introduced by Koepernik _et al._ Koepernik  _et al._ (1997, 1998) is
particularly suitable for the BEB approach. For this reason, we adopt the
notation introduced in Refs. Koepernik  _et al._ , 1997, 1998.

We consider an alloy consisting of $M$ types of atoms (“alloy components”)
denoted by the index $\alpha$. Every lattice site $i$ is uniquely mapped to a
particular component $\alpha$ as expressed by

$i\mapsto\alpha.$ (1)

While this notation corresponds, in principle, to that of Ref. Koepernik  _et
al._ , 1998, we reverse the direction of the arrow to focus on the alloy
components rather than the lattice sites. To address multiple sites we use the
notation
$(i,j\mapsto\alpha,\beta)\coloneqq(i\mapsto\alpha)\land(j\mapsto\beta)$. We
refer to a specific mapping of the $N$ lattice sites to the $M$ components as
a “configuration” ($\operatorname{conf}$) or “disorder realization”, and
denote the set of all possible configurations by
${{\mathbb{C}}}=\\{\operatorname{conf}\\}$. As the specific configuration of a
sample measured in an experiment is unknown, we average over all possible
configurations. Since the concentrations $c^{\alpha}$ of the different
components $\alpha$ are assumed to be known we restrict the average to
configurations with these concentrations, denoted by
${{\mathbb{C}}}|_{\\{c^{\alpha}\\}}$. In the absence of additional information
we further assume that the probability of all physical configurations is the
same:

$P(\operatorname{conf})=1/\absolutevalue{{\mathbb{C}|_{\\{c^{\alpha}\\}}}}\quad\forall\operatorname{conf}\in{\mathbb{C}|_{\\{c^{\alpha}\\}}}.$
(2)

For a random variable $X$, the stochastic average over all physical
configurations is the weighted sum

$\operatorname{\mathbb{E}}(X)=\mkern-12.0mu\sum_{\operatorname{conf}\in{\mathbb{C}|_{\\{c^{\alpha}\\}}}}\mkern-18.0muP(\operatorname{conf})x_{\operatorname{conf}}=\frac{1}{\absolutevalue{{\mathbb{C}|_{\\{c^{\alpha}\\}}}}}\sum_{\operatorname{conf}\in{\mathbb{C}|_{\\{c^{\alpha}\\}}}}\mkern-18.0mux_{\operatorname{conf}}.$
(3)

This situation corresponds to the case of substitutional disorder.

$i,j,\dots$ | site indices | $c^{\alpha}$ | concentration of component $\alpha$
---|---|---|---
$\alpha,\beta,\dots$ | alloy components | $\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha}$ | on-site energy of component $\alpha$
$i\mapsto\alpha$ | mapping of a site to an alloy component | $\boldsymbol{\underline{\mathrm{*}}}{U}^{\alpha}$ | Hubbard parameter of component $\alpha$
$\mathbb{S}^{\alpha}$ | set of sites occupied by component $\alpha$ | $\boldsymbol{\underline{\mathrm{*}}}{t}^{\alpha\beta}(\absolutevalue{\mathbf{r}_{i}-\mathbf{r}_{j}})$ | hopping between component $\alpha$ and $\beta$
${\mathbb{C}|_{\\{c^{\alpha}\\}}}$ | set of configurations (disorder realization) | $\boldsymbol{\underline{\mathrm{*}}}{H}_{ij}^{\alpha\beta}$ | extended Anderson–Hubbard Hamiltonian matrix
$X$ | scalar random variable | ${\boldsymbol{\underline{\mathrm{*}}}{\eta}}^{\alpha}_{ij}$ | indicator tensor, represented as equivalent matrix
$\operatorname{\mathbb{E}}(X)$ | Expectation value defined according to Eq. 3 | $\boldsymbol{\underline{\mathrm{\chi}}}$ | projector onto specific disorder configuration
|  | $\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\beta}$ | dimensionless hopping parameter
|  |  | describing hopping between components $\alpha$ and $\beta$
Table 1: Notations specific to the extended basis and Hamiltonian parameters
used in BEB method.

###  II.2 Multi-component Anderson–Hubbard model

For a specific configuration the Anderson–Hubbard Hamiltonian reads

$\hat{H}=-\sum_{ij\sigma}t_{ij}\hat{c}^{\dagger}_{i\sigma}\hat{c}^{\vphantom{\dagger}}_{j\sigma}+\sum_{i\sigma}(v_{i}-\mu)\hat{n}_{i\sigma}+\sum_{i}U_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow},$
(4)

with the on-site energy $v_{i}$, the local Hubbard interaction $U_{i}$, and
the amplitude $t_{ij}$ for hopping between sites $i$ and $j$. The hopping
parameters are Hermitian $t_{ij}=t^{*}_{ji}$ and off-diagonal, with
$t_{ii}=0$. The Hamiltonian can be written in the compact matrix form

$\hat{H}=\sum_{\sigma}\boldsymbol{\hat{c}^{\vphantom{\dagger}}}^{\dagger}_{\sigma}\boldsymbol{H}_{\sigma}\boldsymbol{\hat{c}^{\vphantom{\dagger}}}_{\sigma}+\boldsymbol{\hat{n}}^{\intercal}_{\uparrow}\boldsymbol{U}\boldsymbol{\hat{n}}_{\downarrow},$
(5)

where we introduced $N\times 1$ matrices to represent the operators. The rows
of the matrix $\boldsymbol{\hat{c}^{\vphantom{\dagger}}_{\sigma}}$ are the
annihilation operators $\hat{c}^{\vphantom{\dagger}}_{i\sigma}$, and the rows
of $\boldsymbol{\hat{n}}_{\sigma}$ are the number operators
$\hat{n}_{i\sigma}$. The one-particle Hamiltonian matrix reads
${(\boldsymbol{H})}_{ij}=-t_{ij}+\delta_{ij}(v_{i}-\mu)$. Here and in the
following we suppress the spin-index $\sigma$ unless explicitly needed, to
simplify the notation. The local interaction is written as a matrix
${(\boldsymbol{U})}_{ij}=\delta_{ij}U_{i}$.

The magnitude of the hopping parameters $t_{ij}$ depends on the alloy
components located on sites $i$ and $j$, respectively, which are referred to
as “terminal points”. In the following we employ the “terminal-point
approximation” Koepernik  _et al._ (1997, 1998) which assumes that parameters
with terminal points $i,j,k,\dots$ depend only on the components located at
$i,j,k,\dots$ and not on the components surrounding these sites. Thus, for a
specific configuration (disorder realization) every parameter
$v_{i},U_{i},t_{ij}$ takes a value depending on the component occupying the
respective site or sites. In the representation of the BEB we denote these
configuration-specific values by an underline and a superscript indicating the
component. For instance, if site $i$ is occupied by component $\alpha$
($i\mapsto\alpha$) the parameter $v_{i}$ takes the value
$\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha}$. For
$i,j\mapsto\alpha,\beta$, we have
$v_{i}=\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha}$,
$U_{i}=\boldsymbol{\underline{\mathrm{*}}}{U}^{\alpha}$, and
$t_{ij}=\boldsymbol{\underline{\mathrm{*}}}{t}^{\alpha\beta}(\absolutevalue{\mathbf{r}_{i}-\mathbf{r}_{j}})$.
This will now be formalized. Denoting the set of sites $i$ occupied by the
component $\alpha$ by

$\mathbb{S}^{\alpha}\coloneqq\\{i|i\mapsto\alpha\\},$ (6)

the terminal point approximation can be expressed conveniently using the
indicator function

$1_{\mathbb{S}^{\alpha}}(i)\coloneqq\begin{cases}1&\mbox{\quad
if\quad}i\in\mathbb{S}^{\alpha},\\\ 0&\mbox{\quad
if\quad}i\not\in\mathbb{S}^{\alpha}.\end{cases}$ (7)

The identity $\sum_{\alpha}1_{\mathbb{S}^{\alpha}}(i)=1$ holds since every
site must be occupied by exactly one component. Thus, the parameters read

$\displaystyle v_{i}$
$\displaystyle=\sum_{\alpha}1_{\mathbb{S}^{\alpha}}(i)\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha},$
(8) $\displaystyle t_{ij}$
$\displaystyle=\sum_{\alpha\beta}1_{\mathbb{S}^{\alpha}}(i)\boldsymbol{\underline{\mathrm{*}}}{t}^{\alpha\beta}(\absolutevalue{\mathbf{r}_{i}-\mathbf{r}_{j}})1_{\mathbb{S}^{\beta}}(j),$
$\displaystyle H_{ij}$
$\displaystyle=\sum_{\alpha\beta}1_{\mathbb{S}^{\alpha}}(i)\boldsymbol{\underline{\mathrm{*}}}{H}^{\alpha\beta}_{ij}1_{\mathbb{S}^{\beta}}(j),$
$\displaystyle U_{i}$
$\displaystyle=\sum_{\alpha}1_{\mathbb{S}^{\alpha}}(i)\boldsymbol{\underline{\mathrm{*}}}{U}^{\alpha},$

with
$\boldsymbol{\underline{\mathrm{*}}}{H}_{ij}^{\alpha\beta}=\delta_{ij}\delta^{\alpha\beta}(\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha}-\mu)-\boldsymbol{\underline{\mathrm{*}}}{t}^{\alpha\beta}(\absolutevalue{\mathbf{r}_{i}-\mathbf{r}_{j}})$.
We further note, that the conditional expectation value of the parameters
equals the underlined component variables:

$\operatorname{\mathbb{E}}(v_{i}|i\mapsto\alpha)=\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha},\quad\operatorname{\mathbb{E}}(t_{ij}|i,j\mapsto\alpha,\beta)=\boldsymbol{\underline{\mathrm{*}}}{t}^{\alpha\beta},\dots$
(9)

The dependence on the components is therefore shifted from the parameters into
the indicator function $1_{\mathbb{S}^{\alpha}}(i)$. Depending on the
component occupying a site, the indicator function selects the corresponding
parameter from a finite set of choices. We note that for elements diagonal in
lattice sites $i$, one has
$1_{\mathbb{S}^{\alpha}}(i)1_{\mathbb{S}^{\beta}}(i)=1_{\mathbb{S}^{\alpha}}(i)\delta^{\alpha\beta}$,
i.e., they are diagonal in the components. In the following we refer to a
quantity with multiple indices, which include both site and component indices,
as a “tensor”.

We introduce the indicator tensor

${\boldsymbol{\underline{\mathrm{*}}}{\eta}}_{ij}^{\alpha}=1_{\mathbb{S}^{\alpha}}(i)\delta_{ij}=\leavevmode\hbox
to56.04pt{\vbox to19.77pt{\pgfpicture\makeatletter\hbox{\hskip
28.61557pt\lower-10.62743pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}}
{{}\pgfsys@rect{-8.53583pt}{-8.53583pt}{17.07166pt}{17.07166pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-2.83888pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\boldsymbol{\underline{\mathrm{\eta}}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {
{}{}{}}{{{}}{{}}}{}{{}}{}
{{{}{}}{}}{}{{}}{}\pgfsys@moveto{-8.73582pt}{4.30554pt}\pgfsys@lineto{-17.27165pt}{4.30554pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{{}{}}}{{}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-25.28256pt}{2.7986pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{\scriptsize$\alpha$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {
{}{}{}}{{{}}{{}}}{}{{}}{}
{{{}{}}{}}{}{{}}{}\pgfsys@moveto{-8.73582pt}{-4.30554pt}\pgfsys@lineto{-17.27165pt}{-4.30554pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{{}{}}}{{}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-23.21625pt}{-6.61388pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{\scriptsize$i$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
{{}{}{}}{{{}}{{}}}{}{{}}{}
{{{}{}}{}}{}{{}}{}\pgfsys@moveto{8.73582pt}{-4.30554pt}\pgfsys@lineto{17.27165pt}{-4.30554pt}\pgfsys@stroke\pgfsys@invoke{
}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{20.80466pt}{-5.93333pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{\scriptsize$j$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys<EMAIL_ADDRESS>(10)

Graphically we represent this tensor as a box with legs as seen on the right-
hand side of Eq. 10. The order of the tensor is given by the number of its
legs, here three. The upper leg carries the alloy component indices $\alpha$,
and the lower legs correspond to the site indices $i,j$. Within our matrix
notation this tensor is equivalent to a $MN\times N$ matrix. We group the left
indices for sites $i$ and components $\alpha$, or in the graphical notation
the legs above each other. In the following we refer to the $MN$-dimensional
vector space of grouped sites and components as “extended space”. Matrix
products in the extended space sum over the grouped $MN$ elements for
component and site indices; they are equivalent to the tensor contraction of
two legs, one for the component and one for the sites. In this matrix notation
the Hamiltonian reads

$\hat{H}=\sum_{\sigma}\boldsymbol{\hat{c}^{\vphantom{\dagger}}}^{\dagger}_{\sigma}\boldsymbol{\underline{\mathrm{\eta}}}^{\intercal}\boldsymbol{\underline{\mathrm{H}}}_{\sigma}\boldsymbol{\underline{\mathrm{\eta}}}\boldsymbol{\hat{c}^{\vphantom{\dagger}}}_{\sigma}+\boldsymbol{\hat{n}}^{\intercal}_{\uparrow}\boldsymbol{\underline{\mathrm{\eta}}}^{\intercal}\boldsymbol{\underline{\mathrm{U}}}\boldsymbol{\underline{\mathrm{\eta}}}\boldsymbol{\hat{n}}_{\downarrow},$
(11)

where we introduced the local interaction tensor
${(\boldsymbol{\underline{\mathrm{U}}})}^{\alpha\beta}_{ij}=\delta_{ij}\boldsymbol{\underline{\mathrm{*}}}{U}^{\alpha}\delta^{\alpha\beta}$.
In the non-interacting case the Hamiltonian matrix of a specific configuration
is the extended matrix sandwiched by the indicator tensors

$\boldsymbol{H}=\boldsymbol{\underline{\mathrm{\eta}}}^{\intercal}\boldsymbol{\underline{\mathrm{H}}}\boldsymbol{\underline{\mathrm{\eta}}}=\leavevmode\hbox
to98.34pt{\vbox to17.47pt{\pgfpicture\makeatletter\hbox{\hskip
17.47165pt\lower-8.73582pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}}
{{}\pgfsys@rect{-8.53583pt}{-8.53583pt}{17.07166pt}{17.07166pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-4.05556pt}{-3.44165pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\boldsymbol{\underline{\mathrm{\eta}}}^{\intercal}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {
{}{}{}}{{{{}{}{}}}}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{
{}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}}
{{}\pgfsys@rect{23.1622pt}{-8.53583pt}{17.07166pt}{17.07166pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{29.19803pt}{-2.83888pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\boldsymbol{\underline{\mathrm{H}}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {
{}{}{}}{{{{}{}{}}}}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{
{}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}}
{{}\pgfsys@rect{54.86023pt}{-8.53583pt}{17.07166pt}{17.07166pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{60.89606pt}{-2.83888pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\boldsymbol{\underline{\mathrm{\eta}}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {
{}{}{}}{{{}}{{}}}{}{{}}{}
{{{}{}}{}}{}{}\pgfsys@moveto{-8.73582pt}{-4.30554pt}\pgfsys@lineto{-17.27165pt}{-4.30554pt}\pgfsys@stroke\pgfsys@invoke{
} {{}{}{}}{{{}}{{}}}{}{{}}{}{ {}{}{}}{{{}}{{}}}
{}{}{{}{}{}}{{{}}{{}}}{}{{}}{}{ {}{}{}}{{{}}{{}}}
{}{}{}\pgfsys@moveto{8.73582pt}{4.30554pt}\pgfsys@lineto{22.9622pt}{4.30554pt}\pgfsys@moveto{40.43385pt}{4.30554pt}\pgfsys@lineto{54.66023pt}{4.30554pt}\pgfsys@stroke\pgfsys@invoke{
} {{}{}{}}{{{}}{{}}}{}{{}}{}{ {}{}{}}{{{}}{{}}}
{}{}{{}{}{}}{{{}}{{}}}{}{{}}{}{ {}{}{}}{{{}}{{}}}
{}{}{}\pgfsys@moveto{8.73582pt}{-4.30554pt}\pgfsys@lineto{22.9622pt}{-4.30554pt}\pgfsys@moveto{40.43385pt}{-4.30554pt}\pgfsys@lineto{54.66023pt}{-4.30554pt}\pgfsys@stroke\pgfsys@invoke{
} {{}{{}{}}{}}{}{{}}{}
{{{}{}}{}}{}{}\pgfsys@moveto{72.13188pt}{-5.04385pt}\pgfsys@lineto{80.66771pt}{-5.04385pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys<EMAIL_ADDRESS>(12)

For every lattice site in the extended representation each component of the
non-interacting Hamiltonian matrix $\boldsymbol{\underline{\mathrm{H}}}$ is
assigned a corresponding element. In this way the non-interacting Hamiltonian
can be generated by $\boldsymbol{\underline{\mathrm{H}}}$ for every disorder
configuration. The matrix product in the algebraic equation equals the tensor
contractions of the internal legs as illustrated by the right hand side of Eq.
12. Appendix A provides an explicit example for a system of $N=3$ sites and
$M=2$ components. In Table 1 we collect the symbols used in our paper.

We note, that the only configuration dependent parts in Eq. 11 are the
matrices $\boldsymbol{\underline{\mathrm{\eta}}}$ and
$\boldsymbol{\underline{\mathrm{\eta}}}^{\intercal}$; the rest is independent
of the specific disorder realization. In other words: Comparing Eq. 11 with
Eq. 5 the configuration dependence of the former equation is moved from the
Hamiltonian matrix to the local indicator tensors Eq. 10. This is the main
point of the BEB algorithm: One can work with a non-random but extended
Hamiltonian matrix $\boldsymbol{\underline{\mathrm{H}}}$, which contains the
parameters for all possible configurations. A specific configuration can be
selected by applying indicator tensors
$\boldsymbol{\underline{\mathrm{\eta}}}$. What remains to be averaged over are
these local indicator tensors.

####  II.2.1 Alloy component Green’s function

In the absence of interaction, $\boldsymbol{\underline{\mathrm{U}}}=0$, the
model can be solved by the generalized CPA introduced by Blackman, Esterling
and Berk Blackman  _et al._ (1971); Esterling (1975) (BEB). Using the
indicator tensor $\boldsymbol{\underline{\mathrm{\eta}}}$ Eq. 10 we define the
projector:

$\boldsymbol{\underline{\mathrm{\chi}}}=\boldsymbol{\underline{\mathrm{\eta\eta}}}^{\intercal}=\leavevmode\hbox
to66.64pt{\vbox to17.47pt{\pgfpicture\makeatletter\hbox{\hskip
17.47165pt\lower-8.73582pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}}
{{}\pgfsys@rect{-8.53583pt}{-8.53583pt}{17.07166pt}{17.07166pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-2.83888pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\boldsymbol{\underline{\mathrm{\eta}}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {
{}{}{}}{{{{}{}{}}}}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{
{}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}}
{{}\pgfsys@rect{23.1622pt}{-8.53583pt}{17.07166pt}{17.07166pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{27.64247pt}{-3.44165pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\boldsymbol{\underline{\mathrm{\eta}}}^{\intercal}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {
{}{}{}}{{{}}{{}}}{}{{}}{}
{{{}{}}{}}{}{}\pgfsys@moveto{-8.73582pt}{-4.30554pt}\pgfsys@lineto{-17.27165pt}{-4.30554pt}\pgfsys@stroke\pgfsys@invoke{
} { {}{}{}}{{{}}{{}}}{}{{}}{}
{{{}{}}{}}{}{}\pgfsys@moveto{-8.73582pt}{4.30554pt}\pgfsys@lineto{-17.27165pt}{4.30554pt}\pgfsys@stroke\pgfsys@invoke{
} {{}{}{}}{{{}}{{}}}{}{{}}{}{ {}{}{}}{{{}}{{}}}
{}{}{}\pgfsys@moveto{8.73582pt}{-4.30554pt}\pgfsys@lineto{22.9622pt}{-4.30554pt}\pgfsys@stroke\pgfsys@invoke{
} {{}{}{}}{{{}}{{}}}{}{{}}{}
{{{}{}}{}}{}{}\pgfsys@moveto{40.43385pt}{-4.30554pt}\pgfsys@lineto{48.96968pt}{-4.30554pt}\pgfsys@stroke\pgfsys@invoke{
} {{}{}{}}{{{}}{{}}}{}{{}}{}
{{{}{}}{}}{}{}\pgfsys@moveto{40.43385pt}{4.30554pt}\pgfsys@lineto{48.96968pt}{4.30554pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}};\qquad\boldsymbol{\underline{\mathrm{\chi}}}^{2}=\boldsymbol{\underline{\mathrm{\chi}}}.$
(13)

It maps a vector in the extended space onto a single configuration; all
elements corresponding to different configurations are set to $0$. The
projector property follows from the indicator identity
$\boldsymbol{\underline{\mathrm{\eta}}}^{\intercal}\boldsymbol{\underline{\mathrm{\eta}}}=\mathds{1}$.
For the non-interacting system, we define the component Green’s function as

$\displaystyle\boldsymbol{\underline{\mathrm{G}}}(z)\coloneqq\boldsymbol{\underline{\mathrm{\eta}}}\boldsymbol{G}(z)\boldsymbol{\underline{\mathrm{\eta}}}^{\intercal}=\leavevmode\hbox
to108.66pt{\vbox to17.47pt{\pgfpicture\makeatletter\hbox{\hskip
17.47165pt\lower-8.73582pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }
{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}}
{{}{{}}}{{}{}}{}{{}{}}
{{}\pgfsys@rect{-8.53583pt}{-8.53583pt}{17.07166pt}{17.07166pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-2.5pt}{-2.83888pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\boldsymbol{\underline{\mathrm{\eta}}}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {
{}{}{}}{{{{}{}{}}}}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{
{}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}}
{{}\pgfsys@rect{23.1622pt}{-8.53583pt}{27.3966pt}{17.07166pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{26.49521pt}{-2.5pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$\boldsymbol{G}(z)$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {
{}{}{}}{{{{}{}{}}}}{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{
{}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}}
{{}\pgfsys@rect{65.18518pt}{-8.53583pt}{17.07166pt}{17.07166pt}\pgfsys@stroke\pgfsys@invoke{
} }{{{{}}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{69.66545pt}{-3.44165pt}\pgfsys@invoke{
}\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\hbox{{$\boldsymbol{\underline{\mathrm{\eta}}}^{\intercal}$}}
}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}
\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} {
{}{}{}}{{{}}{{}}}{}{{}}{}
{{{}{}}{}}{}{}\pgfsys@moveto{-8.73582pt}{-4.30554pt}\pgfsys@lineto{-17.27165pt}{-4.30554pt}\pgfsys@stroke\pgfsys@invoke{
} { {}{}{}}{{{}}{{}}}{}{{}}{}
{{{}{}}{}}{}{}\pgfsys@moveto{-8.73582pt}{4.30554pt}\pgfsys@lineto{-17.27165pt}{4.30554pt}\pgfsys@stroke\pgfsys@invoke{
} {{}{}{}}{{{}}{{}}}{}{{}}{}{ {}{}{}}{{{}}{{}}}
{}{}{}\pgfsys@moveto{8.73582pt}{-4.30554pt}\pgfsys@lineto{22.9622pt}{-4.30554pt}\pgfsys@stroke\pgfsys@invoke{
} {{}{}{}}{{{}}{{}}}{}{{}}{}{ {}{}{}}{{{}}{{}}}
{}{}{}\pgfsys@moveto{50.7588pt}{-4.30554pt}\pgfsys@lineto{64.98518pt}{-4.30554pt}\pgfsys@stroke\pgfsys@invoke{
} {{}{}{}}{{{}}{{}}}{}{{}}{}
{{{}{}}{}}{}{}\pgfsys@moveto{82.45683pt}{-4.30554pt}\pgfsys@lineto{90.99266pt}{-4.30554pt}\pgfsys@stroke\pgfsys@invoke{
} {{}{}{}}{{{}}{{}}}{}{{}}{}
{{{}{}}{}}{}{}\pgfsys@moveto{82.45683pt}{4.30554pt}\pgfsys@lineto{90.99266pt}{4.30554pt}\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{
{}{}{}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}};$ (14)
$\displaystyle\sum_{\alpha\beta}\boldsymbol{\underline{\mathrm{G}}}^{\alpha\beta}(z)=\boldsymbol{G}(z).$

The arrangement of indicator tensors $\boldsymbol{\underline{\mathrm{\eta}}}$
is different compared to Eq. 12: Both the Green’s function $\boldsymbol{G}$
and the component Green’s function $\boldsymbol{\underline{\mathrm{G}}}$ are
configuration dependent. We note that local elements are diagonal in component
space, i.e.,
$\boldsymbol{\underline{\mathrm{*}}}{G}^{\alpha\beta}_{ii}(z)\propto\delta^{\alpha\beta}$.
We sandwich the resolvent for the Green’s function

$\mathds{1}=[\mathds{1}z-\boldsymbol{H}]\boldsymbol{G}(z)$ (15)

by $\boldsymbol{\underline{\mathrm{\eta}}}$ from the left and
$\boldsymbol{\underline{\mathrm{\eta}}}^{\intercal}$ from the right; this
yields the equation for the component Green’s function

$\boldsymbol{\underline{\mathrm{\chi}}}=[\mathds{1}z-\boldsymbol{\underline{\mathrm{\chi}}}\boldsymbol{\underline{\mathrm{H}}}\boldsymbol{\underline{\mathrm{\chi}}}]\boldsymbol{\underline{\mathrm{G}}}(z).$
(16)

The law of total probability Karlin and Taylor (1975) relates the average of
the component Green’s functions and the conditional average of the physical
Green’s function in the following way:

$\operatorname{\mathbb{E}}(\boldsymbol{\underline{\mathrm{*}}}{G}^{\alpha\beta}_{ij})=\begin{cases}{}\phantom{c^{\beta}}c^{\alpha}\operatorname{\mathbb{E}}(G_{ii}(z)|i\mapsto\alpha)\delta^{\alpha\beta}&\mbox{\quad
for\quad}i=j,\\\
c^{\alpha}c^{\beta}\operatorname{\mathbb{E}}(G_{ij}(z)|i,j\mapsto\alpha,\beta)&\mbox{\quad
for\quad}i\neq j.\end{cases}$ (17)

####  II.2.2 Effective medium in the extended space

As in CPA, in the BEB formalism one calculates an effective local Green’s
function $\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}(z)$ from an
effective medium $\boldsymbol{\underline{\mathrm{S}}}(z)$, which approximates
the average local Green’s function
$\operatorname{\mathbb{E}}(\boldsymbol{\underline{\mathrm{G}}}_{ii}(z))$. We
consider only substitutional disorder without structural disorder, i.e., the
lattice structure is assumed to be fixed. Therefore, we decompose the hopping
tensor
$\boldsymbol{\underline{\mathrm{*}}}{t}^{\alpha\beta}(\absolutevalue{\mathbf{r}_{i}-\mathbf{r}_{j}})$
into its component part,
$\boldsymbol{\underline{\mathrm{*}}}{T}^{\alpha\beta}$, and its lattice part,
$t(\absolutevalue{\mathbf{r}_{i}-\mathbf{r}_{j}})$:

$\boldsymbol{\underline{\mathrm{*}}}{t}^{\alpha\beta}(\absolutevalue{\mathbf{r}_{i}-\mathbf{r}_{j}})\eqqcolon\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\beta}\,t(\absolutevalue{\mathbf{r}_{i}-\mathbf{r}_{j}}).$
(18)

Depending on the component of the endpoints, the matrix elements
$\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\beta}$ scale the amplitudes for
hopping on a given lattice structure by a dimensionless factor. In the
following we refer to $\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\beta}$
simply as “dimensionless hopping parameter”. We perform the lattice Fourier
transform for the hopping matrix elements as

$\boldsymbol{\underline{\mathrm{T}}}\frac{1}{N}\sum_{ij}t(\absolutevalue{\mathbf{r}_{i}-\mathbf{r}_{j}})e^{i\mathbf{k}\cdot(\mathbf{r}_{i}-\mathbf{r}_{j})}=\boldsymbol{\underline{\mathrm{T}}}\epsilon_{k}.$
(19)

For a given effective medium $\boldsymbol{\underline{\mathrm{S}}}(z)$, the
effective local Green’s function reads

$\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}(z)=\frac{1}{N}\sum_{k}{[\mathds{1}z-\boldsymbol{\underline{\mathrm{S}}}(z)-\boldsymbol{\underline{\mathrm{T}}}\epsilon_{k}]}^{-1}.$
(20)

The effective medium as well as the effective local Green’s function are
represented by $M\times M$ matrices in the components. Being local quantities
they no longer carry lattice indices. The effective medium
$\boldsymbol{\underline{\mathrm{S}}}(z)$ is determined by demanding that the
averaged $\boldsymbol{\mathsf{t}}$-matrix vanishes:

$\displaystyle\operatorname{\mathbb{E}}(\boldsymbol{\mathsf{t}}(z))\overset{!}{=}0,$
(21)
$\displaystyle\boldsymbol{\mathsf{t}}(z)\coloneqq-{[\boldsymbol{\underline{\mathrm{\chi}}}{[\mathds{1}-\boldsymbol{\underline{\mathrm{\chi}}}(\mathds{1}z-\boldsymbol{\underline{\mathrm{S}}}(z)+\boldsymbol{\underline{\mathrm{v}}})\boldsymbol{\underline{\mathrm{\chi}}}]}^{-1}\boldsymbol{\underline{\mathrm{\chi}}}+\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}(z)]}^{-1}.$
(22)

###  II.3 Inclusion of electronic interactions and the BEB+DMFT self-
consistency loop

We treat the local Hubbard interaction within the dynamical mean-field theory
Metzner and Vollhardt (1989); Georges  _et al._ (1996); Kotliar and Vollhardt
(2004), which assumes a local self-energy
$\Sigma_{ij}(z)=\delta_{ij}\Sigma_{ii}(z)$; this property becomes exact in the
limit of infinite coordination number. The problem of interacting disordered
electrons may equally be viewed as a system of non-interacting particles
moving in an effective local, energy dependent potential $\Sigma_{ii}(z)$; for
details see Ref. Janiš and Vollhardt, 1992; Vollhardt, 2010. The DMFT self-
consistency equations Georges  _et al._ (1996) are equivalent to a fixed-point
problem which can be expressed by a functional $\widehat{\Sigma}$: Given a
self-energy $\Sigma_{ii}$ and the resulting local Green’s function
$G_{ii}(\Sigma_{ii})$ this functional provides a new self-energy
$\widehat{\Sigma}\big{[}G_{ii}(\Sigma_{ii}),\Sigma_{ii}\big{]}$, such that the
DMFT self-energy is determined self-consistently by the fixed-point

$\Sigma_{ii}=\widehat{\Sigma}\big{[}G_{ii}(\Sigma_{ii}),\Sigma_{ii}\big{]}.$
(23)

Within the CPA, the local Green’s function for a given self-energy
$G_{ii}(\Sigma_{ii})$ is replaced by the conditional average
$\operatorname{\mathbb{E}}(G_{ii}(\Sigma_{ii})|i\mapsto\alpha)=\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}*^{\alpha\alpha}(\Sigma_{ii})/c^{\alpha}$,
see Eqs. 17 and 20. Thus, the self-energy
$\widehat{\Sigma}\big{[}\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}*^{\alpha\alpha}(\Sigma_{ii})/c^{\alpha},\Sigma_{ii}\big{]}$
depends on the component $\alpha$. Consequently, the self-energy at the fixed-
point depends on the component $\alpha$, but not on the explicit site $i$:

$\Sigma^{\alpha}=\widehat{\Sigma}\big{[}\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}*^{\alpha\alpha}(\Sigma^{\alpha})/c^{\alpha},\Sigma^{\alpha}\big{]}.$
(24)

This allows one to introduce the BEB+DMFT self-consistency which we will
discuss next.

By merging the BEB formalism with DMFT a two-fold self-consistency arises, one
for the BEB and one for the DMFT corresponding to the fixed-point Eq. 24. The
self-consistency equation of the BEB formalism is pointwise in the frequencies
and is therefore much simpler than the self-consistency condition of the DMFT,
where frequencies mix due to the energy exchange caused by the interaction
between the electrons. We view the former self-consistence as an internal part
of the full self-consistency loop. In the BEB method we calculate an effective
local Green’s function $\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}(z)$,
Eq. 20. The effective medium $\boldsymbol{\underline{\mathrm{S}}}(z)$ and,
therefore, the effective local Green’s function have to be calculate self-
consistently from Eqs. 21 and 22. This condition simplifies to

$\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}^{-1}(z)=\overline{\boldsymbol{\underline{\mathrm{g}}}}^{-1}(z)$
(25)

with the diagonal matrix

${\overline{\boldsymbol{\underline{\mathrm{*}}}{g}}}^{\alpha\beta}(z)=\frac{c^{\alpha}\delta^{\alpha\beta}}{{(\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}^{-1})}^{\alpha\alpha}(z)+\boldsymbol{\underline{\mathrm{*}}}{S}^{\alpha\alpha}(z)+\mu-\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha}-\Sigma^{\alpha}(z)},$
(26)

where $\Sigma^{\alpha}(z)$ is the DMFT self-energy for the component $\alpha$.
The self-consistent Eq. 25 can be solved with standard root-search algorithms
or by simple iteration. In practice, we use an implementation of the BEB
formalism without interactions and merely shift the on-site energy
$\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha}\rightarrow\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha}+\Sigma^{\alpha}(z)$.
An efficient evaluation of the BEB self-consistency equation is discussed in
Appendix B; an implementation is provided in Ref. Weh and Östlin, 2021. To
emphasize the dependence on the self-energy we denote the self-consistently
determined effective local Green’s function for a given self-energy Eq. 25 by
$\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}(z,\boldsymbol{\Sigma}(z))$.

With the BEB self-consistency condition Eq. 25 for the local Green’s function
$\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}(z,\boldsymbol{\Sigma}(z))$,
the combined algorithm corresponds to the conventional DMFT self-consistency
condition Eq. 24, where the local Green’s function, calculated from the
lattice Hilbert transform, is now replaced by the average

$\operatorname{\mathbb{E}}(G_{ii}|i\mapsto\alpha)=\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}*^{\alpha\alpha}(z,\boldsymbol{\Sigma}(z))/c^{\alpha}.$

The reciprocal concentration factor can be avoided by introducing a
renormalized indicator tensor, which leads to a slightly modified BEB self-
consistency as elaborated in Appendix C. We have to solve a separate impurity
problem for every component $\alpha$. Starting from an initial guess for the
DMFT self-energy $\Sigma^{\alpha}(z)$ for every component, the BEB+DMFT scheme
is the following:

  1. 1.

Calculate the effective local Green’s function Eq. 20 using Eqs. 25 and 26,
which yields

$\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}(z,\boldsymbol{\Sigma}(z)),$
(27a)

  2. 2.

calculate the hybridization function

$\Delta^{\alpha}(z)=z+\mu-\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha}-\Sigma^{\alpha}(z)-c^{\alpha}/\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}*^{\alpha\alpha}(z,\boldsymbol{\Sigma}(z)),$
(27b)

for every component $\alpha$,

  3. 3.

solve the impurity problem for the self-energy

$\Sigma^{\alpha}(z)=\Sigma[\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha},\boldsymbol{\underline{\mathrm{*}}}{U}^{\alpha},\Delta^{\alpha}]$
(27c)

for every component $\alpha$,

  4. 4.

repeat from step 1 until self-consistency is reached.

The hybridization function can also be expressed in terms of BEB quantities
using the self-consistency condition Eq. 25:

$\Delta^{\alpha}(z)=z-\boldsymbol{\underline{\mathrm{*}}}{S}^{\alpha\alpha}(z)-{(\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}^{-1})}^{\alpha\alpha}(z,\boldsymbol{\Sigma}(z)).$
(28)

This is different from CPA+DMFT, where only one unique hybridization functions
exists independent of the alloy components. Analogous to the non-disordered
case, an expression for the hybridization function of the Bethe lattice in
terms of the local Green’s function
$\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}(z,\boldsymbol{\Sigma}(z))$
is given in Appendix D.

Central to the DMFT problem is the impurity solver which provides the local
dynamic self-energy Eq. 27c. To this end, we employ a tensor network based
zero temperature solver, the fork tensor-product state (FTPS)
solverBauernfeind  _et al._ (2017). The FTPS impurity solver is a Hamiltonian-
based method which discretizes the hybridization function Eq. 27b using a
large number of bath sites. We use $249$ sites per spin resulting in a median
energy distance of $0.03D$, where the half-bandwidth $D$ sets our energy
scale. We calculate the ground state $\ket{GS}$ of the finite size impurity
problem using the density matrix renormalization group (DMRG)White (1992);
Schollwöck (2005). Subsequently, we perform the time evolution using the time
dependent variational principle (TDVP)Haegeman  _et al._ (2011); Lubich  _et
al._ (2015); Haegeman  _et al._ (2016); Bauernfeind and Aichhorn (2020). To
obtain the retarded time impurity Green’s function ${G}^{\text{ret}}(t)$, the
states $\hat{c}^{\vphantom{\dagger}}_{\sigma}\ket{GS}$,
$\hat{c}^{\dagger}_{\sigma}\ket{GS}$, as well as their adjoint states are
time-evolved, where $\hat{c}^{\vphantom{\dagger}}_{\sigma}$
($\hat{c}^{\dagger}_{\sigma}$) is the annihilation (creation) operator of the
impurity site. For DMRG we chose a truncated weight of $10^{-15}$ and a
maximal bond-dimension of $100$. We perform the TDVP using time-steps of
$0.1/D$ up to a maximal time $t_{\text{max}}=150/D$ with a truncated weight of
$10^{-9}$ and a maximal bond-dimension of $150$. The convergence with respect
to these parameters is checked.

We can calculate the Green’s function ${G}^{\text{ret}}(t)$ only up to a
maximal time, and we have (small) finite size effects due to the
discretization of the bath. Therefore, we cannot evaluate the retarded Green’s
function directly on the real-frequency axis $G(\omega+i0^{+})$. Instead we
calculate it on a parallel contour $G(\omega+i\eta)$ shifted by a fixed finite
$\eta>0$; this corresponds to the Laplace transform:

$G(\omega+i\eta)=\int_{0}^{\infty}dt\operatorname{e}^{i(\omega+i\eta)t}{G}^{\text{ret}}(t)\eqqcolon\operatorname{F}_{\eta}[{G}^{\text{ret}}(t)](\omega).$
(29)

The shift $\eta$ acts as a broadening for the Green’s function $G(z)$ as can
be seen from the Cauchy integral formula

$2\pi
iG(\omega+i\eta)=\oint\\!dz\frac{G(z)}{z-\omega-i\eta}=\int_{-\infty}^{\infty}\\!\\!d\omega^{\prime}\frac{G(\omega^{\prime})}{\omega^{\prime}-\omega-i\eta^{-}},$
(30)

with $\eta^{-}=\eta-0^{+}$, where $0^{+}$ is a positive infinitesimal. We
write the Green’s function on the real axis in terms of the shifted Fourier
transform

$\displaystyle G(\omega+i0^{+})$ $\displaystyle=\lim_{\eta^{\prime}\searrow
0}\operatorname{F}_{{\eta}^{\prime}}[{G}^{\text{ret}}(t)](\omega)=\operatorname{F}_{\eta}[\operatorname{e}^{t\eta^{-}}{G}^{\text{ret}}(t)](\omega)$
(31)
$\displaystyle=\sum_{k}\frac{1}{k!}{(\eta^{-})}^{k}\operatorname{F}_{\eta}[t^{k}{G}^{\text{ret}}(t)](\omega).$

The second equality replaces the limit ${\eta}^{\prime}\searrow 0$ using
$\lim_{{\eta}^{\prime}\searrow 0}\exp(-{\eta}^{\prime}t)=\exp(-\eta
t)\exp(\eta^{-}t)$ introducing a finite variable $\eta$, in the last line we
use the series representation of the exponential function $\exp(\eta^{-}t)$.
The first term $k=0$ is the Green’s function on the shifted contour
$G(\omega+i\eta)$, higher order terms give systematic corrections. In Section
III, we calculate the first order correction

$G(\omega+i0^{+})=G(\omega+i\eta)+\eta\operatorname{F}_{\eta}\left[t{G}^{\text{ret}}(t)\right](\omega)+\order{\eta^{2}}$
(32)

with a typical shift $\eta=0.08$. The self-energy is calculated from the
equation of motion of the impurity model Bulla  _et al._ (1998)

$\displaystyle\Sigma_{\sigma}(z)$
$\displaystyle=UF_{\sigma}(z)/G_{\sigma}(z),$ (33)
$\displaystyle{F}^{\text{ret}}(t)$
$\displaystyle=\expectationvalue{\hat{c}^{\vphantom{\dagger}}_{\sigma}(t)\hat{n}_{-\sigma}(t)\hat{c}^{\dagger}_{\sigma}}{GS},$
(34)

where $F(z)$ is the Laplace transform of ${F}^{\text{ret}}(t)$.

###  II.4 General properties of the BEB formalism

We shortly review some properties of the BEB formalism Gonis and Garland
(1977); Gonis (1992); Koepernik  _et al._ (1998). First, the BEB formalism is
equivalent to the CPA when off-diagonal disorder is absent. This limit was
already proven in the original formulation Blackman  _et al._ (1971);
Esterling (1975). Since the BEB theory includes the off-diagonal disorder in
the single-site approximation, the Herglotz propertyMüller-Hartmann (1973);
Mills and Ratanavararaksa (1978) of the CPA is preserved in the BEB as well
Gonis and Garland (1977).

Second, for a non-interacting tight-binding Hamiltonian the density of states
(DOS) is non-zero only within certain energy ranges, determined by the
Hamiltonian matrix elements. This is also holds for the CPA Velický  _et al._
(1968). Koepernick _et al._ Koepernik  _et al._ (1998) found the same in their
numerical study of one-dimensional chains using the BEB formalism. Likewise,
we find no violations of this property in our numerical results. In the
following subsection we derive exact bounds for the spectral function of a
Bethe lattice with coordination number $Z$ using the BEB formalism in the
limit of independent alloy components, i.e., when there is no hopping between
different components.

###  II.5 Limit of independent components

We consider the limit of vanishing hopping between different components. If
the hopping is diagonal in the components,
$\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\beta}\propto\delta^{\alpha\beta}$,
the BEB effective medium $\boldsymbol{\underline{\mathrm{S}}}(z)$ is also
diagonal in the components and the self-consistency equations Eq. 25 decouple.
In this case, the effective local Green’s function Eq. 20 can be readily
calculated, since the matrix inverse is the reciprocal of the diagonal
elements

$\displaystyle\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}*^{\alpha\beta}(z)$
$\displaystyle=\frac{1}{N}\sum_{k}\frac{\delta^{\alpha\beta}}{z-\boldsymbol{\underline{\mathrm{*}}}{S}^{\alpha\alpha}(z)-\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\alpha}\epsilon_{k}}$
(35)
$\displaystyle=\delta^{\alpha\beta}g_{0}^{\alpha}\big{(}z-\boldsymbol{\underline{\mathrm{*}}}{S}^{\alpha\alpha}(z)\big{)}.$

Here, $g_{0}^{\alpha}$ is the lattice Hilbert transform
$g_{0}(z)=\frac{1}{N}\sum_{k}\frac{1}{z-\epsilon_{k}}$; the superscript
$\alpha$ indicates that the bandwidth is scaled by
$\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\alpha}$:

$g_{0}^{\alpha}(z)=\frac{1}{N}\sum_{k}\frac{1}{z-\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\alpha}\epsilon_{k}}=\frac{1}{\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\alpha}}g_{0}(z/\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\alpha}).$
(36)

For the component $\alpha$, the decoupled self-consistency Eq. 25 reads

$0=\frac{\overline{c}^{\alpha}}{g_{0}^{\alpha}(z-\boldsymbol{\underline{\mathrm{*}}}{S}^{\alpha\alpha})}+\boldsymbol{\underline{\mathrm{*}}}{S}^{\alpha\alpha}-\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha},$
(37)

with the concentration complement $\overline{c}^{\alpha}=1-c^{\alpha}\geq 0$.
For a Bethe lattice with coordination number $Z$ and the lattice Hilbert
transformEconomou (2006)

$\displaystyle
g_{0}(z,Z)=2(Z-2)/\left[z\left(Z-2+Z\sqrt{1-D^{2}/z^{2}}\right)\right],$ (38)

where $D$ is the half-bandwidth, the self-consistency condition is an
algebraic equation and can be solved analytically. The BEB effective medium
reads

$\boldsymbol{\underline{\mathrm{*}}}{S}^{\alpha\alpha}(z,Z)=\frac{(Z-2)\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha}c^{\alpha}+Z\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha}+(Z+2)c^{\alpha}z-Zz-2{(c^{\alpha})}^{2}z-Z\overline{c}^{\alpha}s\sqrt{{(z-\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha})}^{2}-c^{\alpha}{(D^{\alpha})}^{2}\frac{Z-c^{\alpha}}{Z-1}}}{2c^{\alpha}(Z-c^{\alpha})},$
(39)

where $s$ is the sign $s=\operatorname{sign}(\real(z-v^{\alpha}))$, and
$D^{\alpha}$ is the half-bandwidth scaled by
$\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\alpha}$; this is the retarded
solution. A conjugate solution exists with $-s$ and therefore with a plus sign
in front of the square root.

We are interested in the bandwidth of the resulting component spectrum

$A^{\alpha}(\omega)=-\frac{1}{c^{\alpha}\pi}\imaginary
g_{0}^{\alpha}(\omega+i0^{+}-\boldsymbol{\underline{\mathrm{*}}}{S}^{\alpha\alpha}(\omega+i0^{+})).$
(40)

For non-interacting systems, the Gershgorin circle theoremGershgorin (1931)
gives the _maximal_ spectral bounds

$\absolutevalue{z-\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha}}\leq
D^{\alpha}.$ (41)

In the limit
$\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\beta}\propto\delta^{\alpha\beta}$,
we can make a more precise statement and derive exact spectral bounds as will
be discussed below. The spectral function can only vanish when the imaginary
part of the effective medium vanishes. Thus, for non-interacting systems, we
need to check where the argument of the square root is negative. One finds an
imaginary part and therefore spectral weight for

$\absolutevalue{z-\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha}}<\sqrt{c^{\alpha}\frac{Z-c^{\alpha}}{Z-1}}D^{\alpha}.$
(42)

Therefore, for the Bethe lattice with coordination number $Z$ and
$\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\beta}\propto\delta^{\alpha\beta}$,
the bandwidth is reduced due to concentration by a factor
$\sqrt{c^{\alpha}(Z-c^{\alpha})/(Z-1)}$. We then obtain the effective
bandwidth

$D^{\alpha}_{\text{eff}}=\sqrt{c^{\alpha}\frac{Z-c^{\alpha}}{Z-1}}\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\alpha}D.$
(43)

Our numerical results in Section III were obtained for a semicircular DOS,
i.e., the Bethe lattice with infinite coordination number
$Z\rightarrow\infty$. In this limit one finds an effective bandwidth

$D^{\alpha}_{\text{eff}}=\sqrt{c^{\alpha}}\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\alpha}D.$
(44)

The same factor $\sqrt{c}$ was found in Ref. Byczuk  _et al._ , 2004 in the
CPA
($\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{A}}=\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{B}\mathsf{B}}=1$)
in the limit of high disorder strength
$(\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}-\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{A}})/D=\delta\gg\max(1,U/D)$.
While the parameters in these limits are different, both describe the same
physics, namely the decoupling of components. Indeed, the components decouple
not only for vanishing hopping between the components
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=0$, but also in
the case of a large separation in energy ($\delta\gg 1$).

For coordination number $Z=2$ another interesting limit of the Bethe lattice
is obtained; this is the one-dimensional latticeEconomou (2006), where

$g^{1\text{D}}_{0}(z)=g_{0}(z,Z=2)=1/\left[z\sqrt{1-D^{2}/z^{2}}\right].$ (45)

The spectral bounds are given by

$D^{\alpha}_{\text{eff}}=\sqrt{c^{\alpha}(2-c^{\alpha})}\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\alpha}D.$
(46)

Therefore, for the one-dimensional lattice and
$\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\beta}\propto\delta^{\alpha\beta}$
the bandwidth is reduced by the factor $\sqrt{c^{\alpha}(2-c^{\alpha})}$.

##  III Numerical results

The above formalism is now used to study the effect of off-diagonal disorder
in the Anderson–Hubbard model at zero temperature. We employ a Bethe lattice
with infinite coordination number, whose half-bandwidth $D$ sets the energy
scale. Furthermore we consider a discrete binary random alloy distribution
with components $\mathsf{A}$ and $\mathsf{B}$.

In all applications we consider the case of half-filling on
average$\operatorname{\mathbb{E}}(n_{i})=1$; this leads to a fixed chemical
potential which we choose as $\mu=0$. In the following subsection we fix the
alloy component concentration and study the change in the spectral function
starting with the non-interacting case and equal atomic potentials. The alloy
component spectral functions are the concentration-weighted conditional
spectral functions

$\boldsymbol{\underline{\mathrm{*}}}{A}^{\alpha}(\omega)\coloneqq-\frac{1}{\pi}\imaginary\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}*^{\alpha\alpha}(\omega)=-\frac{c^{\alpha}}{\pi}\imaginary\operatorname{\mathbb{E}}(G_{ii}(\omega)|i\mapsto\alpha).$
(47)

The average spectral functions are given by the trace

$\operatorname{\mathbb{E}}(A(\omega))=\sum_{\alpha}\boldsymbol{\underline{\mathrm{*}}}{A}^{\alpha}(\omega)=-\frac{1}{\pi}\imaginary\Tr\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}(\omega).$
(48)

###  III.1 Non-interacting limit

We start with the non-interacting case by setting
$\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{A}}=\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{B}}=0$,
which corresponds to the Anderson disorder model with purely off-diagonal
disorder. Since the non-interacting Green’s function is independent of
temperature the results presented in this section are valid not only for zero,
but also for finite temperatures. We choose the parameters

$\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{A}}=\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}=-U/2=0;\leavevmode\nobreak\
c^{\mathsf{A}}=0.1=1-c^{\mathsf{B}};\leavevmode\nobreak\
\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{A}}=\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{B}\mathsf{B}}=1$

and calculate the average and the alloy component spectral functions for
several values of
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$ at half filling.
The case $\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=1$ is
equivalent to the non-disordered case since
$\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{A}}=\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}$;
in this case the components are indistinguishable. Thus, the average spectral
function is just the spectral function of the non-disordered Bethe lattice,
and the component Green’s function are proportional with a concentration
prefactor.

Figure 1: Non-interacting case: Comparison of spectral functions for different
dimensionless hopping parameters
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$. The parameters
are
$\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{A}}=\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{B}}=0$,
$\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{A}}=\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}=0$,
$c^{\mathsf{A}}=0.1=1-c^{\mathsf{B}}$,
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{A}}=\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{B}\mathsf{B}}=1$.
The solid lines represent the component spectral functions
$\boldsymbol{\underline{\mathrm{*}}}{A}^{\alpha}(z)=c^{\alpha}\operatorname{\mathbb{E}}(A_{i}(z)|i\mapsto\alpha)$,
where $\mathsf{A}$ is red and $\mathsf{B}$ is blue; the dotted yellow line
shows the average spectral function
$\operatorname{\mathbb{E}}(A(z))=\boldsymbol{\underline{\mathrm{*}}}{A}^{\mathsf{A}}(z)+\boldsymbol{\underline{\mathrm{*}}}{A}^{\mathsf{B}}(z)$.
The thin vertical lines show the maximal spectral bounds given by the
Gershgorin circle theoremGershgorin (1931).

Figure 1 shows the spectral function for off-diagonal disorder with
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=0.0,0.5,1.5,5.0$.
The case $\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=0$ was
solved exactly in Section II.5. The panel
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=0$ in Fig. 1
indicates that off-diagonal disorder reduces the bandwidths; according to Eq.
43 the effective bandwidths are given by
$D_{\text{eff}}^{\mathsf{A}}=\sqrt{0.1}D\approx 0.32D$ and
$D_{\text{eff}}^{\mathsf{B}}=\sqrt{0.9}D\approx 0.95D$. For
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}<1=\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\alpha}$
the probabilities for hopping between the alloy components $\mathsf{A}$ and
$\mathsf{B}$ are less than those between the same component $\alpha$. The
spectral functions in the upper half of Fig. 1 correspond to this situation.
In spite of a similar support on the energy axis, the spectral function of the
majority component $\mathsf{B}$ has a larger bandwidth, which encompasses the
effective bandwidth of component $\mathsf{A}$. By contrast, when
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}>1=\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\alpha}$,
$\mathsf{A}$-$\mathsf{B}$ bonds are energetically favorable. The panel
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=1.5$ in Fig. 1
shows that the spectral function of component $\mathsf{A}$ develops shoulders
although both components have similar effective bandwidths. When the value of
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$ is increased
further the shoulders split off from the central band; for the parameters
chosen the split-off is visible for
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}\geq 2.25$.
According to Burdin and Fulde Burdin and Fulde (2007) the split-off upper and
lower bands correspond to bonding and antibonding states, respectively. In the
particle-hole symmetric case, the bonding and antibonding bands in panel
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=5.0$ of Fig. 1
have equal weights. Due to the large value
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=5.0$, the
minority component $\mathsf{A}$ is completely suppressed in the central band.
The components $\mathsf{A}$ and $\mathsf{B}$ contribute roughly equally to the
bonding and antibonding subbands. Therefore the central band for the
$\mathsf{B}$ component is depleted by an amount of $(1-2c_{A})=0.8$ of the
spectral function. The overall results and the spectral weight transfer from
the central band are consistent with those reported in Ref. Burdin and Fulde,
2007.

###  III.2 Alloy components with equal interaction strengths
$\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{A}}=\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{B}}$

In the following, we will discuss the results for the interacting case using
the setup described in Section III.1 at zero temperature ($T=0$). The alloy
components have identical on-site interaction parameters:

$\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{A}}=\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{B}}=U=3D.$

At half-filling, the on-site energies are
$\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{A}}=\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}=-U/2=-1.5D$.

Figure 2: Comparison of spectral functions for different values of the
dimensionless hopping parameter
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$ with
$\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{A}}=\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{B}}=3D$,
$\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{A}}=\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}=-1.5D$,
$c^{\mathsf{A}}=0.1=1-c^{\mathsf{B}}$, and
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{A}}=\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{B}\mathsf{B}}=1$
at $T=0$. For $\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=0$
a shift $\eta=0.12$ had to be used; the other panels were calculated for
$\eta=0.08$.

Figure 2 shows the spectral function for various values of
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$. For
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=0$, when the
self-consistency equations decouple, the spectral function of both components
implies insulating behavior due to the strong interaction $U=3D$. The upper
and lower Hubbard bands centered around $\pm U/2$ are visible. The bandwidth
of the Hubbard bands is effectively reduced according to Eq. 43. Therefore the
gap is wider for component $\mathsf{A}$. Finite values of
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$ lead to wider
spectral functions, and the bandwidth of the minority component $\mathsf{A}$
broadens to the same bandwidth as for $\mathsf{B}$. Although for
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=1.5$ the
imaginary part of the self-energy shows a prominent peak at $\omega=0$, the
spectral function remains finite at the Fermi level. The convergence of the
DMFT computations slows down for this value of
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=1.5$, hinting at
the proximity of a transition. The minority component $\mathsf{A}$ exhibits
shoulders at the band-edge. For
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=1.7$ (not shown)
one observes a pronounced quasi-particle peak at the Fermi level of both
components, indicating that the system is metallic. A further increase to
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=5.0$ leads to an
increased spectral weight at the Fermi level for the majority component
$\mathsf{B}$, while the spectral function of $\mathsf{A}$ has a minimum at the
Fermi level.

For small values of
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$ the spectral gap
results from the local Hubbard physics. Disorder then plays a minor role and
mostly modifies the bandwidth and therefore the gap size. An increase of
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$ leads to a
larger bandwidth compared to that of the CPA+DMFT result for
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=1$. For larger
values of $\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$ the
spectral function of the component $\mathsf{A}$ is seen to open a pseudogap
around the Fermi level which is accompanied by an increase of spectral weight
of the component $\mathsf{B}$. For large
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$, the pseudogap
is a result of the off-diagonal disorder.

###  III.3 Alloy components with different interaction strengths

In a binary alloy the strength of the interaction between electrons may also
depend on the alloy component. Therefore we explore the effect of off-diagonal
disorder in this case. We illustrate the results for an extreme case, namely
for a strong repulsion
$\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{B}}=3D$ of the majority
component $\mathsf{B}$ only, while the minority component remains non-
interacting ($\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{A}}=0$). We
consider half-filling with
$\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{A}}=0$,
$\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}=-\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{B}}/2=-1.5D$
and note that, in spite of the different values
$\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{A}}\neq\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}$,
the effective (diagonal) disorder strength is 0, since the Hartree self-energy
compensates the difference.

Figure 3: Comparison of spectral functions for different values of the
dimensionless hopping parameter
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$ with
$\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{A}}=0$,
$\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{B}}=3D$,
$\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{A}}=0$,
$\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}=-1.5D$,
$c^{\mathsf{A}}=0.1=1-c^{\mathsf{B}}$,
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{A}}=\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{B}\mathsf{B}}=1$
at $T=0$.

Figure 3 shows the evolution of the spectral function for increasing
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$. For
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=0$, the
$\mathsf{A}$ alloy component is metallic, while due to the large
$\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{B}}$ value the $\mathsf{B}$
component is insulating. We note that the two components have different
effective bandwidths due to the different concentrations. The panel with
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=0.5$ shows a
small peak for the $\mathsf{B}$ component, in spite of the large interaction
strength. At $\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=1$,
the $\mathsf{A}$-$\mathsf{A}$, $\mathsf{A}$-$\mathsf{B}$, and
$\mathsf{B}$-$\mathsf{B}$ hopping probabilities are the same, which leads to
the same effective bandwidths, and to the appearance of the metallic state for
both alloy components.

In Fig. 3 the panel with
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=1.5$ shows a
distinct peak for the majority component $\mathsf{B}$ at the Fermi level,
which reduces the spectral function for the $\mathsf{A}$ component at the
Fermi level, leading to a local minimum. Increasing the inter-component
hopping to $\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=5.0$,
the peak of $\mathsf{B}$ becomes even larger, and the spectral weight of
$\mathsf{A}$ almost vanishes at the Fermi level. The panels with
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=5.0$ of Figs. 2
and 3 are seen to be very similar; apparently the interaction of the minority
component has little effect on the spectral function.

Figure 4: Quasi-particle weight $Z$ corresponding to
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=5.0$ in Figs. 2
and 3 for parameters
$\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{A}}=-\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{A}}/2$,
$\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}=-\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{B}}/2$,
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{A}}=\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{B}\mathsf{B}}=1$,
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=5.0$ at $T=0$,
calculated for a shift $\eta=0.12$.

Figure 4 shows the quasi-particle weight

$Z^{\alpha}={\left[1-\evaluated{\partialderivative{\real\Sigma^{\alpha}(\omega+i\eta)}{\omega}}_{\omega=0}\right]}^{-1}$
(49)

corresponding to the spectral functions of panel
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=5.0$ in the Figs.
2 and 3. In spite of the large value of
$\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{A}}=3D$, the quasiparticle
weight $Z^{\mathsf{A}}$ is large, with a magnitude around $0.9$. This gives an
indication why the panels
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=5.0$ of Figs. 2
and 3 are so similar; the large value of
$\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{A}}=3D$ leads only to a small
mass renormalization. Increasing the concentration of the weakly correlated
component $\mathsf{A}$ leads to a significant increasing of the quasiparticle
weight $Z^{\mathsf{B}}$ for both setups. This can be explained by the
increasing number of $\mathsf{A}$-$\mathsf{B}$ bonds, which leads to increased
mobility of the hopping due to the large value of
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=5.0$ compared to
the inter-component hoppings
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{A}}=\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{B}\mathsf{B}}=1$.

###  III.4 Combined effect of diagonal and off-diagonal disorder

In the following, we explore the combined effect of both diagonal and off-
diagonal disorder, and their interplay with interaction. We choose a uniform
interaction strength
$\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{A}}=\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{B}}=U$
and introduce diagonal disorder with on-site potentials
$\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{A}}=-1.5D-U/2$,
$\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}=+1.5D-U/2$. This means
that the scattering strength is of the magnitude
$\delta=(\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}-\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{A}})/D=3$.
We consider components with equal bandwidth
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{A}}=\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{B}\mathsf{B}}=1$
and equal concentration $c^{\mathsf{A}}=c^{\mathsf{B}}=0.5$. Thus, the
components are particle-hole conjugate and fulfill the relation

$\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}^{\mathsf{A}\mathsf{A}}(z)=-{[\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}^{\mathsf{B}\mathsf{B}}(z)]}^{*}.$
(50)

Figure 5: Comparison of the spectral functions for
$\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{A}}=\boldsymbol{\underline{\mathrm{*}}}{U}^{\mathsf{B}}=U$,
$-(\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{A}}+U/2)=\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}+U/2=1.5D$,
$c^{\mathsf{A}}=0.5=c^{\mathsf{B}}$,
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{A}}=\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{B}\mathsf{B}}=1$,
$T=0$ for different values of
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$ and $U$.

Figure 5 shows the average spectral function as well as that for the
individual components for different values of the Hubbard parameter $U$ and
dimensionless inter-component hopping amplitudes
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$. In the case
$U=0$ (first column of Fig. 5) one starts from the split-band limit.

We first investigate the CPA limit
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=1$ (second row of
Fig. 5). In the split band limit, there are no correlation effects: One
component is basically filled $n_{\sigma}^{\mathsf{A}}\approx 1$, the other is
depleted $n_{\sigma}^{\mathsf{B}}\approx 0$. However, the Hartree energy,
$\Sigma^{\text{H}}_{\sigma}=n_{\sigma}U$, decreases the effective disorder
strength

$\displaystyle\delta_{\text{eff}}$
$\displaystyle=\frac{(\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}+n_{\sigma}^{\mathsf{B}}U)-(\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{A}}+n^{\mathsf{A}}U)}{D}$
(51)
$\displaystyle\approx[\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}-\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{A}}-U]/D=\delta-U/D.$

Switching on the interaction $U$ effectively decreases the scattering strength
$\delta$. From $U\approx 2D$ on, the split-band limit at large scattering
strength no longer applies, i.e., there is a combination of disorder and
interaction effects. For $U=4D$, we see the upper and lower Hubbard bands for
each component, as well as a quasiparticle peak at the Fermi level. For even
larger interaction strength ($U=6D$) a Mott insulating phase is observed.
Thus, by increasing $U$ it is possible to tune the system from an alloy-band
insulator, through a metallic phase, to a Mott insulating state. Similar
results were reported by Lombardo _et al._ Lombardo  _et al._ (2006) for
diagonal disorder using CPA+DMFT for somewhat different parameters and a
finite-temperature impurity solver.

The behavior obtained in the CPA limit can now be modified by varying
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$. At $U=2D$ an
off-diagonal hopping
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}<1$ leads to
metallic behavior, while
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}>1$ favors a band
gap. On the other hand, for $U=6D$ a large
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$ favors
metallicity. For
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}\leq 1$ the
spectral function is gapped — similar to the result obtained in the Hubbard-I
approximation Hubbard (1963).

##  IV Summary

We presented a methodological framework for the study of interacting electrons
in the presence of diagonal and off-diagonal disorder. The formalism allows
one to explore a multi-component system of electrons with random on-site
interactions and/or potentials as well as random hopping amplitudes. For this
purpose the Blackman, Esterling, Berk (BEB) formalism for averaging over off-
diagonal disorder was combined with the dynamical mean-field theory (DMFT). We
introduced a tensor notation inspired by Koepernik _et al._ Koepernik  _et
al._ (1997, 1998), by which the randomness of the components is transferred to
local indicator tensors defined in an extended space. In this representation,
the problem can be solved in a single-site approximation, analogous to the
coherent potential approximation (CPA). The computational procedure, including
the two-fold self-consistency and the impurity solverBauernfeind  _et al._
(2017), were discussed in detail.

In the limit of zero intercomponent-hopping an analytic solution for the BEB
approximation of the non-interacting Bethe lattice with the general
coordination $Z$ was provided. We established the exact bounds of the spectral
function, and showed that these lie within the maximal bound given by the
Gershgorin circle theoremGershgorin (1931).

The first application of the BEB+DMFT formalism in tensor formulation was the
computation of spectral functions. In particular, we discussed the changes in
the spectral function for increasing dimensionless hopping amplitude
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$ and interaction
strengths U. For alloy components with the same local interaction, we found a
pseudo-gap in the spectral function of the alloy component with lower
concentration. For larger values of
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$ this is
accompanied by quasiparticles of the dominant alloy component. The pseudo-gap
and the quasiparticle character were inferred by analyzing the respective
self-energies. We found a rather similar behavior in the case of alloy
components with different interaction strengths. For the discussion of both
diagonal and off-diagonal disorder we studied equal alloy concentrations and
equal alloy component bandwidths for a rather large scattering strength
$\delta=3D$. Increasing $U$ while keeping
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}\leq 1.0$ fixed
the electronic system was found to undergo two transitions: One from the non-
correlated alloy insulator to a metallic state, and then from the metal to a
correlated (Mott) insulator. For larger values of
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}$ the formation of
quasiparticles signals the appearance of a metallic state, which eventually
disappears again within the Mott insulating phase for large enough values of
$U$.

Finally, we note that at the CPA+DMFT level our results are in agreement with
previous model-Hamiltonian based publications. As a matter of fact the BEB
formulation is applicable to band structure schemes Koepernik  _et al._ (1997,
1998); Herbig  _et al._ (2017) also for non-orthogonal basis sets. The present
formulation can be naturally extended following the CPA+DMFT Östlin  _et al._
(2018) methodology with similar computational costs. Work along this line is
in progress.

###### Acknowledgements.

Financial support by the Deutsche Forschungsgemeinschaft through TRR80
(project F6) Project number 107745057 is gratefully acknowledged. Y.Z. is
supported in part by the NSF China Grant No. 12004383 and the Fundamental
Research Funds for the Central Universities. H.T. gratefully acknowledges
support from NSF OAC-1931367 and NSF DMR-1944974 grants. K.M.T. is partially
supported by NSF DMR-1728457 and NSF OAC-1931445 grants.

##  Appendix A An explicit example of the extended Hamiltonian matrix

Here we provide an explicit example of the Hamiltonian matrix $\boldsymbol{H}$
and the extended Hamiltonian $\boldsymbol{\underline{\mathrm{H}}}$. We
consider a small system — a chain with $3$ sites ($1,2,3$) and $2$ components
$\mathsf{A}$ and $\mathsf{B}$. For this chain the Hamiltonian matrix reads

$\boldsymbol{H}=\begin{pmatrix}v_{1}&t_{12}&t_{13}\\\ t_{21}&v_{2}&t_{23}\\\
t_{31}&t_{32}&v_{3}\end{pmatrix}-\mathds{1}\mu.$ (52)

This is a random matrix since it depends on the configuration. We consider
$c^{\mathsf{A}}=1/3$ and $c^{\mathsf{B}}=2/3$; then one possible configuration
is $\mathsf{A}-\mathsf{B}-\mathsf{B}$. According to Eq. 8, the Hamiltonian
matrix of this configuration takes the values

$\boldsymbol{H}_{\mathsf{A}\mathsf{B}\mathsf{B}}=\begin{pmatrix}\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{A}}&\boldsymbol{\underline{\mathrm{*}}}{t}^{\mathsf{A}\mathsf{B}}(a)&\boldsymbol{\underline{\mathrm{*}}}{t}^{\mathsf{A}\mathsf{B}}(2a)\\\
\boldsymbol{\underline{\mathrm{*}}}{t}^{\mathsf{B}\mathsf{A}}(a)&\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}&\boldsymbol{\underline{\mathrm{*}}}{t}^{\mathsf{B}\mathsf{B}}(a)\\\
\boldsymbol{\underline{\mathrm{*}}}{t}^{\mathsf{B}\mathsf{A}}(2a)&\boldsymbol{\underline{\mathrm{*}}}{t}^{\mathsf{B}\mathsf{B}}(a)&\boldsymbol{\underline{\mathrm{*}}}{v}^{\mathsf{B}}\end{pmatrix}-\mathds{1}\mu,$
(53)

where $a$ is the distance between neighboring sites. By contrast, the extended
Hamiltonian matrix $\boldsymbol{\underline{\mathrm{H}}}$ does not depend on
the specific configuration. We can choose a matrix representation of the
tensor $\boldsymbol{\underline{\mathrm{H}}}$ shown in Eq. 12, by grouping the
legs $i$ and $\alpha$ on the same side. This is done explicitly defining the
combined index $n=(i,\alpha)$ [$m=(j,\beta)$]. We count $n=2i-1+n_{\alpha}$
with $n_{\mathsf{A}}=0$ and $n_{\mathsf{B}}=1$. Then the extended Hamiltonian
$\boldsymbol{\underline{\mathrm{*}}}{H}^{\alpha\beta}_{ij}=\boldsymbol{\underline{\mathrm{*}}}{H}_{nm}$
reads

$\boldsymbol{\underline{\mathrm{H}}}+\mathds{1}\mu=\\\ \pmqty{undef}.$ (54)

This $MN\times MN=6\times 6$ matrix contains all $M^{N}$ possible
configurations for the $N=3$ site problem with $M=2$ components and is
independent of the concentrations $c^{\alpha}$. A specific configuration can
be selected by applying an appropriate indicator tensor
$\boldsymbol{\underline{\mathrm{\eta}}}$. For the configuration
$\mathsf{A}-\mathsf{B}-\mathsf{B}$,
$\boldsymbol{\underline{\mathrm{*}}}{\eta}^{\alpha}_{ij}=\boldsymbol{\underline{\mathrm{*}}}{\eta}_{nj}$
takes the form

$\boldsymbol{\underline{\mathrm{\eta}}}^{\intercal}_{\mathsf{A}\mathsf{B}\mathsf{B}}=\pmqty{undef},$
(55)

and we obtain the Hamiltonian matrix for this configurations from

$\boldsymbol{H}_{\mathsf{A}\mathsf{B}\mathsf{B}}=\boldsymbol{\underline{\mathrm{\eta}}}^{\intercal}_{\mathsf{A}\mathsf{B}\mathsf{B}}\boldsymbol{\underline{\mathrm{H}}}\boldsymbol{\underline{\mathrm{\eta}}}_{\mathsf{A}\mathsf{B}\mathsf{B}}.$
(56)

While the matrix representation is suitable for performing calculations, the
tensors notation is often clearer. The elements diagonal in components
$\alpha$ of the extended Hamiltonian read

$\boldsymbol{\underline{\mathrm{H}}}^{\alpha\alpha}=\pmqty{undef}-\mathds{1}\mu,$
(57)

the off-diagonal elements $\alpha\neq\beta$ read

$\boldsymbol{\underline{\mathrm{H}}}^{\alpha\beta}=\pmqty{undef},$ (58)

and the indicator tensor for the configuration
$\mathsf{A}-\mathsf{B}-\mathsf{B}$ is given by the elements

$\boldsymbol{\underline{\mathrm{\eta}}}^{\mathsf{A}}_{\mathsf{A}\mathsf{B}\mathsf{B}}=\pmqty{undef}\mbox{\quad
and\quad}\boldsymbol{\underline{\mathrm{\eta}}}^{\mathsf{B}}_{\mathsf{A}\mathsf{B}\mathsf{B}}=\pmqty{undef}.$
(59)

##  Appendix B Efficient evaluation of the local Green’s function and its
inverse

To solve the BEB self-consistency equation, we need to repeatedly evaluate the
effective local Green’s function

$\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}(z)=\frac{1}{N}\sum_{k}{[\boldsymbol{\underline{\mathrm{\xi}}}(z)-\boldsymbol{\underline{\mathrm{T}}}\epsilon_{k}]}^{-1},$
(60)

with
$\boldsymbol{\underline{\mathrm{\xi}}}(z)=\mathds{1}z-\boldsymbol{\underline{\mathrm{S}}}(z)$,
or rather its inverse
${(\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}})}^{-1}$. A naive
evaluation would be computational costly, since one needs to invert a matrix
for every $k$-point and every frequency point. While this is feasible for
small matrices, it has the potential risk of inaccurate $k$-summations (or
integrations), especially for a DOS with singularities as for a one-
dimensional or square lattice. Therefore, we employ an algorithm based on the
_compact_ singular value decomposition (SVD) of the matrix
$\boldsymbol{\underline{\mathrm{T}}}$

$\boldsymbol{\underline{\mathrm{T}}}=\boldsymbol{U\sigma
V}^{\dagger}=\boldsymbol{U\sigma}^{1/2}\,\boldsymbol{\sigma}^{1/2}\boldsymbol{V}^{\dagger}\eqqcolon\tilde{\boldsymbol{U}}\tilde{\boldsymbol{V}}^{\dagger},$
(61)

which we use to split the matrix; here we partitioned the singular values
symmetrically as $\tilde{\boldsymbol{U}}=\boldsymbol{U\sigma}^{1/2}$,
$\tilde{\boldsymbol{V}}^{\dagger}=\sigma^{1/2}\boldsymbol{V}^{\dagger}$. It is
important to use the compact SVD as we will explicitly use the inverse
$\boldsymbol{\sigma}^{-1}$. Numerically, the need to truncate small singular
values arises. We note that for the binary alloy the rank-1 case, where the
SVD has to be truncated, is given for a hopping matrix of the type

$\begin{pmatrix}\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{A}}&\sqrt{\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{A}}\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{B}\mathsf{B}}}\\\
\sqrt{\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{A}}\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{B}\mathsf{B}}}&\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{B}\mathsf{B}}\end{pmatrix}\\!=\\!\begin{pmatrix}\sqrt{\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{A}}}\\\
\sqrt{\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{B}\mathsf{B}}}\end{pmatrix}\\!\\!\\!\begin{pmatrix}\sqrt{\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{A}}}&\sqrt{\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{B}\mathsf{B}}}\end{pmatrix}.$
(62)

This is the structure of the hopping matrix discussed by Shiba Shiba (1971).
Another prominent rank-1 example is the CPA limit with
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{A}}=\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}=\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{B}\mathsf{B}}=1$.
The matrix inverse
${[\boldsymbol{\underline{\mathrm{\xi}}}(z)-\boldsymbol{\underline{\mathrm{T}}}\epsilon_{k}]}^{-1}$
is calculated using the Woodbury matrix identity Higham (2002). Furthermore,
we calculate the eigendecomposition

$\tilde{\boldsymbol{V}}^{\dagger}\boldsymbol{\underline{\mathrm{\xi}}}^{-1}(z)\tilde{\boldsymbol{U}}=\boldsymbol{P}(z)\boldsymbol{d}(z)\boldsymbol{P}^{-1}(z),$
(63)

where $\boldsymbol{d}(z)$ is the diagonal matrix of eigenvalues. The
$k$-dependent Green’s function can be expressed as

$\boldsymbol{\underline{\mathrm{G}}}(z,k)=\boldsymbol{\underline{\mathrm{\xi}}}^{-1}-\boldsymbol{\underline{\mathrm{\xi}}}^{-1}\tilde{\boldsymbol{U}}\boldsymbol{P}{\left[\boldsymbol{d}-\frac{\mathds{1}}{\epsilon_{k}}\right]}^{-1}\boldsymbol{P}^{-1}\tilde{\boldsymbol{V}}^{\dagger}\boldsymbol{\underline{\mathrm{\xi}}}^{-1},$
(64)

where we did not write the $z$-dependence explicitly. We note that only the
term in the square bracket depends on $k$. Since it contains only diagonal
matrices, the matrix inverse only involves the reciprocal matrix elements. We
look at a particular diagonal element with the $z$-dependent eigenvalue
$\boldsymbol{d}_{ii}(z)=\lambda_{i}(z)$

$-{\left[\lambda_{i}(z)-\frac{1}{\epsilon_{k}}\right]}^{-1}=\frac{1}{1/\lambda_{i}(z)-\epsilon_{k}}\frac{1}{\lambda_{i}^{2}(z)}-\frac{1}{\lambda_{i}(z)}.$
(65)

It is straightforward to perform the $k$-summation, as we have the standard
form of the lattice Hilbert transform

$\frac{1}{N}\sum_{k}\frac{1}{z-\epsilon_{k}}\eqqcolon g_{0}(z),$ (66)

evaluated at $1/\lambda_{i}$. For simple lattices like the Bethe lattice we
know the analytic expression for $g_{0}$; a numerical integration can then be
avoided. For the local Green’s function we obtain the lengthy expression

$\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}=\boldsymbol{\underline{\mathrm{\xi}}}^{-1}+\boldsymbol{\underline{\mathrm{\xi}}}^{-1}\tilde{\boldsymbol{U}}\boldsymbol{P}{[\boldsymbol{d}^{-1}g_{0}(\boldsymbol{d}^{-1})-\mathds{1}]}\boldsymbol{d}^{-1}\boldsymbol{P}^{-1}\tilde{\boldsymbol{V}}^{\dagger}\boldsymbol{\underline{\mathrm{\xi}}}^{-1},$
(67)

where again the $z$-dependence was not written out. For the self-consistency
equation of the BEB formalism one only needs the inverse of
$\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}$. The Woodbury matrix
identity yields the simpler expression,

$\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}^{-1}(z)=\\\
\boldsymbol{\underline{\mathrm{\xi}}}(z)+\tilde{\boldsymbol{U}}\boldsymbol{P}(z)\left[\frac{1}{g_{0}\big{(}1/\boldsymbol{d}(z)\big{)}}-1/\boldsymbol{d}(z)\right]\boldsymbol{P}^{-1}(z)\tilde{\boldsymbol{V}}^{\dagger},$
(68)

where we explicitly noted the inverse of the diagonal matrices by reciprocal
matrix elements.

Considering that the main cost of a naive evaluation of
$\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}^{-1}$ arises from the matrix
inversion, this amounts to $N_{z}(N_{k}+1)$ matrix inversions, where $N_{z}$
is the number of frequency points and $N_{k}$ is the number of $k$-points
required for the integration. The alternative algorithm proposed here
requires, on the other hand, $N_{z}$ matrix inversions due to the calculation
of $\boldsymbol{\underline{\mathrm{\xi}}}^{-1}(z)$, and another $N_{z}$ matrix
diagonalizations in the compact space of the SVD. In practice, the
calculations were well behaved, and we encountered no numerical problems
regarding the diagonalization of Eq. 63.

For the common case of full rank $\boldsymbol{\underline{\mathrm{T}}}$, we can
use the unitarity $\boldsymbol{U}^{\dagger}=\boldsymbol{U}^{-1}$ and
$\boldsymbol{V}^{\dagger}=\boldsymbol{V}^{-1}$ to simplify the formulas
further. We obtain the simple formulas for the local Green’s function

$\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}(z)=\boldsymbol{V}\boldsymbol{\sigma}^{-1/2}\boldsymbol{P}(z)g_{0}\big{(}1/\boldsymbol{d}(z)\big{)}\boldsymbol{P}^{-1}(z)\boldsymbol{\sigma}^{-1/2}\boldsymbol{U}^{\dagger}$
(69)

and its inverse

$\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}^{-1}(z)=\tilde{\boldsymbol{U}}\boldsymbol{P}(z)\frac{1}{g_{0}\big{(}1/\boldsymbol{d}(z)\big{)}}\boldsymbol{P}^{-1}(z)\tilde{\boldsymbol{V}}^{\dagger}.$
(70)

Furthermore, we can directly calculate the matrix diagonalization of

$\boldsymbol{\sigma}^{-1/2}\boldsymbol{U}^{\dagger}\boldsymbol{\underline{\mathrm{\xi}}}(z)\boldsymbol{V\sigma}^{-1/2}=\boldsymbol{P}(z)\boldsymbol{d}^{-1}(z)\boldsymbol{P}^{-1}(z),$
(71)

avoiding the need of the $N_{z}$ matrix inversions for
$\boldsymbol{\underline{\mathrm{\xi}}}(z)$.

##  Appendix C BEB self-consistency equation with renormalized indicator
tensors

With $\alpha\neq\beta$, the high-frequency expansion of the effective medium
yields

$\displaystyle\boldsymbol{\underline{\mathrm{*}}}{S}^{\alpha\beta}(z)$
$\displaystyle=-\epsilon_{(1)}\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\beta}+\order{z^{-1}},$
(72) $\displaystyle\boldsymbol{\underline{\mathrm{*}}}{S}^{\alpha\alpha}(z)$
$\displaystyle=\frac{c^{\alpha}-1}{c^{\alpha}}z+\frac{\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha}-\mu+\overline{c}^{\alpha}\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\alpha}\epsilon_{(1)}}{c^{\alpha}}+\order{z^{-1}},$
(73)

where $\epsilon_{(1)}=\int\mathrm{d}\epsilon\rho(\epsilon)\epsilon$ is the
first moment of the DOS, which vanishes for lattices with a symmetric DOS,
$\boldsymbol{\underline{\mathrm{*}}}{v}$ incorporates the static part of the
self-energy $\Sigma^{\alpha}(z)$, and $\overline{c}^{\alpha}=1-c^{\alpha}$ is
the concentration complement. The diagonal of the effective medium,
$\boldsymbol{\underline{\mathrm{*}}}{S}^{\alpha\alpha}(z)$, has a contribution
which grows linearly in $z$, and the on-site energies are multiplied by the
inverse of the concentration. The origin of this peculiar structure is evident
from Eq. 17 and the definition Eq. 20. Unlike the diagonal elements of a one-
particle Green’s function which behave like $1/z$ for large $z$, the effective
local Green’s function $\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}(z)$
behaves like $\boldsymbol{c}/z$. The definition in terms of the effective
medium, however, has the regular form of ${[\mathds{1}z+\dots]}^{-1}$. This
can be resolved by introducing a renormalized version of the component space.
Instead of the indicator tensor $\boldsymbol{\underline{\mathrm{\eta}}}$, Eq.
10, we use the concentration-scaled indicator tensor

${\boldsymbol{\underline{\mathrm{*}}}{\gamma}}^{\alpha}_{ij}=\sqrt{c^{\alpha}}1_{\mathbb{S}^{\alpha}}(i)\delta_{ij}$
(74)

and the Moore–Penrose inverse Moore (1920); Penrose (1955)
$\boldsymbol{\underline{\mathrm{\gamma}}}^{+}$ (which is in this case the
left-inverse, i.e.,
$\boldsymbol{\underline{\mathrm{\gamma}}}^{+}\boldsymbol{\underline{\mathrm{\gamma}}}=\mathds{1}$)
of its equivalent matrix representation. The components of the Moore–Penrose
inverse read

${\boldsymbol{\underline{\mathrm{*}}}{\gamma}^{+}}^{\leavevmode\nobreak\
\alpha}_{ij}=\begin{cases}\frac{1}{\sqrt{c^{\alpha}}}1_{\mathbb{S}^{\alpha}}(i)\delta_{ij}&\mbox{\quad
if\quad}c^{\alpha}>0,\\\ 0&\mbox{\quad if\quad}c^{\alpha}=0.\end{cases}$ (75)

We can express the projector Eq. 13 with the
$\boldsymbol{\underline{\mathrm{\gamma}}}$ tensor:

$\boldsymbol{\underline{\mathrm{\chi}}}=\boldsymbol{\underline{\mathrm{\gamma\gamma}}}^{+}.$
(76)

For the renormalized BEB formalism, we define the component Green’s function
and the Hamiltonian matrix in terms of
$\boldsymbol{\underline{\mathrm{\gamma}}}$ and the inverse
$\boldsymbol{\underline{\mathrm{\gamma}}}^{+}$ as

$\displaystyle\boldsymbol{\underline{\mathrm{\tilde{G}}}}(z)$
$\displaystyle\coloneqq{(\boldsymbol{\underline{\mathrm{\gamma}}}^{+})}^{\intercal}\boldsymbol{G}(z)\boldsymbol{\underline{\mathrm{\gamma}}}^{+},$
(77) $\displaystyle\boldsymbol{H}$
$\displaystyle\eqqcolon{\boldsymbol{\underline{\mathrm{\gamma}}}}\boldsymbol{\underline{\mathrm{\tilde{H}}}}{\boldsymbol{\underline{\mathrm{\gamma}}}}^{\intercal}.$

Compared to the definitions in Section II.2, the Green’s function and the
Hamiltonian are scaled by the concentration. This can be conveniently
demonstrated in the locator expansion

$\boldsymbol{G}(z)=\boldsymbol{g}(z)+\boldsymbol{g}(z)\boldsymbol{T}\boldsymbol{G}(z),$
(78)

where $\boldsymbol{g}(z)={[\mathds{1}z-\boldsymbol{v}]}^{-1}$ is the locator
and ${(\boldsymbol{T})}_{ij}=t_{ij}$. Sandwiching this equation by
${(\boldsymbol{\underline{\mathrm{\gamma}}}^{+})}^{\intercal}$ and
$\boldsymbol{\underline{\mathrm{\gamma}}}^{+}$, we obtain

${(\boldsymbol{\underline{\mathrm{\gamma}}}^{+})}^{\intercal}\boldsymbol{G}(z)\boldsymbol{\underline{\mathrm{\gamma}}}^{+}={(\boldsymbol{\underline{\mathrm{\gamma}}}^{+})}^{\intercal}\boldsymbol{g}(z)\boldsymbol{\underline{\mathrm{\gamma}}}^{+}\\\
+{(\boldsymbol{\underline{\mathrm{\gamma}}}^{+})}^{\intercal}\boldsymbol{g}(z)\boldsymbol{\underline{\mathrm{\gamma}}}^{+}{\boldsymbol{\underline{\mathrm{\gamma}}}}\boldsymbol{T}\boldsymbol{\underline{\mathrm{\gamma}}}^{\intercal}{(\boldsymbol{\underline{\mathrm{\gamma}}}^{+})}^{\intercal}\boldsymbol{G}(z)\boldsymbol{\underline{\mathrm{\gamma}}}^{+},$
(79)

where we inserted the identity
$\boldsymbol{\underline{\mathrm{\gamma}}}^{+}\boldsymbol{\underline{\mathrm{\gamma}}}=\mathds{1}$.
This can be written in terms of the renormalized component quantities

$\boldsymbol{\underline{\mathrm{\tilde{G}}}}(z)=\boldsymbol{\underline{\mathrm{\tilde{g}}}}(z)+\boldsymbol{\underline{\mathrm{\tilde{g}}}}(z)\boldsymbol{\underline{\mathrm{\tilde{T}}}}\,\boldsymbol{\underline{\mathrm{\tilde{G}}}}(z).$
(80)

Compared to the regular BEB formalism the Green’s functions are scaled with
the reciprocal concentration, and the hopping matrix with the concentration
$\boldsymbol{\underline{\mathrm{\tilde{*}}}}{T}^{\alpha\beta}_{ij}=\sqrt{c^{\alpha}}\boldsymbol{\underline{\mathrm{*}}}{t}^{\alpha\beta}(\absolutevalue{\mathbf{r}_{i}-\mathbf{r}_{j}})\sqrt{c^{\beta}}$.
The renormalized component Green’s function relates now to the one-particle
Green’s function in the following way:

$\operatorname{\mathbb{E}}(\boldsymbol{\underline{\mathrm{\tilde{*}}}}{G}^{\alpha\beta}_{ij})=\begin{cases}{}\operatorname{\mathbb{E}}(G_{ii}(z)|i\mapsto\alpha)\delta^{\alpha\beta}&\mbox{\quad
if\quad}i=j,\\\
\sqrt{c^{\alpha}c^{\beta}}\operatorname{\mathbb{E}}(G_{ii}(z)|i,j\mapsto\alpha,\beta).&\mbox{\quad
if\quad}i\neq j\end{cases}$ (81)

There is no more concentration prefactor for the local Green’s function, which
is the central quantity of the BEB formalism, since it is a local theory. The
renormalized version of the self-consistency equation Eqs. 25 and 26 reads

$0=\boldsymbol{\underline{\mathrm{\tilde{g}}}}_{\text{loc}}^{-1}(z)-\overline{\boldsymbol{\underline{\mathrm{\tilde{g}}}}}^{-1}(z),$
(82)

with the diagonal matrix

$\overline{\boldsymbol{\underline{\mathrm{\tilde{*}}}}{g}}^{\alpha\beta}(z)=\frac{c^{\alpha}\delta^{\alpha\beta}}{{(\boldsymbol{\underline{\mathrm{\tilde{g}}}}_{\text{loc}}^{-1})}^{\alpha\alpha}+\boldsymbol{\underline{\mathrm{\tilde{*}}}}{S}^{\alpha\alpha}-\overline{c}^{\alpha}z+c^{\alpha}(\mu-\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha}-\Sigma^{\alpha})}.$
(83)

With $\alpha\neq\beta$, the high-frequency expansion of the renormalized
effective medium yields

$\displaystyle\boldsymbol{\underline{\mathrm{\tilde{*}}}}{S}^{\alpha\beta}(z)$
$\displaystyle=-\epsilon_{(1)}\sqrt{c^{\alpha}}\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\beta}\sqrt{c^{\beta}}+\order{z^{-1}}$
(84)
$\displaystyle\boldsymbol{\underline{\mathrm{\tilde{*}}}}{S}^{\alpha\alpha}(z)$
$\displaystyle=\boldsymbol{\underline{\mathrm{*}}}{v}^{\alpha}-\mu+\overline{c}^{\alpha}\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\alpha}\epsilon_{(1)}+\order{z^{-1}}.$
(85)

The scaling removes the contribution proportional to $z$, and for a symmetric
DOS the static part is simply the on-site energy of the components.
Furthermore, the static part remains finite with vanishing concentration.

##  Appendix D Hybridization of the Bethe lattice

Specific to the BEB+DMFT scheme is an alloy-component dependent hybridization
function $\Delta^{\alpha}(z)$ computed according to Eq. 28. The hybridization
function describes the hopping process into and out-of the impurity site. The
special form of the lattice Hilbert transform $g_{0}(z)$ for the Bethe lattice

$z-1/g_{0}(z)={({D}/{2})}^{2}g_{0}(z)$ (86)

gives a direct relation between the hybridization function
$\Delta^{\alpha}(z)$ and the effective local Green’s function
$\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}(z)$ Eq. 20. We promote the
hybridization function Eq. 28 to a full hybridization matrix

$\boldsymbol{\underline{\mathrm{\Delta}}}(z)=\mathds{1}z-\boldsymbol{\underline{\mathrm{S}}}(z)-\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}^{-1}(z)\eqqcolon\boldsymbol{\underline{\mathrm{\xi}}}(z)-\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}^{-1}(z),$
(87)

whose diagonal elements are the physical hybridization function
$\boldsymbol{\underline{\mathrm{*}}}{\Delta}^{\alpha\alpha}=\Delta^{\alpha}$.
Using the representation Eq. 68 given in Appendix B, we can apply the identity
Eq. 86:

$\boldsymbol{\underline{\mathrm{\Delta}}}(z)={(D/2)}^{2}\tilde{\boldsymbol{U}}\boldsymbol{P}(z)g_{0}\big{(}1/\boldsymbol{d}(z)\big{)}\boldsymbol{P}^{-1}(z)\tilde{\boldsymbol{V}}^{\dagger}.$
(88)

Likewise, the matrix
$\boldsymbol{\underline{\mathrm{T}}}\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}(z)\boldsymbol{\underline{\mathrm{T}}}$
can be expressed using Eqs. 61, 63 and 67:

$\boldsymbol{\underline{\mathrm{T}}}\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}(z)\boldsymbol{\underline{\mathrm{T}}}=\tilde{\boldsymbol{U}}\boldsymbol{P}(z)g_{0}\big{(}1/\boldsymbol{d}(z)\big{)}\boldsymbol{P}^{-1}(z)\tilde{\boldsymbol{V}}^{\dagger}.$
(89)

Thus, comparing Eqs. 88 and 89, we identify the relation for the Bethe
lattice:

$\boldsymbol{\underline{\mathrm{\Delta}}}(z)={({D}/{2})}^{2}\boldsymbol{\underline{\mathrm{T}}}\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}(z)\boldsymbol{\underline{\mathrm{T}}}.$
(90)

The diagonal elements are the hybridization function, which reads

$\Delta^{\alpha}(z)={(D/2)}^{2}\sum_{\beta}\absolutevalue{\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\beta}}^{2}\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}^{\beta\beta}(z).$
(91)

In the CPA case $\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\beta}=1$ we
evidently recover the component independent hybridization:

$\Delta^{\alpha}(z)={(D/2)}^{2}\sum_{\alpha}\boldsymbol{\underline{\mathrm{g}}}_{\text{loc}}^{\alpha\alpha}(z)={(D/2)}^{2}\operatorname{\mathbb{E}}\big{(}G(z)\big{)}.$
(92)

For the case of a binary alloy shown by Fig. 5, we see that for
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}<1=\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\alpha}$
the hybridization for $\mathsf{A}$ stems mostly from $\mathsf{A}$ sites, while
for
$\boldsymbol{\underline{\mathrm{T}}}*^{\mathsf{A}\mathsf{B}}>1=\boldsymbol{\underline{\mathrm{T}}}*^{\alpha\alpha}$
the main contribution to the hybridization of $\mathsf{A}$ comes from the
$\mathsf{B}$ sites.

## References

  * Lee and Ramakrishnan (1985) P. A. Lee and T. V. Ramakrishnan, Rev. Mod. Phys. 57, 287 (1985).
  * Mott (1990) N. F. Mott,  _Metal-insulator transitions_ , 2nd ed. (Taylor & Francis, London; New York, 1990).
  * Belitz and Kirkpatrick (1994) D. Belitz and T. R. Kirkpatrick, Rev. Mod. Phys. 66, 261 (1994).
  * Mott (1949) N. F. Mott, Proc. Phys. Soc. A 62, 416 (1949).
  * Imada  _et al._ (1998) M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).
  * Anderson (1958) P. W. Anderson, Phys. Rev. 109, 1492 (1958).
  * Evers and Mirlin (2008) F. Evers and A. D. Mirlin, Rev. Mod. Phys. 80, 1355 (2008).
  * Abrahams (2010) E. Abrahams, ed.,  _50 Years Of Anderson Localization_ (World Scientific, Singapore, 2010).
  * Abrahams  _et al._ (2001) E. Abrahams, S. V. Kravchenko, and M. P. Sarachik, Rev. Mod. Phys. 73, 251 (2001).
  * Dobrosavljevic  _et al._ (2012) V. Dobrosavljevic, N. Trivedi, and J. M. Valles Jr., eds.,  _Conductor-Insulator Quantum Phase Transitions_ (Oxford University Press, Oxford, 2012).
  * Nandkishore and Huse (2015) R. Nandkishore and D. A. Huse, Annu. Rev. Condens. Matter Phys. 6, 15 (2015).
  * Abanin  _et al._ (2019) D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).
  * Ziman (1979) J. M. Ziman,  _Models of Disorder: The Theoretical Physics of Homogeneously Disordered System_ (Cambridge University Press, 1979).
  * Soven (1967) P. Soven, Phys. Rev. 156, 809 (1967).
  * Taylor (1967) D. W. Taylor, Phys. Rev. 156, 1017 (1967).
  * Velický  _et al._ (1968) B. Velický, S. Kirkpatrick, and H. Ehrenreich, Phys. Rev. 175, 747 (1968).
  * Elliott  _et al._ (1974) R. J. Elliott, J. A. Krumhansl, and P. L. Leath, Rev. Mod. Phys. 46, 465 (1974).
  * Economou (2006) E. N. Economou,  _Green’s functions in quantum physics_ , Vol. 7 (Springer Science & Business Media, 2006).
  * Gonis (1992) A. Gonis,  _Green functions for ordered and disordered systems_ (North-Holland, Amsterdam, 1992).
  * Metzner and Vollhardt (1989) W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 (1989).
  * Georges and Krauth (1992) A. Georges and W. Krauth, Phys. Rev. Lett. 69, 1240 (1992).
  * Georges  _et al._ (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
  * Kotliar and Vollhardt (2004) G. Kotliar and D. Vollhardt, Phys. Today 57, 53 (2004).
  * Vlaming and Vollhardt (1992) R. Vlaming and D. Vollhardt, Phys. Rev. B 45, 4637 (1992).
  * Janiš and Vollhardt (1992) V. Janiš and D. Vollhardt, Phys. Rev. B 46, 15712 (1992).
  * Milovanović  _et al._ (1989) M. Milovanović, S. Sachdev, and R. N. Bhatt, Phys. Rev. Lett. 63, 82 (1989).
  * Dobrosavljević and Kotliar (1993) V. Dobrosavljević and G. Kotliar, Phys. Rev. Lett. 71, 3218 (1993).
  * Dobrosavljević and Kotliar (1994) V. Dobrosavljević and G. Kotliar, Phys. Rev. B 50, 1430 (1994).
  * Ulmke and Scalettar (1997) M. Ulmke and R. T. Scalettar, Phys. Rev. B 55, 4149 (1997).
  * Denteneer  _et al._ (2001) P. J. H. Denteneer, R. T. Scalettar, and N. Trivedi, Phys. Rev. Lett. 87, 146401 (2001).
  * Miranda and Dobrosavljevic (2012) E. Miranda and V. Dobrosavljevic, in  _Conductor-Insulator Quantum Phase Transitions_, edited by V. Dobrosavljević, N. Trivedi, and J. M. Valles Jr. (Oxford University Press, 2012) p. 161.
  * Blackman  _et al._ (1971) J. A. Blackman, D. M. Esterling, and N. F. Berk, Phys. Rev. B 4, 2412 (1971).
  * Esterling (1975) D. M. Esterling, Phys. Rev. B 12, 1596 (1975).
  * Terletska  _et al._ (2014) H. Terletska, C. E. Ekuma, C. Moore, K.-M. Tam, J. Moreno, and M. Jarrell, Phys. Rev. B 90, 094208 (2014).
  * Zhang  _et al._ (2016) Y. Zhang, R. Nelson, E. Siddiqui, K.-M. Tam, U. Yu, T. Berlijn, W. Ku, N. S. Vidhyadhiraja, J. Moreno, and M. Jarrell, Phys. Rev. B 94, 224208 (2016).
  * Koepernik  _et al._ (1997) K. Koepernik, B. Velický, R. Hayn, and H. Eschrig, Phys. Rev. B 55, 5717 (1997).
  * Koepernik  _et al._ (1998) K. Koepernik, B. Velický, R. Hayn, and H. Eschrig, Phys. Rev. B 58, 6944 (1998).
  * Bauernfeind  _et al._ (2017) D. Bauernfeind, M. Zingl, R. Triebl, M. Aichhorn, and H. G. Evertz, Phys. Rev. X 7, 031013 (2017).
  * Ulmke  _et al._ (1995) M. Ulmke, V. Janiš, and D. Vollhardt, Phys. Rev. B 51, 10411 (1995).
  * Lombardo  _et al._ (2006) P. Lombardo, R. Hayn, and G. I. Japaridze, Phys. Rev. B 74, 085116 (2006).
  * Gonis and Garland (1977) A. Gonis and J. W. Garland, Phys. Rev. B 16, 1495 (1977).
  * Papaconstantopoulos  _et al._ (1989) D. A. Papaconstantopoulos, A. Gonis, and P. M. Laufer, Phys. Rev. B 40, 12196 (1989).
  * Herbig  _et al._ (2017) A. Herbig, R. Heid, and R. Eder, Phys. Rev. B 96, 205128 (2017).
  * Shvaika (2003) A. M. Shvaika, Phys. Rev. B 67, 075101 (2003).
  * Burdin and Fulde (2007) S. Burdin and P. Fulde, Phys. Rev. B 76, 104425 (2007).
  * Karlin and Taylor (1975) S. Karlin and H. E. Taylor, (1975).
  * Vollhardt (2010) D. Vollhardt, in  _Lectures on the Physics of Strongly Correlated Systems XIV_, Vol. 1297, edited by A. Avella and F. Mancini (American Institute of Physics, Melville, 2010) p. 339.
  * Weh and Östlin (2021) A. Weh and A. Östlin, “Derweh/gftools: Gftool 0.9.0 release,” Zenodo http://doi.org/10.5281/zenodo.4744546 (2021).
  * White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * Schollwöck (2005) U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
  * Haegeman  _et al._ (2011) J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Phys. Rev. Lett. 107, 070601 (2011).
  * Lubich  _et al._ (2015) C. Lubich, I. V. Oseledets, and B. Vandereycken, SIAM J. Numer. Anal. 53, 917 (2015).
  * Haegeman  _et al._ (2016) J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Phys. Rev. B 94, 165116 (2016).
  * Bauernfeind and Aichhorn (2020) D. Bauernfeind and M. Aichhorn, SciPost Phys. 8, 24 (2020).
  * Bulla  _et al._ (1998) R. Bulla, A. C. Hewson, and T. Pruschke, J. Phys.: Condens. Matter 10, 8365 (1998).
  * Müller-Hartmann (1973) E. Müller-Hartmann, Solid State Commun. 12, 1269 (1973).
  * Mills and Ratanavararaksa (1978) R. Mills and P. Ratanavararaksa, Phys. Rev. B 18, 5291 (1978).
  * Velický  _et al._ (1968) B. Velický, S. Kirkpatrick, and H. Ehrenreich, Phys. Rev. 175, 747 (1968).
  * Gershgorin (1931) S. A. Gershgorin, Bull. Acad. Sci. URSS 1931, 749 (1931).
  * Byczuk  _et al._ (2004) K. Byczuk, W. Hofstetter, and D. Vollhardt, Phys. Rev. B 69, 045112 (2004).
  * Hubbard (1963) J. Hubbard, Proc. R. Soc. Lond. A 276, 238 (1963).
  * Östlin  _et al._ (2018) A. Östlin, L. Vitos, and L. Chioncel, Phys. Rev. B 98, 235135 (2018).
  * Shiba (1971) H. Shiba, Prog. Theor. Phys. 46, 77 (1971).
  * Higham (2002) N. J. Higham,  _Accuracy and Stability of Numerical Algorithms_, Other Titles in Applied Mathematics (Society for Industrial and Applied Mathematics, 2002).
  * Moore (1920) E. H. Moore, Bull. Am. Math. Soc. 26, 394 (1920).
  * Penrose (1955) R. Penrose, Math. Proc. Camb. Phil. Soc. 51, 406 (1955).