Electric-field-driven resistive switching in dissipative Hubbard model

Electric-field-driven resistive switching in dissipative Hubbard model Jiajun Lil,Camille Aron23,Gabriel Kotliar2 and Jong E.Han! Department of Physics,Sta University of New ersity,New Jersey 08455,USA We study how st correlated electrons on a dissipative lattice evolve from equilibrium under actions,by means of a non-ed dyna n than th -Da Fo arge b st but experi 8 metal into Mott insu lator. We pre non-monotoni field due to an interpla neous metal-insulator mixed state PACS numbers:7127 4a 71.30.+h.7220.Ht anding of 1d driven out equilib ium by lated system.The derived en for the are within the feasible tal ro We study the Hubbard model in a constant and homo sistive transitions. Multiple studies of this phenomenor 10.m ew- toa dramatic drop of resistivity up to5orders of mag anced by coupling the system to a thermostat which The relatively accessible threshold fields E cess of onergy transfer [14,15,21- 包keth ide TB)erm ries.A Landau-Zener hanism [11]see the Coulomb ga ulikely as it predicts a threshold field on the order trostatic pote ential -E imposed on the (-th TB site of 10 V/m. (E= 6 倒wa ated fermion bath [15 Yet The ng Hath ant mod In organic charge -transfer Ho=->(dfladto +H.c.)- ∑daCtaz+Ho ple a the electro-c emica ch nary oxides such as Nio 7]and VO28-10:the electric +god-2+h0 driven current locally heats up the sample which whered are the tight-binding electron creation opera- ng. tors at the -th site with spin o-↑or↓， and the electro ope rators attach librium under an external field,and how we describe the rsion relatiot defined with res ect to the electro quilibrium steadys tates that onsequently emerge static potential-E.g is the overlap between the TB we velope ded ba microscopic theory o cha and the rese ervoirs of length which will be sent to in ty, ical efforts [11-271 we tify in a canonical model of strongly interacting elec will extend this chain into higer dimensional lattice.The resistive switch electric field does not act within each res voirs whos ng ta emonstrat Jo ne ho ating e role is to cha eom th the line me loul vidth)for the e the di ing parameter asr)We work with

Electric-field-driven resistive switching in dissipative Hubbard model Jiajun Li1 , Camille Aron2,3 , Gabriel Kotliar2 and Jong E. Han1 1 Department of Physics, State University of New York at Buffalo, Buffalo, New York 14260, USA 2 Department of Physics, Rutgers University, New Jersey 08854, USA 3 Department of Electrical Engineering, Princeton University, New Jersey 08455, USA. (Dated: May 18, 2015) We study how strongly correlated electrons on a dissipative lattice evolve from equilibrium under a constant electric field, focusing on the extent of the linear regime and hysteretic non-linear effects at higher fields. We access the non-equilibrium steady states, non-perturbatively in both the field and the electronic interactions, by means of a non-equilibrium dynamical mean-field theory in the Coulomb gauge. The linear response regime, limited by Joule heating, breaks down at fields much smaller than the quasi-particle energy scale. For large electronic interactions, strong but experi￾mentally accessible electric fields can induce a resistive switching by driving the strongly correlated metal into a Mott insulator. We predict a non-monotonic upper switching field due to an interplay of particle renormalization and the field-driven temperature. Hysteretic I-V curves suggest that the non-equilibrium current is carried through a spatially inhomogeneous metal-insulator mixed state. PACS numbers: 71.27.+a, 71.30.+h, 72.20.Ht Understanding of solids driven out of equilibrium by external fields [1, 2] has been one of the central goals in condensed matter physics for the past century and is very relevant to nanotechnology applications such as re￾sistive transitions. Multiple studies of this phenomenon have been performed in semiconductors and oxides [3– 10]. In oxides, the application of an electric field can lead to a dramatic drop of resistivity up to 5 orders of mag￾nitude. The relatively accessible threshold fields Eth ∼ 104−6 V/m and the hysteretic I-V curves make them good candidates for the fabrication of novel electronic memories. A Landau-Zener type of mechanism [11] seems unlikely as it predicts a threshold field on the order of 108−9 V/m. In narrow gap chalcogenide Mott in￾sulators, an avalanche breakdown was suggested with Eth ∼ E2.5 gap [3]. Yet, the resistive switchings in other classes of correlated materials do not seem to involve solely electronic mechanisms. In organic charge-transfer complexes, it is believed to occur via the electro-chemical migration of ions [4, 5]. Finally, there are strong indica￾tions that a Joule heating mechanism occurs in some bi￾nary oxides such as NiO [7] and VO2 [8–10]: the electric- field-driven current locally heats up the sample which experiences a temperature-driven resistive switching. These experiments raise basic questions of how a strongly correlated state continuously evolves out of equi￾librium under an external field, and how we describe the non-equilibrium steady states that consequently emerge. We develope a much needed basic microscopic theory of the driven metal-insulator transition. Building on earlier theoretical efforts [11–27] we iden￾tify in a canonical model of strongly interacting elec￾trons a region where electric-field-driven resistive switch￾ing takes place. We demonstrate how Joule heating ef￾fects modify the linear response regime and how, away from the linear regime, the same Joule physics leads to the hysteretic resistive transitions of the strongly corre￾lated system. The derived energy scales for resistive tran￾sitions are orders of magnitude smaller than bare model parameters, within the feasible experimental range. We study the Hubbard model in a constant and homo￾geneous electric field E which induces electric current J. After a transient regime, a non-equilibrium steady state establishes if the power injected in the system, J·E, is bal￾anced by coupling the system to a thermostat which can absorb the excess of energy via heat transfer [14, 15, 21– 24]. The thermostat is modeled by identical fermion reservoirs attached to each tight-binding (TB) sites. In the Coulomb gauge, the electric field amounts in an elec￾trostatic potential −`E imposed on the `-th TB site (` = −∞, · · · , ∞) and on its associated fermion bath [15]. The model is fully consistent with gauge-covariant mod￾els [23]. The non-interacting Hamiltonian reads, Hˆ 0 = −γ X `σ (d † `+1,σd`σ + H.c.) − g √ V X `ασ (d † `σc`ασ + H.c) + X `ασ αc † `ασc`ασ− X `σ `E(d † `σd`σ+ X α c † `ασc`ασ), (1) where d † `σ are the tight-binding electron creation opera￾tors at the `-th site with spin σ =↑ or ↓, and c † `ασ are the corresponding reservoir electron operators attached. α is a continuum index corresponding to the reservoir dispersion relation α defined with respect to the electro￾static potential −`E. g is the overlap between the TB chain and the reservoirs of length V which will be sent to infinity, assuming furthermore that the reservoirs re￾main in equilibrium at bath temperature Tb. Later we will extend this chain into higer dimensional lattice. The electric field does not act within each reservoirs whose role is to extract energy but not electric charge from the system [15]. We use a flat density of states (infinite band￾width) for the reservoir spectra α, and define the damp￾ing parameter as Γ = V −1πg2 P α δ(α). We work with arXiv:1410.0626v2 [cond-mat.str-el] 15 May 2015

2 h-e-kB-a-1 in which e is the electronic charge and a is the latti W)and r solution of the non-interacting model in Eq.(1)has been shown [14.15]to reprodu .The the conventional Boltzmann elerc-d. with the on-site Coulomb interaction parameter as FIG =∑(-）)(du-) tric field E. (2) =4.The lin e in the syn ric limi (DMFT JIG and int ting (mode ty deviates from t 28])to treat the many-body interaction via a self istent local approximation of the self nergies.Note the self-emergy ha oth the idth W the d case.All following energis are in unit of une (w)andS(）=2ifE w)+ (w)with the Fermi-Dirac(FD)distribution fer exp(w/T川 the local retarde means and lo ions(GFs). p a homoe eous non in the Coulomb y(t) sites are quivalent. In the Coulomb gauge, this lead U2(t)2(t).The GFs are updated with this self Gu( the energy using the above Dyson's equations and the proce dure is repeated until convergence is achieved formation from the temporal gauge Wiwe gen to hig her 1 Below,we present the implementation of our DMFT scheme the Coul ts in gauge directly g ou pendicular direction and the above construction of the Dyson's equation can be carried out independen tly pe semi-infinite dissipative Hul ailed divec e S bard chains and its own reservoir)with a sel results of the model in one and three dimension E The ng We first discuss the linear response regime. Withir treated hy means of an imnurity solve the DMFT,the DC in the limit of zere For given self-energy w)r(w+E),the on site Green's functions obey the following Dyson equations the Ku G(w)-1_ w-)-2Fo() )LfED() the spectral function at a given wave-vector k p() k+订 ong a in which 2E< Thi imi to the left and lar to the one used by prange F(w +E)+F and for the (w)is the on-sit I th ain（=1whic electron-phonon interaction.Recent calculat nilar dy e FI(w)-1=w-Eiot(w)-Y2FI(w+E), (5 on is in dependent of the interaction strength i 261 both in which can b solved recursively after (a)one and (b)three-dimension The linear behav and can be 0.00 i (a). or of the impurity GFs,,are constructed using ing pa the W+-W05 mith th g)-1=w+i-2F) (6 factor=l-Re话/a g(w)Ig"(w)22irfeD(w)+2F(w)].(7) With increasing E-field,the contribution atE=U/2is

2 ~ = e = kB = a = 1 in which e is the electronic charge and a is the lattice constant. In the rest of this Letter, we measure energies in units of the full TB bandwidth W = 4γ = 1 (1-d) and W = 12γ = 1 (3-d). The exact solution of the non-interacting model in Eq. (1) has been shown [14, 15] to reproduce the conventional Boltzmann transport theory despite the lack of momentum transfer scattering. The Hubbard model Hˆ = Hˆ 0 + Hˆ 1 is defined with the on-site Coulomb interaction parameter U as Hˆ 1 = U X `  d † `↑ d`↑ − 1 2  d † `↓ d`↓ − 1 2  . (2) Our calculations are in the particle-hole symmetric limit. We use the dynamical mean-field theory (DMFT [16, 28]) to treat the many-body interaction via a self￾consistent local approximation of the self-energies. Note that the self-energy has contributions from both the many-body interaction Hˆ 1 and the coupling to the reser￾voirs: Σr tot(ω) = −iΓ+Σr U (ω) and Σ< tot(ω) = 2iΓfFD(ω)+ Σ < U (ω) with the Fermi-Dirac (FD) distribution fFD(ω) ≡ [1 + exp(ω/Tb)]−1 . Once the local retarded and lesser self-energies are computed, one can access the full re￾tarded and lesser Green’s functions (GFs). Note that in a homogeneous non-equilibrium steady state, all the TB sites are equivalent. In the Coulomb gauge, this leads to G r,< ``0 (ω) = G r,< `+k,`0+k (ω + kE) and similarly for the self-energies [15, 25], as can be derived via a gauge trans￾formation from the temporal gauge. Below, we present the implementation of our DMFT scheme in the Coulomb gauge directly in the steady states. It consists in singling out one TB site – say ` = 0 – (often referred as impurity) and replacing its direct environment (i.e. semi-infinite dissipative Hub￾bard chains and its own reservoir) with a self-consistently determined non-interacting environment (often referred as Weiss “fields”). The local electronic problem is then treated by means of an impurity solver. For given self-energy [Σr,< ` (ω) ≡ Σ r,< U (ω+`E)], the on￾site Green’s functions obey the following Dyson equations G r (ω) −1 = ω − Σ r tot(ω) − γ 2F r tot(ω), (3) G <(ω) = |G r (ω)| 2 [Σ< tot(ω) + γ 2F < tot(ω)], (4) in which γ 2F r,< tot are the total hybridization functions to the left and right semi-infinite chains, F r,< tot (ω) = F r,< + (ω + E) + F r,< − (ω − E). F+(ω) is the on-site re￾tarded GF at the end of the RHS-chain (` = 1) which obeys the self-similar Dyson equation F r +(ω) −1 = ω − Σ r tot(ω) − γ 2F r +(ω + E), (5) which can be solved recursively after more than 500 iter￾ations. F−(ω) corresponds to the GF of the LHS-chain and can be obtained similarly. The non-interacting parts of the impurity GFs, G, are constructed using G r (ω) −1 = ω + iΓ − γ 2F r tot(ω) (6) G <(ω) = |Gr (ω)| 2 [2iΓfFD(ω) + γ 2F < tot(ω)]. (7) 0 0.5 1 1.5 electric-field, E/W 0 0.02 0.04 0.06 0.08 0.1 current, J U = 0 U = 1.5 1-d tight-binding chain, K/W = 0.0625, Tbath/W = 0.00125 0 0.005 0.01 E/W 0 0.005 0.01 J linear-response U/2 - peak U - peak (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 electric-field, E/W 0 0.005 0.01 0.015 current, J U = 0 U = 0.33 U = 0.50 U = 0.67 U = 0.83 3-d tight-binding lattice, Γ/W = 0.0083, Tbath = 0.0 0 0.008 0.016 E/W 0 0.002 0.004 0.006 J (b) FIG. 1: (color online) Electric current (per spin) J vs. elec￾tric field E. (a) 1-d chain with damping Γ = 0.0625W and fermion bath temperature Tb = 0.00125W with the 1-d TB bandwidth W = 4γ. The linear conductance in the small field limit (magnified in the inset) is the same for non-interacting (U = 0) and interacting (U = 1.5W) models. After the con￾ductivity deviates from the linear response behavior, inelastic contributions appear at E = U/2 and E = U. (b) 3-d lat￾tice with Γ = 0.0083W and Tb = 0.00042W with the 3-d TB bandwidth W = 12γ. The main features remain similar to the 1-d case. All following energies are in unit of W, unless otherwise mentioned. The local self-energies are obtained by means of the iterative-perturbation theory (IPT) up to the second-order in the Coulomb parameter U: Σ≷ U (t) = U 2 [G ≷(t)]2G ≶(t). The GFs are updated with this self￾energy using the above Dyson’s equations and the proce￾dure is repeated until convergence is achieved. We generalize the above method to higher dimensions. With the electric-field along the principal axis direction, E = Exˆ, the lattice is translation invariant in the per￾pendicular direction and the above construction of the Dyson’s equation can be carried out independently per each perpendicular momentum vector. See Supplemen￾tary Material for a detailed discussion. Below, we present results of the model in one and three dimensions. We first discuss the linear response regime. Within the DMFT, the DC conductivity in the limit of zero temperature and zero electric field can be obtained via the Kubo formula as σDC ∝ limω→0 P k R dνρk(ν)ρk(ν + ω)[fFD(ν) − fFD(ν + ω)]/ω = P k R dν[ρk(ν)]2 δ(ν) with the spectral function at a given wave-vector k ρk(ν) = −π −1 Im[ν − k + iΓ − Σ r U (ν)]−1 . Therefore, as long as Σ r U (ν) → 0 as ν → 0, T → 0, the DC conductivity is independent of the interaction. This argument is simi￾lar to the one used by Prange and Kadanoff [29] for the electron-phonon interaction. Recent calculations did not have access to the linear response regime [21, 23, 24]. FIG. 1 confirms the validity of the linear response analysis. The initial slope of the J − E relation is in￾dependent of the interaction strength U [26] both in (a) one and (b) three-dimension. The linear behav￾ior deviates at the field Elin ≈ 0.003 in (a), orders of magnitude smaller than the renormalized bandwidth W∗ = zW ≈ 0.5 with the equilibrium renormalization factor z = [1 − Re∂Σ r U (ω)/∂ω] −1 ω=E=Tb=0. With increasing E-field, the contribution at E = U/2 is

m-15 IG.2:(col of Nott-n 12 0.06 with the inte te at ze d lines at U=1 and field becomes ins gat electric fieldof vith thee pr r the electric-field is inc the perat that ence (b)Spectral functi and distribution functio ight rapidly disa near the MITdriven b with qualitative agreement beyond the linear response limit. onotonic cold-hot-cold temperature evo a two-step resonant process which can be viewed as acon r the the lowe o limit in th h27 e of interaction.T he in at the Fom then deviates this behavior outside the nar row linear regime,as discussed below. Hub The scattering rate can be directly related to the elec 1 c() tatCeariowhichtheCoiombicitoracion interacting limit.the linear ductivity can be written the dominant scattering process and is rapidly modified as 0o.DC =22/(VT2+4)[15.nFIG.2(b),w ormula with he ing in Fig.2(a)that the ttering rate from the Coulomb interaction, =-m(w=0), a wide range of the E-field,well beyond the linear regime for different sets of the I collapse onto a scaling e as a function of (E proximate as brium retarded self-en eparture from the linear (83).This estimate is valid away from and the metal-insulator limit,and agrees well hile nega eriodic structures due to the Bloch the dashed lines (U=0)in Fig.1.the NDR here come from strong non-linear scattering enhanced by the Joule of weak dissipation and str 1 comes more dramatic.With the effective temperature Eq.(8),having a singular limit as →0，the ectro dashed lines in Fig.2(a of th DOS at=0.The robust agreement in the self-energies enormalized coherent energy scales causes the system leaves no doubt that the electron scattering is dominated to immediately deviate from the linear response regime, by the Joule heating with T given with Eq.(8)in the preventing itself from overheating.This mechanism,in

3 0 0.005 0.01 0.015 0.02 (E/K) 2 0 0.001 0.002 0.003 0.004 0.005 o U -1 = - Im Y U r(0) U/W = 1.5 U/W = 1 (a) 0 0.05 0.1 0.15 0.2 0.25 electric-field, E/W 0 0.01 0.02 0.03 current, J U = 1 U = 1.5 K/W = 0.025 (b) FIG. 2: (color online) (a) Interacting scattering rate, τ −1 U = −ImΣr U (ω = 0), plotted against (E/Γ)2 . Different colors de￾note different damping Γ = 0.0125, · · · , 0.06 with the interval of 0.0025. For small (E/Γ), the numerical results on the 1-d chain collapse on well-defined lines at U = 1 and 1.5. The dashed lines are predictions based on the equilibrium self￾energy with the temperature replaced by the non-interacting effective temperature Teff given in Eq. (8). The remarkable agreement proves that Joule heating controls the scattering in the small field limit. (b) Comparison of the current and the Drude formula estimate with the total scattering rate Γ+τ −1 U , with qualitative agreement beyond the linear response limit. a two-step resonant process which can be viewed as a con￾sequence of the energy overlap between the lower/upper Hubbard bands of the left/right neighboring sites with the in-gap states present at the Fermi level [27]. The current peak at E = U is due to the direct overlap of the Hubbard bands on neighboring sites [18, 27]. The immediate departure from the linear conductivity at very small fields can be well understood with a Joule heating scenario in which the Coulombic interaction is the dominant scattering process and is rapidly modified by an increasing effective temperature as the field is in￾creased. We first demonstrate this effective temperature effect by showing in Fig. 2(a) that the scattering rates from the Coulomb interaction, τ −1 U = −ImΣr U (ω = 0), for different sets of the damping Γ collapse onto a scaling curve as a function of (E/Γ)2 for small E. This scal￾ing is clearly evocative of the well known T 2 behavior of equilibrium retarded self-energies. In the non-interacting 1-d chain with Tb = 0, the ef￾fective temperature has been obtained in the small field limit as [15, 17] Teff = √ 6 π γ E Γ . (8) Inserting this Teff into the equilibrium perturbative self￾energy [30], we obtain in the weak-U limit τ −1 U = −ImΣr eq(ω = 0, Teff) ≈ π 3 2 A0(0)3 U 2 T 2 eff, (9) which is represented by the dashed lines in Fig. 2(a). Here A0(0) = (π p Γ2 + 4γ 2) −1 is the non-interacting DOS at ω = 0. The robust agreement in the self-energies leaves no doubt that the electron scattering is dominated by the Joule heating with Teff given with Eq. (8) in the 0 0.0002 0.0004 0.0006 0 0.001 0.002 0.003 0.004 0.005 current, J electric-field, E/W U = 1.225 (a) increasing E decreasing E metal-insulator coexistence 0 1 2 -1 -0.5 0 0.5 1 spectral function, A(ω) frequency, ω/W (b) E = 0.0 0.0017 0.0050 -0.02 0 0.02 0 0.5 1 ω/W floc(ω) FIG. 3: (color online) (a) Electric-field driven metal-to￾insulator transition (MIT) in the vicinity of a Mott-insulator at U = 1.225, Γ = 0.00167 and Tb = 0.0025 in a 3-dimensional cubic lattice with electric field in x-direction. The metallic state at zero field becomes insulating at electric field of magni￾tude orders of magnitude smaller than bare energy scales. De￾pending on whether the electric-field is increased or decreased, metal-insulator hysteresis occurs with a window for phase￾coexistence. (b) Spectral function and distribution function floc(ω) with increasing electric-field. The quasi-particle (QP) spectral weight rapidly disappears near the MIT driven by the electric-field, opening an insulating gap. The non-equilibrium energy distribution function indicates that the system under￾goes a highly non-monotonic cold-hot-cold temperature evo￾lution near the MIT. linear response limit in the presence of interaction. Teff then deviates strongly from this behavior outside the nar￾row linear regime, as discussed below. The scattering rate can be directly related to the elec￾tric current via the Drude conductivity J(E) = σDC(E)E with the non-linear DC conductivity σDC(E). In the non￾interacting limit, the linear conductivity can be written as σ0,DC = 2γ 2/(πΓ p Γ2 + 4γ 2) [15]. In FIG. 2(b), we plot the Drude formula with the scattering rate Γ re￾placed by the total scattering Γtot = Γ +τ −1 U . The quali￾tative agreement with the numerical results extends over a wide range of the E-field, well beyond the linear regime. Using Eq. (9), the current at small field can be ap￾proximated as J = σ0,DC E/(1 + E2/E2 lin) with the departure from the linear behavior occuring around (from the condition Γ ≈ τ −1 U at E = Elin), Elin ≈ (8π 2/3)1/2γ 1/2Γ 3/2/U. This estimate is valid away from U = 0 and the metal-insulator limit, and agrees well with FIG. 2(b) [31]. We emphasize that, while negative￾differential-resistance (NDR) behaviors occur typically in periodic structures due to the Bloch oscillations [32] as the dashed lines (U = 0) in Fig. 1, the NDR here comes from strong non-linear scattering enhanced by the Joule heating. In the presence of weak dissipation and strong elec￾tronic interactions, the non-equilibrium evolution be￾comes more dramatic. With the effective temperature, Eq. (8), having a singular limit as Γ → 0, the electron temperature tends to rise very sharply as the field is ap￾plied. This effect, together with a small value of the renormalized coherent energy scales, causes the system to immediately deviate from the linear response regime, preventing itself from overheating. This mechanism, in

U =1.225 can be converted to EMur =107-10 V/m rment 35 implics that the critical field decreases with damping as EMrr VT.Therefore,accounting for the rang of expormcnta el th e mo litative fea tures of the resistive switching phenomenon but a more quantitative ana ysis calls for a better modelling of the ve m for the RS of Fig.4(a) erally reflects that of the) the upturn of the upper critical E-field (black line)in en the metal- ulator transition (black line)with in 上1g 4(a)wit ineeasiseomtcrinthtitie e)with doc rT.P=0.00167t canng reg eMIT.Tm即 ng E,wit white line (dashed line)at about Uer 139 For smal U<U ing field. the QP bandwidth W is larger than Tf and h de en away Q rticle s ected from incoherent sp <T istical property comes strongly no therm eakly [35 as seen in Fig.4(c).This slow increase of T allows a vicinity of a quantum phase transition. larger critical field and leads to modify the state of a system.Indeed,we will show that experimen in Fig.4(d)for U show the QP states spec on-equili trally disconnected incoherent electrons, and a strong in the other to an insulator. non- ior even at E/W th Te Dee In Fig.3(a),we start from a metallic state at U 025 otheelectric-fel We use Even though the calculations performed here are on ho mogeneous lattices,the phase coexistence suggests that. has an extremely narrow linear response window with followed by an NDR behavior. As the have inhomo plex thermodynamic states.The hot metallic regions will oi-belio be oriented in the direction of the held,forming experi nanges in the spe ig BS.The energy distribution function defined as dered films 36].Our calculations of the coexistence of -ImG<(w)/ImG(w),evolves from the FD function at two distinct osteady-tateoth ld h the systems where metal-to-insulator transitions occur with he l emphasize that the energy scale hierarchy increasing temperature Our calculation ignores long nge an nagne an Ein《EMIT《IW (10 extensions to cluster-DMFT would allow a realistic trea observed above differs markedly from that in the quan- ment of the electronic structure and could successfully address the case of VO2 quantum th Fig.4(ab)show the metalinsulato.Our jay Singh.This work has been supported by the Na estimate of the threshold electric field Er0.004 at tional Science Foundation with the Grants No.DMR

4 0 0.002 0.004 0.006 Electric field, E 0.005 0.01 0.015 1 1.1 1.2 1.3 1.4 1.5 Temperature, Tb U/W effective temperature 1 1.1 1.2 1.3 1.4 1.5 U/W 0.005 0.01 0.015 0 1 2 -0.04 -0.02 0 0.02 0.04 ω/W U = 1.38 E = 0.0023 A(ω) floc(ω) T b = 0.0025 (a) Teff , increasing E (b) (c) (d) E = 0.0 metal coexistence region insulator metal insulator metal ins. FIG. 4: (color online) Phase diagram of metal-insulator tran￾sition in a cubic lattice driven by (a) electric field and (b) temperature. The metal-insulator coexistent phase exists be￾tween the metal-to-insulator transition (black line) with in￾creasing E or Tb, and the insulator-to-metal transition (red line) with decreasing E or Tb. Γ = 0.00167. (c) Effective tem￾perature Teff map with increasing E, with the white line for the MIT. The white dashed line becomes the phase boundary with decreasing field. (d) Spectral and distribution functions for strong U beyond the crossover line [black dashed in (c)]. Quasi-particle states are disconnected from incoherent spectra and their statistical property becomes strongly non-thermal. a vicinity of a quantum phase transition, can strongly modify the state of a system. Indeed, we will show that there is a region of the parameters U and E for which the non-equilibrium Dyson’s equations have two distinct solutions, one corresponding to an incoherent metal and the other to an insulator. In Fig. 3(a), we start from a metallic state at U = 1.225, and increase the electric-field from zero. We use the self-consistent solution at a certain E-field as an in￾put to the next E run. As discussed above, the system has an extremely narrow linear response window with Elin ∼ 10−4 , followed by an NDR behavior. As the electric-field is further increased, an electric-field-driven metal-to-insulator RS occurs at EMIT ≈ 0.004. Simi￾lar strong non-linear I-V behavior followed by a resistive transition has been observed in NiO [7]. After gradual changes in the spectral functions in Fig. 3(b), a finite insulating gap opens abruptly after the RS. The local energy distribution function floc(ω), defined as floc(ω) = − 1 2 ImG Ucross with W∗ . Teff, Teff increases with E much weakly [35], as seen in Fig. 4(c). This slow increase of Teff allows a larger critical field and leads to the maximum EMIT(U) near U = Ucross – a prediction which can be experimen￾tally verified. The spectral and distribution functions in Fig. 4(d) for U > Ucross, show the QP states spec￾trally disconnected incoherent electrons, and a strong non-thermal behavior even at E/W∗ ∼ 0.1. To evalu￾ate Teff, fit to a Fermi-Dirac function with Teff has been performed on data satisfying |floc(ω) − 0.5| < 0.25. Even though the calculations performed here are on ho￾mogeneous lattices, the phase coexistence suggests that, under a uniform field, the system can be spatially seg￾regated into metal and insulator regions which in turn have inhomogeneous temperature distribution with com￾plex thermodynamic states. The hot metallic regions will be oriented in the direction of the field, forming experi￾mentally observed current-carrying filaments. The Joule heating scenario has been previously in￾voked in the literature for resistive switching in disor￾dered films [36]. Our calculations of the coexistence of two distinct non-equilibrium steady-state solutions in the framework of a relatively simple quantum mechanical model could be applicable to NiO [7] and CrxV2−xO3 [37] systems where metal-to-insulator transitions occur with increasing temperature. Our calculation ignores long￾range anti-ferromagnetic correlations and does not ad￾dress switching from ordered insulating phases. Further extensions to cluster-DMFT would allow a realistic treat￾ment of the electronic structure and could successfully address the case of VO2. The authors are grateful for helpful discussions with Satoshi Okamoto, Sambandamurthy Ganapathy and Su￾jay Singh. This work has been supported by the Na￾tional Science Foundation with the Grants No. DMR-

oda and Naoto Nagaosa. [21]N. gToaadHAoiPhws.Re.B78,23512 回4em女 and P.Prelovsek ☒G.D.Maha Mny-Particle Physicsrd Ed.Chap. ro G.Kotiar,and C.Weber,Py Rev.Let.10 2 C.Weber,M.Capone,and C.Kotliar (201 Tokur.Science 24.165 T Graf T M.G. R.G.S e,S.H.Chang,J.S.Lee et al.,Appl B网 266 (1970) 3 M.Di Ventra,and Appl.Phys 95.043503(2009) Rev. Lett.91 ces the Joule heat 34D. [1]v andJ.K.Freericks,Phys.Rev.B 71. 51042005 010100) Science 281 I V.Lerne and I.L. Ra次gno4风 .86,779(20140 3(2009). 171 37 W.F 2008 B nann et al 50 20131. ka,and P.Werner,Phys.Rev.Lett

5 0907150, DMR-115181, DMR-1308141, PHYS-1066293 and the hospitality of the Aspen Center for Physics. [1] Leo P. Kadanoff and Gordon Baym, Quantum Statistical Mechanics, Westview Press (1994). [2] G. D. Mahan, Many-Particle Physics 3rd Ed., Chap. 8, Kluwer Academic (2000). [3] V. Guiot, L. Cario, E. Janod, B. Corraze, V. Ta Phuoc, M. Rozenberg, P. Stoliar, T. Cren, and D. Roditchev, Nat Commun 4, 1722 (2013); P. Stoliar, L. Cario, E. Janod, B. Corraze, C. Guillot-Deudon, S. Salmon-Bourmand, V. Guiot, J. Tranchant, and M. Rozenberg, Advanced Ma￾terials 25, 3222 (2013). [4] R. Kumai, Y. Okimoto, Y. Tokura, Science 284, 1645 (1999). [5] J. Jeong, N. Aetukuri, T. Graf, T.D. Schladt, M.G. Samant, and S.S.P. Parkin, Science 339, 1402 (2013). [6] S. Lee, A. Fursina, J.T. Mayo, C.T. Yavuz, V.L. Colvin, R.G. Sumesh Sofin, I.V. Shvets, and D. Natelson, Nat. Mat. 7, 130 (2007). [7] S. B. Lee, S. C. Chae, S. H. Chang, J. S. Lee et al., Appl. Phys. Lett. 93, 252102 (2008). [8] J. Duchene, M. Terraillon, P. Pailly, and G. Adam, Appl. Phys. Lett. 19, 115 (1971). [9] T. Driscoll, H.-T. Kim, B.-G. Chae, M. Di Ventra, and D. N. Basov, Appl. Phys. Lett. 95. 043503 (2009). [10] A. Zimmers, L. Aigouy, M. Mortier, A. Sharoni, Siming Wang, K.G. West, J.G. Ramirez, and I.K. Schuller, Phys. Rev. Lett. 110, 056601 (2013). [11] T. Oka, R. Arita, and H. Aoki, Phys. Rev. Lett. 91, 066406 (2003); T. Oka and H. Aoki, Phys. Rev. B 81, 033103 (2010); T. Oka, Phys. Rev. B 86, 075148 (2012). [12] V. Turkowski and J. K. Freericks, Phys. Rev. B 71, 085104 (2005). [13] J. K. Freericks, Phys. Rev. B 77, 075109 (2008). [14] J. E. Han, Phys. Rev. B 87, 085119 (2013). [15] J. E. Han and J. Li, Phys. Rev. B 88, 075113 (2013). [16] H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and . Werner, Rev. Mod. Phys. 86, 779 (2014). [17] A. Mitra and A. Millis, Phys. Rev. B 77, 220404(R) (2008). [18] A. V. Joura, J. K. Freericks and Th. Pruschke, Phys. Rev. Lett. 101, 196401 (2008). [19] M. Eckstein, T. Oka, and P. Werner, Phys. Rev. Lett. 105, 146404 (2010). [20] Naoyuki Sugimoto, Shigeki Onoda and Naoto Nagaosa, Phys. Rev. B 78, 155104 (2008). [21] N. Tsuji, T. Oka, and H. Aoki, Phys. Rev. B 78, 235124 (2008). [22] M. Mierzejewski, L. Vidmar, J. Bonca, and P. Prelovsek, Phys. Rev. Lett. 106, 196401 (2011); L. Vidmar, J. Bonca, T. Tohyama, and S. Maekawa, ibid. 107, 246404 (2011). [23] C. Aron, G. Kotliar, and C. Weber, Phys. Rev. Lett. 108, 086401 (2012). [24] A. Amaricci, C. Weber, M. Capone, and G. Kotliar, Phys. Rev. B 86, 085110 (2012). [25] Satoshi Okamoto, Phys. Rev. Lett. 101, 116807 (2008). [26] G. Mazza, A. Amaricci, M. Capone, M. Fabrizio, arXiv:1412.6415 (2014). [27] C. Aron, Phys. Rev. B 86, 085127 (2012). [28] A. Georges et al, Rev. Mod. Phys. 68, 13 (1996). [29] R. E. Prange and L. P. Kadanoff, Phys. Rev. 134, A566 (1964). [30] K. Yamada, Prog. Theor. Phys. 54, 316 (1975). [31] In a more realistic model with impurity scattering which becomes more effective than the dissipation at small fields, the critical field has different behavior Elin ∝ τ −1/2 imp Γ/U. [32] Paul A. Lebwohl and Raphael Tsu, J. Appl. Phys. 41, 2664 (1970). [33] At the metal-to-insulator RS, the Teff cools as far as the insulating state is allowed as a meta-stable solution in the equilibrium phase diagram, Fig. 4(b). Therefore the Teff at the E-field immediately after the upper switching field EMIT maps to the red line in Fig. 4(b). After the RS to insulator, the current, while reduced by orders of mag￾nitude, self-consistently produces the Joule heat enough to support the insulating solution. Also see Fig. 4(c). [34] D. Goldhaber-Gordon, Hadas Shtrikman, D. Mahalu, David Abusch-Magder, U. Meirav, and M. A. Kastner, Nature (London) 391, 158 (1998); S. M. Cronenwett, T. H. Oosterkamp, and L. P. Kouwenhoven, Science 281, 540 (1998). [35] See Supplementary Material. [36] B. L. Altshuler, V. E. Kravtsov, I. V. Lerner, and I. L. Aleiner, Phys. Rev. Lett. 102, 176803 (2009). [37] D. B. McWhan, A. Menth, J. P. Remeika, W. F. Brinkman, and T. M. Rice, Phys. Rev. B 7, 1920 (1973); P. Hansmann et al., Phys. Status Solidi B 250, 1251 (2013)

Supplementary Material:Metal-to-insulator phase transition in field-driven electron lattice coupled to dissipative baths Jiaiun Li.Camille Aron.Gabriel Kotliar and Jong E.Han C reservoir n the Coulomb ital a +1 diates theenergyy the o heating hain which I.FORMULATION OF ONE-DIMENSIONAL CHAIN UNDER ELECTRIC-FIELD We model a one-dimensional Hubbard model in the Coulomb gauge,as shown in Fig.S1, with the Hamiltonian H=-∑di1oda+he)+u∑(t-2)(nu-）+∑tathoeCoa 是Esa+he)-EE(4d+∑we 1) where all orbitals on the (-th TB site,and their chemical potential,are shifted by (E.Here we use the unit e=a=1.The fermionic chain reservoirs drain the excess energy of the excited electrons on the main tight-binding (TB)lattice.This Hamiltonian in the Coulomb gauge [2] is equivalent to the temporal gauge [4]and the gauge-covariant form [5].In a long-time limit, we assume we have already reached a well-defined nonequilibrium steady state in the presence of reservoirs.With respect to d,we have two sources for the electronic self-energy,one from the fermion reservoirs and the other from the Coulomb interaction,which we denote as Er and Eu,respectively.Here we make a dynamical mean-field theory (DMFT)assumption that the self-energies are local and identical,except for the energy shift due to the voltage drop along the TB lattice [1]. Ghi(w)=Glic (w+E)and r(w)=(w+(E), (2) oo.oo denotes the lattice site 121.Here the subscript'loc'refers to the local quantity at the central site =0.The Dyson's equation in the steady state for the full retarded Green's 1

Supplementary Material: Metal-to-insulator phase transition in field-driven electron lattice coupled to dissipative baths Jiajun Li, Camille Aron, Gabriel Kotliar and Jong E. Han energy d` d`+1 E reservoir Fig. S1: One-dimensional tight-binding chain driven by a uniform electric-field within the Coulomb gauge. Each orbital along the transport chain is connected to a semi-infinite fermionic chain which dissipates the excess energy accumulated by the Joule heating. I. FORMULATION OF ONE-DIMENSIONAL CHAIN UNDER ELECTRIC-FIELD We model a one-dimensional Hubbard model in the Coulomb gauge, as shown in Fig. S1, with the Hamiltonian H = −γ X `σ (d † `+1,σd`σ + h.c.) + U X `  n`↑ − 1 2  n`↓ − 1 2  + X `ασ αc † `ασc`ασ − g √ L X `ασ (d † `σc`ασ + h.c.) − X `σ `E d † `σd`σ + X α c † `ασc`ασ! , (1) where all orbitals on the `-th TB site, and their chemical potential, are shifted by `E. Here we use the unit e = a = 1. The fermionic chain reservoirs drain the excess energy of the excited electrons on the main tight-binding (TB) lattice. This Hamiltonian in the Coulomb gauge [2] is equivalent to the temporal gauge [4] and the gauge-covariant form [5]. In a long-time limit, we assume we have already reached a well-defined nonequilibrium steady state in the presence of reservoirs. With respect to d` , we have two sources for the electronic self-energy, one from the fermion reservoirs and the other from the Coulomb interaction, which we denote as ΣΓ and ΣU , respectively. Here we make a dynamical mean-field theory (DMFT) assumption that the self-energies are local and identical, except for the energy shift due to the voltage drop along the TB lattice [1], G r,< `` (ω) = G r,< loc (ω + `E) and Σr,< `` (ω) = Σr,< loc (ω + `E), (2) ` = −∞, ∞ denotes the lattice site [2]. Here the subscript ‘loc’ refers to the local quantity at the central site ` = 0. The Dyson’s equation in the steady state for the full retarded Green’s 1 arXiv:1410.0626v2 [cond-mat.str-el] 15 May 2015

function can be expressed in the familiar form as w-E+ir-光o(w-E) 0 Gr(w)-1= 0 w+ir-the刨 w+E+ir-o(+E) =w+eE+ir-EtLoe(w+lE)]6ue +youe-el1. 3 energy is omitted the central site, [g"(w)-I]u=[G"(w)-I]u+EUe(w)608o=[G"(w)-1+Uoelu, (4) with2%oe=diag…,0,0,boc(u),0,0,小 The inversion of the above infinite matrix can be achieved efficiently by a recursive method. We divide the lattice into three parts with the central site=0,the left(=-1,-2..) and right (=1,2,..,)semi-infinite chains.We denote the retarded GF matrixF on the RHS semi-infinite chain as [F()-lle=+E+i-oc(w+EEe+de-, (5) with t,=1,2,...,oo.The local GF at the end of the chain (1)[F(+E)F(w)] can be expressed as a continued fraction 1 F(w+E)= w+E+ir-EUloc(w+E)- w+2E+ir-EB.loe(w+2E)- kw+E+ir-Ei.loc(w+E)-7FI(w+2E)]-1, (6) or FI(w+E)-1=w+E+ir-EUloe(w+E)-7F(w+2E), (7) from the self-similarity of the semi-infinite chain.The recursive relation Eq.(7)is solved mumerically with iteration number M over 500.Practically,we start from an initial GFF(w+ ME)=w+ME +ir-Uloc(w+ME)]-and by Eq.(7)we generate F(w+(M-1)E). We repeat the process Eq.(7)until we reach F(w+E).The LHS GF,F(w-E),can be similarly obtained through F(w -E)1=w-E+ir-EUloc(w -E)-Y2F(w-2E) (⑧) Once we obtain fully convergent GFsFE),the full local GF for the infinite chain can be constructed as Gioc(w)=w+ir-Etloc(w)-[F4(w+E)+F(w-E)]. (9) g"(w)-1=w+ir-72[F(w+E)+F(w-E)]=Gloe(w)-1+Loe(w).(10) 2

function can be expressed in the familiar form as Gr (ω) −1 =          . . . ω − E + iΓ − Σ r U,loc(ω − E) γ 0 γ ω + iΓ − Σ r U,loc(ω) γ 0 γ ω + E + iΓ − Σ r U,loc(ω + E) . . .          = [ω + `E + iΓ − Σ r U,loc(ω + `E)]δ``0 + γδ|`−` 0 |,1 . (3) The Weiss-field Green function G can be expressed similarly except that the interacting self￾energy is omitted at the central site, [G r (ω) −1 ]``0 = [Gr (ω) −1 ]``0 + Σr U,loc(ω)δ`0δ` 00 ≡ [Gr (ω) −1 + Σr U,loc]``0, (4) with Σr U,loc = diag[· · · , 0, 0, Σ r U,loc(ω), 0, 0, · · ·]. The inversion of the above infinite matrix can be achieved efficiently by a recursive method. We divide the lattice into three parts with the central site ` = 0, the left (` = −1, −2, · · · , −∞) and right (` = 1, 2, · · · , ∞) semi-infinite chains. We denote the retarded GF matrix F r + on the RHS semi-infinite chain as [F r +(ω) −1 ]``0 = [ω + `E + iΓ − Σ r U,loc(ω + `E)]δ``0 + γδ|`−` 0 |,1 , (5) with `, `0 = 1, 2, · · · , ∞. The local GF at the end of the chain (` = 1) [F r +(ω + E) ≡ Fr +(ω)11] can be expressed as a continued fraction F r +(ω + E) = 1 ω + E + iΓ − Σr U,loc(ω + E) − γ 2 ω+2E+iΓ−Σr U,loc(ω+2E)− γ2 ··· = [ω + E + iΓ − Σ r U,loc(ω + E) − γ 2F r +(ω + 2E)]−1 , (6) or F r +(ω + E) −1 = ω + E + iΓ − Σ r U,loc(ω + E) − γ 2F r +(ω + 2E), (7) from the self-similarity of the semi-infinite chain. The recursive relation Eq. (7) is solved numerically with iteration number M over 500. Practically, we start from an initial GF F r +(ω+ ME) = [ω + ME + iΓ − Σ r U,loc(ω + ME)]−1 and by Eq. (7) we generate F r +(ω + (M − 1)E). We repeat the process Eq. (7) until we reach F r +(ω + E). The LHS GF, F r −(ω − E), can be similarly obtained through F r −(ω − E) −1 = ω − E + iΓ − Σ r U,loc(ω − E) − γ 2F r −(ω − 2E). (8) Once we obtain fully convergent GFs F r ±(ω ± E), the full local GF for the infinite chain can be constructed as G r loc(ω) −1 = ω + iΓ − Σ r U,loc(ω) − γ 2 [F r +(ω + E) + F r −(ω − E)]. (9) The Weiss-field GF G r (ω), omits the interacting self-energy only on the central site (` = 0) and we have G r (ω) −1 = ω + iΓ − γ 2 [F r +(ω + E) + F r −(ω − E)] = G r loc(ω) −1 + Σr U,loc(ω). (10) 2

Now,we turn to the Dyson's equation for lesser GFs.When the lattice of de is connected to the reservoirs and with finite interaction,its steady-state dynamics is governed by the transport equation G(u)=∑Gp(u)∑tot(a)Gge(u), (11) sum of contributions from teractio For th =0,we use a G()Gi.()i.+()) )Cpoww). (12) For the RHS summation (p 0),one can write the Dyson's equation Go(w)= Go))(w)and similarly for the advanced FGs,and therefore we have the third term as 2Gie()2F(w)(w)F(w). (13) B>0 The summed expression is nothing but the local lesser GF Fw+E)=Fw)within the LHS semi-infinite chain.and we obtain Gfe(w)=IGie(w){()+F(+E)+F5(w-E)]}. (14) Fw+E)can be obtained from (w)following similar steps e+月=Aoo回 =lF(u+E26.:u)+21Fu+E)∑F子.2p(u)6,p(a)年p2(w), =2 where the tilde denotes that the GFs are on the semi-infinite chains of p=2,3,...,oo.Using the self-similarity of the chains=1,.,oo and =2...,o,we have F坠w±E)=F(w±E)t(w±E)+F(w±2E) (15) The lesser Weiss-field GF can be written as G(o)=I9(2{fe(w)+Y2F吓w+E)+F匹w-EJ, (16) with the damping part of the self-energy Er.Using Eq.(14), gro-ar(意品-6 (17) In the main paper,the subscript loc'has been omitted for brevity.Electric current per spin is calculated as Jd doo dba di +dh d-do do)=Re[Gi (0)-G(0)] Re)- (18) =-Re (G()[F(+)-P2(-)]+G"(w)[F(w+-F5(w-]). 3

Now, we turn to the Dyson’s equation for lesser GFs. When the lattice of d` is connected to the reservoirs and with finite interaction, its steady-state dynamics is governed by the transport equation G 0 G r 0p (ω)Σ 0), one can write the Dyson’s equation Gr 0p (ω) = Gr loc(ω)(−γ)F r +,1p (ω) and similarly for the advanced FGs, and therefore we have the third term as γ 2 |G r loc(ω)| 2 X p>0 F r +,1p (ω)Σ< tot,p(ω)F a +,p1 (ω). (13) The summed expression is nothing but the local lesser GF F < + (ω + E) = F < +,11(ω) within the LHS semi-infinite chain, and we obtain G < loc(ω) = |G r loc(ω)| 2 n Σ < tot,loc(ω) + γ 2 [F < + (ω + E) + F < − (ω − E)]o . (14) F < ± (ω ± E) can be obtained from Σ< loc(ω) following similar steps. F < + (ω + E) = X∞ p=1 F r +,1p (ω)Σ< tot,p(ω)F a +,p1 (ω) = |F r +(ω + E)| 2Σ < tot,1 (ω) + γ 2 |F r +(ω + E)| 2 X∞ p=2 F˜r +,2p (ω)Σ< tot,p(ω)F˜a +,p2 (ω), where the tilde denotes that the GFs are on the semi-infinite chains of p = 2, 3, · · · , ∞. Using the self-similarity of the chains ` = 1, · · · , ∞ and ` = 2, · · · , ∞, we have F < ± (ω ± E) = |F r ±(ω ± E)| 2 h Σ < tot,loc(ω ± E) + γ 2F < ± (ω ± 2E) i . (15) The lesser Weiss-field GF can be written as G <(ω) = |Gr (ω)| 2 n Σ < Γ,loc(ω) + γ 2 [F < + (ω + E) + F < − (ω − E)]o , (16) with the damping part of the self-energy ΣΓ. Using Eq. (14), G <(ω) = |Gr (ω)| 2 G< loc(ω) |Gr loc(ω)| 2 − Σ < U,loc(ω) ! . (17) In the main paper, the subscript ‘loc’ has been omitted for brevity. Electric current per spin is calculated as J = i 2 γhd † 1σ d0σ − d † 0σ d1σ + d † 0σ d−1σ − d † −1σ d0σi = γRe[G < 01(t = 0) − G < 0−1 (t = 0)] = γRe Z dω 2π [G < 01(ω) − G < 0−1 (ω)] (18) = −γ 2Re Z dω 2π  G <(ω) F a +(ω + Ω) − F a −(ω − Ω) + G r (ω) F < + (ω + Ω) − F < − (ω − Ω) . 3

To summarize,given the local self-energies o(w),GFs for the semi-infinite chains F are calculated via Eqs.(7,8,15).The retarded GFs are obtained via Eqs.(9,10),and finally the lesser GFs follow via Eqs.(14,17).This procedure,formulated on real-space,corresponds to the k-summation of the impurity GF in equailibrium DMFT formalism. Multi-dimensional lattice under electric-field:For higher dimensional cubic lattice with the field along an axial direction(E=E&),the lattice has translational invariance perpendic- ular to the field,and the problem is block-diagonalized with the transverse wave-vector k We solve the Dyson's equation as above with the k-space (the self-energy o(w)does not have k dependence),and then sum over k to get the local GF.For hypercubic TB lattice the dispersion is ek =-2ycos()+e(k).Then adding e(k)to the on-site energy of the 1-d tight-binding chain and carrying out the 1-d Dyson's equation in the previous section,we obtain the GF Gs(w).By summing over k in the d-1 dimensional Brillouin zone,we get the full local GFs r dd-1k ()=()=deD-(e)G() (19) with the d-1 dimensional DoS Da-1(eL).The Weiss-field GFs are obtained via Eqs.(10,17). II.CROSSOVER OF EFFECTIVE TEMPERATURE NEAR U1.32w The effective temperature Te in Fig.4(c)shows different behavior around Ueross1.32W. where UU the increase of Ter with the E-field is rapid in a similar fashion as that well away from the coexisitence region with smaller U,whereas for U>Ueross the increase of Te is much slower after the linear response limit.Therefore it leads to a large electric-field to reach the metal-insulator transition near the crossover,resulting in a maximum of the EMrr(U) curve in Fig.4(a)and(c).Here,we give a sketch of the different scaling behaviors in the two regimes. From the balance of the Joule heating and the dissipation of energy into the fermion baths, we arrive at the rigorous relation Eq.(39)of Han and Li [2], JE=2rwA(w)[floc(w)-fo(w)]dw, (20) where A(w)is is the ith the effecti (i)In the gime of U Ucroes with W>T To:we can apply the Sommerfeld (21) es with the phenor uation of Altshuler et al [3l with on given by the fermion baths.Aw which ag the dissin response limit,and E2xT Te/W2. (22) 4

To summarize, given the local self-energies Σr, Ucross the increase of Teff is much slower after the linear response limit. Therefore it leads to a large electric-field to reach the metal-insulator transition near the crossover, resulting in a maximum of the EMIT(U) curve in Fig. 4(a) and (c). Here, we give a sketch of the different scaling behaviors in the two regimes. From the balance of the Joule heating and the dissipation of energy into the fermion baths, we arrive at the rigorous relation Eq. (39) of Han and Li [2], JE = 2Γ Z ωA(ω)[floc(ω) − fb(ω)]dω, (20) where A(ω) is the on-site spectral function of the tight-binding lattice, fb(ω) the Fermi-Dirac function of the bath. We approximate the local distribution function floc(ω) as a Fermi-Dirac function with the effective temperature Teff. (i) In the regime of U Teff  Tb, we can apply the Sommerfeld expansion and obtain JE ≈ π 2 3 ΓA(0)(T 2 eff − T 2 b ), (21) which agrees with the phenomenological energy balance equation of Altshuler et al [3] with the dissipation given by the fermion baths. Away from the linear response limit, τ −1 U  Γ and J ∝ γE/τ −1 U and we have E 2 ∝ Γτ −1 U T 2 eff/W2 . (22) 4

From the perturbative self-energy for the scattering rate.as used in the main text. =-1mu=0,Tm)≈交An(03U2Ta (23) we arrive to the scaling relation, E2TU2T/W5,or Ten/W(E/U)2 (24) with To increasing as vE beyond the linear response regime.This scaling agrees well with the numerical calculations. (ii)In the regime of U>Ueroes with Tet W,the scaling relation is quite different Although somewhat exaggerated,we assume TW in the following argument for the sake of simplicity.In such limit,the Sommerfeld expansion is not applicable to both the Eqs.(21) and (23),and the half-QP-bandwidth W*/2 replaces the role of temperature Tef,leading to the approximate relations JE (r/W)w*2 and U2w+2/W3, (25) which de ennence tive switching,and the maximum behavior of 0.005 U-1225,160.002 0.00 Fig.S2:Scaling relation of the critical field for the resistive switching vs.the damping parameter I as the electric-field is increased.The relation relation Err(T)is consistent with the thermal scenario. The above scaling relations can be used to derive EMrr(U)'s dependence on I.At a given U,the nonequilibrium MIT occurs when the effective temperature matches the equilibrium transition temperature Tefr=TgMIT(U).Then Eq.(24),for U<s leads to ErU)/W5,and Er(U)T. (26) The numerical calculations for U=1.225 shown in Fig.S2 confirms the relation. moto, 3 B.L.Altshuler,V.E.Kravtsov,I.V.Lerner,and I.L.Aleiner,Phys.Rev.Lett.102.176803(2009). eber,Phs.Rev.Lett.108,086401(2012). 5

From the perturbative self-energy for the scattering rate, as used in the main text, τ −1 U = −ImΣr eq(ω = 0, Teff) ≈ π 3 2 A0(0)3 U 2 T 2 eff, (23) we arrive to the scaling relation, E 2 ∝ ΓU 2T 4 eff/W5 , or Teff/W ∝ (E/U) 1/2 , (24) with Teff increasing as √ E beyond the linear response regime. This scaling agrees well with the numerical calculations. (ii) In the regime of U > Ucross with Teff Ucross, away from the linear response regime. This slow increase of Teff leads to the enhenced upper switching field EMIT(U) to reach the resistive switching, and the maximum behavior of EMIT(U) results near the crossover value U ≈ Ucross. 0.005 0 damping, Γ 0.004 0.003 0.002 0.001 0 0.001 0.002 MIT electric field, E U = 1.225, Tb = 0.0025 EMIT ~ √Γ Fig. S2: Scaling relation of the critical field for the resistive switching vs. the damping parameter Γ, as the electric-field is increased. The relation relation EMIT(Γ) ∝ √ Γ is consistent with the thermal scenario. The above scaling relations can be used to derive EMIT(U)’s dependence on Γ. At a given U, the nonequilibrium MIT occurs when the effective temperature matches the equilibrium transition temperature Teff = Teq,MIT(U). Then Eq. (24), for U < Ucross leads to EMIT(U) 2 ∝ ΓU 2T 4 eq,MIT/W5 , and EMIT(U) ∝ √ Γ. (26) The numerical calculations for U = 1.225 shown in Fig. S2 confirms the relation. [1] Satoshi Okamoto, Phys. Rev. Lett. 101, 116807 (2008). [2] Jong E. Han and Jiajun Li, Phys. Rev. B 88, 075113 (2013). [3] B. L. Altshuler, V. E. Kravtsov, I. V. Lerner, and I. L. Aleiner, Phys. Rev. Lett. 102, 176803 (2009). [4] V. Turkowski and J. K. Freericks, Phys. Rev. B 71, 085104 (2005). [5] Camille Aron, Gabriel Kotliar, and Cedric Weber, Phys. Rev. Lett. 108, 086401 (2012). 5