Electric-field-driven resistive switching in dissipative Hubbard model.pdf_大学文库

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 experimentally 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 resistive 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 magnitude. 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 insulators, 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 indications that a Joule heating mechanism occurs in some binary 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 equilibrium 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 identify in a canonical model of strongly interacting electrons a region where electric-field-driven resistive switching takes place. We demonstrate how Joule heating effects 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 correlated system. The derived energy scales for resistive transitions are orders of magnitude smaller than bare model parameters, within the feasible experimental range. We study the Hubbard model in a constant and homogeneous 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 balanced 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 electrostatic potential −`E imposed on the `-th TB site (` = −∞, · · · , ∞) and on its associated fermion bath [15]. The model is fully consistent with gauge-covariant models [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 operators 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 electrostatic 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 remain 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 bandwidth) for the reservoir spectra α, and define the damping 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 selfconsistent 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 reservoirs: Σ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 retarded 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 transformation 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 Hubbard 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 onsite 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 retarded 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 iterations. 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. electric 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 conductivity deviates from the linear response behavior, inelastic contributions appear at E = U/2 and E = U. (b) 3-d lattice 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 selfenergy using the above Dyson’s equations and the procedure 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 perpendicular direction and the above construction of the Dyson’s equation can be carried out independently per each perpendicular momentum vector. See Supplementary 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 similar 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 independent of the interaction strength U [26] both in (a) one and (b) three-dimension. The linear behavior 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

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 denote 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 selfenergy 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 consequence 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 increased. 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 scaling 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 effective temperature has been obtained in the small field limit as [15, 17] Teff = √ 6 π γ E Γ . (8) Inserting this Teff into the equilibrium perturbative selfenergy [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-toinsulator 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 magnitude orders of magnitude smaller than bare energy scales. Depending on whether the electric-field is increased or decreased, metal-insulator hysteresis occurs with a window for phasecoexistence. (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 undergoes a highly non-monotonic cold-hot-cold temperature evolution near the MIT. linear response limit in the presence of interaction. Teff then deviates strongly from this behavior outside the narrow linear regime, as discussed below. The scattering rate can be directly related to the electric current via the Drude conductivity J(E) = σDC(E)E with the non-linear DC conductivity σDC(E). In the noninteracting 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 Γ replaced by the total scattering Γtot = Γ +τ −1 U . The qualitative 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 approximated 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 negativedifferential-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 electronic interactions, the non-equilibrium evolution becomes 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 applied. 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 transition in a cubic lattice driven by (a) electric field and (b) temperature. The metal-insulator coexistent phase exists between the metal-to-insulator transition (black line) with increasing E or Tb, and the insulator-to-metal transition (red line) with decreasing E or Tb. Γ = 0.00167. (c) Effective temperature 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 input 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. Similar 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 experimentally verified. The spectral and distribution functions in Fig. 4(d) for U > Ucross, show the QP states spectrally disconnected incoherent electrons, and a strong non-thermal behavior even at E/W∗ ∼ 0.1. To evaluate 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 homogeneous lattices, the phase coexistence suggests that, under a uniform field, the system can be spatially segregated into metal and insulator regions which in turn have inhomogeneous temperature distribution with complex thermodynamic states. The hot metallic regions will be oriented in the direction of the field, forming experimentally observed current-carrying filaments. The Joule heating scenario has been previously invoked in the literature for resistive switching in disordered 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 longrange anti-ferromagnetic correlations and does not address switching from ordered insulating phases. Further extensions to cluster-DMFT would allow a realistic treatment 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 Sujay Singh. This work has been supported by the National Science Foundation with the Grants No. DMR-

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 Materials 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 magnitude, 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 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 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 selfenergy 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

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