Gapless regime in the charge density wave phase of the finite dimensional Falicov-Kimball model

TL;DR

This paper investigates the gapless regime in the charge density wave phase of the finite-dimensional Falicov-Kimball model, revealing how localized states and subgap density of states influence charge transport.

Contribution

It demonstrates the existence of a stable gapless phase in the charge density wave regime and links it to localized single-particle excitations and subgap states.

Findings

Presence of a stable gapless phase with high density of states at the Fermi level.

Strong localization indicated by inverse participation ratio analysis.

Enhanced charge transport in heterostructures with the Falicov-Kimball system.

Abstract

The ground-state density of states of the half-filled Falicov-Kimball model contains a charge-density-wave gap. At finite temperature, this gap is not immediately closed, but is rather filled in by subgap states. For a specific combination of parameters, this leads to a stable phase where the system is in an ordered charge-density-wave phase, but there is high density of states at the Fermi level. We show that this property can be, in finite dimensions, traced to a crossing of sharp states resulting from the single particle excitations of the localized subsystem. The analysis of the inverse participation ratio points to a strong localization in the discussed regime. However, the pronounced subgap density of states can still lead to a notable increase of charge transport through a finite size system. We show this by focusing on the transmission in heterostructures where a Falicov-Kimball system is sandwiched between two metallic leads.

Equations

$$H_{\mathrm{FK}}$$

$$H_w=\sum _{\mathbf{l},{\mathbf{l}}^{\prime }}{h}_{{\mathbf{l}\mathbf{l}}^{\prime }}{d}_{\mathbf{l}}^{†}{d}_{{\mathbf{l}}^{\prime }}^{}-\frac{U}{2}{N}_f=\sum _{\alpha }{\epsilon }_{\alpha }{\tilde{d}}_{\alpha }^{†}{\tilde{d}}_{\alpha }^{}-\frac{U}{2}{N}_f.$$

$$⟨\hat{O}⟩={\mathrm{Tr}}_w{⟨\hat{O}⟩}_d\equiv \frac{1}{Z}\sum _w{e}^{-\beta F\left(w\right)}{⟨\hat{O}⟩}_d,$$

$$F\left(w\right)=-\frac{1}{\beta }\sum _{\alpha }\ln \left[1+e^{-\beta {\epsilon }_{\alpha }}\right]-\frac{U}{2}{N}_f,$$

$$\mathrm{DOS}(\epsilon )=\frac{1}{L}{\mathrm{Tr}}_w\sum _{\alpha }\delta \left(\epsilon -{\epsilon }_{\alpha }\right)$$

$$\mathrm{IPR}(\epsilon )={\mathrm{Tr}}_w\frac{{\sum }_i{\mathrm{DOS}}_i(\epsilon ,w{)}^2}{\mathrm{DOS}(\epsilon ,w{)}^2},$$

$$H_{\mathrm{lead}}^l$$

$$H_{\mathrm{hyb}}^l$$

$${\rho }_l\left(\epsilon \right)=\frac{2}{\pi B^2}\sqrt{B^2-{\left(\epsilon -{ϵ}_l\right)}^2},$$

$${\mathbf{G}}^{r,a}\left(\epsilon \right)$$

$${\mathbf{G}}^{<}\left(\epsilon \right)$$

$$g_{\alpha \beta }^{r,a}\left(\epsilon \right)=\frac{{\delta }_{\alpha \beta }}{\epsilon -{\epsilon }_{\alpha }\pm i0}.$$

$${\Sigma }_{l,\alpha \beta }^{r,a}(\epsilon )$$

$${\Sigma }_{l,\alpha \beta }^{<}(\epsilon )$$

$${\Gamma }_{l,\alpha \beta }(\epsilon )$$

$${\Lambda }_{l,\alpha \beta }(\epsilon )$$

$${\mathcal{U}}_{\alpha \beta }^{\{s^l\}}$$

$${\Theta }^w\left(\epsilon \right)=\mathrm{T}{\mathrm{r}}_{\mathrm{d}}\left[{\mathbf{\Gamma }}_L\left(\epsilon \right){\mathbf{G}}^r\left(\epsilon \right){\mathbf{\Gamma }}_R\left(\epsilon \right){\mathbf{G}}^a\left(\epsilon \right)\right],$$

$$\mathrm{LDOS}{\mathrm{h}}_i(\epsilon ,w)$$

$$\mathrm{gIPR}(E)={\mathrm{Tr}}_w\frac{{\sum }_i\mathrm{LDOS}{\mathrm{h}}_i^2(\epsilon ,w)}{{\mathrm{DOSh}}^2(\epsilon ,w)},$$

$$H_{\mathrm{FK}}^{\mathrm{CDW}}$$

$$d_{\mathbf{k}}^{†}=\frac{1}{\sqrt{L}}\sum _{\mathbf{l}}{e}^{i\mathbf{k}\cdot \mathbf{l}}{d}_{\mathbf{l}}^{†},$$

$$H_{\mathrm{FK}}^{\mathrm{CDW}}=\sum _{\mathbf{k}}{ϵ}_{\mathbf{k}}{d}_{\mathbf{k}}^{†}{d}_{\mathbf{k}}^{}+\frac{U}{2}\sum _{\mathbf{k}}{d}_{\mathbf{k}}^{†}{d}_{\mathbf{k}+\mathbf{\pi }}^{},$$

$${ϵ}_{\mathbf{k}}=-2t\sum _{i=1}^D\cos \left(k_i\right).$$

$$H_{\mathrm{FK}}^{CDW}$$

$${\epsilon }_{\mathbf{k}}^{\pm }=\pm \sqrt{{ϵ}_{\mathbf{k}}^2+{\left(\frac{U}{2}\right)}^2}.$$

$$Q_{\boldsymbol{k}}=\left(\begin{array}[]{cc}\frac{\epsilon_{\mathbf{\boldsymbol{k}}}+\varepsilon_{\mathbf{\boldsymbol{k}}}^{+}}{\sqrt{\left(\frac{U}{2}\right)^{2}+\left(\epsilon_{\mathbf{\boldsymbol{k}}}+\varepsilon_{\mathbf{\boldsymbol{k}}}^{+}\right)^{2}}},&\frac{\epsilon_{\mathbf{\boldsymbol{k}}}+\varepsilon_{\mathbf{\boldsymbol{k}}}^{-}}{\sqrt{\left(\frac{U}{2}\right)^{2}+\left(\epsilon_{\mathbf{\boldsymbol{k}}}+\varepsilon_{\mathbf{\boldsymbol{k}}}^{-}\right)^{2}}}\\ \frac{U/2}{\sqrt{\left(\frac{U}{2}\right)^{2}+\left(\epsilon_{\mathbf{\boldsymbol{k}}}+\varepsilon_{\mathbf{\boldsymbol{k}}}^{+}\right)^{2}}},&\frac{U/2}{\sqrt{\left(\frac{U}{2}\right)^{2}+\left(\epsilon_{\mathbf{\boldsymbol{k}}}+\varepsilon_{\mathbf{\boldsymbol{k}}}^{-}\right)^{2}}}\end{array}\right).$$

$$H_{\mathrm{ex}}=-U{e}^{i\mathbf{\pi }\mathbf{m}}{d}_{\mathbf{m}}^{†}{d}_{\mathbf{m}}^{}.$$

$${\tilde{d}}_{\alpha }^{†}$$

$${\tilde{d}}_{\alpha }^{}$$

$$H=\sum _{\alpha }{\epsilon }_{\alpha }{\tilde{d}}_{\alpha }^{†}{\tilde{d}}_{\alpha }^{}-U{e}^{i\mathbf{\pi }\mathbf{m}}\sum _{\alpha ,\beta }{\psi }_{\alpha \mathbf{m}}{\psi }_{\mathbf{m}\beta }^{†}{\tilde{d}}_{\alpha }^{†}{\tilde{d}}_{\beta }^{}.$$

$$G_{\alpha ,\beta }(\tau )=-i\Theta (\tau )⟨{\left[{\tilde{d}}_{\alpha }^{}(\tau ),{\tilde{d}}_{\beta }^{†}(0)\right]}_{+}⟩$$

$$\left(i\frac{d}{d\tau }-{\epsilon }_{\alpha }\right)G_{\alpha ,\beta }(\tau )$$

$$\left(\epsilon -{\epsilon }_{\alpha }+i0\right)G_{\alpha ,\beta }(\epsilon )={\delta }_{\alpha ,\beta }-U{e}^{i\mathbf{\pi m}}\sum _{{\alpha }^{\prime }}{\psi }_{\alpha \mathbf{m}}{\psi }_{\mathbf{m}{\alpha }^{\prime }}^{†}{G}_{{\alpha }^{\prime },\beta }(\epsilon ),$$

Full text

I Introduction

The Falicov-Kimball model (FKM) Falicov and Kimball (1969) is one of the simplest models for the description of correlated electrons on a lattice. Over time, it has become a standard tool for the investigation of various phenomena including crystallization Kennedy and Lieb (1986); Jedrzejewski et al. (1989); Gruber et al. (1994), metal-insulator and valence transitions Plischke (1972); Michielsen (1994); Farkašovský (1995a, b); Portengen et al. (1996); Czycholl (1999); Byczuk (2005); Lemański et al. (2017); Haldar et al. (2017a); Farkašovský (2019); Nasu et al. (2019); Haldar et al. (2019), inhomogeneous charge and spin orderings  Lemberger (1992); Freericks et al. (1999); Lemański et al. (2002); Freericks et al. (2002); Lemanski et al. (2004); Tran (2006); Čenčariková and Farkašovský (2011); Žonda (2012); Debski (2016); Kapcia et al. (2017), nonlocal correlations Schiller (1999); Hettler et al. (2000); Ribic et al. (2016, 2017), ferroelectricity Batista (2002); Farkašovský (2002, 2008); Schneider and Czycholl (2008); Golosov (2013), mixtures of heavy and light cold atoms in optical lattices Maśka et al. (2008); Iskin and Freericks (2009); Maśka et al. (2011); Hu et al. (2015), transport through layered systems Freericks et al. (2001); Freericks (2004, 2006); Hale and Freericks (2012); Kaneko et al. (2013); Žonda and Thoss (2019) or different nonequilibrium phenomena Freericks et al. (2006); Turkowski and Freericks (2007); Eckstein and Kollar (2008); Eckstein et al. (2009); Herrmann et al. (2016); Haldar et al. (2016, 2017b); Smith et al. (2017); Qin and Hofstetter (2018); Herrmann et al. (2018).

Its biggest advantage over the paradigmatic Hubbard model J. (1963) lies in the fact that it is accessible by exact methods. It is exactly solvable in the limit of infinite dimensions (infinite coordination number) by means of dynamical mean field theory (DMFT) Brandt and Mielsch (1989, 1990); Janiš (1991); Freericks and Zlatic (2003) and it can be addressed by an exact, sign-problem-free Monte Carlo (MC) method Maśka and Czajka (2006); Žonda et al. (2009); Žonda (2012); Huang and Wang (2017) in finite dimensions. Both methods take advantage of the fact that the FKM combines quantum and classical degrees of freedom.

Despite the simplicity and accessibility of this model, its research continues to offer new and often surprising results. This is true even for its simplest spin-less version at the particle-hole symmetric point. For example, nonequilibrium DMFT and cluster approximation method studies showed that its quantum subsystem does not thermalize after an interaction quench Eckstein and Kollar (2008); Herrmann et al. (2016, 2018); simple generalizations of the FKM can be used to study the interplay of topology and interactions at finite temperatures Gonçalves et al. (2019) or fractionalization of particles into charge and spin objects Tran (2019); it was used to derive universal features of the critical metal-insulator transition that are transferable to other Hubbard-like models Janiš and Pokorný (2014); Haldar et al. (2017a), utilized in studies of different quasiparticles Prosko et al. (2017); Kauch et al. (2019); and, as discussed below, the phase diagram of the model is still in question as well.

The Hamiltonian of the spinless FKM at half filling reads

[TABLE]

where the first term describes nearest-neighbor hopping of spinless electrons on a lattice. The second term represents a Coulomb-like local interaction between the localized $f$ particle and itinerant $d$ electron on the same lattice site. The terms with factor $-\frac{1}{2}$ secure the half-filling conditions $N_{f}+N_{d}=L/2$ for chemical potential $\mu=0$. Here $N_{f(d)}$ is the total number of $f$ ($d$) particles and $L$ is the total number of lattice points.

The phase diagram of this model in finite as well as infinite dimensions contains three main phases (see Fig. 1): an ordered charge density wave (CDW) phase (OP) that exists at low temperatures Brandt and Mielsch (1989, 1990, 1991); Chen et al. (2003); Žonda et al. (2009), a gapless disordered phase for weak interaction $U$ and high temperatures (DPw), and a gapped disordered phase for strong interaction $U$ and high temperatures (DPs) Gruber and Macris (1996); Freericks and Zlatic (2003).

However, this is not a complete picture. A recent study of the model on a two-dimensional ($D=2$) lattice Antipov et al. (2016) showed that in the thermodynamic limit DPw exhibits Anderson localization which destabilizes the metallic-like phase reported in older works focused on relatively small lattice sizes Maśka and Czajka (2006); Žonda et al. (2009). Therefore, all three phases are insulating in the thermodynamics limit. Nevertheless, for any finite system size, there is a crossover from an Anderson localized insulating phase at intermediate $U$ through a weakly localized regime, with the above mentioned metallic-like character, to a Fermi gas at $U=0$. In addition, a series of papers proved that in infinite dimensions ($D\rightarrow\infty$) there is a stable ordered CDW phase without a gap at the Fermi level Hassan and Krishnamurthy (2007); Matveev et al. (2008); Lemański and Ziegler (2014); Lemański (2016); Kapcia et al. (2019).

The gapless CDW phase in the infinite dimensional FKM is related to the existence of narrow bands in the density of states (DOS) that are formed inside the CDW gap at finite temperatures. These subgap bands come in pairs placed symmetrically around the center of the gap and merge for a finite range of interaction strengths $U$ and temperatures. The resulting merged single subgap band is centered around the Fermi level and, consequently, there is no gap at the Fermi level. Lemański argued that this merging is related to the inversion of the subgap bands belonging to two sublattices of a bipartite lattice at critical interaction $U_{c}$ Lemański (2016). Here a sublattice A constitutes such lattice points that all their nearest neighbors belong to the sublattice B and vice versa (they are alternating). The DOS calculated for each sublattice separately contains both subgap bands placed symmetrically around the Fermi energy. However, they differ in width and height. This property is related to the CDW ordering because the average $f$-electron occupancy differs for the sublattices as was discussed in various studies Hassan and Krishnamurthy (2007); Matveev et al. (2008); Lemański and Ziegler (2014). What is interesting is that at some critical $U_{c}$ the qualitatively different subgap bands belonging to different sublattices flip positions.

The open question is if there is an analog of such a band crossing in finite dimensions as well. The subgap states, respective bands, had already been discussed in both $D=2$ and $D=3$ Maśka and Czajka (2006); Tran (2006); Žonda et al. (2009). However, a systematic study focused on the region around $U_{c}$ is missing. The present paper has the aim to fill this gap. We show that there is indeed a clear crossing of distinct subgap states in finite dimensions. Moreover, we reveal the underlying mechanism of the crossing by focusing on the single particle excitations from the CDW ground state. We show that these excitations significantly influence the density of states up to surprisingly high temperatures approaching $T_{c}$ of the order-disorder transition. We also demonstrate that, despite the presence of strong localization, the crossing can support charge transmission through a finite-sized system in the gapped regime. This is done by addressing a heterostructure where the finite system modeled by the FKM is sandwiched between two metallic leads. We mostly focus on the $D=2$ case, but $D=1$ and $D=3$ are addressed as well.

The rest of the paper is organized as follows. In Sec. II we outline the methodology for addressing the thermodynamic properties of the FKM and introduce the model of the heterostructure as well as the method for studying charge transport in it. The main results are presented in Sec. III, where we first analyze the origin and properties of the sharp features of the subgap DOS in Sec. III.1 and then show in Sec. III.2 how these affect the charge transport through finite FKM coupled to metallic leads. Section IV concludes with a summary. In the Appendices we show some analytical results on the existence of the gap and the positions of sharp subgap states.

II Methods

II.1 Thermodynamic properties

The $f$ particles in the FKM represent classical degrees of freedom. This can be seen from the fact that the $f$-particle occupation numbers $f_{\boldsymbol{l}}^{\dagger}f_{\boldsymbol{l}}^{\phantom{\dagger}}$ commute with the entire Hamiltonian in Eq. (1). This allows us to replace any operator $f_{\boldsymbol{l}}^{\dagger}f_{\boldsymbol{l}}^{\phantom{\dagger}}$ by its eigenvalues $w_{\boldsymbol{l}}=0$ or $w_{\boldsymbol{l}}=1$ and write a partial Hamiltonian for a particular classical configuration $w$. After neglecting the constant energy shift $-U/4$, the Hamiltonian in Eq. (1) reads for a chosen configuration $w$

[TABLE]

Thereby, $\varepsilon_{j}$ are the eigenvalues of the matrix with elements $h_{\boldsymbol{l\,l}^{\prime}}=Uw_{\boldsymbol{l}}\delta_{\boldsymbol{l\,l}^{\prime}}-t\delta_{\left\langle\boldsymbol{l,l}^{\prime}\right\rangle}$, where $\delta_{\left\langle\boldsymbol{l,l}^{\prime}\right\rangle}=1$ when the lattice positions represented by vectors $\boldsymbol{l}$ and $\boldsymbol{l}^{\prime}$ are the nearest neighbors and zero otherwise. The ground-state configuration for any bipartite lattice at the particle-hole symmetric point is the checkerboard ordering of $f$ particles Kennedy and Lieb (1986); Gruber and Macris (1996). The corresponding configuration $w$ can be written as $w_{\mathbf{l}}=(1+e^{i\boldsymbol{\pi}\boldsymbol{l}})/2$ for any hypercubic lattice. It is easy to show (see Appendix A) that such an effective potential opens a gap of the width $U$ in the band structure which is centered around the Fermi energy.

An advantage of the Falicov-Kimball model is that the mean values of any $d$-electron operator $\hat{O}$ can be written in the form

[TABLE]

where

[TABLE]

with $Z=\sum_{w}e^{-\beta F\left(w\right)}$ being the partition function (we assume $\mu=0$). Here $\left\langle.\right\rangle_{d}$ is the trace over the $d$-electron subsystem for fixed $w$ Maśka and Czajka (2006). As this is a single-particle problem, the trace can be calculated effectively using exact numerical diagonalization. The sum over configurations $w$ can then be calculated using a Metropolis algorithm based Monte Carlo method Maśka and Czajka (2006); Žonda et al. (2009, 2009); Žonda (2012); Huang and Wang (2017).

The quantities that we are mostly interested in are the normalized DOS defined as

[TABLE]

where $\mathrm{Tr}_{w}\equiv\frac{1}{Z}\sum_{w}e^{-\beta F\left(w\right)}$ and averaged inverse participation ratio (IPR)

[TABLE]

where $\mathrm{DOS}_{i}(\varepsilon,w)=\sum_{\alpha}\delta\left(\varepsilon-\varepsilon_{\alpha}\right)\mathcal{U}_{i\alpha}\mathcal{U}_{\alpha i}^{\dagger}/L$ is the local DOS and the matrix $\mathcal{U}$ consists of the eigenvectors belonging to eigenvalues $\varepsilon_{\alpha}$ of the matrix $\boldsymbol{h}$ from Eq. (2) calculated for a particular configuration $w$ and arranged in columns. The matrix $\mathcal{U}$ can be evaluated numerically for a finite system and it can be chosen to be real.

The finite size scaling of the IPR can be used for studying localization of the itinerant electrons in the system Evers and Mirlin (2008); Murphy et al. (2011); Perera and Wortis (2015). The $\mathrm{IPR}$ scales as $1/L$ for a completely itinerant system states and converges to a finite value with increasing $L$ for strongly localized states. In the case of perfect spacial localization to a single lattice point, the IPR converges to one. Note, that because of the finite size of our lattices, we use a Gaussian broadening of the $\delta$-functions, $\delta\left(\varepsilon-\varepsilon_{\alpha}\right)\approx\exp[-(\varepsilon-\varepsilon_{\alpha})^{2}/(2b^{2})]/(b\sqrt{2\pi})$. In most of presented cases, we set the broadening constant to $b=0.002t$. This small value is a compromise between the effort to suppress the influence of the artificial broadening on our results (especially IPR) and the preservation of the stability of the calculations of the IPR for a broad range of temperatures and lattice sizes. We discuss the influence of the Gaussian broadening on our results in detail below.

II.2 Heterostructure

Besides studying an isolated FKM, we also address a heterostructure $H_{\mathrm{h}}=H_{\mathrm{FK}}+\sum_{l=L,R}H_{\mathrm{lead}}^{l}+H_{\mathrm{hyb}}^{l}$, where a two dimensional FK system is sandwiched between two metallic leads as illustrated in Fig. 2. The central system $H_{\mathrm{FK}}$ is finite in $x$ but in principle infinite in the $y$ direction (modeled by periodic boundary conditions). The leads and hybridization terms read

[TABLE]

where $t_{l}$ is the hopping for lead $l=L,R$ , $\epsilon$ represents an energy shift of the lead, and $\gamma_{l}$ is hopping parameter between the system and the lead $l$.

We have two main motivations for addressing this more complex setup. First, the broadening of the system states is in the case of the heterostructure provided naturally by the coupling $\gamma_{l}$ to the semi-infinite leads, which allows to test the results obtained for isolated system potentially influenced by an artificial broadening of the $\delta$ functions. Second, we want to reveal how the existence of the finite DOS in the CDW gap influences the transport properties of the model.

The properties of the heterostructure can be effectively addressed by a combination of a sign-problem free Monte-Carlo with the nonequilibrium Green’s function technique approach Žonda and Thoss (2019). We assume in our analysis that the central FK system was in the distant past decoupled from the leads and that both system and leads had been in thermal equilibrium. The occupation numbers of the $f$ electrons are integrals of motion, therefore, their distribution can be calculated for the isolated system as it will not change after coupling to the leads. Here, in contrast to the previous subsection, we assume open boundary conditions of the central system at the system-lead interface. Test calculations show, that if system is large enough ($L_{x}>10$) the influence of the boundary conditions on the $f$-electron distribution is negligible. We further assume, that the semi-infinite leads are unaffected by the system and are modeled by parallel chains coupled to the central system as shown Fig. 2. Therefore, the leads can be characterized by their surface density of states

[TABLE]

with band half-width $B=2t_{l}$ centered around the band energy shift $\epsilon_{l}$ from Eqs. (7). This allow us to write an exact formal solution for the Green’s function of the heterostructure for a particular configuration $w$ (for details, see Refs. Žonda and Thoss (2019); Jauho et al. (1994)):

[TABLE]

Here, $\mathbf{G}^{r\,(a)}$ is the retarded (advanced) Green’s function of the coupled system, $\mathbf{G}^{<}$ is the lesser Green’s function of the coupled system, and $\boldsymbol{g}^{r\,(a)}\left(\varepsilon\right)$ is the retarded (advance) Green’s function of the bare system with components:

[TABLE]

The total tunneling self-energies $\mathbf{\Sigma}^{r,a,<}=\mathbf{\Sigma}_{L}^{r,a,<}+\mathbf{\Sigma}_{R}^{r,a,<}$ have the components

[TABLE]

where $\{s^{L,R}\}$ are the sets of system lattice positions at the left and right interfaces, $f_{l}(\varepsilon)$ is the Fermi function, and $\mu_{l=L,R}$ is the chemical potential of the leads.

The transmission function, which contains the most detailed information on charge transport, has for a specific $w$ a compact form Haug and Jauho (2008),

[TABLE]

where the trace goes over the $d$-electron subsystem. Its total mean value is obtained by a trace over the $f$-electron subsystem which is done by the MC method. Similarly, the local density of states of a coupled system (heterostructure) for a specific $w$ can be calculated as

[TABLE]

and it allows us to define an averaged generalized inverse participation ratio (gIPR)

[TABLE]

where $\mathrm{DOSh}(\varepsilon,w)=\mathrm{Tr}_{d}i\left[\mathbf{G}^{r}\left(\varepsilon\right)-\mathbf{G}^{a}\left(\varepsilon\right)\right]/2\pi L$.

In the following, we set $\mu_{L}=\epsilon_{L}$ and $\mu_{R}=\epsilon_{R}$, which corresponds to half-filled lead bands and we introduce a finite voltage drop as $V=\mu_{L}-\mu_{R}$ with antisymmetric condition $\mu_{L}=-\mu_{R}$. In equilibrium, we set $\mu_{L}=\mu_{R}=0$ and assume that the electrostatic potential (set to zero) of the central system is uninfluenced by the leads. This fixes the half-filling condition for the central system. We focus on equivalent leads $\gamma=\gamma_{L}=\gamma_{R}$ with a broad band half-width $B=10t$.

III Results

III.1 Subgap density of states

The typical subgap bands, calculated for $D=2$ at temperatures in the vicinity of the CDW phase transition Maśka and Czajka (2006); Žonda et al. (2009); Antipov et al. (2016), resemble the exact DMFT results calculated for infinite dimensions Hassan and Krishnamurthy (2007); Matveev et al. (2008); Lemański and Ziegler (2014); Lemański (2016). However, in finite dimensions the subgap bands reduce with the decreasing temperature into sharp features pinned mostly to few distinct energies. This is illustrated in Fig. 3, where we show the subgap ($\left|\varepsilon\right|<U/2$) low-temperature DOS calculated for $L=24\times 24$, three values of $U$, and various temperatures below $T_{c}$.

The character of the subgap DOS changes with $U$. For $U=1t$, there are two pronounced maxima placed symmetrically at $\varepsilon\sim\pm 0.25U$; for $U=2.5t$, a sharp maximum is centered around $\varepsilon=0$ and is accompanied by two sharp features of approximately a third of its height, which are placed at $\varepsilon\sim\pm 0.15U$; for $U=4t$, four sharp local maxima of comparable heights exist at $\varepsilon\sim\pm 0.13U$ and $\varepsilon\sim\pm 0.21U$. Considering their positions, the subgap maxima have the same qualitative dependence on $U$ that was described for the subgap bands for $D\rightarrow\infty$ Hassan and Krishnamurthy (2007); Matveev et al. (2008); Lemański and Ziegler (2014); Lemański (2016). Namely, the maxima approach each other with increasing $U$ until they merge at some critical $U_{c}$ and then, above it, draw apart.

The position of these most pronounced distinct local maxima does not depend on the temperature and, in contrast to other states, their weight is not completely negligible even for a very low temperatures.

Because the $d$-electron DOS is dictated by the distribution of $f$ particles, one can expect that the origin of these distinct maxima must be in some states reflecting the low energy disruptions of the ground-state checkerboard ordering. To analyze this conjecture we focus on three typical single $f$-electron excitations from the checkerboard ordering illustrated in the top panels of Fig. 4. They represent an addition of one $f$ electron to the otherwise perfect checkerboard ordering; a removal of one $f$ electron and, finally, a displacement of a single $f$ electron to the nearest unoccupied lattice point. The actual spacial position of these three defects does not play a role, as we are assuming periodic boundary conditions. Note that because $f$ electrons can be also interpreted as ions, the studied defects can be seen as an equivalent of typical lattice defects such as vacancies (Schottky defects) or interstitial defects (Frenkel defects).

Because the above disruptions of the checkerboard pattern play a role of classical single impurities, they lead to sharp bound states inside of the CDW gap Economou (2006). We plot the normalized position of these states (i.e., their calculated eigenenergies) in Fig. 4 as a function of $U$ in all three realistic dimensions (for details and analytical results, see Appendix B). The positions of the the main peaks in the finite temperature subgap DOS plotted in Fig. 3 coincide with the eigenenergies at the same $U$ and $D=2$ shown in Fig. 4(b) (they are also signaled by the arrows in Fig. 3). The highest peaks in Fig. 3 reflect an addition or removal of an $f$ electron and the second in order reflect a single $f$-electron displacement. Therefore we can conclude that these simple $f$-electron excitations are the main underlaying mechanism for the stable finite-temperature subgap anomalies.

The results show that displacing an $f$ electron leads to two subgap bound states placed symmetrically around $\varepsilon=0$. Adding or removing an $f$ electron leads to a single subgap state. These are related by $\varepsilon_{\mathrm{add}}=-\varepsilon_{\mathrm{rem}}$ and cross each other at critical $U_{c}$ which depends on the dimension. For the one-dimensional case we get $U_{c}=4/\sqrt{3}$, for $D=2$ it is $U_{c}\cong 2.5$, and $D=3$ leads to $U_{c}\cong 3.18t$ (see Appendix B)111Note, that this crossing resembles the situation in other gaped systems with bound states, e.g. superconducting quantum dots Balatsky et al. (2006); Meden (2019).. This means that similarly to the infinite dimensional case Lemański (2016), the main subgap features exchange roles at $U_{c}$.

Another qualitative change that takes place at $U_{c}$ can be seen from the difference between the minimal energy of our disrupted configurations ($E_{e}$) and the real ground-state energy of the checkerboard ordering ($E_{g}$) plotted in Fig. 4(d). There is a kink in the $E_{e}-E_{g}$ dependence on $U$ (solid red line) exactly at $U_{c}$ for and added as well as removed $f$ electron. This is because below $U_{c}$ the minimal energy requirement always sets $N_{d}=L/2$, which leads to $N_{d}+N_{f}=L+1$ for a configuration with an additional $f$ electron and $N_{d}+N_{f}=L-1$ for a removed one. This may seem strange, considering that $\mu=0$ should ensure the half-filling, but it is a straightforward consequence of the finite lattice size. The number of eigenenergies $\varepsilon_{\alpha}$ equals $L$. These are, for a perfect checkerboard pattern, divided equally into two bands (Appendix A). For a disrupted configuration, the subgap state $\varepsilon_{\mathrm{add}}$ is pulled out from the lower band. Consequently, as it is energetically advantageous to occupy this state if $\varepsilon_{\mathrm{add}}<0$, which is the case below $U_{c}$, this leads to $N_{d}=L/2$. Above $U_{c}$, we have $\varepsilon_{\mathrm{add}}>0$, which leaves the state unoccupied and therefore $N_{d}=L/2-1$. The opposite is true for $\varepsilon_{\mathrm{rem}}$. Therefore, the half-filling for any single of these excitations is restored only above $U_{c}$ and we need a combination of them to fulfill this requirement below $U_{c}$. The energy profile for forced condition $N_{f}+N_{d}=L$ is plotted in Figs. 4(d,e) using a dotted line. Note that the situation for $D=1$ is somewhat more complicated as here the displacement of a single $f$ particle can have a lower energy than adding or removing a localized particle at both weak and strong interaction limits [see Fig. 4(e)].

The above discussed single $f$-electron excitations have a profound effect on the DOS even at surprisingly high temperatures. This is already illustrated by the sharp subgap features in the finite temperature DOS plotted in Fig. 3. Nevertheless, we are especially interested in the ordered phase with high DOS at the Fermi level analogous to the one observed in infinite dimensions Lemański and Ziegler (2014). Therefore, we show in Fig. 5(a) the dependence of DOS$(\varepsilon=0)$ on $U$ at various temperatures. The highest temperature $T=1t$ represents the disordered phase, $T=0.12t$ is just above the critical temperature for $U=2.5t$ and the rest is below it. Figure 5(a) is a direct $D=2$ analog of the infinite dimensional case shown in Fig. 7 of the work by Lemański and Ziegler Lemański and Ziegler (2014). Although both cases share some qualitative characteristics, like the increase of the DOS around $U_{c}$ for small temperatures, there is one distinct difference. The increase of the DOS around $U_{c}$ is not only much sharper, but for $T\lesssim T_{c}$ it also leads to values of DOS which are higher than the high temperature limit where the gap is completely closed. This is again illustrated in Fig. 6(a), where we show the temperature dependence of the DOS$(\varepsilon=0)$ for $D=2$. The position of the maximum of DOS ($T_{m}\sim 0.095t$) is clearly below the critical temperature ($T_{c}\sim 0.12t$). Moreover, the DOS profile is stable with respect to the lattice size as it is illustrated in Fig. 6(a), Fig. 5(b), and in the inset of Fig. 3(b), where we show the detail of the DOS around the Fermi level.

The above-discussed DOS were calculated with an artificial Gaussian broadening of the $\delta$-functions with $b=0.002t$. The effect of the Gaussian broadening on the DOS$(\varepsilon=0)$ is completely negligible in the disordered phase as illustrated in Fig. 5(c) by the red circles. The situation in the ordered CDW phase is more complicated, especially for critical $U_{c}$ and low temperatures. Figure 5(c) shows that the DOS$(\varepsilon=0)$ calculated at $T=0.095t$ (black circles) increases with the decreasing broadening. Also in Fig. 6(a), one can see by following the dashed line, which represents the DOS$(\varepsilon=0)$ calculated for $L=20\times 20$ and broader broadening of $b=0.004t$, that bigger artificial broadening leads to a decrease of the calculated DOS$(\varepsilon=0)$ maximum.

The reason is twofold. For the $D=2$, the band around $\varepsilon=0$ has a fine structure, as is shown in the inset of Fig. 3(b), and has a clear maximum at $\varepsilon=0$. A wide artificial broadening smooths this structures, which significantly lowers the DOS at the Fermi level. From our data it is not clear if this fine structure will disappear in the thermodynamic limit. It seems not to be the case, as the detail of the DOS shown in the inset of Fig. 3(b) depends only weakly on the lattice size, but a study on much bigger lattices is required to confirm this.

Even more important is that despite the broadening into a band provided by multi-particle excitations, the states at the Fermi level stay very sharp even for large lattices, as is shown in the inset of Fig. 6(b). This effect becomes even more crucial with decreasing temperature as here the configurations with just one $f$ particle added or removed from a perfect checkerboard ordering can have a very high weight. This is shown in Fig. 6(b), where we plotted the total MC weight of these configuration $w\sim\left(e^{-\beta F(w_{\mathrm{add}})}+e^{-\beta F(w_{\mathrm{rem}})}\right)/Z$ as a function of temperature for various lattice sizes (solid lines) and compare it with the analogous weight for checkerboard ordering (dashed line). There is a clear maximum at which the combined weight of the above excited states is almost one-third of the total one. This can explain the extremely sharp subgap states for low temperatures in the finite system shown in Fig. 3 222Some caution is in place regarding this statement. The used MC is of a single-flip character which could in principle distort the statistics if there would be some low energy $f$-electron configurations, which are hard to reach by a series of single particle updates from e.g. the checkerboard ordering. Nevertheless, we have tested this by using a big set of random initial conditions and this does not seem to be the case..

The enhanced DOS at the Fermi level caused by the crossing of the sharp maxima might raise a question if the phase has still an insulating character. However, knowing that the prevailing mechanism behind this effect is the impurity like single-particle excitation of the localized subsystem, one can expect a strong localization of the $d$ electrons. We show by studying the scaling of the averaged IPR($\varepsilon=0$) depicted in Fig. 7(a) that this is really the case. For localized states, the IPR should saturate with increasing lattice size to a finite value. Note that IPR$\rightarrow 1$ would point to a complete localization of the $d$ electrons to a single lattice point. Therefore, the saturation of IPR to $\sim 0.16$ for $T=0.095t$ shown in Fig. 7(a) (red circles) still points to a strong localization. On the other hand, the slow decline of the IPR for $T=1t$ (blue circles) points to a weak localization as expected in this regime for a finite size system Antipov et al. (2016); Žonda and Thoss (2019).

As shown in Fig. 7(b), the IPR is, in contrast to the DOS results discussed above, relatively stable for a broad range of broadening parameter $b$ (note the logarithmic scale) even below the critical temperature. Therefore, the conclusions of strong localization of the crossed states is not affected by the artificial broadening of the $\delta$-functions.

To further support these findings, we provide an alternative test of the above conclusions in the next subsection. Instead of artificial broadening, we consider a heterostructure where the system is coupled to two semi-infinite metallic leads. These provide a different kind of broadening of the system states. It can be argued that this broadening is more natural, but it is also uneven, because the broadening of the LDOS decreases with increasing distance from the system-leads interface Freericks (2004); Žonda and Thoss (2019). In addition, this setup allows us to study the transport thought a finite system.

III.2 Transport properties of a heterostructure

We have shown recently Žonda and Thoss (2019) that the typical phases of the FKM have different influence on the transport properties of a heterostructure. However, that study did not focus on the particular case where the system is in CDW phase but has a large density of states at the Fermi level as is the case in Fig. 6(a). Here we study the effect of the finite DOS in this regime on the transport properties for a finite dimensional central system. Because of that, we first have to readdress the problem of DOS and localization for the heterostructure.

In contrast to the isolated FKM studied in the previous section, in the heterostructure the broadening of the system states results naturally from the coupling to the semi-infinite leads. In Fig. 8, we show the averaged DOSh$(\varepsilon=0)$ for $L=20\times 20$ and three values of $\gamma$. The DOSh$(\varepsilon=0)$ above the critical temperature, signaled by vertical dotted line, is identical to the one in Fig. 6(a) and does not depend on the lattice size as can be seen in the inset of Fig. 8 (blue line). The positions of the maxima are in compliance as well ($T\sim 0.095t$) and the maximum of DOSh$(\varepsilon=0)$ is well above its value in the disordered phase. In addition, the DOSh$(\varepsilon=0)$ shown in Fig. 8 depends only weakly on the chosen values of $\gamma$’s. This supports the conclusion that the crossing of the major subgap bands at critical $U$ can lead to a DOS at the Fermi level, which exceeds its values in the gapless disordered phase. On the other hand, the DOSh$(\varepsilon=0)$ calculated for $T=0.095t$ depends much stronger on the lattice size (red circles in the inset of Fig. 8) than the equivalent DOS calculated for isolated FKM. This, as well as the increasing error bars, is a direct consequence of the fact that the broadening coming from the leads is not homogeneous in the central system as its effect on the LDOSh decreases with the distance from the system-leads interfaces Žonda and Thoss (2019). Consequently, the broadening in the central part of the system coming from the leads might vanish fast with the increased lattice size.

The qualitative differences between the natural and artificial broadening allow us to perform an alternative investigation of the localization based on the gIPR defined in Eq. (16). A direct comparison of the finite size scaling of the gIPR calculated for $\gamma=2t$ with the IPR is shown Fig. 7(a). The gIPR for $T=0.095t$ (black squares) has the same profile as IPR and it convergences to a similar finite value, which confirms strong localization even for the coupled system. Similarly, the scaling of the gIPR in the disordered phase represented by $T=1t$ (green squares) points to a weakly localized central system at best (see also Ref. Žonda and Thoss (2019)). The comparison with the IPR also reveals that the coupling to the leads can suppress the localization in the finite system in this regime.

We analyze the effect of coupling to the leads on the localization in the finite system ($L=20\times 20$) in more detail in Fig. 9. Both gIPR curves calculated at $T=0.095t$ and $T=1t$ are saturated at low values of the coupling $\gamma$. The saturated values are in good agreement with the ones calculated for the isolated system [Fig. 7(b)]. This can also be interpreted as evidence that the coupling to the leads provides a good independent method for studying the problem of the localization.

However, as we increase $\gamma$ (note the logarithmic scale in Fig. 9) the gIPR significantly decreases for $\gamma\apprge 0.3t$ at $T=1t$ and $\gamma\apprge 1t$ for $T=0.095t$. This effect is stronger in the disordered phase, where at strong coupling the already small gIPR drops to half its weak coupling value. Nevertheless, the localization is weakened in the ordered phase as well. It is therefore worth it to examine how the coupling to the leads affects the transport properties in a finite system.

We focus on the transmission function as this provides the most detail information on the charge transport. Figure 10(a) shows the transmission function at $\varepsilon=0$ as a function of temperature for $U=2.5t$, central system size $L=20\times 20$ and for three values of system lead hopping $\gamma$. We focus on the equilibrium situation ($V=0$) because a voltage that is smaller than the CDW gap ($V<U/2$) has only a small effect on the transmission function at $\varepsilon=0$. This is shown in Fig. 10(b) for $U=2.5t$ and $T=0.095t$, where voltage values are spread on a logarithmic scale.

The transmission function in Fig. 10(a) is negligible for low temperatures, but a clear “hump” for $\gamma=2t$ and local maximum for $\gamma<2t$ (signaled by an arrow) are present close to the temperature where the DOS of the central system has its maximum [see Fig. 8]. The transmission function at Fermi level is still several magnitudes smaller than its values outside the CDW gap, as seen in the inset of Fig. 10(a), and it is negligible for any energy within the gap with the exception of the close vicinity of $\varepsilon=0$. Nevertheless, its clear increase close to $T\sim 0.1t$ shows that, despite the relatively strong localization, the crossing of the subgap states can influence the charge transport through a finite system. Interestingly, it actually enables a finite charge transmission otherwise blocked by the CDW gap. Still, the transmission function significantly increases above the critical temperature of the order-disorder transition and, despite much lower DOS in this phase, the transmission function can be several times higher than at $T\sim 0.1t$.

This clearly reflects the difference of the localization in the ordered and disordered phases already shown in Fig. 7(a). This effect can also be seen in the system size scaling of the transmission functions plotted in Fig. 10(c). To highlight the difference, we scaled the transmission function by the linear size of the system. The weak localization in the disordered phase ($T=1t$), further lowered by the coupling to the leads, results in a situation where the scaled transmission function increases with $L$. On the other hand, the scaled transmission for $T=0.095t$ , which even for small system lattices is approximately ten times smaller than for $T=1t$, rapidly decreases with growing lattice size. The strong localization in this regime clearly overrules even the increasing DOSh($\varepsilon=0$) shown in the inset of Fig. 8.

We can therefore conclude that the high density of states in the CDW phase observed in the vicinity of $U_{c}$ can lead to a notable increase of charge transport through a finite system, but this effect rapidly vanishes with increasing system size. The reason is the relatively strong localization of the states at the Fermi level.

IV Conclusions

We have studied in realistic dimensions the structure of the subgap states of the FKM in the ordered CDW phase. We have shown that, similar to exact results in infinite dimensions, there are pronounced maxima in the subgap DOS placed symmetrically around the Fermi level which merge around some critical value of $U$ and exchange their positions above it. The position of these maxima does not depend on the temperature, because they mainly reflect the underlying sharp states resulting from single-particle excitations of the localized $f$-electron subsystem. The most pronounced of these states belong to an addition and removal of a single $f$ electron from a perfect checkerboard ordering.

The crossing of the most pronounced subgap maxima leads to a rise of the DOS at the Fermi level which can exceed even the DOS in the disordered gapless phase. However, the states in this regime are strongly localized. We have confirmed this by studying DOS and IPR for two different setups. First, we addressed an isolated FKM, where we used an artificial broadening of the states. Second, we studied a heterostructure, where the states had been broadened by coupling the system to simple semi-infinite leads. Although the coupling to the leads can lower the localization in the disordered phase, both studies have lead to the same qualitative conclusions.

In the case of the heterostructure, we have also shown that the significant increase of the DOS at critical $U$ can boost the charge transport trough a small finite system in the ordered phase. It increases the charge transmission which is otherwise suppressed in the whole range of energies within the CDW gap. Here, the strength of the coupling to the leads plays an important role. Nevertheless, with increasing lattice size this effect is quickly suppressed by the strong localization.

Acknowledgements.

The authors acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through Grant No. INST 40/467-1 FUGG. J. O. acknowledges support from Georg H. Endress Foundation. We thank Romuald Lemański and James K. Freericks for introducing us to the problem of the crossing of the subgap states in FKM and Tomáš Novotný for very helpful discussions.

Appendix A: CDW gap

Here we present the ground-state solution of the spinless FKM on a hypercubic lattice. We show that the checkerboard ordering of the $f$ electrons opens a gap in the DOS and derive general formulas for eigenvalues and eigenvectors.

The Hamiltonian Eq. (1) for the CDW ordering of $f$ electrons in hyper cubic lattice reads

[TABLE]

where $\boldsymbol{\delta}$ indices’s the relevant nearest neighbors (e.g., $1\boldsymbol{i}_{x}+0\boldsymbol{i}_{y}+0\boldsymbol{i}_{z}$). Fourier transformation, given by

[TABLE]

leads to

[TABLE]

where

[TABLE]

By replacing of the summation by integration and careful reshaping of the Brillouin zone, the Hamiltonian in Eq. (19) can be written in a matrix form

[TABLE]

Therefore, the states are divided into two bands separated by CDW gap with width $U$

[TABLE]

The related eigenvectors stored as a columns in matrix $Q_{\boldsymbol{k}}$ read

[TABLE]

Appendix B: Subgap states

In this Appendix, some analytical results for the positions of the subgap states resulting from single $f$-electron excitations are derived.

Let us excite the groundstate Hamiltonian in Eq. (19) by adding or removing a single $f$-electron on position $\boldsymbol{m}$

[TABLE]

Our aim here is to find the position of the sharp states in the spectra resulting from this single impurity. The Hamiltonian (17) can be formally diagonalized by unitary transformation,

[TABLE]

after which the total system Hamiltonian including $H_{\mathrm{ex}}$ reads

[TABLE]

We define the Green’s function,

[TABLE]

and use the equation of motion technique to get

[TABLE]

After the transformation into energy domain this reads

[TABLE]

or written in matrix form

[TABLE]

where $\mathbf{E}$ is diagonal matrix with elements $E_{\alpha,\alpha}=\varepsilon-\varepsilon_{\alpha}+i0$. We are primarily interested in the position of the bound states inside the CDW gap which can be calculated from the zeros of the determinant of the inverse Green’s function. Because $\mathbf{E}$ is diagonal and $\psi_{\alpha\alpha^{\prime}}^{\boldsymbol{m}}=\psi_{\alpha\boldsymbol{m}}\psi_{\boldsymbol{m}\alpha^{\prime}}^{\dagger}$ we can use Sylvester’s determinant theorem [$\det\left(X+cr\right)=\det(X)(1+rX^{-1}c)$]:

[TABLE]

The position of the sharp states in the CDW gap can therefore be obtained by solving

[TABLE]

We evaluate this equation by using the transformation from Eq. (18) together with the decomposition $M=Q\Lambda Q^{-1}$, where $M$ is the inner matrix from Eq. (26) and elements of the diagonal matrix $\Lambda$ are $\varepsilon_{\boldsymbol{k}}^{+}$ and $\varepsilon_{\boldsymbol{k}}^{-}$. The product in Eq. (29) then reads

[TABLE]

and the sum over index $\alpha$ changes to

[TABLE]

Using the new notation, Eq. (29) can be rewritten to a simple form

[TABLE]

In $D=1$ this equation can be easily evaluated and reads

[TABLE]

Therefore, the solution for $U$ as a function of the energy of bound states is

[TABLE]

where $\tilde{\varepsilon}=\varepsilon/U$. At the crossing of the states, we have $\varepsilon=0$ and therefore the critical interaction in $D=1$ is $U_{c}=4/\sqrt{3}$.

The evaluation of Eq. (30) in $D=2$ leads to

[TABLE]

where $u=U^{2}(0.25-\tilde{\varepsilon}^{2})$ and $\mathcal{K}(x)$ is the elliptic integral of the first kind. The critical interaction in $D=2$ is $U_{c}\doteq 2.5t$ (with rounding at the third decimal place). The $D=3$ case can be evaluated numerically and has the critical interaction $U_{c}\cong 3.18t$.

The displacement of the $f$ electron to a neighboring position can be seen as adding and removing an $f$ electron at adjoining positions and can be represented by an additional term,

[TABLE]

where $i_{n}$ choses the direction of the $f$-electron shift. This leads to a more complicated formula for the determinant

[TABLE]

The sums in this formula can be evaluated analogously to the previous case, but their evaluation leads to much more complicated expressions, which we therefore omit here.