# Electrically driven insulator-to-metal transition in a correlated insulator:
Electronic mechanism and thermal description

Manuel I. Díaz<EMAIL_ADDRESS>Laboratoire de Physique, École Normale
Supérieure, CNRS, Université PSL, Sorbonne Université, Université de Paris,
75005 Paris, France  Jong E. Han<EMAIL_ADDRESS>Department of Physics,
State University of New York at Buffalo, Buffalo, New York 14260, USA  Camille
Aron<EMAIL_ADDRESS>Laboratoire de Physique, École Normale Supérieure, CNRS,
Université PSL, Sorbonne Université, Université de Paris, 75005 Paris, France
Institute of Physics, École Polytechnique Fédérale de Lausanne (EPFL), CH-1015
Lausanne, Switzerland

###### Abstract

Motivated by the resistive switchings in transition-metal oxides (TMOs)
induced by a voltage bias, we study the far-from-equilibrium dynamics of an
electric-field-driven strongly-correlated model featuring a first-order
insulator-to-metal transition at equilibrium, namely the dimer-Hubbard model.
We use a non-equilibrium implementation of the dynamical cluster approximation
to access the steady-state spectral and transport properties. We show that the
electric field can drive both metal-to-insulator and insulator-to-metal
transitions. While they proceed by quite distinct mechanisms, specifically
simple heating of the metal versus non-equilibrium effects in the correlated
charge gap, we show that both of these non-equilibrium transitions can be
unified in a single framework once the excitations are accounted for in terms
of an effective temperature. This conceptual advance brings together the two
sides of the long-lasting debate over the origins of the electrically-driven
resistive switching in TMOs.

##  I Introduction

Insulator-to-metal transitions (IMTs) in transition-metal compounds are among
the most abrupt first-order phase transitions in nature: temperature
variations of a few Kelvins can cause dramatic resistivity drops of three to
six orders of magnitude, on timescales of only tens of nanoseconds. For over
half a century, theories and experiments have examined the different
mechanisms behind these transitions, in particular the role of strong
correlation among the electrons [1, 2, 3]. These resistive switchings can also
be triggered electrically, at ambient temperature and pressure, with
relatively small electric fields, $E_{\rm IMT}\sim 10^{2}$ – $10^{3}$ kV/cm
[4, 5]. This makes them excellent candidates for modern electronic switches
with reduced response time and power consumption suitable, in particular, to
neuromorphic computing [6, 7, 8, 9]. They proceed with the creation of
metallic filaments; these heterogeneous dynamics are a manifestation of a
bistable phase where both the insulator and the metal can coexist [10, 11,
12]. This bistability also accounts for the hysteretic behavior by which the
metal-to-insulator transition (MIT) occurs at lower threshold fields $E_{\rm
MIT}\ll E_{\rm IMT}$.

The colossal energy mismatch between the electron-Volt (eV) electronic scales
and the meV scales corresponding to the measured threshold fields has fuelled
a long-standing debate on whether electrically-driven resistive switchings
follow the same mechanisms as in equilibrium, with a temperature increase in
the filament driven by Joule heating effects, or whether intrinsically non-
equilibrium electronic effects are at play [13, 14, 15]. The case of vanadium
oxides has been particularly scrutinized as their thermally-driven IMT is
often presented as archetypal of the physics of strongly-correlated electrons.
Experimental evidence seems to indicate that the resistive switching in
pristine Vanadium dioxide (VO2) follows the Joule-heating scenario, whereas
the one in Vanadium sesquioxide (V2O3) is thought to result from non-thermal
effects [16, 17, 18, 19, 20, 21, 22, 23]. However, the controversy is still
fierce and elements of both scenarios could participate towards the
transition.

From a theoretical perspective, it is a great challenge to address the
electrically-driven resistive switchings in these correlated insulators. On
the one hand, we often start with a partial understanding of the _equilibrium_
IMT, and on the other hand, it requires solving far-from-equilibrium dynamics
of complex open quantum many-body systems.

Conceptual progress has been made in the context of elementary correlated band
insulators treated by means of static mean-field techniques, allowing access
to analytical solutions [24, 25]. The ubiquity of non-equilibrium first-order
switchings was unveiled by showing that even when the equilibrium system
features a continuous phase transition, such as in V3O5, a finite electric
field can open a bistable regime between the metal and the insulator, and turn
the resistive switching into a non-equilibrium first-order transition.
Moreover, it was proposed to unify the two competing scenarios of resistive
switching in a single framework. The pivotal concept that has been put forward
is the notion of effective temperature, $T_{\rm eff}$, which quantifies the
number of electronic excitations irrespective of their thermal or purely non-
equilibrium origin [26]. In that view, IMT and MIT may proceed by distinct
mechanisms, at much different threshold fields, but one should recover the
equilibrium phase diagram once the non-equilibrium phase diagram is
parameterized in terms of $T_{\rm eff}$.

The last decade has also seen methodological progress with the development of
non-equilibrium formulations of Dynamical Mean-Field Theory (DMFT) allowing to
treat both the finite electronic interaction in driven-dissipative correlated
lattices and their distance to equilibrium in a non-perturbative fashion [27,
28, 29, 30]. In particular, non-equilibrium steady-state (NESS)
implementations offer to bypass time-dependent transient regimes that are
otherwise numerically intensive to resolve [31, 32, 33, 34, 35].

Several scenarios of non-equilibrium resistive switching have been studied in
the context of the Hubbard model, such as the dielectric breakdown of the Mott
insulator by an intense electric field [36, 37], the role of a finite-sized
sample [38], the role of filament formation and heterogeneities [39], the role
of the nature of the dissipative environment [40, 41, 42], the role of Hund’s
coupling [43], photo-induced resistive switching [44], etc. Notably, all the
DMFT studies were performed with models featuring a thermally-driven MIT at
equilibrium rather than the more common IMT, excluding e.g. the case of VO2.

In this manuscript, we harness those techniques to study the electric-field-
driven resistive switching of a correlated insulator that features a
thermally-driven first-order IMT. Our particular model is inspired by the
physics of VO2, with strong intra-orbital correlation and dimerization. We
devise a non-equilibrium version of a two-site cluster-DMFT approach to solve
for its NESS. We study the non-equilibrium spectral and transport properties
as the electric field is varied. We identify the mechanisms behind both MITs
and IMTs, discussing in particular how the self-heating effects are strongly
suppressed in the insulator, giving way to truly non-equilibrium electronic
mechanisms. Importantly, we establish that the effective temperature is a
fruitful concept that serves as a Rosetta stone between the equilibrium phase
diagram and its non-equilibrium counterpart.

The paper is organized as follows. The driven-dissipative model is introduced
in Sect. II, and its equilibrium phase diagram is discussed along Fig. 1.
Section III is devoted to describing the non-equilibrium steady-state DMFT
methodology we use to solve the Schwinger-Keldysh Green’s functions. The
results are presented in Sect. IV. In particular, the non-equilibrium phase
diagram parameterized in terms of the effective temperature $T_{\rm eff}$ is
displayed in Fig. 10. We discuss the outlooks and conclude in Section V.

##  II Model

$t$ | $t^{\perp}$ | $U$ | $\Gamma$ |  $T_{\rm IMT}$(VO2) |  $E_{\rm IMT}$(VO2)
---|---|---|---|---|---
0.25 | 0.3 | 2.5 | $10^{-3}$ | $2.9\times 10^{-2}$ |  $10^{-3}$ – $10^{-2}$
Table 1:  Typical values (in units of eV) of the parameters of the Dimer-
Hubbard model (DHM) in Eq. (1), which are pertinent to a description of
vanadium dioxide (VO2) [45]. $E_{\rm IMT}$ values were extracted from Refs.
[17, 20] and converted to eV using the potential drop between neighboring
sites, $|q|Ea$, and the interatomic distance $a=4.5$ Å [46].

We consider a correlated insulator driven out of equilibrium by a DC electric
field and coupled to a heat sink. The total Hamiltonian of the electronic
many-body system, its non-equilibrium drive, and its dissipative environment
reads

$\displaystyle H=H_{\mathrm{DHM}}+H_{E}+H_{\Gamma}\,.$ (1)

The electronic system is inspired by the physics of VO2. It is given by the
dimer-Hubbard model (DHM) which is a variation on the standard Hubbard model
incorporating two ingredients that are now firmly established to be key to the
equilibrium IMT in VO2: strong electronic interaction and the dimerization of
the vanadium atoms [47, 48]. The DHM was first proposed in the context of the
thermally-driven IMT in VO2 in Ref. [49] and was recently studied in Refs.
[50, 45]. We should note that a quantitative description of the thermally-
driven IMT in VO2 requires more sophisticated modeling [47, 51, 52, 53, 54]
whose out-of-equilibrium treatment is simply out of reach for current state-
of-the-art non-equilibrium methodologies. The DHM Hamiltonian is given by (we
set $\hbar=1$)

$\displaystyle H_{\mathrm{DHM}}=$ $\displaystyle-t\\!\\!\sum_{\langle
ij\rangle
a\sigma}\\!(c_{ia\sigma}^{\dagger}c_{ja\sigma}+\mathrm{H.c.})-\frac{U}{2}\\!\sum_{ia\sigma}c^{\dagger}_{ia\sigma}c_{ia\sigma}$
(2)
$\displaystyle\\!\\!+U\\!\sum_{ia}c^{\dagger}_{ia\uparrow}c_{ia\uparrow}c^{\dagger}_{ia\downarrow}c_{ia\downarrow}-t^{\perp}\\!\sum_{i\sigma}(c^{\dagger}_{i1\sigma}c_{i2\sigma}+\mathrm{H.c.}).$

The $c_{ia\sigma}^{\dagger}$ ($c_{ia\sigma}$) are the creation (annihilation)
operators of electrons at the site $i$ of a two-dimensional square lattice, in
the orbital $a=1,2$, with spin $\sigma=\;\uparrow,\downarrow$. In two
dimensions, the DHM is also referred to as the bilayer-Hubbard model [55, 56,
57, 58]. $t>0$ sets the hopping amplitude between nearest neighboring sites.
$U>0$ sets the strength of the local electronic interaction originating from
the intra-orbital Coulombic repulsion. $t^{\perp}>0$ couples the two orbitals
within each dimer site. At $t^{\perp}=0$, one recovers two uncoupled copies of
the standard single-orbital Hubbard model (SOHM). A finite $t^{\perp}$ tends
to dimerize the overall system by favoring the formation of local spin
singlets at each site. This model is symmetric under orbital and spin
permutations, $1\leftrightarrow 2$ and $\uparrow\;\leftrightarrow\;\downarrow$
respectively. Furthermore, the presence of the quadratic term in $U/2$ ensures
the particle-hole symmetry of the DHM. Below, we assume that these symmetries
are not spontaneously broken.

The electrons are driven out of equilibrium by a constant and uniform external
electric field $\bm{E}$ which we choose to be aligned along the $x$-axis of
the square lattice: $\bm{E}=E\bm{u}_{x}$. After a transient regime, this is
expected to generate a steady electric current $\bm{J}=J\bm{u}_{x}$. We work
in the Coulomb gauge where the electric field enters the problem as a ramp
potential:

$\displaystyle
H_{E}=-|q|E\sum_{ia\sigma}x_{i}c^{\dagger}_{ia\sigma}c_{ia\sigma}\,,$ (3)

where $x_{i}$ is the spatial coordinate along the $x$-axis of the site $i$ and
$-|q|$ is the charge of the electron. We set the interatomic distance $a=1$
and $|q|=1$.

Figure 1:  Equilibrium phase diagram in the $U-T$ plane of (a) the
paramagnetic single-orbital Hubbard model (SOHM, $t^{\perp}=0$) ; and (b) the
Dimer-Hubbard model (DHM, $t^{\perp}=0.3$) in Eq. (1). The solid first-order
transition lines delimit the extent of the coexistence region between the
metal and the insulator. The dashed lines are the spinodal curves. The SOHM
paramagnet features a coexistence region whose left-leaning tilt supports
temperature-driven metal-to-insulator transitions (MIT). The opposite right-
leaning tilt of the DHM brings temperature-driven insulator-to-metal
transitions (IMT). ($\Gamma=7.5\times 10^{-4}$).

To allow for non-trivial steady states, it is crucial to include an energy
dissipation channel. Otherwise, the work performed on the electronic system
per unit of time and volume, $\mathcal{W}=\bm{J}\cdot\bm{E}>0$, would lead to
a trivial infinite-temperature steady state. In practice, we couple the
orbitals to independent reservoirs of non-interacting electrons [59, 60, 33,
61, 62, 63, 64, 42]

$\displaystyle
H_{\Gamma}=\gamma\sum_{ia\sigma\,l}c_{ia\sigma}^{\dagger}b_{ia\sigma
l}+\sum_{ia\sigma\,l}(\epsilon_{l}-|q|Ex_{i})b_{ia\sigma
l}^{\dagger}b_{ia\sigma l}\,,$ (4)

where $\gamma$ sets the hopping amplitude to the reservoirs and $\epsilon_{l}$
are the many energy levels of the reservoirs. We impose that the fermionic
reservoirs act as a good thermal bath with no back-action from the system:
independently of the state of the system, they remain in equilibrium at the
temperature $T$ and chemical potential $\mu-|q|Ex_{i}$. In practice, we take
reservoirs with a flat density of states,
$\sum_{l}\delta(\omega-\epsilon_{l})=\rho_{\rm b}$, providing dissipation
channels at all energies. The electronic reservoirs are controlling the
electronic filling of the DHM: we set $\mu=0$ to work at half-filling, i.e.
with one electron per orbital on average. Owing to their non-interacting
nature, the reservoir degrees of freedom can be explicitly integrated out.
Eventually, this simple form of dissipation enters the problem via two energy
scales: the bath temperature $T$ (we set $k_{\rm B}=1$) and
$\Gamma=\pi\gamma^{2}\rho_{\rm b}$ which sets the rate at which electrons are
exchanged with the environment. We stress that this choice of a fermionic
environment should be seen as a simple heuristic way to single-handedly
account for the different channels of energy dissipation present in actual
physical systems.

We work in units of $4t$, which corresponds to the half-bandwidth of the
system at $t^{\perp}=U=0$. Incidentally, this roughly corresponds to working
in units of eV. In Table 1 we collect the typical values of the model
parameters that are experimentally pertinent and that we use throughout this
manuscript. The dissipative rate $\Gamma$ was estimated using experimental
data, see App. A for more details.

Figure 2:  Metal/insulator coexistence: equilibrium local density of states
(DOS) of the DHM in Eq. (1) for a set of parameters corresponding to the
coexistence region. In red: the metallic solution. In black: the insulating
solution, characterized by a Mott pseudogap controlled by the electronic
interaction $U$ which features a cleaner gap $\Delta$ controlled by
$t^{\perp}$, see the close up in Fig. 3 (a). The parameters are $U=3$,
$t^{\perp}=0.3$, $T=2.5\times 10^{-3}$, and $\Gamma=7.5\times 10^{-4}$.

The equilibrium phase diagram (at $E=0$ and $\Gamma=0$) of the DHM has been
mapped out in Refs. [65, 66, 57, 50, 45] by means of dynamical mean-field
theory methods. At $U=0$, the non-interacting Hamiltonian can be diagonalized
exactly. For small values of $t^{\perp}$, the non-interacting system is a
metal. For large $t^{\perp}>4t$, the system experiences a Peierls transition
from a metal to an insulator with the opening of band gap
$\Delta^{\perp}=2t^{\perp}-8t$. In this work, following the typical parameter
values collected in Table 1, we stay below this transition.

In Figure 1, we report the $U-T$ equilibrium phase diagrams of the SOHM
($t^{\perp}=0$) and of the DHM ($0<t^{\perp}<4t$) computed by means of single-
site DMFT utilizing a simple impurity solver (namely Iterated Perturbation
Theory), where short-range correlations are neglected (see details below).
Directly relevant to resistive switching is the metal-to-insulator transition
(MIT) driven by a finite electronic repulsion $U$. The corresponding Mott
transition proceeds as a first-order phase transition at low temperatures,
with the opening of a correlated gap in the density of states. $U_{c1}$ marks
the insulator-to-metal transition (IMT) below which the system is a stable
metal, while $U_{c2}>U_{c1}$ marks the MIT transition above which the system
is a stable insulator. In the intermediate region, $U_{c1}<U<U_{c2}$, both the
metal and the insulator may coexist. This is illustrated in Fig. 2 where we
display the local density of states of the coexisting metal and insulator.
This coexistence region subsists until a critical temperature $T^{*}$, where
$U_{c1}=U_{c2}=U^{*}$. For temperatures $T>T^{*}$, the electronic interaction
drives a smooth crossover between a (bad) metal and a (dirty) insulator.

Importantly, within single-site DMFT, the shape and the nature of the
coexistence region depends crucially on whether $t^{\perp}=0$ or $t^{\perp}$
is finite. The resulting qualitative picture has been confirmed using a more
sophisticated single-site solver [50]. At $t^{\perp}=0$, the overall
triangular shape of the coexistence region in the $U$-$T$ plane allows for
metal-to-insulator transitions: a metal prepared in the metastable region can
be heated up to experience an MIT. At finite $t^{\perp}$, the overall tilt of
the coexistence region is reversed: an insulator prepared in the metastable
region can be heated up so as to cross the IMT line. This presence of a
thermally-driven first-order IMT corresponds to the equilibrium physics of
VO2.

The difference in the overall tilt of the coexistence region between the
(half-filled) paramagnetic SOHM at $t^{\perp}=0$ and the DHM at $t^{\perp}>0$
is related to the different nature of their respective insulating ground state
111This difference in the nature of the insulating state is analogous to the
difference between the insulating states of the SOHM computed by means of
single-site DMFT versus cluster-DMFT. In the latter, short-range correlations
are kept such as to allow for a lower entropy insulating ground state. [68].
In the SOHM, neglecting short-range correlations, the paramagnetic insulator
consists of a collection of independent spins at every site. This is
responsible for a macroscopic zero-temperature residual entropy (of value
$\log{2}$ per site). Comparing the free energies as the temperature $T$ is
increased, this favors the insulating solution against the metallic solution
for which the Fermi-liquid entropy only grows linearly in $T$. In the DHM, the
insulating ground state corresponds to a collection of pairs of spins locked
in spin singlets at every site. This corresponds to a zero entropy state, as
for the metal. However, as the temperature $T$ is increased, the entropy of
the Mott insulator increases slower than that of the metal because excitations
require tunneling through the charge gap. Hence, the free energy of the metal
is favored and the DHM undergoes a thermally-driven IMT.

Figure 3:  (a) Close-up of the gapped local density of states (DOS) of the
SOHM ($t^{\perp}=0$) and of the DHM ($t^{\perp}=0.3$) at equilibrium for the
same $U=3.5$. The level of the local DOS at the Fermi energy of the DHM is 25
times smaller than for the SOHM. (b) Width of the gap $\Delta$ as a function
of $t^{\perp}$. The solid lines are linear fits. The parameters are
$T=2.5\times 10^{-3}$ and $\Gamma=7.5\times 10^{-4}$.

This difference in the nature of the insulating state between the paramagnetic
SOHM and the DHM has also important consequences on their spectral features.
While both insulators exhibit a residual in-gap density of states due to the
dissipation $\Gamma>0$, the one of the DHM is smaller by at least one order of
magnitude. This difference can be understood as follows. The in-gap states of
the SOHM are predominantly caused by the dissipative leaking of the Hubbard
lobes inside the Mott gap: $\rho(\omega=0)\sim\Gamma/U^{2}$ [37]. In the DHM
case, the local spin singlets are bound tighter as $t^{\perp}$ increases, and
this results in the opening of an additional gap of width $\Delta$ within the
correlated gap. This gap is illustrated in Fig. 3 (a), and we report the
linear dependence of $\Delta$ with respect to $t^{\perp}$ in Fig. 3 (b). The
residual in-gap density of states is now controlled by the local coupling of
the singlets to the dissipative reservoirs,
$\rho(\omega=0)\sim{t^{\perp}}^{2}\Gamma/U^{4}$, which makes the gap much
cleaner than in the SOHM case. See App. C for more details.

In turn, the extremely low density of states of the DHM insulator makes it
impervious to perturbations around the Fermi energy, such as the micro-
currents driven by an external electric field. Joule heating is therefore
expected to be much less effective than in the case of the SOHM insulator.

##  III Methods

We solve the electric-field-driven many-body problem directly in the non-
equilibrium steady state (NESS), bypassing the transient dynamics. In
practice, the non-equilibrium Green’s functions are obtained by solving a
Schwinger-Keldysh formulation of the Schwinger-Dyson’s equations, assuming
that the dynamics reach a well-defined NESS with time and space translation
symmetries. Such a functional approach allows us to tackle non-equilibrium
regimes far from the linear-response theory and to properly treat the quantum
fluctuations of the system and its dissipative environment.

We account for the finite electronic interaction $U$ by means of non-
equilibrium DMFT. This mean-field approximation simplifies the task of solving
the original extended lattice model by mapping it to a self-consistently
determined local impurity model. In our case, the impurity model consists of a
dimer site (with two orbitals) coupled to an out-of-equilibrium fermionic
bath. This constitutes a non-equilibrium implementation of cluster-DMFT [69,
70, 71] with a cluster of size two.

Given the permutation symmetry between the orbitals 1 and 2 of the DHM
Hamiltonian in Eq. (2), it is simpler to work in the bonding (B) - antibonding
(A) basis which diagonalizes the non-interacting problem. In terms of the
original orbital degrees of freedom, this amounts to working with the creation
operators $c^{\dagger}_{i{\rm
B}\sigma}=(c_{i1\sigma}^{\dagger}+c_{i2\sigma}^{\dagger})/\sqrt{2}$ and
$c^{\dagger}_{i{\rm
A}\sigma}=(c_{i1\sigma}^{\dagger}-c_{i2\sigma}^{\dagger})/\sqrt{2}$. The real-
time retarded and Keldysh Green’s functions are defined as

$\displaystyle\begin{array}[]{rl}G^{\rm R}_{ij\,m}(t,t^{\prime})&=-\frac{{\rm
i}}{2}\langle[c_{im\sigma}(t),c^{\dagger}_{jm\sigma}(t^{\prime})]\rangle\Theta(t-t^{\prime})\,,\\\
\vphantom{\biggl{(}}G^{\rm K}_{ij\,m}(t,t^{\prime})&=-\frac{{\rm
i}}{2}\langle\\{c_{im\sigma}(t),c^{\dagger}_{jm\sigma}(t^{\prime})\\}\rangle\,,\end{array}$
(7)

respectively, where $m={\rm B},{\rm A}$ is the band index and $\Theta(t)$ is
the Heaviside step function. Assuming paramagnetic solutions, and to simplify
notations, we dropped the spin index.

The DHM Hamiltonian is invariant under translations in the directions
perpendicular to the electric field $\bm{E}=E\bm{u}_{x}$. This ensures that
the spatial extent of the problem in the $y$-direction can be simply accounted
for by the quantum number $k_{y}$. However, our choice to work with the
Coulomb gauge implies that the translational invariance is formally broken
along $\bm{u}_{x}$. This prompts us to work with Green’s functions that are
evaluated at equal coordinates along the electric field direction: $G^{\rm
R/K}(\omega,k_{y};x_{i}=x_{j})$. Because we are targeting spatially
homogeneous steady-states, all sites are assumed _physically_ equivalent and,
for convenience, we pick $x_{i}=x_{j}=0$. These Green’s functions obey the
following Dyson’s equations on the lattice [29, 35]

$\displaystyle\begin{array}[]{l}G^{\rm
R}_{m}(\omega,k_{y})^{-1}\\!=\omega\\!-\\!\epsilon^{\perp}_{m}(k_{y})\\!-\\!\Sigma^{\rm
R}_{m}(\omega)\\!-\\!t^{2}F_{m}^{\rm R}(\omega,k_{y}),\\\
\vphantom{\biggl{(}}G^{\rm K}_{m}(\omega,k_{y})=\lvert G^{\rm
R}_{m}(\omega)\rvert^{2}\\!\left[\Sigma^{\rm K}_{m}(\omega)\\!+\\!t^{2}F^{\rm
K}_{m}(\omega,k_{y})\right].\end{array}$ (10)

$\epsilon^{\perp}_{\rm A}(k_{y})=-2t\cos(k_{y})+t^{\perp}$ and
$\epsilon^{\perp}_{\rm B}(k_{y})=-2t\cos(k_{y})-t^{\perp}$ are the band
dispersion relations in the sublattice corresponding to the directions
perpendicular to the electric field. We introduced the quantities $F^{\rm
R/K}_{m}(\omega,k_{y})=F^{\rm R/K}_{m+}(\omega+E,k_{y})+F^{\rm
R/K}_{m-}(\omega-E,k_{y})$ which stem from the hybridization of a given site
to the semi-infinite chains of its neighbors along $\pm\bm{u}_{x}$. They obey
the self-consistent equations

$\displaystyle F^{\rm
R}_{m\pm}(\omega,k_{y})^{-1}\\!=\\!\omega\\!-\\!\epsilon^{\perp}_{m}(k_{y})\\!-\\!\Sigma^{\rm
R}_{m}(\omega)\\!-\\!t^{2}F^{\rm R}_{m\pm}(\omega\\!\pm\\!E,k_{y}),$
$\displaystyle F^{\rm K}_{m\pm}(\omega,k_{y})\\!=\\!\lvert F^{\rm
R}_{m\pm}(\omega,k_{y})\rvert^{2}\big{[}\Sigma^{\rm
K}_{m}(\omega)\\!+\\!t^{2}F^{\rm K}_{m\pm}(\omega\\!\pm\\!E,k_{y})\big{]}.$
(11)

The self-energies have contributions from the dissipative environment and from
the Hubbard $U$ interaction: $\Sigma=\Sigma_{\Gamma}+\Sigma_{U}$ where

$\displaystyle\begin{array}[]{rl}\vphantom{\Bigl{(}}\Sigma^{\rm
R}_{\Gamma}(\omega)&=-{\rm i}\Gamma\,,\\\ \Sigma^{\rm
K}_{\Gamma}(\omega)&=-2{\rm
i}\Gamma\tanh{\left(\frac{\omega}{2T}\right)}\,,\end{array}$ (14)

and $\Sigma_{U}$ will be discussed below. The expression of $\Sigma^{\rm
K}_{\Gamma}$ is dictated by the fluctuation-dissipation theorem (FDT),
$\Sigma^{\rm K}_{\Gamma}(\omega)=2{\rm
i}\tanh{\left(\frac{\omega}{2T}\right)}\mathrm{Im}\,\Sigma^{\rm
R}_{\Gamma}(\omega)$, which is applicable since the dissipative environment is
assumed to remain in equilibrium. In Eqs. (10) and (11) above, we already
implemented the DMFT approximation which consists in assuming that
$\Sigma_{U}$ is local, namely that it does not depend on the momentum $k_{y}$
but only on the frequency $\omega$.

In the spirit of the DMFT, the local Green’s functions $G_{m}^{\rm
R/K}(\omega)=\sum_{k_{y}}G_{m}^{\rm R/K}(\omega,k_{y})$ are identified to
those of a quantum impurity problem consisting of a single dimer site coupled
to an ad hoc non-interacting non-equilibrium environment. The non-interacting
Green’s functions $\mathcal{G}^{\rm R/K}(\omega)$ of the impurity problem are
determined by the following Schwinger-Dyson’s equations,

$\displaystyle\begin{array}[]{rl}&\mathcal{G}^{\rm
R}_{m}(\omega)^{-1}\\!=G^{\rm R}_{m}(\omega)^{-1}+\Sigma^{\rm
R}_{Um}(\omega)\;,\\\ \vphantom{\biggl{(}}&\mathcal{G}^{\rm
K}_{m}(\omega)=\lvert\mathcal{G}^{\rm
R}_{m}(\omega)\rvert^{2}\left[\frac{G^{\rm K}_{m}(\omega)}{\lvert G^{\rm
R}_{m}(\omega)\rvert^{2}}-\Sigma^{\rm K}_{Um}(\omega)\right].\end{array}$ (17)

Note that the Green’s functions and the self-energies are diagonal in both the
spin and band index due to the original spin and orbital permutation
symmetries of the driven-dissipative DHM.

The local self-energy $\Sigma_{U}$ is computed by means of iterated
perturbation theory (IPT) [72, 27] which treats the electronic interaction $U$
to second order. IPT has already been used in the context of the DHM in
thermal equilibrium [65, 45]. For a compact expression in the non-equilibrium
steady state, it is convenient to go back to the original Kadanoff-Baym-
Keldysh contour ($\alpha,\beta=+,-$), to real time, and to the orbital basis
($a,b=1,2$). At half-filling, it reads

$\displaystyle\Sigma^{\alpha\beta}_{U\,ab}(\tau)=-\alpha\beta\,U^{2}\,\mathcal{G}^{\alpha\beta}_{ab}(\tau)^{2}\mathcal{G}^{\beta\alpha}_{ba}(-\tau)\,.$
(18)

The corresponding expression in the Keldysh basis is given in App. B.

The IPT expressions are valid in the weakly-interacting regime $U/t\rightarrow
0$, and they have also been shown to be exact in the dimer limit
$t/U\rightarrow 0$ at half-filling in zero-temperature equilibrium. For
intermediate, finite, values of $U/t$, the IPT provides a crude approximation
to the exact impurity self-energy which has already proven to be extremely
effective at capturing the qualitative aspects of the Mott equilibrium phase
transition of the paramagnetic SOHM [73, 27]. It has also been successfully
used in similar non-equilibrium contexts [33, 37, 35]. Notably, the electric-
field induced dimensional crossover at strong fields also justifies the use
IPT out of equilibrium since it correctly reproduces the limit $E/U\to\infty$
[33]. At finite temperatures, the IPT produces spurious in-gap states in the
insulating DOS of the DHM. However, we argue in App. B that they have a very
limited impact on the equilibrium and non-equilibrium phase transitions
discussed in this manuscript.

Our NESS DMFT algorithm proceeds as follows. We start from an educated guess
for the self-energy kernels. In the metallic side, we start from the non-
interacting limit $\Sigma^{\rm R}_{U}(\omega)=0$, and in the insulating side
from the non-dissipative dimer limit at zero-temperature, $\Sigma^{\rm
R}_{U\,{\rm B/A}}(\omega)=\frac{U^{2}}{4}\frac{1}{\omega\mp 3t^{\perp}+{\rm
i}0^{+}}$, see App. B. ($i$) The $F^{\rm R/K}(\omega,k_{y})$ are determined by
solving self-consistently the set of equations (11). ($ii$) The lattice Greens
functions $G^{\rm R/K}(\omega,k_{y})$ are computed via Dyson’s Eqs. (10).
($iii$) The impurity non-interacting Green’s functions $\mathcal{G}^{\rm
R/K}(\omega)$ are determined via Dyson’s Eqs. (17). ($iv$) The self-energies
$\Sigma^{\rm R/K}_{U}(\omega)$ are updated using the IPT expression in Eq.
(18). Steps ($i$) to ($iv$) are repeated until convergence is achieved. When
varying the electric field, we use the previously converged solution as the
educated guess for the self-energy kernels.

##  IV Results

Let us now present the fate of the driven-dissipative DHM when the electric
field is varied by studying both its non-equilibrium spectral and transport
properties.

###  IV.1 Non-equilibrium spectral features

Figure 4:  Close-up of the gapped local density of states in the non-
equilibrium DHM insulator (note the semi-log scale). At equilibrium ($E=0$,
dashed line), $t^{\perp}$ opens a clean gap of width $\Delta\approx 0.5$
around the Fermi energy. In the presence of the electric field ($E=0.1$, solid
line), the edges of the gap are replicated into Bloch-Zener side-bands located
at $\pm\Delta/2+n|q|Ea$, $n\in\mathbb{Z}$. The parameters are $t^{\perp}=0.3$,
$U=2.5$, $\Gamma=7.5\times 10^{-4}$, $T=2.5\times 10^{-3}$.

We first discuss the local density of states of an orbital, $\rho_{\rm
loc}(\omega)=-\frac{1}{\pi}\sum_{k_{y}}\frac{1}{2}\sum_{m}\mathrm{Im}\,G_{m}^{\rm
R}(\omega,k_{y})$, when the electric field $E$ is turned on. The electric
field has little impact on the density of states of the metal. However, it
induces distinctive features in the density of states of the insulator: in-gap
states are created within the correlated gap, sometimes referred to as Bloch-
Zener or Wannier-Stark side-bands. This is illustrated in Fig. 4. The origin
of these in-gap structures is easily understood in the context of an electric-
field-driven non-interacting one-dimensional band insulator: distinct regions
of finite density of states appear beyond the edges of the bands. They are
distant in energy by multiples of the potential drop $|q|Ea$ and are
exponentially suppressed in the gap on a scale $t^{1/3}(|q|Ea)^{2/3}$ [74].
Notably, their spectral content is inherited from the spectral content of the
band insulator at equilibrium but in lower dimensions: the dimensions
perpendicular to the electric field. A similar scenario was already identified
and described in the context of a correlated insulator (SOHM) in Refs. [37,
64, 42]. In the case of the DHM with a finite $t^{\perp}$, the edges of the
gap located at $\pm\Delta/2$ are responsible for Bloch-Zener side-bands
located at $\pm\Delta/2+n|q|Ea$ ($n\in\mathbb{Z}$), see Fig. 4. Features
corresponding to values of $|n|>2$ are strongly suppressed.

###  IV.2 Energy distribution function

Figure 5:  Local distribution function $f_{\rm loc}(\omega)$ of the electric-
field-driven DHM, defined in Eq. (19). (a) In the metallic phase, $f_{\rm
loc}(\omega)$ (solid line) cannot be distinguished from a Fermi-Dirac
distribution (dashed line) at the temperature $T_{\rm eff}=0.38$ determined
via Eq. (21). (b) In the insulating phase, $f_{\rm loc}(\omega)$ corresponds
to $T_{\rm eff}=9.3\times 10^{-3}$ but it markedly departs from a Fermi-Dirac
distribution. Both examples were taken at the same bath temperature
$T=2.5\times 10^{-3}$, and parameters $E=0.1$, $t^{\perp}=0.3$, $U=2.5$, and
$\Gamma=7.5\times 10^{-4}$.

Let us now discuss the local distribution function of an electronic orbital

$\displaystyle f_{\rm
loc}(\omega)=\frac{1}{2}\left(1-\frac{1}{2}\frac{\sum_{m}\mathrm{Im}\,G^{\rm
K}_{m}(\omega)}{\sum_{m}\mathrm{Im}\,G^{\rm R}_{m}(\omega)}\right),$ (19)

which informs on the non-equilibrium energy fluctuations. Similarly to the
DOS, the electric field has a mild impact on the energy distribution function
of the metal: it perfectly matches the equilibrium Fermi-Dirac function albeit
at a higher temperature; see Fig. 5 and the discussion on effective
temperature below. On the other hand, the energy distribution function of the
insulator significantly departs from equilibrium: it features structures
equally spaced by the energy $|q|Ea$. These structures have already been
reported in studies regarding the SOHM. They can be seen as a consequence of
the coupling of the impurity site to its neighboring sites in the direction of
the field, $\pm\bm{u}_{x}$. These sites act as baths with chemical potentials
shifted by the potential drop $\pm|q|Ea$ [75].

###  IV.3 Non-equilibrium phase diagram

Figure 6:  Non-equilibrium phase diagram of the electric-field-driven Dimer-
Hubbard model (DHM) in the $U-E$ plane. The IMT (MIT) line is in blue (green).
The inset shows a larger region of the coexistence region. The parameters are
$t^{\perp}=0.3$, $T=2.5\times 10^{-3}$, and $\Gamma=7.5\times 10^{-4}$.

Let us now present the non-equilibrium phase diagram of the electric-field-
driven DHM. It is computed in Fig. 6 in the $U-E$ plane for a fixed value of
$t^{\perp}$. The different phases are determined depending on the presence of
a gap at the Fermi level in the local density of states, or the lack thereof.
In practice, we track the value of $\rho_{\rm loc}(\omega=0)$, which presents
a sharp variation of at least 2 orders of magnitude between phases, as is
shown in Fig. 7. Similarly to the equilibrium phase diagram presented in Fig.
1, the solid lines delimit a region (shaded blue) where both the metal and the
insulator coexist. However, while the IMT and MIT lines (blue and green,
respectively) have general trends that resemble their equilibrium
counterparts, the electric-field scales on which they develop are drastically
smaller on the MIT side than on the IMT side: $E_{\rm IMT}\gg E_{\rm MIT}$.
This results in a very wide coexistence region that extends to very large
values of $U$ and $E$. For very large electric fields, outside of the ranges
of Fig. 1, the IMT and MIT become smooth crossovers between bad metals and
dirty insulators.

This is quite in contrast to the case of the electric-field-driven SOHM whose
corresponding non-equilibrium phase diagram studied in Ref. [35] for
$T=5\times 10^{-3}$ and $\Gamma=3.4\times 10^{-3}$ was found to sport a closed
coexistence region qualitatively similar to its equilibrium phase diagram in
the $U-T$ plane, c.f. the $t^{\perp}=0$ case in Fig. 1.

Figure 7:  Local density of states at the Fermi level as a function of the
electric field. Starting from an insulating solution at $E=0$, the sharp leap
of two orders of magnitude is used to locate the insulator-to-metal transition
(IMT). The parameters are $U=2.5$, $t^{\perp}=0.3$, $T=2.5\times 10^{-3}$, and
$\Gamma=7.5\times 10^{-4}$.

###  IV.4 Transport properties

Let us now discuss the far-from-equilibrium transport properties of the
electric-field-driven DHM. The electric field generates an electric current
per site and per spin $\bm{J}=J\bm{u}_{x}$ which is computed as,

$\displaystyle J$ $\displaystyle=\frac{t^{2}}{2}\sum_{m}\int\frac{{\rm
d}\omega}{2\pi}\sum_{k_{y}}\left\\{\mathrm{Im}G^{\rm
R}_{m}(\omega,k_{y})\mathrm{Im}F^{\rm K}_{m}(\omega,k_{y})\right.$
$\displaystyle\left.\qquad\qquad\qquad-\mathrm{Im}G^{\rm
K}_{m}(\omega,k_{y})\mathrm{Im}F^{\rm R}_{m}(\omega,k_{y})\right\\},$ (20)

where the first summation is performed over the bonding ($m={\rm B}$) and
antibonding ($m={\rm A}$) bands.

In Fig. 8, we illustrate the $I-V$ characteristics of the DHM starting from
low-temperature equilibrium states in the coexistence region and ramping up,
and then down, the electric field. The $I-V$ curves display hysteresis both
when starting from the metal or the insulator. When starting from a metallic
state in the coexistence region, one first expects a linear regime
$J=\sigma_{\rm dc}E$, where $\sigma_{\rm dc}$ is the DC conductivity. The
latter quickly renormalizes as the temperature of the sample increases,
relegating the linear regime down to extremely small values of electric field
which cannot be seen on the scale used here [35]. At $E=E_{\rm MIT}$, the DHM
experiences an MIT where the current density drops by 9 orders of magnitude.
The DHM remains an insulator after a subsequent decrease in the electric
field.

When starting from an insulating state in the coexistence region, the applied
electric field has little effect on the slowly increasing current. Indeed, the
very clean nature of the gap opened by $t^{\perp}$, as discussed above, makes
$J$ much suppressed compared to the case of the SOHM. The Bloch-Zener
structures described above have an impact on the current at finite electric
fields: the current is enhanced whenever they provide a pathway through the
gap. Indeed, states at the lower edge of the gap can be excited by a sequence
of electric-field driven transitions of energy $|q|Ea$. This happens at values
$|q|Ea=\Delta/m$, where $m$ is a small positive integer. In Fig. 8, this is
seen for $m=4$ (corresponding to $n=2$ in Fig. 4). Smaller values of $m\leq 3$
correspond to field strengths that lay outside the range of the figure and
after the IMT. Larger values of $m\geq 5$ correspond to exponentially
suppressed BZ in-gap structures and to smaller features in the $I-V$ which are
hard to isolate. Note that these Bloch-Zener effects cannot be seen in the MIT
case given the much smaller values of the electric field involved. As the
electric field is furthermore increased, the DHM experiences an IMT at
$E=E_{\rm IMT}\gg E_{\rm MIT}$. Note the violent drop of resistivity by 4
orders of magnitude. The DHM remains metallic after a subsequent decrease in
the electric field.

Figure 8:  Hysteretic $I-V$ characteristics of the DMH (note the semi-log
scale). (a) Starting from a metallic state in the coexistence region at
$U=3.25$, the electric field is increased until the DHM experiences an MIT at
$E_{\rm MIT}$. Later, the electric field is decreased and the DHM remains
insulating. (b) Starting from an insulating state in the coexistence region at
$U=2.5$, the electric field is increased until the DMH experiences an IMT at
$E_{\rm IMT}\gg E_{\rm MIT}$. The parameters are $t^{\perp}=0.3$ and
$\Gamma=7.5\times 10^{-4}$.

###  IV.5 Effective temperature

Figure 9:  Effective temperature $T_{\rm eff}$ computed using Eq. (21) as a
function of the increasing and decreasing electric field (note the semi-log
scale). (a) Starting from a metal with the same parameters as in Fig. 8(a).
(b) Starting from an insulator with the same parameters as in Fig. 8(b). Both
MIT and IMT take place once $T_{\mathrm{eff}}$ reaches the corresponding
equilibrium transition temperature, $T_{\rm MIT}^{\rm eq}(U=3.25)=0.017$ and
$T_{\rm IMT}^{\rm eq}(U=2.5)=0.023$, marked by horizontal dashed lines.

Let us now discuss the thermodynamic properties of the electric-field-driven
DHM. The concept of temperature is a priori ill-defined for a non-equilibrium
steady state. Indeed, a single scalar quantity cannot in principle inform on
the full energy and momentum-dependent non-equilibrium distribution function.
In practice, the notion of effective temperature has however proven useful. In
the long-lasting debate on the microscopic origin of the resistive switching,
a strong advocate for a thermally-driven scenario has been a series of
experiments where the effective temperature of the sample was monitored in
situ via the emission of a fluorescent particle acting as a local temperature
probe [18]. Moreover, we have already shown that the non-equilibrium local
distribution function of the metallic state is a Fermi-Dirac function with a
renormalized temperature (see Fig. 5). In that spirit, we adopt a Sommerfeld-
like definition for the local effective temperature [39]:

$\displaystyle T_{\mathrm{eff}}^{2}=\frac{6}{\pi^{2}}\int{\rm
d}\omega\,\omega\left[f_{\rm loc}(\omega)-\Theta(-\omega)\right]\,,$ (21)

where the local orbital energy distribution function $f_{\rm loc}(\omega)$ has
been defined in Eq. (19). The above definition of $T_{\rm eff}$ can be seen as
a simple measure of the ability of the system to create and sustain
excitations above the Fermi level. Naturally, $T_{\rm eff}$ boils down to the
thermodynamic temperature in equilibrium, when $f_{\rm loc}(\omega)$ is the
Fermi-Dirac distribution $f_{\rm FD}(\omega)=[1+\exp(\omega/T)]^{-1}$. Note
that other choices could have been made to extract a $T_{\rm eff}$ from
$f_{\rm loc}(\omega)$, such as fitting it to a Fermi-Dirac distribution or
using its slope at the Fermi energy,
$-\left(\partial_{\omega}f(\omega)\lvert_{\omega=0}\right)^{-1}/4$. We checked
that our results are qualitatively robust with respect to these alternative
choices. The strength of our definition lies in the fact that it does not
require any fitting parameter and it does not require $f_{\rm loc}(\omega)$ to
be close to a Fermi-Dirac distribution.

In Fig. 9, we monitor the effective temperature of the DHM starting from a
low-temperature equilibrium state and ramping up, and then down, the electric
field. Similarly to Fig. 8, we start from both a metallic state and an
insulating state in the coexistence region. In both cases, $T_{\mathrm{eff}}$
grows monotonously with $E$, starting from the bath temperature $T$ at $E=0$.
Naturally, the stronger currents produced in the metallic states are
responsible for higher effective temperatures than in the insulating states.
This is reflected in the huge discontinuities of $T_{\mathrm{eff}}$ at the
transitions. Moreover, the high $T_{\mathrm{eff}}$ after the IMT induces
states known as ‘bad metals’, characterized by a low (but not gapped) density
of states around the Fermi energy. A first outcome of Fig. 9 is the notion of
state-dependent effective temperature: similarly to the hysteretic $I-V$
curves, the effective temperature has a region of bistability tied to the
metastability of the state. In other words, for the same system parameters and
the same electric field, the effective temperature greatly depends on the
state the system, whether insulating or metallic. A second outcome of Fig. 9
is the fact that both non-equilibrium transitions, IMT and MIT, occur whenever
the effective temperature matches the corresponding equilibrium transition
temperature (represented with dotted lines).

In Fig. 10, we repeat the previous analysis on the entire non-equilibrium
phase diagram. We prepare low-temperature equilibrium metallic and insulating
samples at all values of $U$ in the coexistence region, we increase the
electric field, and we report the effective temperatures $T^{\rm eq}_{\rm
MIT}$ and $T^{\rm eq}_{\rm IMT}$ measured exactly at MIT and the IMT,
respectively. Remarkably, we confirm the previous observations made in Fig. 9:

$\displaystyle\begin{array}[]{rl}T_{\rm eff}(E_{\rm
IMT};\mathrm{insulator})&\approx T^{\rm eq}_{\rm IMT}\,,\\\
\vphantom{\biggl{(}}T_{\rm eff}(E_{\rm MIT};\mathrm{metal})&\approx T^{\rm
eq}_{\rm MIT}\,.\end{array}$ (24)

This matching is one of the main results of this manuscript. Equations (24)
account for the large difference in magnitude of the threshold fields $E_{\rm
IMT}\gg E_{\rm MIT}$ even though $T^{\rm eq}_{\rm IMT}$ and $T^{\rm eq}_{\rm
IMT}$ are on the same order of magnitude. Notably, this is in agreement with
the experimental findings of Refs. [18, 22] where the resistivity of VO2
driven by a DC bias was found to match the equilibrium resistivity, once the
voltage bias is parameterized in terms of the local effective temperature of
the sample. It was already demonstrated numerically for the MIT occurring in
the electric-field-driven dissipative SOHM [35], and it was given some
analytical insight in the context of a simplified one-dimensional model with a
correlated gap treated within mean-field theory [25]. This result articulates
both the thermal and the electronic scenarios for resistive switching in a
unified framework: once the non-equilibrium electronic mechanisms are
accounted for by an effective temperature, the non-equilibrium phase
transition reduces to the thermally-driven equilibrium phase transition.

Figure 10:  Phase diagram of the electric-field-driven DHM in the $U-T_{\rm
eff}$ plane. The effective temperatures $T_{\rm eff}$ were computed using Eq.
(21). The blue (green) dots correspond to $T_{\rm eff}$ computed just at the
IMT (MIT) induced by increasing the electric field. The magenta crosses
indicate instances of $T_{\rm eff}$ at the transition computed by fitting
$f_{\mathrm{loc}}(\omega)$ to a Fermi-Dirac distribution. The solid lines
correspond to the equilibrium IMT and MIT driven by increasing the
temperature. The parameters are $t^{\perp}=0.3$, $T=2.5\times 10^{-3}$, and
$\Gamma=7.5\times 10^{-4}$.

##  V Discussion and conclusion

In this work, we investigated a driven-dissipative version of the DHM in which
both first-order MIT and IMT can be driven by a DC electric field rather than
temperature.

The field-driven MIT was found to follow a simple Joule-heating scenario in
which the driven metal essentially remains in local thermal equilibrium
despite the finite electronic current and heat flow. The balance between the
power injected by the electric field and the heat released to the environment
brings the sample to its equilibrium transition temperature, triggering the
MIT at threshold fields controlled by the dissipation rate, $E_{\rm
MIT}\sim\sqrt{\Gamma}$ [35], and corresponding to meV energy scales.

The field-driven IMT was found to follow a different scenario in which both
the inter-orbital hopping and Bloch-Zener (BZ) effect play a crucial role.
Indeed, $t^{\perp}$ is responsible for opening a clean gap $\Delta$ within the
pseudogap of the Mott insulator. Such a clean gap provides a very good
insulating ground state where Joule-heating effects are strongly suppressed,
auguring very large threshold fields to overcome the gap $\Delta$. However, we
found a substantial contribution of the BZ effect to Joule-heating: the
electric field is responsible for the formation of in-gap states, located at
energies given by multiple values of the electric field, which can form a non-
equilibrium pathway to bridge the gap. Altogether, we found threshold fields
controlled by fractions of the gap, $E_{\rm MIT}\sim\Delta/4$. This large 100
meV scale is at odds with the available experimental data on threshold fields.
This calls for further investigation in the modeling of these driven-
dissipative correlated materials. A promising avenue is to explore the
influence of a phononic dissipative environment rather than the simple
electronic dissipation used here [40, 42]. Notably, it was recently reported
that the interplay of the electric field with such a bosonic bath could create
a strong non-equilibrium pathway through the gap resulting in electronic
avalanches with fields at sub meV scales [41].

Importantly, the study of the electric-field-driven SOHM and DHM has shown
that there are two key factors in building up toward the electric-field-driven
transition: the nature of the insulating state, and its hybridization with a
dissipative environment. The MIT in both models proceeds with a simple Joule
heating mechanism involving excitations at the Fermi level. The IMT in the
SOHM proceeds similarly but the in-gap excitations are produced by the
hybridization with the dissipative bath. For the IMT in the DHM, the low-
energy insulating state is impervious to the bath hybridization and the
destabilization of the insulator comes from a higher energy mechanism, namely,
the Bloch-Zener effect. The fact that these three distinct non-equilibrium
scenarios can be unified in a single equilibrium framework once their effects
are measured in terms of effective temperature is a non-trivial observation.
We believe this observation is not limited to the particular model at hand and
generalizes to other cases of RS in correlated materials with possibly other
microscopic mechanisms at play.

Most experiments show IMTs with spatially heterogeneous solutions in the form
of metallic filaments which were not discussed in this manuscript. Since the
insulating phase is spatially homogeneous, it is reasonable to think that both
the mechanism driving the onset of the switching and the value threshold field
do not involve heterogeneous dynamics. More generally, our homogeneous
solutions can be seen as mesoscopic building blocks for a larger sample with a
heterogeneous landscape. Relaxing the assumption of homogeneity would require
the use of inhomogeneous DMFT [76, 77, 29, 78, 79] which is computationally
challenging. Let us rather highlight the steps for developing an effective
field theory description for the electric-field-driven RS which leverages the
teachings of our microscopic analysis and operates at a coarse-grained level.
With our results in mind, this can simply be achieved by starting from the
existing equilibrium field theories and by replacing the thermodynamic
temperature by the state-dependent effective temperature $T_{\rm eff}$. In the
context of the Hubbard model, we can therefore start from the results of Refs.
[80, 81, 82, 83] and propose an effective free-energy functional of the form

$\displaystyle\mathcal{F}_{\rm NESS}[\phi]=$ $\displaystyle\int{\rm
d}^{d}x\frac{1}{2}(\nabla_{x}\phi)^{2}+\frac{1}{2}\vec{\mu}\cdot\left(\vec{K}(\phi)-\vec{K}^{*}\right)\phi^{2}$
$\displaystyle\qquad+\frac{1}{4}\lambda\phi^{4}-\vec{h}\cdot\left(\vec{K}(\phi)-\vec{K}^{*}\right)\phi\;,$
(25)

where the field $\phi$ is the scalar order parameter, $d$ is the dimension of
the space, $\vec{K}(\phi)-\vec{K}^{*}=\left(T_{\rm
eff}(\phi;E)-T^{*},U-U^{*}\right)$ is a two-component vector which measures
the distance to criticality, $\vec{\mu}$ and $\vec{h}$ are two-component
parameters, and $\lambda>0$. In practice, the free energy has to be
supplemented with kinematic and geometric constraints on the fields such as
the Maxwell and Kirchhoff laws. Such a low-energy description, to be presented
in Ref. [84], promises to bridge the gap between the microscopic computations
discussed in this manuscript and the heuristic models of resistor networks
that have successfully been used to account for the filament formation
observed in experiments [85, 86].

#### Acknowledgements.

MD and CA acknowledge the support from the French ANR “MoMA” project
ANR-19-CE30-0020. CA is grateful for the support from project 6004-1 of the
Indo-French Centre for the Promotion of Advanced Research (IFCPAR). We warmly
thank Marcelo Rozenberg and Marcello Civelli for inspiring discussions. MD
wishes to thank Lorenzo Fratino and Soumen Bag for many helpful conversations
and advice.

##  Appendix A Estimation of $\Gamma$ from experimental data

Heat transfer to the electronic baths, controlled by the single rate $\Gamma$,
is the only dissipative mechanism included in the driven-dissipative DHM given
in Eq. (1). In reality, the dissipative environment of the active electronic
degrees of freedom is more complex with, in particular, the coupling to
lattice phonons. Consequently, the heuristic parameter $\Gamma$ is to be
understood as an effective dissipative rate that aims to encompass the
different dissipative processes that are not explicitly included in our model.
In this Section, we discuss how we estimated from experimental data the value
$\Gamma=7.5\times 10^{-4}$ eV used throughout this manuscript.

In a steady state, the electric field is constantly pumping energy into the
system, which is dissipated into the reservoirs. The rate of heat exchange
with the reservoir per unit volume is computed as [63]

$\mathcal{Q}=2\Gamma\sum_{m}\int{\rm
d}\omega\rho_{m}(\omega)\left[f_{m}(\omega)-f_{\mathrm{th}}(\omega)\right]\omega\,.$
(26)

Index $m$ labels the bonding and antibonding bands ${\rm B}$ and ${\rm A}$,
and $\rho_{m}(\omega)$, $f_{m}(\omega)$ are the band-specific DOS and local
energy distribution functions. $f_{\textrm{th}}(\omega)$ is the energy
distribution function of the bath, which is a Fermi-Dirac function at
temperature $T$.

According to the first law of thermodynamics, the heat and work rates
compensate in a steady state, $\mathcal{W+Q}=0$, and we therefore can
determine $\mathcal{W}$ through Eq. (26).

In a metallic state, we can reasonably approximate $f_{m}(\omega)\approx
f_{\rm loc}(\omega)$, and $\rho_{m}(\omega)\approx 1/8t$ for $\omega\sim 0$,
where the integrand in the r.h.s. of (26) is nonzero. Furthermore, using the
definition of $T_{\rm eff}$ in Eq. (21), and reinstating the universal
constants, we obtain an expression for the dissipative rate $\Gamma$ that can
be used to match to experimental data:

$\Gamma=\frac{12}{\pi^{2}}\frac{\hbar}{k_{\rm
B}^{2}}\frac{t\mathcal{W}}{T_{\textrm{eff}}^{2}-T^{2}}.$ (27)

We extract the corresponding experimental values from Ref. [19]. In this
experiment, a thin film of VO2 of volume $6,0\times 10^{-21}$
$\mathrm{m^{3}}$, placed on a TiO2 substrate (which plays the role of the
bath), is driven through the transition by an external electric field. In the
presence of this substrate, the critical temperature is brought down to
$T^{\rm eq}_{\rm IMT}=270$ K. $T_{\rm eff}$ is shown to reach this value at
the electrically-driven transition, which is where we choose to evaluate the
r.h.s. of Eq. (27). The total power needed to bring the system to the
transition is measured to be $\mathcal{W}\times N_{\rm s}=6$ mW for the bath
temperature $T=230$ K. $N_{\rm s}$ denotes the number of vanadium atoms in the
sample, which we estimated from the volume of the sample and the lattice
constant $a=4.5$ Å(caption, Table 1).

Replacing the experimental values above in Eq. (27), we compute an estimate
for the dissipative rate $\Gamma\sim 10^{-3}\textrm{eV}$. In the main text, we
use a slightly smaller value of $7.5\times 10^{-4}$ eV. This ensures that
$\Gamma$ remains the smallest energy scale in the system (below the bath
temperature $T$) and we stay within the bounds of the IPT approximation (see
App. B).

##  Appendix B Iterated Perturbation Theory (IPT) for the DHM

###  B.1 Practical expressions

Figure 11: Second-order Feynman diagrams contributing to the IPT self-energy
$\Sigma_{U\,m}^{\alpha\beta}$ in the bonding - antibonding basis ($m=$ A, B ;
$\overline{\mathrm{A}}=$ B, $\overline{\mathrm{B}}=$ A) and in the Kadanoff-
Baym-Keldysh contour basis ($\alpha,\beta=+,-$). The lines are the non-
interacting impurity Green’s function $\mathcal{G}^{\alpha\beta}_{m}$. Spin
indices were dropped for the sake of compactness.

In this Section, we provide an alternative expression for the IPT self-energy
of the DHM, useful for its numerical implementation.

Figure 11 displays the Feynman diagrams contributing to the IPT self-energy.
In Eq. (18), we reported its expression at half-filling in the basis of the
original Kadanoff-Baym-Keldysh contour ($\alpha,\beta=+,-$) and in the real
time domain. In practice, we employ a more practical formulation, and we work
with the retarded and Keldysh components. Moreover, it is most efficient to
compute only the imaginary parts, since they are compactly supported functions
in Fourier space (frequency domain). The real part of the retarded component
can eventually be recovered through the Kramers-Kronig relations, and the
Keldysh component as defined in Eq. (7) is purely imaginary. Applying the
particle-hole relations at half-filling, given by

$\displaystyle\mathrm{Im}\mathcal{G}^{\rm R}_{ab}(\omega)$
$\displaystyle=\;\;\,(2\delta_{ab}-1)\mathrm{Im}\mathcal{G}^{\rm
R}_{ab}(-\omega)\;,$ $\displaystyle\mathrm{Im}\mathcal{G}^{\rm
K}_{ab}(\omega)$ $\displaystyle=-(2\delta_{ab}-1)\mathrm{Im}\mathcal{G}^{\rm
K}_{ab}(-\omega)\;,$ (28)

we come to

$\begin{split}\mathrm{Im}\Sigma_{U\,ab}^{\rm
R}(\omega)=\frac{U^{2}}{4}\\{-3\mathcal{G}_{ab}^{\rm
K}(\omega)\circ\mathcal{G}_{ab}^{\rm
K}(\omega)\circ\mathrm{Im}\mathcal{G}_{ab}^{\rm R}(\omega)&\\\
+4\mathrm{Im}\mathcal{G}_{ab}^{\rm
R}(\omega)\circ\mathrm{Im}\mathcal{G}_{ab}^{\rm
R}(\omega)\circ\mathrm{Im}\mathcal{G}_{ab}^{\rm
R}(\omega)\\}(2\delta_{ab}-1)\;,&\end{split}$ (29)

and

$\begin{split}\Sigma_{U\,ab}^{\rm
K}(\omega)=\frac{U^{2}}{4}\\{12\mathcal{G}_{ab}^{\rm
K}(\omega)\circ\mathrm{Im}\mathcal{G}_{ab}^{\rm
R}(\omega)\circ\mathrm{Im}\mathcal{G}_{ab}^{\rm R}(\omega)&\\\
-\mathcal{G}_{ab}^{\rm K}(\omega)\circ\mathcal{G}_{ab}^{\rm
K}(\omega)\circ\mathcal{G}_{ab}^{\rm
K}(\omega)\\}(2\delta_{ab}-1)\;,&\end{split}$ (30)

where we have denoted the convolution between two functions as $(f\circ
g)(\omega)\equiv\int\frac{{\rm
d}\omega^{\prime}}{2\pi}f(\omega-\omega^{\prime})g(\omega^{\prime})$, and
$a,b=1,2$ are the orbital basis indices. Rotation to the bonding - antibonding
basis ($m=\mathrm{B}$, $\mathrm{A}$) used throughout the manuscript can be
made via the following dictionary:

$\displaystyle G_{\rm B}=G_{11}+G_{12}\;,$ (31) $\displaystyle G_{\rm
A}=G_{11}-G_{12}\;,$

which applies as well for the self-energies.

###  B.2 Dimer limit

Energy level | Eigenstate | $N$
---|---|---
$-\frac{U}{2}\left(1-\sqrt{1+\left(\frac{4t^{\perp}}{U}\right)^{2}}\right)$ | $\left(1-\sqrt{1+\left(\frac{4t^{\perp}}{U}\right)^{2}}\right)(\ket{\uparrow,\downarrow}-\ket{\downarrow,\uparrow})+\frac{4t^{\perp}}{U}\left(\ket{\uparrow\downarrow,0}+\ket{0,\uparrow\downarrow}\right)$ | $2$
| $\ket{\uparrow\downarrow,0}-\ket{0,\uparrow\downarrow}$ | $2$
$0$ | $\ket{\uparrow\downarrow,\uparrow\downarrow}$ | 4
| $\ket{0}$ | 0
$t^{\perp}-U/2$ | $\ket{\uparrow,\uparrow\downarrow}+\ket{\uparrow\downarrow,\uparrow},\;\;\ket{\downarrow,\uparrow\downarrow}+\ket{\uparrow\downarrow,\downarrow}$ | $3$
| $\ket{\uparrow,0}-\ket{0,\uparrow},\;\;\ket{\downarrow,0}-\ket{0,\downarrow}$ | $1$
$-t^{\perp}-U/2$ | $\ket{\uparrow,\uparrow\downarrow}-\ket{\uparrow\downarrow,\uparrow},\;\;\ket{\downarrow,\uparrow\downarrow}-\ket{\uparrow\downarrow,\downarrow}$ | $3$
| $\ket{\uparrow,0}+\ket{0,\uparrow},\;\;\ket{\downarrow,0}+\ket{0,\downarrow}$ | $1$
$-U$ | $\ket{\uparrow,\uparrow},\;\;\ket{\downarrow,\downarrow},\;\;\ket{\uparrow,\downarrow}+\ket{\downarrow,\uparrow}$ | $2$
$-\frac{U}{2}\left(1+\sqrt{1+\left(\frac{4t^{\perp}}{U}\right)^{2}}\right)$ | $\ket{\phi_{\mathrm{GS}}}=\left(1+\sqrt{1+\left(\frac{4t^{\perp}}{U}\right)^{2}}\right)(\ket{\uparrow,\downarrow}-\ket{\downarrow,\uparrow})+\frac{4t^{\perp}}{U}\left(\ket{\uparrow\downarrow,0}+\ket{0,\uparrow\downarrow}\right)$ | $2$
Table 2:  Eigensystem of the DHM Hamiltonian in the dimer limit. We use the
shorthand notation $\ket{\psi,\phi}\equiv
c^{\dagger}_{1\psi}c^{\dagger}_{2\phi}\ket{0}$, with
$\psi,\phi\in\\{\uparrow,\downarrow,\uparrow\downarrow\\}$,
$c^{\dagger}_{a\uparrow\downarrow}\equiv
c^{\dagger}_{a\uparrow}c^{\dagger}_{a\downarrow}$. The states are ordered from
the highest excited state (top) to the ground state (bottom), and they are not
normalized. $N$ denotes the number of particles in each state, which is a
well-defined quantum number for all eigenstates since $[H,N]=0$. The other
good quantum numbers are the total spin $\bm{S}^{2}$, the spin projection
$S_{z}$, and the orbital permutation parity $P_{12}$.

Let us now discuss the range of validity of the IPT approximation. To this
end, we work in a limit where the system is diagonalizable and we can derive
exact expressions, in order to compare them to the IPT scheme.

We consider the equilibrium DHM at half-filling, in the dimer limit
($t/U\rightarrow 0$, keeping $t^{\perp}$ finite) and in the absence of
dissipation ($\Gamma=0$). This defines the two-orbital, single-site
Hamiltonian

$\displaystyle H_{\mathrm{dim}}=$
$\displaystyle-t^{\perp}\sum_{\sigma}(c^{\dagger}_{1\sigma}c_{2\sigma}+c^{\dagger}_{2\sigma}c_{1\sigma})+U(n_{1\uparrow}n_{1\downarrow}+n_{2\uparrow}n_{2\downarrow})$
$\displaystyle\qquad-\frac{U}{2}\sum_{\sigma}(n_{1\sigma}+n_{2\sigma})\,,$
(32)

which can be diagonalized exactly. This limit is usually referred to as the
atomic limit when $t^{\perp}=0$ (SOHM). The energy spectrum and eigenstates
are shown in Table 2. We can subsequently write the bonding retarded Green’s
function of the dimer, $G^{\rm R\,({\rm dim})}_{\rm B}(\omega)$, in the
spectral representation:

$\displaystyle G^{\rm R\,({\rm dim})}_{\rm B}(\omega)=\sum_{\sigma,j,k}$
$\displaystyle\frac{|\bra{\phi_{j}}c^{\dagger}_{{\rm
B}\sigma}\ket{\phi_{k}}|^{2}}{\omega-(E_{j}-E_{k})+{\rm i}0^{+}}$ (33)
$\displaystyle\times(\mathrm{e}^{(E_{j}-E_{\mathrm{GS}})/T}+\mathrm{e}^{(E_{k}-E_{\mathrm{GS}})/T})$
$\displaystyle\stackrel{{\scriptstyle\left(T\to
0\right)}}{{=}}\sum_{\sigma,j}$
$\displaystyle\frac{|\bra{\phi_{j}}c^{\dagger}_{{\rm
B}\sigma}\ket{\phi_{\mathrm{GS}}}|^{2}}{\omega-(E_{j}-E_{\mathrm{GS}})+{\rm
i}0^{+}}$ $\displaystyle+$
$\displaystyle\frac{|\bra{\phi_{\mathrm{GS}}}c_{{\rm
B}\sigma}\ket{\phi_{j}}|^{2}}{\omega-(E_{\mathrm{GS}}-E_{j})+{\rm i}0^{+}}$
(34) $\displaystyle=\;\;\;\;\;\;\;$
$\displaystyle\frac{1}{\omega+t^{\perp}-\frac{U^{2}}{4}\frac{1}{\omega-3t^{\perp}+{\rm
i}0^{+}}}.$ (35)

$\ket{\phi_{j}}$ are the eigenstates of $H_{\rm at}$, $E_{j}$ are the
corresponding energies, and ‘GS’ labels the ground state of the system. From
the last line we extract the bonding self-energy at zero temperature,

$\Sigma^{\rm R\,({\rm
dim})}_{U\,\rm{B}}(\omega)=\frac{U^{2}}{4}\frac{1}{\omega-3t^{\perp}+{\rm
i}0^{+}}\;.$ (36)

Let us now compare this expression to the self-energy computed with IPT in the
dimer limit. In this setting, the non-interacting Green’s functions read

$\displaystyle\mathcal{G}^{\rm R}_{\rm B/A}(\omega)=\frac{1}{\omega\pm
t^{\perp}+{\rm i}0^{+}}\;,$ (37) $\displaystyle\mathcal{G}^{{\rm K}}_{\rm
B/A}=2{\rm i}\tanh{\left(\frac{\omega}{2T}\right)}\mathrm{Im}\mathcal{G}^{\rm
R}_{\rm B/A}(\omega)\;.$

Replacing them in equation (LABEL:eq:IPT_R), we obtain

$\displaystyle\Sigma^{\mathrm{R}}_{U\,\rm{B}}(\omega)=$
$\displaystyle\frac{U^{2}}{16}\biggl{(}\left[1+3\mathrm{tanh}^{2}\left(\frac{t^{\perp}}{2T}\right)\right]\frac{1}{\omega-3t^{\perp}+{\rm
i}0^{+}}$ $\displaystyle+$ $\displaystyle
3\left[1-\mathrm{tanh}^{2}\left(\frac{t^{\perp}}{2T}\right)\right]\frac{1}{\omega+t^{\perp}+{\rm
i}0^{+}}\biggr{)}$ (38) $\displaystyle\stackrel{{\scriptstyle\left(T\to
0\right)}}{{=}}$
$\displaystyle\;\frac{U^{2}}{4}\frac{1}{\omega-3t^{\perp}+{\rm i}0^{+}}\;.$
(39)

Note that the terms with a pole at $3t^{\perp}$ originate from the third
Feynman diagram in Fig. 11. The bonding retarded Green’s function is

$G^{\mathrm{R}}_{\rm
B}(\omega)=\frac{1}{\omega+t^{\perp}-\Sigma^{\mathrm{R}}_{U\,{\rm
B}}(\omega)}\;.$ (40)

Expression (39) perfectly matches $\Sigma^{\rm R\,({\rm dim})}_{U\,\rm B}$ in
Eq. (36), which implies that the IPT approximation at $T=0$ is exact in the
limit $t/U\to 0$. This is remarkable for a perturbative approach, which is
expected to be valid a priori only around the opposite limit $U/t\to 0$.

It is important to remark that the comparisons in this Section are performed
in the absence of a dissipative environment (meaning $\Gamma=0$). Once the
bath degrees of freedom are integrated out (see Section III), it is not
possible to perform exact spectral diagonalization. By keeping $\Gamma$ as the
smallest scale of the system, we minimize any effect it may have on the
validity of the approximation.

The picture is different at finite temperature, where the second term in Eq.
(38) becomes finite. The presence of this term produces a second root in the
denominator of $G^{\mathrm{R}}_{\rm B}$ in Eq. (40), yielding a secondary peak
in the DOS located at $\omega\sim-t^{\perp}$. Comparing to the exact
expression in Eq. (33), there is no peak located at a similar frequency:
inspecting Table 2 shows that $|E_{j}-E_{k}|\neq~{}t^{\perp}$
$\forall\,j,\,k$. Analogous reasoning applies for the antibonding component,
producing a secondary peak at $\omega\sim t^{\perp}$. We note that these
spurious peaks are relatively small compared to the primary structures. In
order to keep this spurious effect under control, when $t^{\perp}$ is finite
we make sure to work in regimes $T\ll t^{\perp}$, where the second term is
suppressed.

This range of validity extends to all temperatures when $t^{\perp}$ is zero or
very small. Indeed, in the lattice model, the width of the spectral features
is controlled by $t$. If $t^{\perp}\ll t$, the spurious peaks described above
cannot be distinguished from the other structures. This means that in this
case, the approximation is still valid at high temperatures. Furthermore, in
the case of the SOHM ($t^{\perp}=0$) the matching is exact for all
temperatures: Eq. (38) has only one term centered around $\omega\sim 0$, and
extracting a self-energy from Eq. (33) yields the same result.

##  Appendix C Estimation of the in-gap density of states

In the following, we derive semi-analytical expressions for the in-gap density
of states of the SOHM and the DHM. The objective is to account for the large
difference observed in Fig. 3.

The density of states of each band $m=$ B, A can be written, in terms of the
non-interacting Green’s functions of the impurity $\mathcal{G}_{m}(\omega)$,
as

$\rho_{m}(\omega)=\frac{1}{\pi}\frac{\mathrm{Im}({\mathcal{G}_{m}^{\rm
R}}^{-1})(\omega)-\mathrm{Im}\Sigma^{\rm
R}_{U\,m}(\omega)}{\left[\mathrm{Re}({\mathcal{G}_{m}^{\rm
R}}^{-1})(\omega)-\mathrm{Re}\Sigma^{\rm
R}_{U\,m}(\omega)\right]^{2}+\left[(\mathrm{Im}{\mathcal{G}_{m}^{\rm
R}}^{-1})(\omega)-\mathrm{Im}\Sigma^{\rm R}_{U\,m}(\omega)\right]^{2}}\,,$
(41)

and the local DOS as $\rho_{\rm
loc}(\omega)=\frac{1}{2}\sum_{m}\rho_{m}(\omega)$.

In the SOHM, this simplifies as
$\rho(\omega)=\rho_{B}(\omega)=\rho_{A}(\omega)$, which holds in the DHM at
$\omega=0$ due to the orbital permutation and particle-hole symmetries in Eq.
(B.1).

In the dimer limit and at $T=0$, we can estimate the density of states at
$\omega=0$ by working with the expression (41). In this limit, we can
explicitly write $\mathcal{G}_{m}^{\rm R}$ and $\Sigma^{\rm R}_{U\,m}$, using
the expressions given in App. B. The smallest energy scale of the system is
always $\Gamma$ (see the typical values in Table 1), and we therefore include
it as the regularizing term (replacing $\rm{i}0^{+}\to\rm{i}\Gamma$). In the
SOHM (atomic limit), we have:

$\displaystyle\mathrm{Re}\Sigma^{\rm R}_{U\,m}|_{\omega=0}=\mathrm{Re}($
$\displaystyle\mathcal{G}_{m}^{\rm R\,-1})|_{\omega=0}=0\,,$ (42)
$\displaystyle|\mathrm{Im}\Sigma^{\rm
R}_{U\,m}|_{\omega=0}\sim\frac{U^{2}}{\Gamma}\gg\mathrm{Im}($
$\displaystyle\mathcal{G}_{m}^{\rm R\,-1})|_{\omega=0}=\Gamma\,,$

and therefore the gap density is given by

$\displaystyle\rho(\omega=0)\sim\Gamma/U^{2}.$ (43)

For the DHM, the typical scales given in Table 1 lead us to work in the regime
$U\gg t^{\perp}\gg\Gamma$. In the dimer limit,
$\mathrm{Re}({\mathcal{G}_{B}^{\rm R}}^{-1})|_{\omega=0}=t^{\perp}$ and
$\mathrm{Im}({\mathcal{G}_{B}^{\rm R}}^{-1})|_{\omega=0}=\Gamma$. The real
part of the self-energy is given by $\mathrm{Re}\Sigma^{\rm
R}_{U\,B}|_{\omega=0}=-\frac{U^{2}}{4}\frac{3t^{\perp}}{(3t^{\perp})^{2}+\Gamma^{2}}\sim~{}-U^{2}/12t^{\perp}$,
and numerical computations of $\mathrm{Im}\Sigma^{\rm R}_{U\,B}|_{\omega=0}$
at finite $\Gamma$ show that $\mathrm{Im}\Sigma^{\rm
R}_{U\,B}|_{\omega=0}\ll\Gamma$. The latter term can therefore be neglected
from Eq. (41). Consequently, the density of states in the gap is given by

$\displaystyle\rho(\omega=0)\sim\frac{{t^{\perp}}^{2}\Gamma}{U^{4}}.$ (44)

Comparing $\rho(\omega=0)$ in the two models, we see that they differ by a
factor of $(U/t^{\perp})^{2}$, which accounts for the difference observed in
Fig. 3.

## References

  * Cox [2010] P. A. Cox, Transition metal oxides: an introduction to their electronic structure and properties, Oxford University Press, Oxford, UK (2010).
  * Maekawa  _et al._ [2004] S. Maekawa, T. Tohyama, S. E. Barnes, S. Ishihara, W. Koshibae, and G. Khaliullin, Physics of transition metal oxides, Springer Science & Business Media, Berlin, Germany (2004).
  * Khomskii [2014] D. Khomskii, Transition metal compounds, Cambridge University Press, Cambridge, UK (2014).
  * Sawa [2008] A. Sawa, Resistive switching in transition metal oxides, Mater. Today 11, 28 (2008).
  * Lee  _et al._ [2015] J. S. Lee, S. Lee, and T. W. Noh, Resistive switching phenomena: A review of statistical physics approaches, Appl. Phys. Rev. 2, 031303 (2015).
  * Cope and Penn [1968] R. Cope and A. Penn, High-speed solid-state thermal switches based on vanadium dioxide, J. Phys. D: Appl. Phys. 1, 161 (1968).
  * Lin  _et al._ [2018] J. Lin, S. Ramanathan, and S. Guha, Electrically driven insulator–metal transition-based devices—part ii: Transient characteristics, IEEE Trans. Electron Devices 65, 3989 (2018).
  * del Valle  _et al._ [2018] J. del Valle, J. G. Ramírez, M. J. Rozenberg, and I. K. Schuller, Challenges in materials and devices for resistive-switching-based neuromorphic computing, J. Appl. Phys. 124, 211101 (2018).
  * Hoffmann  _et al._ [2022] A. Hoffmann, S. Ramanathan, J. Grollier, A. D. Kent, M. J. Rozenberg, I. K. Schuller, O. G. Shpyrko, R. C. Dynes, Y. Fainman, A. Frano, _et al._ , Quantum materials for energy-efficient neuromorphic computing: Opportunities and challenges, APL Mater. 10, 070904 (2022).
  * Ridley [1963] B. K. Ridley, Specific negative resistance in solids, Proc. Phys. Soc. 82, 954 (1963).
  * Volkov and Kogan [1969] A. F. Volkov and S. M. Kogan, Physical phenomena in semiconductors with negative differential conductivity, Phys.-Usp. 11, 881 (1969).
  * Duchene  _et al._ [1971] J. Duchene, M. Terraillon, P. Pailly, and G. Adam, Filamentary conduction in vo2 coplanar thin‐film devices, Appl. Phys. Lett. 19, 115 (1971).
  * Stefanovich  _et al._ [2000] G. Stefanovich, A. Pergament, and D. Stefanovich, Electrical switching and Mott transition in VO2, J. Phys.: Condens. Matter 12, 8837 (2000).
  * Gopalakrishnan  _et al._ [2009] G. Gopalakrishnan, D. Ruzmetov, and S. Ramanathan, On the triggering mechanism for the metal–insulator transition in thin film vo2 devices: electric field versus thermal effects, J. Mater. Sci. 44, 5345 (2009).
  * Janod  _et al._ [2015] E. Janod, J. Tranchant, B. Corraze, M. Querré, P. Stoliar, M. Rozenberg, T. Cren, D. Roditchev, V. T. Phuoc, M.-P. Besland, and L. Cario, Resistive switching in Mott insulators and correlated systems, Adv. Funct. Mater. 25, 6287 (2015).
  * Higgins  _et al._ [1977] J. Higgins, B. Temple, and J. Lewis, Electrical properties of vanadate-glass threshold switches, J. Non-Cryst. Solids 23, 187 (1977).
  * Wu  _et al._ [2011] B. Wu, A. Zimmers, H. Aubin, R. Ghosh, Y. Liu, and R. Lopez, Electric-field-driven phase transition in vanadium dioxide, Phys. Rev. B 84, 241410 (2011).
  * Zimmers  _et al._ [2013] A. Zimmers, L. Aigouy, M. Mortier, A. Sharoni, S. Wang, K. G. West, J. G. Ramirez, and I. K. Schuller, Role of thermal heating on the voltage induced insulator-metal transition in VO2, Phys. Rev. Lett. 110, 056601 (2013).
  * Li  _et al._ [2016] D. Li, A. A. Sharma, D. K. Gala, N. Shukla, H. Paik, S. Datta, D. G. Schlom, J. A. Bain, and M. Skowronski, Joule heating-induced metal–insulator transition in epitaxial VO2/TiO2 devices, ACS Appl. Mater. Interfaces 8, 12908 (2016), pMID: 27136956\.
  * del Valle  _et al._ [2017] J. del Valle, Y. Kalcheim, J. Trastoy, A. Charnukha, D. N. Basov, and I. K. Schuller, Electrically induced multiple metal-insulator transitions in oxide nanodevices, Phys. Rev. Applied 8, 054041 (2017).
  * Polozov  _et al._ [2020] V. I. Polozov, S. S. Maklakov, A. L. Rakhmanov, S. A. Maklakov, and V. N. Kisel, Blow-up overheating instability in vanadium dioxide thin films, Phys. Rev. B 101, 214310 (2020).
  * Kalcheim  _et al._ [2020] Y. Kalcheim, A. Camjayi, J. del Valle, P. Salev, M. Rozenberg, and I. K. Schuller, Non-thermal resistive switching in mott insulator nanowires, Nature Communications 11, 2985 (2020).
  * del Valle  _et al._ [2021] J. del Valle, N. M. Vargas, R. Rocco, P. Salev, Y. Kalcheim, P. N. Lapa, C. Adda, M.-H. Lee, P. Y. Wang, L. Fratino, M. J. Rozenberg, and I. K. Schuller, Spatiotemporal characterization of the field-induced insulator-to-metal transition, Science 373, 907 (2021).
  * Sugimoto  _et al._ [2008] N. Sugimoto, S. Onoda, and N. Nagaosa, Field-induced metal-insulator transition and switching phenomenon in correlated insulators, Phys. Rev. B 78, 155104 (2008).
  * Han  _et al._ [2018] J. E. Han, J. Li, C. Aron, and G. Kotliar, Nonequilibrium mean-field theory of resistive phase transitions, Phys. Rev. B 98, 035145 (2018).
  * Mitra and Millis [2008] A. Mitra and A. J. Millis, Current-driven quantum criticality in itinerant electron ferromagnets, Phys. Rev. B 77, 220404 (2008).
  * Georges  _et al._ [1996] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996).
  * Freericks  _et al._ [2006] J. K. Freericks, V. M. Turkowski, and V. Zlatić, Nonequilibrium dynamical mean-field theory, Phys. Rev. Lett. 97, 266408 (2006).
  * Okamoto [2008] S. Okamoto, Nonlinear transport through strongly correlated two-terminal heterostructures: A dynamical mean-field approach, Phys. Rev. Lett. 101, 116807 (2008).
  * Aoki  _et al._ [2014] H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Nonequilibrium dynamical mean-field theory and its applications, Rev. Mod. Phys. 86, 779 (2014).
  * Joura  _et al._ [2008] A. V. Joura, J. K. Freericks, and T. Pruschke, Steady-state nonequilibrium density of states of driven strongly correlated lattice models in infinite dimensions, Phys. Rev. Lett. 101, 196401 (2008).
  * Tsuji  _et al._ [2008] N. Tsuji, T. Oka, and H. Aoki, Correlated electron systems periodically driven out of equilibrium: $\text{Floquet}+\text{DMFT}$ formalism, Phys. Rev. B 78, 235124 (2008).
  * Aron  _et al._ [2012] C. Aron, G. Kotliar, and C. Weber, Dimensional crossover driven by an electric field, Phys. Rev. Lett. 108, 086401 (2012).
  * Arrigoni  _et al._ [2013] E. Arrigoni, M. Knap, and W. von der Linden, Nonequilibrium dynamical mean-field theory: An auxiliary quantum master equation approach, Phys. Rev. Lett. 110, 086403 (2013).
  * Li  _et al._ [2015] J. Li, C. Aron, G. Kotliar, and J. E. Han, Electric-field-driven resistive switching in the dissipative Hubbard model, Phys. Rev. Lett. 114, 226403 (2015).
  * Eckstein  _et al._ [2010] M. Eckstein, T. Oka, and P. Werner, Dielectric breakdown of Mott insulators in dynamical mean-field theory, Phys. Rev. Lett. 105, 146404 (2010).
  * Aron [2012] C. Aron, Dielectric breakdown of a mott insulator, Phys. Rev. B 86, 085127 (2012).
  * Mazza  _et al._ [2015] G. Mazza, A. Amaricci, M. Capone, and M. Fabrizio, Electronic transport and dynamics in correlated heterostructures, Phys. Rev. B 91, 195124 (2015).
  * Li  _et al._ [2017] J. Li, C. Aron, G. Kotliar, and J. E. Han, Microscopic theory of resistive switching in ordered insulators: Electronic versus thermal mechanisms, Nano Lett. 17, 2994 (2017), pMID: 28394624\.
  * Chiriacò and Millis [2018] G. Chiriacò and A. J. Millis, Voltage-induced metal-insulator transition in a one-dimensional charge density wave, Phys. Rev. B 98, 205152 (2018).
  * Han  _et al._ [2023] J. E. Han, C. Aron, X. Chen, I. Mansaray, J.-H. Han, K.-S. Kim, M. Randle, and J. P. Bird, Correlated insulator collapse due to quantum avalanche via in-gap ladder states, Nature Communications 14, 2936 (2023).
  * Mazzocchi  _et al._ [2022] T. M. Mazzocchi, P. Gazzaneo, J. Lotze, and E. Arrigoni, Correlated Mott insulators in strong electric fields: Role of phonons in heat dissipation (2022).
  * Dasari  _et al._ [2020] N. Dasari, J. Li, P. Werner, and M. Eckstein, Revealing Hund’s multiplets in Mott insulators under strong electric fields, Phys. Rev. B 101, 161107 (2020).
  * Grandi and Eckstein [2021] F. Grandi and M. Eckstein, Ultrafast metal-to-insulator switching in a strongly correlated system, arXiv:2104.03644 (2021).
  * Nájera  _et al._ [2018] O. Nájera, M. Civelli, V. Dobrosavljević, and M. J. Rozenberg, Multiple crossovers and coherent states in a mott-peierls insulator, Phys. Rev. B 97, 045108 (2018).
  * Lu and Robertson [2019] H. Lu and J. Robertson, Density functional theory studies of the metal–insulator transition in vanadium dioxide alloys, Phys. Status Solidi B 256 (2019).
  * Biermann  _et al._ [2005] S. Biermann, A. Poteryaev, A. I. Lichtenstein, and A. Georges, Dynamical singlets and correlation-assisted peierls transition in ${\mathrm{v}\mathrm{o}}_{2}$, Phys. Rev. Lett. 94, 026404 (2005).
  * Shao  _et al._ [2018] Z. Shao, X. Cao, H. Luo, and P. Jin, Recent progress in the phase-transition mechanism and modulation of vanadium dioxide materials, NPG Asia Mater. 10, 581 (2018).
  * Sommers and Doniach [1978] C. Sommers and S. Doniach, First principles calculation of the intra-atomic correlation energy in vo2, Solid State Commun. 28, 133 (1978).
  * Nájera  _et al._ [2017] O. Nájera, M. Civelli, V. Dobrosavljević, and M. J. Rozenberg, Resolving the VO2 controversy: Mott mechanism dominates the insulator-to-metal transition, Phys. Rev. B 95, 035113 (2017).
  * Tomczak  _et al._ [2008] J. M. Tomczak, F. Aryasetiawan, and S. Biermann, Effective bandstructure in the insulating phase versus strong dynamical correlations in metallic VO2, Phys. Rev. B 78, 115103 (2008).
  * Eyert [2011] V. Eyert, VO2: A novel view from band theory, Phys. Rev. Lett. 107, 016401 (2011).
  * Weber  _et al._ [2012] C. Weber, D. D. O’Regan, N. D. M. Hine, M. C. Payne, G. Kotliar, and P. B. Littlewood, Vanadium dioxide: A Peierls-Mott insulator stable against disorder, Phys. Rev. Lett. 108, 256402 (2012).
  * Brito  _et al._ [2016] W. H. Brito, M. C. O. Aguiar, K. Haule, and G. Kotliar, Metal-insulator transition in VO2: A $\mathrm{DFT}+\mathrm{DMFT}$ perspective, Phys. Rev. Lett. 117, 056402 (2016).
  * dos Santos [1995] R. R. dos Santos, Magnetism and pairing in Hubbard bilayers, Phys. Rev. B 51, 15540 (1995).
  * Kancharla and Okamoto [2007] S. S. Kancharla and S. Okamoto, Band insulator to Mott insulator transition in a bilayer Hubbard model, Phys. Rev. B 75, 193103 (2007).
  * Rüger  _et al._ [2014] R. Rüger, L. F. Tocchio, R. Valentí, and C. Gros, The phase diagram of the square lattice bilayer Hubbard model: a variational Monte Carlo study, New J. Phys. 16, 033010 (2014).
  * Xu  _et al._ [2022] Y. Xu, K. Kang, K. Watanabe, T. Taniguchi, K. F. Mak, and J. Shan, A tunable bilayer Hubbard model in twisted WSe2, Nat. Nanotechnol. , 1 (2022).
  * Büttiker [1985] M. Büttiker, Small normal-metal loop coupled to an electron reservoir, Phys. Rev. B 32, 1846 (1985).
  * Tsuji  _et al._ [2009] N. Tsuji, T. Oka, and H. Aoki, Nonequilibrium steady state of photoexcited correlated electrons in the presence of dissipation, Phys. Rev. Lett. 103, 047403 (2009).
  * Amaricci  _et al._ [2012] A. Amaricci, C. Weber, M. Capone, and G. Kotliar, Approach to a stationary state in a driven Hubbard model coupled to a thermostat, Phys. Rev. B 86, 085110 (2012).
  * Han [2013] J. E. Han, Solution of electric-field-driven tight-binding lattice coupled to fermion reservoirs, Phys. Rev. B 87, 085119 (2013).
  * Han and Li [2013] J. E. Han and J. Li, Energy dissipation in a dc-field-driven electron lattice coupled to fermion baths, Phys. Rev. B 88, 075113 (2013).
  * Murakami and Werner [2018] Y. Murakami and P. Werner, Nonequilibrium steady states of electric field driven Mott insulators, Phys. Rev. B 98, 075102 (2018).
  * Moeller  _et al._ [1999] G. Moeller, V. Dobrosavljević, and A. E. Ruckenstein, Rkky interactions and the Mott transition, Phys. Rev. B 59, 6846 (1999).
  * Monien  _et al._ [1997] H. Monien, N. Elstner, and A. J. Millis, Possible explanation for the absence of bilayer splitting in ybco (1997).
  * Note [1] This difference in the nature of the insulating state is analogous to the difference between the insulating states of the SOHM computed by means of single-site DMFT versus cluster-DMFT. In the latter, short-range correlations are kept such as to allow for a lower entropy insulating ground state.
  * Park  _et al._ [2008] H. Park, K. Haule, and G. Kotliar, Cluster dynamical mean-field theory of the Mott transition, Phys. Rev. Lett. 101, 186403 (2008).
  * Biroli  _et al._ [2004] G. Biroli, O. Parcollet, and G. Kotliar, Cluster dynamical mean-field theories: Causality and classical limit, Phys. Rev. B 69, 205108 (2004).
  * Maier  _et al._ [2005] T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Quantum cluster theories, Rev. Mod. Phys. 77, 1027 (2005).
  * Tsuji  _et al._ [2014] N. Tsuji, P. Barmettler, H. Aoki, and P. Werner, Nonequilibrium dynamical cluster theory, Phys. Rev. B 90, 075117 (2014).
  * Georges and Kotliar [1992] A. Georges and G. Kotliar, Hubbard model in infinite dimensions, Phys. Rev. B 45, 6479 (1992).
  * Zhang  _et al._ [1993] X. Y. Zhang, M. J. Rozenberg, and G. Kotliar, Mott transition in the $d$=$\infty$ Hubbard model at zero temperature, Phys. Rev. Lett. 70, 1666 (1993).
  * Davies and Wilkins [1988] J. H. Davies and J. W. Wilkins, Narrow electronic bands in high electric fields: Static properties, Phys. Rev. B 38, 1667 (1988).
  * Aron  _et al._ [2013] C. Aron, C. Weber, and G. Kotliar, Impurity model for non-equilibrium steady states, Phys. Rev. B 87, 125113 (2013).
  * Bakalov  _et al._ [2016] P. Bakalov, D. Nasr Esfahani, L. Covaci, F. M. Peeters, J. Tempere, and J.-P. Locquet, Electric-field-driven mott metal-insulator transition in correlated thin films: An inhomogeneous dynamical mean-field theory approach, Phys. Rev. B 93, 165112 (2016).
  * Okamoto [2007] S. Okamoto, Nonequilibrium transport and optical properties of model metal–mott-insulator–metal heterostructures, Phys. Rev. B 76, 035105 (2007).
  * Potthoff and Nolting [1999a] M. Potthoff and W. Nolting, Metallic surface of a mott insulator–mott insulating surface of a metal, Phys. Rev. B 60, 7834 (1999a).
  * Potthoff and Nolting [1999b] M. Potthoff and W. Nolting, Surface metal-insulator transition in the hubbard model, Phys. Rev. B 59, 2549 (1999b).
  * Kotliar, G. [1999] Kotliar, G., Landau theory of the Mott transition in the fully frustrated Hubbard model in infinite dimensions, Eur. Phys. J. B 11, 27 (1999).
  * Kotliar  _et al._ [2000] G. Kotliar, E. Lange, and M. J. Rozenberg, Landau theory of the finite temperature Mott transition, Phys. Rev. Lett. 84, 5180 (2000).
  * Chitra and Kotliar [2001] R. Chitra and G. Kotliar, Effective-action approach to strongly correlated fermion systems, Phys. Rev. B 63, 115110 (2001).
  * Kotliar [2021] G. Kotliar, 14 – The Mott Transition, in  _More is Different, Fifty Years of Condensed Matter Physics_, edited by N.-P. Ong and R. Bhatt (Princeton University Press, Princeton, 2021) pp. 211–236.
  * Díaz and Aron [2022] M. I. Díaz and C. Aron, Effective field theory of electric field-driven resistive switching (2022), in preparation.
  * Stoliar  _et al._ [2013] P. Stoliar, L. Cario, E. Janod, B. Corraze, C. Guillot-Deudon, S. Salmon-Bourmand, V. Guiot, J. Tranchant, and M. Rozenberg, Universal electric-field-driven resistive transition in narrow-gap Mott insulators, Adv. Mater. 25, 3222 (2013).
  * Adda  _et al._ [2022] C. Adda, M.-H. Lee, Y. Kalcheim, P. Salev, R. Rocco, N. M. Vargas, N. Ghazikhanian, C.-P. Li, G. Albright, M. Rozenberg, and I. K. Schuller, Direct observation of the electrically triggered insulator-metal transition in ${\mathrm{v}}_{3}{\mathrm{o}}_{5}$ far below the transition temperature, Phys. Rev. X 12, 011025 (2022).