Iterated Perturbation Theory for Mott Insulators in a Static Electric Field with Optical Phonons

This manuscript aims to compare the so-called iterated perturbation theory (IPT) and auxiliary master equation approach (AMEA) impurity solvers for a Mott insulating system driven out of equilibrium by a static electric field. Electronic heat bath and optical phonons are employed as dissipation mechanism of the current-induced Joule heat that the excited electrons of the lattice experience as the result of the field's driving. Despite its simplicity, the IPT approach yields results which qualitatively are in good agreement with those obtained within the AMEA impurity solver, although fails to reproduce some correlation effects.


Introduction
2][3][4][5][6] Due to this property, they could be used to model the so-called insulator-to-metal transition, which has been investigated theoretically in the seminal works 7,8 and then observed experimentally. 9,10 e to their strongly interacting nature, Mott insulating systems require a nonperturbative method to be dealt with.Nowadays, the most well-established approach is given by the dynamical mean-field theory (DMFT) [11][12][13][14][15][16] which holds under both equilibrium and non-equilibrium conditions.The DMFT relies on impurity solvers to address the non-equilibrium steady-state of the system, which often are the bottleneck of the approach as they can be computationally costly.9][20][21][22] The IPT solver yields results which are in quite good agreement with those obtained by employing AMEA for the very same setup investigated in a previous work from the authors. 23The rest of the manuscript is organized as follows: In Section 2 we introduced the model at hand while Section 3 will be devoted to a short recap of the Green's function formalism and the IPT impurity solver.Results are discussed in Section 4 while Section 5 is left for conclusions and further considerations.The details of the physical setup under investigation and the derivation of all the relevant observables of interest can be found in the previous work from the authors. 23

Model Hamiltonian
We study the single-band Hubbard model in the presence of a constant electric field in the temporal gauge, 23 the Hamiltonian of which is given by Ĥ(t) = ĤU (t) + Ĥbath + Ĥe-ph + Ĥph . ( The units employed in this manuscript are such that the lattice spacing a, Planck constant ℏ and the electron charge q are chosen as a = ℏ = −q = 1, hence the Hubbard Hamiltonian ĤU (t) in Equation ( 1) is given by ĤU arXiv:2311.00617v2 [cond-mat.str-el]25 Mar 2024 where f † iσ ( fiσ ) is the creation (annihilation) operator of an electron of spin σ = {↑, ↓} at the i-th lattice site and nf iσ ≡ f † iσ fiσ the corresponding density operator.Sums over nearest neighbor sites are denoted by (i, j).The electrons' onsite energy is chosen as ε c ≡ −U/2 for the system to be fulfil particle-hole symmetry and t c is the bare hopping amplitude.The homogeneous vector potential A(t) = −F t is chosen such that the static electric field is constant and oriented along the body diagonal of a hypercubic lattice e 0 = (1, 1, . . ., 1) and is given by F = −∂ t A(t).We take the infinite-dimension limit, 6,24 i.e. d → ∞, with the usual rescaling of the hopping t c = t * /(2 √ d), which allows to perform summations over the electron crystal momentum using the joint density of states 25,26 ρ(ϵ, ϵ) = 1/(πt To stabilize the DMFT loop we include electronic heat baths, consisting of a collection of noninteracting fermionic degrees of freedom, coupled locally to each lattice site which are described by the Hamiltonian Ĥbath , the details of which will be specified in Section 3.1.1,see Equation ( 9).

Electron Dyson equation
The interacting electron GF obeys the Dyson equation with G 0 denoting the GF corresponding to the noninteracting part of the Hamiltonian in Equation (2).By Σ bath , Σ and Σ e-ph we denote the fermionic heat bath, electron and electron-phonon (e-ph) self-energy (SE), respectively.By means of the DMFT and Migdal 35,36 approximations both Σ and Σ e-ph in Equation (4) are local.
Every quantity X with an underline denotes the so-called Keldysh structure, namely with X R,A,K being the retarded, advanced and Keldysh components where X K ≡ X > + X < and X ≷ being the greater and lesser components.39] However, given the time-translation invariant character of the problem at hand, 23,25,26,38 only the diagonal components and especially the time-averaged element X 00 are nonvanishing 1 .For this reason we will omit the Floquet indices in the rest of the manuscript.Details concerning the computation of the e-ph SE Σ e-ph within the Migdal approximation can be found in our previous work. 23,26  recall the definition of the electron spectral function (SF) where the local electron GF is given by being the inverse time-averaged retarded component of the GF in Equation ( 4).
In terms of the contour-times z, z ′ , and in the Migdal approximation, 23,26,35,36 the e-ph SE reads and corresponds to the lowest-order expansion in the phonon GF D ph , the form of which will be discussed below.G loc (z, z ′ ) is the contour-times local GF, the retarded component of which obeys Equation (7).The retarded and Keldysh components of Equation ( 8) can be found in previous work from the authors. 26n this manuscript we will use the wide-band limit approximation 40 for the heat bath described by Ĥbath where Γ e is the so-called electronic dephasing rate. 23The Keldysh component Σ K bath is obtained by the fluctuation-dissipation theorem for fermions, i.e.Σ K bath (ω ] with β the inverse temperature and µ the chemical potential of the bath.Other important observables are the current J flowing in the direction of the applied field and the kinetic energy E kin of the electrons of the lattice: the derivation of both these quantities can be found in previous works from the authors. 23,26

Phonon Dyson equation
The optical phonon branch consists of Einstein phonons coupled to an ohmic bath, 23,36 the Dyson equation of which reads with the non-interacting retarded component of the Einstein phonon GF given by The Einstein phonon is coupled to an ohmic bath Ĥph,ohm , the real retarded GF of which is obtained from the Kramers-Krönig relations, 36 while the Keldysh component is given by the fluctuation-dissipation theorem for bosons, i.e.
We choose the following form for the ohmic bath density of states (DOS) in ( 12) where In Equation ( 13) ω c denotes the ohmic bath cutoff frequency and v c the hybridization strength to the ohmic bath: we only stress that the DOS in Equation ( 13) is linear for According to the DMFT approximation, the polarization diagram Π e-ph only depends on the local electron GF.Within the Migdal approximation, the contour times polarization diagram 35,36 in Equation ( 10) reads where the factor 2 accounts for spin degeneracy.We denote the scheme in which Π e-ph (z, z ′ ) is set to zero as non-self-consistent (NSC) while within the self-consistent (SC) treatment the phonon SE in Equation ( 14) is non-vanishing.The real time components of Equation ( 14) have been derived in previous work from the authors. 23

IPT impurity solver and DMFT loop
The IPT impurity solver is based on the perturbative expansion of the electron SE in terms of the Weiss field, 17 namely g 0 being the non-interacting GF of the isolated site, the retarded component of which (in frequency domain) reads g R 0 (ω) = (ω − ε c ) −1 .2][13] The real-time components of the electron SE can then be written as The retarded and Keldysh components of the electron SE can be obtained from Equation ( 16) as At the steady-state all quantities in Equation ( 17) are dependent on the difference t − t ′ alone, hence their Fourier transform in frequency domain is straightforward.The retarded and Keldysh components of the Weiss field then read We stress that at half-filling the Hartree term U/2 must be explicitly added to Σ R (ω) before computing the quantities in Equation (18).The main steps of DMFT employing the IPT as impurity solver are: i. Guess Σ(ω), Σ e-ph (ω) and Π e-ph (ω) ii. Compute G loc (ω) and D ph (ω) as in Equation ( 7) and (10)   iii.Extract e-ph (t, t ′ ) vi. Update Σ(ω), Σ e-ph (ω) and Π e-ph (ω) via back Fourier-transform.
The steps ii. to vi. are then repeated until convergence.
g ω E 8 -4 0 0.05 0.6 0.055 0.4 0.6 Table 1: Default values of the main parameters used in this manuscript.All values are in units of the hopping t * .

Results
In this section we benchmark the results yielded by the IPT approach against those obtained within the AMEA impurity solver for the system at hand.

SC and NSC phonons
We set off by analyzing the steady-state current and kinetic energy in both the SC and NSC schemes obtained within the IPT impurity solver.We then compare these results with those obtained within the AMEA only at a later time.

Steady-state current and kinetic energy
The current J and and kinetic energy E kin as functions of the applied field F for both the SC and NSC schemes obtained employing the IPT impurity solver are shown in Figure 1.We observe that the SC  treatment has no effect for field strengths F < U/2 as the SC and NSC curves for both J and E kin lie on top of each other.On the other hand, the SC scheme does contribute corrections, albeit tiny, once F gets past field strengths of the order of the resonance U/2, for both J and E kin .In particular we observe a suppression of the peak values of the current J in the SC approach at F ≈ U with respect to the NSC treatment.This is accompanied by an overall smearing of the J-F curves around the maximum for all the values of Γ e employed in this manuscript, see Figure 1(a), 1(b) and 1(c).The kinetic energy is also affected by the SC treatment as its maximum value is suppressed, while its minimum is raised with respect to the NSC scheme regardless of the value of Γ e , see Figure 1(d), 1(e) and 1(f).This behavior is qualitatively in agreement with the results obtained within the AMEA impurity solver presented in a recent work from the authors. 23As a matter of fact, using AMEA the drop in the current J at the resonance F ≈ U observed in the SC scheme is way more pronounced than in the IPT case as one can see by comparing Figure 1 and 2 for corresponding values of the dephasing rate Γ e .On the other hand, the differences between the impurity solvers are less pronounced when it comes to the kinetic energy E kin , as it can be observed by comparing the corresponding curves in Figure 1 and 2 for the same Γ e 's.
We want to mention that the orders of magnitude of both J and E kin in the IPT scheme quantitatively agrees with those obtained within the AMEA impurity solver.However, using the IPT approach the two main resonances at F ≈ U/2 and F ≈ U are shifted towards smaller field strengths and the differences between the SC and NSC schemes are not as pronounced as in the AMEA scheme [see again Figure 2].In addition, the two impurity solvers differ in that the IPT cannot capture correlation effects like the tiny resonance in the current J at F ∼ U/3 ≈ 2.62 , which instead can be distinguished using the AMEA solver as soon as Γ e is small enough, see Figure 2(a), 2(d) and in particular 2(g).

The effect of the Hubbard U
This section is devoted to the analysis of the role of the Hubbard U on the insulating phase: once again we focus on the current J and kinetic energy E kin in the SC and NSC approaches to then compare them to the corresponding quantities obtained within the AMEA impurity solver.
In Figure 3 is displayed the behavior of J and E kin as function of the applied field for selected values of U , see panels (a) and (b) for the NSC scheme and (c) and (d) for SC treatment.We observe that in the SC case both the J-F and E kin -F curves are broadened, see panels (b) and (d).
In particular we observe an increase in the current J at F ≈ U/2 as U is decreased for both NSC and SC schemes.Also, while the values of the two peaks at U/2 and U stay approximately the same in height in the NSC and SC schemes, we find that the latter treatment enhances J for field strengths that lie in between the two main resonances, compare Figure 3(a) and 3(c).When it comes to the kinetic energy, on the other hand, let alone an overall smearing of the curves the SC treatment does not affect E kin significantly, compare Figure 3(b) and 3(d).As already discussed in Reference 1 both the J-F and E kin -F curves collapse on one another when they are plotted as function of the difference F − U (not shown), signalling that the breakdown of the insulating phase depends on the value of the Hubbard U alone.Once again we stress the qualitative agreement between the results of Figure 3 and those presented in the previous work 23 obtained within AMEA that we reproduce in Figure 4.However, the two solvers do differ in that, as it can be shown with the help of Figure 1(c) and 2(g), the IPT approach fails to capture the resonance in J at F ∼ U/33 .We see that such resonance is missing within the IPT solver in both the SC and NSC schemes regardless of the value of U , see Figure 3(a) and 3(c).On the other hand, when using the AMEA one can appreciate it, even though it tends to be smeared out by the SC treatment, especially for small values of U , see Figure 4(a) and 4(d).Furthermore, the enhancement of the current J in between the two main resonances F = U/2 and U occurring in the SC treatment within the IPT approach [see again Figure 3(a) and 3(c)] is even more pronounced when the AMEA impurity solver is used, especially for small values of U , as one can see by comparing the curves for the current J in Figure 4(a) and 4(d).
The main difference in the kinetic energy between the ITP and AMEA impurity solvers lies in the overall sharper E kin -F curves in the latter approach, as one can see by comparing the results shown in Figure 3 and 4. It is worth mentioning that the SC treatment smears out the E kin -F curves more when the AMEA solver is employed especially for field strengths F ∼ U/2.

Equilibrium spectral features
In this section we briefly compare the electronic spectral features at equilibrium, i.e.F = 0.
Figure 5 shows the electron SF at F = 0 within the (a) NSC and (b) SC schemes for both the IPT and AMEA impurity solvers.The IPT-resolved SF shows a much more pronounced quasi-particle peak (QPP) at around ω = 0 26 as U is reduced, together with an underestimation of the width of the Hubbard bands within both the NSC and SC schemes, as compared to the results obtained with the AMEA impurity solver.On the other hand, at U = 8t * the height of the QPP in the IPT and AMEA approaches is in quite good agreement, even though the Hubbard bands are still narrower in the IPT scheme in both the NSC and SC approaches.
The larger amount of in-gap states due to the QPP at ω ≈ 0 in the IPT with respect to AMEA is common to all electric field strengths used in this manuscript (not shown) and reveals the systematic underestimation of the band gap committed by the IPT solver.This underestimation, in turn, explains the shifting of the main resonances (at F = U/2 and U ) in the current J toward smaller values of the electric field discussed in Section 4.1.1.

The threshold field
In models of Mott insulating systems one expects the J-F curve to display the threshold behavior ?, 32 J ∝ F exp(−F th /F ), (19)   due to the opening of a gap that prevents the free motion of electrons from the lower to the upper band.Equation (19) shows that only when the applied field F gets past the threshold F th can electrons cross the band gap and thus give rise to a steady-state current.By determining the threshold F th one can then infer the magnitude of the effective4 band gap of the model.However, as it has been shown in previous works investigating Mott insulators with a large5 band gap, 6,26 correlation effects are responsible for resonances in the current J at F ≈ U/n, see also Section 4.1 [Figure 1 and 2], so we may expect the existence of at least two threshold fields.
It is worth recalling that these resonances are determined by emergence of the Wannier-Stark6, 26 states in the local electron SF which effectively allow electron tunnelling to the upper band by bridging the band gap and are sometimes referred to as Landau-Zener excitations in other models. 1 However, in this manuscript, F th identifies the field strength necessary for the current J to reach nonnegligible values for the first time, so Equation ( 19) should hold for F < U/2 at least for the IPT solver as it cannot capture the whole of electronic correlation 6 , i.e. there is no other resonance before field strengths of the order of half of the band gap.This section is then devoted to the analysis of these aspects limited to the comparison of the IPT and AMEA solvers.In fact, given the simple model at hand the study of the dependence of the effective band gap on the threshold F th goes beyond the purpose of this manuscript: for a detailed study can be found in the previous work.?, ?
We benchmark the results for the J-F curve obtained within the IPT and AMEA impurity solvers for selected values of U against a linear regression fit according to Equation (19).

Conclusion
In this manuscript a Mott insulating system has been characterized in terms of its conducting properties when subject to an external static electric.Optical phonons and electronic heat bath provide the relaxation pathways for the extra energy injected by the field, so that the electron of the lattice can relax back to the valence band and a steady-state current be established.The iterated perturbation theory (IPT) approach has been used as impurity solver to address the steadystate of the system.The corresponding results have been benchmarked against those obtained within a much more computationally costly impurity solver developed by the authors, the so-called auxiliary master equation approach (AMEA).It has been shown that the results obtained employing the IPT qualitatively agree with those of the AMEA impurity solver even though the former approach does not capture the resonances in the current characteristics which are directly related to the correlated nature of the electrons of the lattice.Being computationally cheaper, the IPT solver could be used to span the parameter space when investigating novel setups and gather information about the interesting regions to be addressed by a more reliable and computationally costly impurity solver.The SC treatment has no effect on F th , which is the reason why panel (c) does not distinguish between the two cases.Default parameters are specified in Table 1.(Here Γ e = 0.12t * .)

Energy
Schematic representation of the current flowing within the lattice.When the field is applied electrons can tunnel through the former gap of an adjacent site, creating the necessary states to populate the conduction band of a site which is twice lattice spacings apart.
* 2 ) exp[−(ϵ 2 + ϵ 2 )/t * 2 ] with ϵ = −2t c d i=1 cos k i and ϵ = −2t c d i=1 sin k i .An optical phonon branch is attached to each lattice site by means of the Hamiltonian Ĥe-ph = g iσ nf iσ xi (3) with xi ≡ ( b † i + bi )/ √ 2, where b † i ( bi ) creates (annihilates) an optical phonon with energy ω E at the lattice site i.The optical phonon Hamiltonian consists of an Einstein phonon Ĥph,E = ω E i nb i with nb i = b † i bi the phonon density, coupled to a noninteracting, ohmic bath Ĥph,ohm with spectral density given in Equation (13).

Figure 1 :
Figure 1: Current J as function of the applied field F for (a) Γ e = 0.20t * (b) Γ e = 0.16t * and (c) Γ e = 0.12t * for both SC and NSC schemes.(d), (e) and (f) show the kinetic energy E kin as function of F for the same coupling strengths, respectively.Default parameters are specified in Table 1.Results have been obtained using the IPT impurity solver.(Here U = 8t * .)

Figure 6 (
a) shows the ratio J/F as function of 1/F in both SC and NSC schemes obtained using the IPT: the values of the inverse field used for the fit can be deduced by the extent of the black line.As one can see by direct inspection the curves for all values of U exhibit a linear behavior for a wide range of values of inverse field strengths up to F −1 ∼ 0.25t * −1 , which correspond to the resonance at F = U/2.The results obtained with the AMEA impurity solver, instead, are shown in Figure6(b), see corresponding inset for a close-up of the region around F −1 ∼ 0.4t * −1 , corresponding to F ≈ U/3.The most evident feature is the small kink occurring at F −1 ∼ 0.4t * −1 , clearly visible especially for U > 5t * , which is absent in the IPT approach as it can be seen by comparing Figure6(a) and 6(b).Due to this additional resonance, the linear regression fit within the AMEA impurity solver can be performed over a far smaller range of values of F −1 and indeed there are two regions where linearity holds: the first occurs for F −1 < 0.4t * −1 and the second for inverse field smaller than 0.25t * −1 , see again Figure6(b).The threshold fields obtained with the linear regression fit are shown in Figure6(c): the SC treatment basically does not affect the results within the numerical accuracy, thus leaving the threshold fields unaltered in both the IPT and AMEA approaches.It should be noted that by performing the linear regression fit for F −1 < 0.4t * −1 as in this case, the IPT-and AMEA-resolved curves are basically on top of one another.The AMEA scheme, however, yields a slightly larger threshold F th than the IPT one for U = 8t * , see again Figure6(c).As expected, a larger F th is required to compensate for a larger band gap (and hence a larger U ) and promote particles across it.However, due to the extension of the Hubbard bands the effective gap, and thus the threshold field, is way smaller than the Hubbard U so that a naive relation of the form F th ∼ U does not hold (see Figure6(c)), as already argued in previous work ? .

Figure 2 :
Figure 2: (a) Current, (b) double occupation and (c) kinetic energy as function of the applied field F at Γ e = 0.20t * for both SC and NSC schemes.(d), (e) and (f) show the same quantities for Γ = 0.16t* , while (g), (h) and (i) refer to Γ e = 0.12t * .Default parameters are specified in Table 1.These results have been obtained using the AMEA impurity solver and originally published by the authors in https://link.aps.org/doi/10.1103/PhysRevB.107.155103.(Here U = 8t * .)

Figure 3 :
Figure 3: Current J within the (a) NSC and (c) SC scheme for selected values of the Hubbard U as function of the applied field F .(b) and (d) show the kinetic energy E kin as function of F in the NSC and SC scheme, respectively, for the same values of the Hubbard U .Default parameters are specified in Table1.Results have been obtained using the IPT impurity solver.(Here Γ e = 0.12t * .)

Figure 4 :Figure 5 :Figure 6 :
Figure 4: Current J within (a) NSC and (d) SC scheme for selected values of the Hubbard U as function of the applied field F .Black arrows in (a) highlight the progressive merge of the resonances at F ∼ U/3 and F = U/2 as U is lowered which is enhanced by the SC treatment.(b) and (e) show the double occupation d as function of F for the NSC and SC scheme, respectively, while the kinetic energy E kin is shown in (c) and (f) for the same values of U .Default parameters are specified in Table1.These results have been obtained using the AMEA impurity solver and originally published by the authors in https://link.aps.org/doi/10.1103/PhysRevB.107.155103.(Here Γ e = 0.12t * .) . Default parameters are specified in Table1.These results have been obtained using the AMEA impurity solver and originally published by the authors in https://link.aps.org/doi/10.1103/PhysRevB.107.155103.(Here U = 8t * .) * *

Table of Contents t
* t * site i site (i + 1) site (i + 2)