Non Fermi liquid behavior and continuously tunable resistivity exponents in the Anderson-Hubbard model at finite temperature

TL;DR

This study uses advanced computational methods to explore the finite-temperature behavior of the Anderson-Hubbard model, revealing non-Fermi liquid resistivity scaling and a tunable exponent influenced by disorder and interactions.

Contribution

First application of a new computational technique to analyze the Anderson-Hubbard model at finite temperature across various interaction and disorder strengths.

Findings

Resistivity scales as $T^{\alpha}$ with a continuously tunable exponent.

Non-Fermi liquid behavior observed in resistivity scaling.

Disorder and interactions systematically influence the effective disorder.

Abstract

We employ a recently developed computational many-body technique to study for the first time the half-filled Anderson-Hubbard model at finite temperature and arbitrary correlation ($U$) and disorder ($V$) strengths. Interestingly, the narrow zero temperature metallic range induced by disorder from the Mott insulator expands with increasing temperature in a manner resembling a quantum critical point. Our study of the resistivity temperature scaling $T^{\alpha}$ for this metal reveals non Fermi liquid characteristics. Moreover, a continuous dependence of $\alpha$ on $U$ and $V$ from linear to nearly quadratic was observed. We argue that these exotic results arise from a systematic change with $U$ and $V$ of the "effective" disorder, a combination of quenched disorder and intrinsic localized spins.

Equations

$$H=-t\sum _{⟨i,j⟩\sigma }{c}_{i,\sigma }^{†}{c}_{i,\sigma }^{}+\sum _{i}U{n}_{i\uparrow }{n}_{i\downarrow }+\sum _i(V_i-\mu )(n_{i\uparrow }+n_{i\downarrow }),$$

$$\begin{array}{rl}H& =-t\sum _{⟨i,j⟩\sigma }{c}_{i,\sigma }^{†}{c}_{i,\sigma }^{}+\sum _{i}U{n}_{i\uparrow }{n}_{i\downarrow }\\ & +\sum _i(V_i-\mu )(n_{i\uparrow }+n_{i\downarrow }),\end{array}$$

$$n_{i,\uparrow }{n}_{i,\downarrow }$$

$$\begin{array}{rl}{H}_{eff}& =-t\sum _{⟨i,j⟩\sigma }{c}_{i,\sigma }^{†}{c}_{i,\sigma }^{}+\frac{U}{2}\sum _i(⟨n_i⟩n_i-{\text{m}}_i.{\sigma }_i)\\ & +\frac{U}{4}\sum _i({\text{m}}_i^2-⟨n_i{⟩}^2)+\sum _i(V_i-\mu )(n_{i\uparrow }+n_{i\downarrow }).\end{array}$$

$$S(q)=\frac{1}{N^2}\sum _{i,j}{e}^{iq\cdot (r_i-r_j)}⟨{\mathbb{S}}_i\cdot {\mathbb{S}}_j⟩,$$

$$C(|r|)=\frac{1}{P}\sum _{|r|=|i-j|,a}(-1{)}^{|i-j|}⟨S_i^a{S}_j^a⟩.$$

$$\sigma (\omega )=\frac{\pi e^2}{N\hslash a_0}\sum _{\alpha ,\beta }(n_{\alpha }-n_{\beta })\frac{|f_{\alpha \beta }{|}^2}{{ϵ}_{\beta }-{ϵ}_{\alpha }}\delta (\omega -({ϵ}_{\beta }-{ϵ}_{\alpha })).$$

Full text

A hallmark of a conventional Fermi liquid (FL) in good metals is the $T^{2}$ scaling of the resistivity ($\rho$) with temperature $T$. However, deviations from this behavior have been reported in several correlated electronic materials such as heavy fermions Gegenwart et al. (1998); Trovarelli et al. (2000); Julian et al. (1996); Coleman et al. (2001), rare earth nickelates Jaramillo et al. (2014), layered dichalcogenides Lahoud et al. (2014), and cuprates Damascelli et al. (2003); Hussey et al. (2003); Keimer et al. (2015). Various ideas for explaining non Fermi liquid (NFL) states have been proposed. For instance, a $T=0$ quantum critical point (QCP) could induce the linear $\rho\sim T$ scaling in the cuprates Anderson (2006); Keimer et al. (2015); Casey and Anderson (2011). In the NFL observed in the two-dimensional electron gas (2DEG) Bogdanovich and Popović (2002); Jaroszyński et al. (2002); Kravchenko and Sarachik (2004), charge or spin glassy metallic states could provide an alternative Arrachea et al. (2004); Dalidovich and Dobrosavljević (2002); Sachdev and Read (1996); Miranda and Dobrosavljević (2005). In spite of these important efforts, the understanding of NFLs in correlated systems still eludes theorists. Moreover, in heavy fermion experiments a puzzling continuous variation of the $\rho$ vs. $T$ scaling exponent $\alpha$ between 1 and 1.6 was found Gegenwart et al. (1998); Trovarelli et al. (2000); Julian et al. (1996); Coleman et al. (2001). Considering that the microscopic physics of the several NFL material families are quite different, it is a challenge to find a global understanding of NFL states in correlated systems. In particular, we need to identify concrete model Hamiltonian systems that not only support NFL states but also, within a single framework, capture various NFL systematics observed across different material families.

To address these issues, here we study the temperature characteristics of the unconventional metal known to develop at $T=0$ from the competition between strong electron interactions and disorder in the half-filled Anderson-Hubbard model on a square lattice. In the clean limit, the ground state is a Mott insulator (MI) and correlated metals arising from doping MI’s Deng et al. (2013) violate the $\rho\sim T^{2}$ scaling. In the other limit where quenched disorder dominates, single particle states are localized in two dimensions and these disorder-induced Anderson insulators often display variable range hopping behavior Imry (2002).

The surprising $T=0$ intermediate metallic state that results from the combination of correlations and disorder has been studied theoretically using Dynamical Mean Field Theory (DMFT) Dobrosavljevic et al. (2012); Song et al. (2008); Semmler et al. (2010); Tanasković et al. (2003); Aguiar et al. (2008), Quantum Monte Carlo Otsuka and Hatsugai (2000); Srinivasan et al. (2003); Chakraborty et al. (2007), Exact Diagonalization Kotlyar and Das Sarma (2001), and Hartree-Fock Heidarian and Trivedi (2004). Experimental results Kravchenko et al. (1995); Sarma et al. (1998); Nakatsuji et al. (2004); Sanchez-Palencia and Lewenstein (2010); Lahoud et al. (2014) are compatible with the zero temperature calculations. However, the finite temperature understanding of this exotic metal and its scaling is limited and several questions remain. How does a metal that arises from competing Mott and Anderson insulators behave at finite temperatures? What temperature scaling does the resistivity of the ensuing metal display? Is there a dependence of the exponent $\alpha$ on disorder and interaction strengths that can be tuned? Are spin or charge cluster states Miranda and Dobrosavljević (2005) responsible for such scaling behavior? Answers to these open questions are of relevance for experiments and theory alike.

In this publication, we study the half-filled Anderson-Hubbard model at finite temperature using the recently developed Mean Field-Monte Carlo (MF-MC) technique. This approach properly incorporates thermal fluctuations in a mean field theory mfm . Details and benchmarks are in Supplementary Material Sec. I. Using MF-MC, here we establish the disorder-interaction-temperature ($V-U-T$) phase diagram. In particular, we observed a disorder-induced continuous evolution from the Mott to the Anderson insulators with a strange metal in between. Our temperature analysis of this region unveils an intriguing quasi QCP behavior, with a narrow metallic region increasing in width with increasing temperature resembling a quantum phase transition. Through optical conductivity and resistivity calculations, we uncover a striking behavior: by changing $U$ and $V$, $\alpha$ can be tuned (akin to heavy fermions) from linear $\alpha=1$, as in cuprates, to near quadratic. Then disorder and interactions can be used to modify the scaling $\rho\sim T^{\alpha}$.

The model is:

[TABLE]

where the first term is the kinetic energy and the second the standard Hubbard repulsion. $c^{\phantom{\dagger}}_{i\sigma}$ ( $c^{\dagger}_{i\sigma}$) annihilates (creates) an electron at site $i$ with spin $\sigma$. The number operator is $n_{i\sigma}=c^{\dagger}_{i\sigma}c^{\phantom{\dagger}}_{i\sigma}$. The disorder $V_{i}$ at each site is chosen randomly in the interval $[-V,V]$ with uniform probability. The chemical potential $\mu$ is adjusted to achieve half filling globally. In MF-MC, we first Hubbard-Stratonovich decouple the interaction term, by introducing a vector auxiliary $\bm{m_{i}}$ and a scalar $\phi_{i}$ field at every site. The former couples to spin and latter to charge. Dropping the time dependence of the auxiliary fields (Aux. F.) a model with “spin fermion” characteristics arises. The Aux. F. are treated by classical MC that admits thermal fluctuations, and the fermionic sector is solved using Exact Diagonalization. Details of the considerable numerical effort involved are discussed in Supplementary Material Sec I.

1. Phase diagram. Consider the phase diagram shown in Fig. 1 at the representative value $U/t=4$. Various indicators, such as the $(\pi,\pi)$ static magnetic structure factor, density of states (DOS), and optical conductivity $\sigma(\omega)$ were used (Fig. 2). At $T/t=0.005$ the antiferromagnetic (AFM) order is progressively reduced increasing $V/t$ as shown in Fig. 2 (a), and for $V\geq 2.6t$ the signal becomes negligible. At $V=0$, the magnetic order starts at $T/t=0.10$ upon cooling but the system remains insulating above this temperature, as expected, and a paramagnetic insulator (PM-I), is deduced based on the optical conductivity behavior discussed below. Increasing $V/t$, $T_{N}$ initially slightly increases and then reduces with increasing disorder upt . The metal insulator boundary decreases roughly linearly with $V/t$, collapsing to zero at $V/t=2.6$. Panel (b) shows the low-$T$ DOS, $N(\omega)$, for various disorder strengths. We find that the clean-limit Mott gap evolves into a pseudogap gap at $V/t\sim 2$. This pseudogap persists up to $V/t=2.8$ and flattens out for larger $V/t$, with the weight around $\omega=0$ decreasing gradually with disorder as the DOS spreads over a larger energy range due to increasing scattering. Thus, the AF-I, PM-M and CA-I phases in Fig. 1 are gapless.

Figure 2 (d) shows $\omega\sigma(\omega)$ at $T=0$ for different disorder values covering the Mott-Insulator, the strange metal, and the large disorder (CA-I) phases. For $V/t\leq 2.0$, $\sigma(\omega)$ is clearly gapped. In the narrow range $2.0<V/t<2.6$, $\omega\sigma(\omega)$ tends to zero as $\omega\rightarrow 0$ in a nonlinear manner. For $2.6\leq V/t\leq 3.8$, $\omega\sigma(\omega)\rightarrow 0$ linearly i.e. $\sigma(\omega)$ is constant at small $\omega$ indicating a metal. For the CA-I, variable range hopping is expected to provide a $\omega^{3}ln^{3}(I/\omega)$ behavior Imry (2002), where $I$ is a typical energy scale depending on the localization length. Supplementary Material Sec II Fig. 2 shows that the same behavior holds across the finite $T$ insulator to metal transitions as well. Our $T=0$ phase boundaries in Fig. 1 are in excellent agreement with earlier literature Heidarian and Trivedi (2004); Heidarian (2006) and the T = 0 metal is robust against finite size scaling (Supplementary Material, Fig. 3). We now shift the focus to our main contributions at finite $T$.

2. Non Fermi liquid metal. The resistivity extracted from $1/\sigma(\omega)$ at small $\omega$ (see supplementary) is in Fig. 3 (a) for $U/t=4$. There are several important features: (i) $d\rho/dT$ becomes positive at large $T$ for all values of $V/t$; (ii) For the PM-I regime at $V/t=0$ and 1.6, $d\rho/dT$ becomes negative with eventual divergence at the critical AFM temperature. For the metallic ($V/t=2.8$ and $3.4$) and the CA-I ($V/t=4.4$ and $6.6$) phases, $\rho(T)$ saturates at the lowest $T$’s investigated fin . (iii) There are resistivity minima at finite $T$ which coincide with the corresponding location of the peaks in $N(0)$ in Fig. 2 (c) at, e.g., $V/t=3.4$. The NFL nature of the disorder-induced metallic state can be inferred from panel (b) which shows that $\sigma(\omega)$ has a non-Drude form with a peak at finite frequency. This peak is further pushed to higher frequency with increasing $T$, except at very low $T$ when the peak converges to $\omega/t\sim 1$.

Both $\rho(T)$ here, and specific heat ($C_{v}$) in Supplementary Material Fig. 4, show low $T$ deviations from FL behavior consistent with literature on disorder induced NFL’s Miranda and Dobrosavljević (2005), justifying our ′NFL metal′ nomenclature.

Consider now the $\rho$ vs. $T$ behavior for the NFL state. In the metallic regime ($V/t=3.4$), the resistivity minimum occurs at $T/t\sim 0.1$. From Fig. 2 (c) at $V/t=3.4$, the location of $\rho_{min}$ coincides with the peak in $N(0)$ at $T/t\sim 0.055$. This non-monotonic dependence of $N(0)$ on $T$ agrees with DQMC studies of the Hubbard model Paiva et al. (2011) for $V=0$. The initial increase of $N(0)$ is due to thermally induced fluctuations that enhance the DOS weight at $\omega=0$. At high $T$, the scattering of fermions from the Aux. F.’s suppress $N(0)$, and this non-monotonicity is reflected in the metallic-like thermal behavior of $\rho(T)$. In summary, at low $T$ the initial increase of the DOS at the Fermi level forces $d\rho/dT$ to be negative, while at high $T$ this DOS is suppressed again because of the localized spins and $d\rho/dT$ changes sign.

3. Scaling of resistivity. In Fig. 3 (c) we show $\rho(T)$ for combinations of $U/t$ and $V/t$ where the system is a metal over a wide temperature range. The full map of the low temperature metallic region in the $U/t-V/t$ plane is in Fig. 3 (d). The resistivity data is fitted to $\rho(0)+AT^{\alpha}$ for each case to extract $\alpha$ fit (a). For small/intermediate values of $(U/t,~{}V/t)$ (open circles), $\rho(T)$ grows linearly with $T$ in the range analyzed. For larger $U/t$ (and corresponding $V/t$) $\alpha$ increases from $\sim$1.0 to $1.7$ for $U/t=2,~{}V/t=3.4$. As shown in Fig. 3 (d), the metallic window at $T/t=0.005$ occurs roughly around the line $U\sim 1.25V$ RS- (a). The dashed line guides the eye and it envelops the metallic region. The metallic-regime color scale indicates the value of $\alpha$ in the temperature fit of $\rho(T)$. For up to $U/t=2.5$, $\alpha\sim 1$ growing slowly with $U/t$ (the smallest values checked are $U/t=0.5,~{}V/t=0.5$). For larger $U/t,V/t>2.5$, $\alpha$ grows reaching a maximum value close to 2 for $U/t=5$ fit (b). This $\alpha\sim 2$ does not imply a FL but we believe it is just one of the possible transport exponents that occurs in our system in its slow evolution. For even larger $U/t$, within our precision the metallic region closes.

4. Discussion. To better understand the combined effect of disorder and interaction in the metallic phase, in Fig. 4 (a) we show the variance of the local density $\{n_{i}\}$, defined as $\delta n(U,~{}V,~{}T)=\langle\sqrt{\langle n\rangle^{2}-\langle n^{2}\rangle}\rangle$. The outer angular brackets imply averaging over MC samples at fixed $T$, while the inner ones are the quantum average within a single MC sample. There are two sources of disorder that control the variance of $\delta n(U,~{}V,~{}T)$. First is the static disorder ($V$) with variance $\delta V$ and the second are the Aux. F.’s. These Aux. F.’s directly couple to the fermions and indirectly to the disorder through the local fermion occupations (Sec. I of the Supplementary Material). Then, the Aux. F. also provide an inhomogeneous background that at low temperatures follows the presence of intrinsic static disorder. But at temperatures where $k_{B}T\geq\delta V$, the intrinsic disorder ($V$) is unimportant and the Aux. F. configurations become homogeneous in average. Thus, $\delta n(U,~{}V,~{}T)$ remains non zero at low $T$ while it tends to zero at large $T$, on MC sample averaging. This behavior is observed in the inset of Fig. 4 (a). In the main panel we show the same data between $T/t=0.01$ and $1$. We find a systematic increase in $\delta n(U,~{}V,~{}T)$ with increasing $U/t$ and $V/t$. This variance manifests as real-space charge clusters as shown in Fig. 4 (b) to (d) that contains maps of $|\psi(\boldmath{r})|^{2}$ at fixed $T$, MC sample averaged, for states at the Fermi energy. With increasing $U/t$ and $V/t$, the charge clustering and the charge fluctuations magnitude increase systematically following the increase in $\delta n(U,~{}V,~{}T)$. This provides a controlled enhancement in spatially inhomogeneous background from which the fermions scatter RS- (b).

It is known that fermions coupled to classical variables such as disorder Kumar and Majumdar (2005a), adiabatic phonons Blawid et al. (2003) etc. can exhibit charge clustering, metallic glasses, and NFL behavior. In our case, not only the quenched disorder but also the Aux. F. fluctuations play the role of the classical scatterer that give rise to NFL scaling. Such charge clusters and NFL behavior has been experimentally observed in the 2DEG near a $T=0$ quantum critical point Bogdanovich and Popović (2002); Jaroszyński et al. (2002); Kravchenko and Sarachik (2004). Here, within MF-MC, we have found such a “charge cluster metal” in the half-filled Anderson Hubbard model and also observed that their deviation from FL theory can be tuned. This tunability allows us to show that, in a single model Hamiltonian, the resistivity scaling with $T$ can vary between linear to near quadratic, features observed in real NFLs like cuprates and heavy fermions. Our results thus represent progress towards identifying a single model system with NFL behavior transcending many material families.

Acknowledgments. A. Mukherjee acknowledges useful discussions with P. Majumdar, P. Chakraborty, and H. R. Krishnamurthy. N.P. and N.K. were supported by the National Science Foundation, under Grant No. DMR-1404375. E.D. and A. Moreo were supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Materials Sciences and Engineering Division.

I I.   Model and Method

In this section we discuss the Mean field-Monte Carlo (MF-MC) method. The approach properly incorporates thermal fluctuations in a mean field theory and has been shown to reproduce accurately the properties of the clean-limit Hubbard model when compared with Determinantal Quantum Monte Carlo (DQMC) results Mukherjee et al. (2014). MF-MC allows for efficient simulations on large lattices Mukherjee et al. (2015) and has also been used to study frustrated systems Tiwari and Majumdar (2014); Swain et al. (2016), the BCS-BEC crossover Karmakar and Majumdar (2016a); Tarat and Majumdar (2015), imbalanced Fermi gases Karmakar and Majumdar (2016b), and the multiorbital Hubbard model Mukherjee et al. (2016); Swain and Majumdar (2016). Below we present detailed account of the formalism, its numerical implementation, benchmarks against clean system DQMC results, observables calculated and typical parameters (system size, broadening etc.) used in the present calculation.

*1. Effective spin-fermion from the Hubbard model: * As mentioned in the main paper, for our study we use the Anderson-Hubbard model on a square lattice at half filling ($n=1$). The Hamiltonian is:

[TABLE]

where the first term is the kinetic energy and the second term is the standard Hubbard repulsion. $c^{\phantom{\dagger}}_{i\sigma}$ ( $c^{\dagger}_{i\sigma}$) annihilates (creates) an electron at site $i$ with spin $\sigma$. The number operator is $n_{i\sigma}=c^{\dagger}_{i\sigma}c^{\phantom{\dagger}}_{i\sigma}$. The disorder $V_{i}$ at each site is chosen randomly in the interval $[-V,V]$ with a uniform probability. We perform a rotationally invariant decoupling of the interaction term and use the Hubbard-Stratonovich (HS) decomposition to obtain an effective Hamiltonian under the approximations described below. Briefly, the rotational invariant form for the interaction term gives

[TABLE]

Here, the spin operator is ${\mathbb{S}_{i}}=\frac{\hbar}{2}\sum_{\alpha,\beta}c^{\dagger}_{i,\alpha}{\mathbb{\sigma}}^{\phantom{\dagger}}_{\alpha,\beta}c^{\phantom{\dagger}}_{i,\beta}$, $\hbar=1$, $\{\sigma^{x},\sigma^{y},\sigma^{z}\}$ are the Pauli matrices, and $\hat{\Omega}$ is an arbitrary unit vector. To decouple this term we introduce two auxiliary fields $\bm{m}_{i}(\tau)$ and $\phi_{i}(\tau)$ at every site $i$ and every imaginary time $\tau$, with each auxiliary field coupled to the spin and the charge degrees of freedom, respectively. In our recent work Mukherjee et al. (2014) we showed that for $V_{i}=0$, if we consider time independent auxiliary fields and replace the $\phi_{i}$ auxiliary field by its saddle point value, $\phi_{i}=-i\frac{U}{2}\langle n_{i}\rangle$, we obtain the first three terms of the effective Hamiltonian in Eq. 3. Since the disorder term is a one-body interaction it remains unchanged from Eq. 1 to Eq. 3:

[TABLE]

Equation 3 is, thus, a one-body Hamiltonian that belongs to the spin fermion class as it is a functional of the classical auxiliary fields $\{\bm{m}_{i}\}$. We note that although there is no direct coupling between the auxiliary fields $\{\bm{m}_{i}\}$ and the disorder, there is an indirect coupling through the fermion density which is sensitive to the disorder profile and are in turn coupled to the auxiliary fields $\{\bm{m}_{i}\}$.

At $T=0$, $H_{eff}$ reduces to a mean field Hamiltonian while at finite $T$, admitting thermal fluctuations (via employing a Monte Carlo scheme as discussed below) renders the calculation much better than a simple finite $T$ mean-field theory. These conclusions are based on results shown below and in previous literature. Because of the mixing of mean field theory and Monte Carlo, we call the method the “Mean Field-Monte Carlo” (MF-MC) approach.

*2. Treatment of the spin-fermion model: *

Our strategy for solving the resulting spin-fermion model is as follows. We aim to generate equilibrium configurations of $\{\bm{m}_{i}\}$ at a given temperature that minimizes the free energy of the system. The free energy of course contains a quantum component for which we diagonalize $H_{eff}$ for a given configuration $\{\bm{m}_{i}\}$. Since we have broken translation invariance by the quenched disorder, we cannot set $\langle n_{i}\rangle=1$ as could be done in the half filling clean problem. Thus, the values of $\langle n_{i}\rangle$ need to be determined at each iteration. Only then we can diagonalize $H_{eff}$. For this task, we start with an average value of $\langle n_{i}\rangle=1$ at the beginning of our calculation at high $T$. We then choose a random value set of the Aux.F. and diagonalize $H_{eff}$ and extract the associated eigenvalues. We then add the purely classical term $\frac{U}{4}\sum_{i}\textbf{m}_{i}^{2}$ to the fermionic energy. The final free energy is then used within a Metropolis algorithm to update the $\{\bm{m}_{i}\}$ variables. By this repeated process the static auxiliary fields $\{\bm{m}_{i}\}$ are annealed with a classical Monte-Carlo (MC). We repeat the classical MC system sweep 50 times. We then stop updating the Aux.F. and run a self consistency loop for the $\langle n_{i}\rangle$ in the fixed Aux. F. background. Once the $\langle n_{i}\rangle$ have converged for all the sites we stop the self consistency and again run the classical MC for the Aux. F. now for the fixed $\langle n_{i}\rangle$ configuration. The $\langle n_{i}\rangle$ self consistency is repeated after every 50 steps of the classical MC system sweeps during the thermalization portion of the classical MC. This usually consists of 2000 MC system sweeps. During the measurement steps, namely another 2000 MC system sweeps at the same temperature, the density is no longer adjusted, as we found that observables are not sensitive to the self consistency after the first 2000 MC system sweeps. The process then is repeated at successive lower $T$, where we use the $\langle n_{i}\rangle$ and Aux. F. configurations of the last system sweep of the higher $T$ as the initial condition for the lower temperature. This is a computational demanding procedure that requires hundreds of nodes because of the fine grid in $U$, $V$, and $T$ required plus the average over quenched disorder configurations.

To summarize, the static auxiliary fields $\{\bm{m}_{i}\}$ are annealed with a classical Monte-Carlo (MC) while the fermion problem is solved by ED, with some self consistency steps to obtain $\langle n_{i}\rangle$. The chemical potential $\mu$ is freely adjusted to fix the global density to half filling. We performed a total of 4000 MC system sweeps at every temperature from which 200 of the last 2000 were used to calculate the observables discussed below and in the main paper. We begin MC simulation at high temperature, usually ($T/t=1.0$) with a random configuration of auxiliary variables $\{\bm{m}_{i}\}$ and then cool down to a lowest temperature of ($T/t\sim 5\times 10^{-3}$), where up to $100$ intermediate temperatures are used to allow for a gradual evolution in the configurations allowing us to avoid potential metastable states. The observables are averaged over 10 to 20 quenched disorder realizations.

To study the effects of random disorder on a strongly interacting system, it is important to consider large lattices. We use a variant of the above mentioned ED+MC technique known as the traveling cluster approximation (TCA) Kumar and Majumdar (2006) that can be easily parallelized to access large system sizes Mukherjee et al. (2015). We mainly present results for $32^{2}$ systems although some of our results were cross checked on lattices up to $56^{2}$. We used the thermalized configurations of the auxiliary fields and the corresponding eigenvalues and eigenvectors to calculate various observables: density of states $N(\omega$), magnetic structure factor in real and momentum space $S(\textbf{q})$ and $C(|{\mathbb{r}}|)$, respectively, real-space maps of charge densities ($|\psi(\bm{r_{i}},E_{F})|^{2}$)) for states at the Fermi energy, inverse participation ratio ($|\psi(E_{F})|^{4}$) at Fermi energy and optical conductivity ($\sigma(\omega)$) using the Kubo-Greenwood formula. We also extract the resistivity ($\rho$) vs. temperature. We discuss these observables in detail in the next section.

The MF-MC method can be used to study the Hubbard model at arbitrary coupling and temperature and is free of the sign problem. The MF-MC approach has been extensively benchmarked against Determinantal Quantum Monte-Carlo for the Hubbard model Mukherjee et al. (2014) and its $t-t^{\prime}$ variant. Also the method can be parallelized to allow access to very large lattices Mukherjee et al. (2015) up to $250^{2}$ in two dimensions for the one-band Hubbard model.

Finally, establishing the phase diagram on a fine $V$, $U$ and $T$ grid, involving also the averaging of the data over many disorder realizations, demanded a substantial computational effort. Each disorder configuration, for a single set of $U$ and $V$ on a 322 system, requires $\sim 1500$ cpu hours for running the simulation with 100 temperature points. This is multiplied by a factor of 20 for averaging over disorder realizations. Finally, for the full $U$-$V$ phase diagram 400 ($U$, $V$) combinations were used. In addition, calculation of disorder averaged optical conductivity, that scales as $L^{3}$, $L=32^{2}$, for every $(U,V)$ set, also added substantial numerical cost. This was possible due to access to the large cluster computing facility (Newton) at the University of Tennessee, Knoxville.

I.1 Observables

*(i) Static magnetic structure factor: *Information regarding the Néel antiferromagnet order is found by computing the static magnetic structure factor. Here we focus on the wavevector $\bm{q}=(\pi,\pi)$ as that is the expected order at half-filling for any finite value of $U$. The spins ${\mathbb{S}_{i}}$ are constructed from the eigenvectors of the equilibrated configurations. Below the magnetic ordering temperature, the electronic spin degree of freedom are enslaved to magnetic auxiliary variables $\bm{m}_{i}$. Thus, it is interesting to note that even a magnetic structure factor constructed only out of the magnetic auxiliary fields closely resembles the full spin structure factor. The static magnetic structure factor is defined as:

[TABLE]

where ${\mathbb{q}}=\{\pi,\pi\}$ is the wavevector of interest. The normalization $N^{2}$ is equal to $L^{4}$, where $L$ is the linear dimension of a $L^{2}$ system.

(ii) Density of states: We calculate the density of states (DOS), $N(\omega)=\sum_{m}\delta(\omega-\omega_{m})$, where $\omega_{m}$ are the eigenvalues of the fermionic sector and the summation runs up to $2L^{2}$, i.e. the total number of eigenvalues of a $L^{2}$ system with spin. $N(\omega)$ is calculated by implementing the usual Lorentzian representation of the $\delta$ function. The broadening needed to obtain $N(\omega)$ from the Lorentzians is $\sim BW/2L^{2}$, where $BW$ is the fermionic bandwidth at $U=0$. Numerically for the 322 system, the broadening is about $0.004t$. Hundred $N(\omega)$ samples are obtained from the 4000 measurement system sweeps at every temperature. The 200 $N(\omega)$ samples are used to obtain the thermally averaged $\langle N(\omega)\rangle_{T}$ at a fixed temperature. These are further averaged over data obtained from 10-20 independent runs with different random number seeds. We also plot $N(\omega-\mu=0)$ vs temperature and disorder.

(iii) Real space $q=(\pi,\pi)$ magnetic correlations: We also calculate the real space correlation function between the ${\mathbb{S}_{i}}$ vectors. This correlation function is defined as,

[TABLE]

In $C(|{\mathbb{r}}|)$ the summation runs over all P pairs of sites at a distance $|{\mathbb{r}}|$ and is normalized accordingly. The sum over $a$ runs over the three directions $x$, $y$, and $z$. This indicator is used to investigate short range magnetic order in the system in real space.

(iv)Transport calculations: The d.c conductivity $\sigma_{dc}$ is estimated by the Kubo-Greenwood expression Mahan (1990) for the optical conductivity. In a one-electron model system:

[TABLE]

The $f_{\alpha\beta}$ are the matrix elements of the current operator, e.g., $\langle\psi_{\alpha}|j_{x}|\psi_{\beta}\rangle$, and the current operator itself (in the tight-binding model) is given by $j_{x}=ita_{0}e\sum_{i,\sigma}(c^{\dagger}_{{i+a_{0}\hat{x}},\sigma}c^{\phantom{\dagger}}_{i,\sigma}-h.c)$. The $\psi_{\alpha}$ are single-particle eigenstates, and $\epsilon_{\alpha}$ are the corresponding eigenvalues. The $n_{\alpha}=f(\mu-\epsilon_{\alpha})$ are Fermi factors.

We can compute the low-frequency average, $\sigma_{av}(\mu,\Delta\omega,N)=(\Delta\omega)^{-1}\int_{0}^{\Delta\omega}\sigma(\mu,\omega,N)d\omega$, using periodic boundary conditions in all directions. The averaging interval is reduced with increasing $N$, with $\Delta\omega\sim B/N$. Here the constant $B$ is fixed by setting $\Delta\omega=0.008t$ at $N=32^{2}$. Ideally, the d.c. conductivity is finally obtained as $\sigma_{dc}(\mu)={\lim}_{L\rightarrow\infty}{\sigma}_{av}(\mu,B/L,L)$. However, given the extensive numerical cost of our calculation, we simply use the result of $32^{2}$ system as our $\sigma_{dc}(\mu)$. The chemical potential is set to target the required electron density $n$. This approach to d.c. transport calculations has been benchmarked in a previous work Kumar and Majumdar (2005b).

I.2 Half-filled Hubbard model: Results $\&$ benchmarking against DQMC

In this section, we briefly reproduce some pertinent results from earlier publications on the finite temperature properties of the half-filled Hubbard model at zero disorder Mukherjee et al. (2014). We wish to establish a comparison with DQMC for the finite $T$ properties and to discuss the nature of auxiliary field fluctuations with $T$.

1. Zero disorder phase diagram: Figure 1 (a) shows the $U-T$ phase diagram at half filling on a three-dimensional $4^{3}$ lattice from our earlier work Mukherjee et al. (2014). The shaded region (below) solid square symbols is a $(\pi,\pi,\pi)$ antiferromagnetic Mott insulator. The solid squares are obtained using our approach and the open square symbols are DQMC results Staudt et al. (2000). The dashed line with triangles are the $T_{N}$ obtained within finite temperature Hartree-Fock formalism. Clearly, our approach captures the correct non monotonic behavior of $T_{N}$ with $U$ and has semiquantitative agreement with DQMC. This is much better than simple finite temperature Hartree-Fock that gives qualitatively incorrect results. In addition the blue shaded region corresponds to preformed local moments, another feature that a simplistic finite $T$ mean-field theory cannot capture. The gray region to the left of the preformed local moments area is a PM-M. Also as shown in our previous work Mukherjee et al. (2014), the MF-MC approach can capture the two peak specific heat related to moment formation and moment ordering.

2. Auxiliary field fluctuations at zero disorder: As discussed above, the auxiliary fields $\{\bm{m}_{i}\}$ couple linearly to the fermions and they have a quadratic classical cost as well (See Eq. 3). In Fig. 1 (b) we show the system averaged (auxiliary fields) $\{\bm{m}_{i}\}$ as a function of temperature up to $T=100t$ (in the inset) and up to $T=5t$ in the main panel. The data is shown for different $U$’s as indicated in panel (b). For $U=4t$ at low $T$, the $\{\bm{m}_{i}\}$ fields are enslaved to the magnetic moment and thus have a fixed value. With increasing $T$, the magnetic moments weakened and so does $\{\bm{m}_{i}\}$. At $T_{N}=0.12t$, the system is in a paramagnetic metallic state (the gray region in Fig. 1). At higher temperatures with negligible local moments, the value of $\langle|\textbf{m}|\rangle$ is governed by thermal fluctuations. At these temperatures the auxiliary fields behave as harmonic oscillators with a mean amplitude proportional to $\sqrt{T/U}$. Thus, $\langle|\textbf{m}|\rangle$ grows with increasing temperature. With increasing $U$, the dominance of thermal fluctuations in governing $\langle|\textbf{m}|\rangle$ is pushed to progressively higher temperatures that increase with $U$.

This unrestricted growth of the Aux. F. at very high $T$ is unavoidable in the continuous Hubbard-Stratonovich based MF-MC. However, to compare with DQMC, where the magnitude of the Aux. F are fixed to 1, we restrict the comparison to a temperature range for which $\langle|\textbf{m}|\rangle\leq 1$. For all data of focus in our publication, particularly with regards to the temperature scaling of $\rho(T)$, we are in a regime where $\langle|\textbf{m}|\rangle\leq 1$.

In panel (c) we show the clean limit $N(0)$ vs $T$. We find that the temperature where $N(0)$ has a peak coincides with both the $T_{MIT}$ in Fig.1 (a) and the temperature where $\langle|\textbf{m}|\rangle$ start growing as $\sqrt{(T/U)}$. This suggests that auxiliary field amplitude fluctuations are the main scattering effect in the zero disorder problem at high temperature. This non monotonic dependence of $N(0)$ has been found in DQMC calculations as well. In the main text, we show that the same qualitative behavior holds even in the presence of disorder.

II II.  Anderson Hubbard model: determining metal insulator boundary

As mentioned in the main paper, the slope of the resistivity vs temperature is not a good criterion to decide between metals and insulators in strongly correlated systems with disorder. Thus, we have used the behavior of the optical conductivity vs frequency ($\omega/t$) at small frequency as an indicator of metal insulator transitions. As discussed in the main paper, we look for linear behavior of $\omega\sigma(\omega)$ with $\omega$. This ensures a constant conductivity close to $\omega/t=0$. Previously at $T=0$ this criterion has been often employed to decide between metal and insulator Heidarian and Trivedi (2004). In Fig. 2 (d) of the main paper we have shown the corresponding data from our calculations to determine the low temperature (in our case $T/t=0.005$) metal insulator boundary. In this Supplementary Material Fig. 2 (a) to (c), we show typical supporting data that we have used to determine the finite-temperature metal insulator boundary for $U/t=4$. In (a), the $V=0$ case is shown. The straight (black) line is a guide to the eye. We see a clear departure from linear dependence for $T/t=0.08$, while for $T/t=0.2$ the linear behavior is recovered. In (b) we show the data for $V/t=3.4$ which is in the metallic regime. Here we see that, as in the main text Fig 2 (d), for $V/t=3.4$ at $T/t=0.005$ the linear behavior continues at all explored temperatures (we have checked this up to $T/t=2$). The highest $T$ data shown here is $T/t=0.5$. In (c) we show that the deviation from linearity reduces with increasing $T$ for the CA-I phase as well. In the data shown, $T/t=1.2$ is a metal while the lower $T$ system is insulating.

At $T/t=0.005$, we also calculated the inverse participation ratio (IPR) for states at the Fermi energy ($E_{F}$). In practice we averaged over a few states in a small window around the $E_{F}$. This is done for two reasons. Firstly, it allows an alternative independent way to cross check the metal insulator boundaries at low $T$. Secondly, numerically this calculation is much less expensive than calculating $\sigma(\omega)$. Thus, at least at low $T$ we could perform standard finite size scaling to show that IPR extrapolates to zero for the metallic phase and converges to a finite value for both the M-I and the CA-I phases. The data is shown in Fig. 3.

III III.  Specific heat vs temperature

Figure 4 shows the evolution of the specific heat ($C_{v}$) with temperature for the disorder-induced metallic phase. As seen in the main panel, $C_{v}/T$ has a $Log(T/T_{0})$ dependence on temperature over a wide temperature range. This behavior holds true for all the three cases of disorder-induced metal shown in the figure. The specific heat at $T/t<0.3$ diverges at a faster rate. Thus, our calculation also captures another hallmark of a non Fermi liquid, e.g. the divergence of the specific heat with decreasing temperature, in contrast to FL behavior.

Details of the variations in $\gamma_{0}$ and $T_{0}$, in the range of logarithmic scaling of $C_{v}/T$, will be reported elsewhere.