Thermoelectric properties of a strongly correlated layer

TL;DR

This study explores how strong electronic interactions influence the thermoelectric properties of a correlated layer, revealing that such correlations can enhance thermopower and providing insights into optimizing thermoelectric efficiency.

Contribution

The paper applies nonequilibrium dynamical mean-field theory with an auxiliary master equation solver to analyze thermoelectric effects in a strongly correlated layer beyond linear response.

Findings

Strong correlations increase thermopower.

Thermoelectric efficiency depends on model parameters.

Theoretical criteria for efficient thermoelectric devices.

Abstract

In this paper we investigate the effect of strong electronic interactions on the thermoelectric properties of a simple generic system, consisting of a single correlated layer sandwiched between two metallic leads. Results will be given for the linear response regime as well as beyond, for which a full nonequilibrium many-body calculation is performed, based on nonequilibrium dynamical mean-field theory (DMFT). As impurity solver we use the auxiliary master equation approach, which addresses the impurity problem within a finite auxiliary system consisting of a correlated impurity, a small number of uncorrelated bath sites and Markovian environments. For the linear response regime, results will be presented for the Seebeck coefficient, the electrical conductance and the electronic contribution to the thermal conductance. Beyond linear response, i.e. for finite differences in the temperatures and/or the bias voltages in the leads, we study the dependence of the current on various model parameters, such as gate voltage and Hubbard interaction of the central layer. We will give a detailed parameter study as far as the thermoelectric efficiency is concerned. We find that strong correlations can indeed increase the thermopower of the device. In addition, some general theoretical requirements for an efficient thermoelectric device will be given.

Figures

Equations

$$ZT$$

$$L$$

$$\mathcal{H}=-\sum _{z,⟨{\mathbf{r}}_{\perp }^{},{\mathbf{r}}_{\perp }^{\prime }⟩,\sigma }{t}_z{c}_{z,{\mathbf{r}}_{\perp },\sigma }^{†}{c}_{z,{\mathbf{r}}_{\perp }^{\prime },\sigma }^{}-\sum _{⟨z,z^{\prime }⟩,{\mathbf{r}}_{\perp },\sigma }{t}_{z{z}^{\prime }}{c}_{z,{\mathbf{r}}_{\perp },\sigma }^{†}{c}_{z^{\prime },{\mathbf{r}}_{\perp },\sigma }^{}$$

$$\displaystyle+U\sum_{{\bf r}_{\perp}}\big{(}n_{0,{\bf r}_{\perp},\uparrow}-\frac{1}{2}\big{)}\big{(}n_{0,{\bf r}_{\perp},\downarrow}-\frac{1}{2}\big{)}+\sum_{z,{\bf r}_{\perp},\sigma}\varepsilon_{z}n_{z,{\bf r}_{\perp},\sigma}\,.$$

$$j_0$$

$$\underline{G}=\left(\begin{array}[]{cc}G^{R}&G^{K}\\ 0&G^{A}\end{array}\right)\,,$$

$${\underline{G}}_c^{-1}(\omega ,{\mathbf{k}}_{\perp })={\underline{g}}_c^{-1}(\omega ,{\mathbf{k}}_{\perp })-\sum _{\alpha =l,r}{v}_{\alpha }^2{\underline{g}}_{\alpha }(\omega ,{\mathbf{k}}_{\perp })-\underline{\Sigma }(\omega ),$$

$$g_c^R(\omega ,{\mathbf{k}}_{\perp })$$

$$g_{\alpha =l,r}^R(\omega ,{\mathbf{k}}_{\perp })$$

$$g_{\alpha =l,r}^K(\omega ,{\mathbf{k}}_{\perp })=2i(1-2f_{\alpha }(\omega ))\mathrm{Im}{g}_{\alpha }^R(\omega ,{\mathbf{k}}_{\perp }),$$

$${\underline{\Delta }}_{\mathrm{ph}}(\omega )={\underline{g}}_0^{-1}(\omega )-{\underline{G}}_{\mathrm{loc}}^{-1}(\omega )-\underline{\Sigma }(\omega ),$$

$${\underline{G}}_{\mathrm{loc}}(\omega )=\underset{\mathrm{BZ}}{\int }\frac{d{\mathbf{k}}_{\perp }}{(2\pi {)}^2}{\underline{G}}_c(\omega ,{\mathbf{k}}_{\perp }).$$

$$J_{\text{c/en}}=-2\underset{-\infty }{\overset{\infty }{\int }}\frac{d\omega }{2\pi }\mathcal{T}(\omega )\left(f_l(\omega )-f_r(\omega )\right){\zeta }_{\text{c/en}}(\omega ).$$

$$\mathcal{T}(\omega )$$

$$\mathcal{T}(\omega ,{\mathbf{k}}_{\perp })$$

$$A(\omega ,{\mathbf{k}}_{\perp })=-\frac{1}{\pi }\mathrm{Im}{G}_c^R(\omega ,{\mathbf{k}}_{\perp }),$$

$${\gamma }_{\alpha }(\omega ,{\mathbf{k}}_{\perp })=-2v_{\alpha }^2\mathrm{Im}{g}_{\alpha }^R(\omega ,{\mathbf{k}}_{\perp }).$$

$$\epsilon ({\mathbf{k}}_{\perp })$$

$$\mathcal{T}(\omega )$$

$${\rho }_{2D}(\epsilon )$$

$$J_{\text{el}}$$

$$J_{\text{en}}$$

$$\sigma ={\frac{\partial J}{\partial \Phi }|}_{\Delta T=0}=I_0,$$

$$S=-{\frac{\partial \Phi }{\partial \Delta T}|}_{J_{\text{el}}=0}=-\frac{I_1}{I_0}\beta ,$$

$${\kappa }_e={\frac{\partial J_{\text{en}}}{\partial \Delta T}|}_{J_{\text{el}}=0}=\sigma \left[\frac{I_2}{I_0}-\frac{I_1^2}{I_0^2}\right]\beta .$$

$$I_m=2\underset{-\infty }{\overset{\infty }{\int }}\frac{d\omega }{2\pi }\mathcal{T}(\omega )\left|f^{\prime }(\omega )\right|{\omega }^m,$$

$$L$$

$$L=L_0=\frac{{\pi }^2}{3}$$

$$|f^{\prime }(\omega )|=\frac{\beta }{4{\cosh }^2(\frac{\beta \omega }{2})}.$$

$$⟨{\omega }^m⟩=\int d\omega \rho (\omega ){\omega }^m.$$

$$\rho (\omega )$$

$$S$$

$$L$$

Full text

I Introduction

Thermoelectric devices on the nano scale are an interesting topic and research object, not only from the pure scientific point of view, but much more so for technological reasons. All aspects of thermoelectric devices are nicely reviewed by Benenti et al.Benenti et al. (2013, 2017). In particular, the efficient harvesting of the thermal energy is of great socio-economical importance. It depends on the physical properties of the given material and in the linear response regime it can be characterized by the dimensionless quantity, which is called thermoelectric figure of merit

[TABLE]

which is expressed by the Seebeck coefficient $S$ and the Lorenz number $L$. Here $T$ is the temperature, $\sigma$ is the electrical conductance and $\kappa$ is the thermal conductance. Any inelastic scattering will contribute to the thermal conductivity. However, for simplicity in the present work we consider only the electrical contribution $\kappa=\kappa_{e}$. A high value of the figure of merit is required to obtain high efficiency.Benenti et al. (2013); Kirchner, Zamani, and Muñoz (2013) In spite of the many theoretical and experimental effortsHarman et al. (2002); Chowdhury et al. (2009); Zhang et al. (2011); Majumdar (2004); Poudel et al. (2008); Venkatasubramanian et al. (2001); Arita et al. (2008); Hsu et al. (2004); Beyer et al. (2002); Sales, Mandrus, and Williams (1996); Pardo, Botana, and Baldomir (2013, 2012); Snyder and Toberer (2008); Dresselhaus et al. (2007); Dubi and Di Ventra (2011); Vineis et al. (2010); Sootsman, Chung, and Kanatzidis (2009); Shakouri (2011); Dorda et al. (2016); Freericks and Zlatić (2007) most efficient actual devices still operate at $ZT\lesssim 1$.Benenti et al. (2013) In this paper, we will study the impact of strong correlations on the thermoelectric properties of layered systems. It will turn out that the Lorenz ratio in these systems is close to the universal Wiedemann-Franz value $L_{0}$. Therefore, the figure of merit is solely governed by the Seebeck coefficient. The system we studied numerically, is composed of a single correlated layer sandwiched between two metallic leads. We will present results in and beyond the linear response regime. In the linear response regime we study the behavior of the Seebeck coefficient, electrical conductance and the electronic contribution to the thermal conductance, while in the non-linear regime we will address the electrical current $J$ and the corresponding power $P=J\cdot\Phi$, with $\Phi$ being the applied voltage.

To investigate the behavior of a many-body system out of equilibrium we use dynamical mean-field theory (DMFT)Georges et al. (1996); Vollhardt (2010); Metzner and Vollhardt (1989) combined with the nonequilibrium Green’s function approach originally proposed by KuboKubo (1957), SchwingerSchwinger (1961), Kadanoff, BaymBaym and Kadanoff (1961); Kadanoff and Baym (1962) and KeldyshKeldysh (1965). DMFT is a comprehensive, thermodynamically consistent and non-perturbative scheme. It is one of the powerful methods to investigate high-dimensional correlated systems and becomes exact in infinite dimensions. The sole approximation is the assumption of a local self-energy. In order to determine it, the original lattice problem is mapped onto a single impurity Anderson model (SIAM)Anderson (1961), which implies a self-consistency cycle. To solve the impurity problem, particularly in the nonequilibrium case, we used a recently developed auxiliary master equation approach (AMEA)Arrigoni, Knap, and von der Linden (2013); Dorda et al. (2014); Titvinidze et al. (2015); Dorda et al. (2017, 2016). The latter directly treats the nonequilibrium steady-state situation. However, AMEA is not restricted to the nonequilibrium situation, it can very efficiently be used in the equilibrium case as well. AMEA treats the impurity problem within an auxiliary system consisting of the correlated impurity, a small number of uncorrelated bath sites and two Markovian environments, which are encoded in a Lindblad equation.

The paper is organized as follows: In Sec. II.1 we introduce the Hamiltonian of the system. In Sec. II.2 we briefly discuss the nonequilibrium DMFT approach together with AMEA, while in Sec. II.3 we define important thermoelectric properties which we use in our calculations. Afterwards in Sec. III we present our results. Sec. III.1 is dedicated to the linear response results, while in Sec. III.2 we present our results in presence of the large temperature difference between the leads. Finally, in Sec. IV the results are summarized.

II Model and Method

II.1 Model

We consider a model system consisting of a single correlated infinite mono-layer ($z=0$) sandwiched between two metallic leads ($z<0$ and $z>0$), which are semi-infinite (see Fig. 1). The system is translationally invariant in the $xy$ plane (parallel to the correlated layer). The Hamiltonian reads

[TABLE]

Here $c_{z,{\bf r}_{\perp},\sigma}^{\dagger}$ creates an electron at site ${\bf r}_{\perp}$ of layer $z$ with spin ${\sigma}$ and $n_{z,{\bf r}_{\perp},\sigma}=c_{z,{\bf r}_{\perp},\sigma}^{\dagger}c_{z,{\bf r}_{\perp},\sigma}^{\phantom{\dagger}}$ denotes the corresponding occupation-number operator. $\langle z,z^{\prime}\rangle$ stands for neighboring layers and $\langle{\bf r}^{\phantom{\dagger}}_{\perp},{\bf r}^{\prime}_{\perp}\rangle$ stands for neighboring sites within a layer.

The first two terms of the Hamiltonian (3) describe nearest-neighbor intra-layer and inter-layer hoppings, with hopping amplitudes $t_{z}$ and $t_{zz^{\prime}}$, respectively. We consider identical nearest-neighbor hopping parameters within the leads. $t_{z}=t_{zz^{\prime}}=t_{l}$ for the left lead ($z,z^{\prime}<0$) and $t_{z}=t_{zz^{\prime}}=t_{r}$ for the right lead ($z,z^{\prime}>0$). The hopping amplitude inside the correlated layer is $t_{0}=t_{c}$, while the hybridization between left (right) lead and correlated layer $t_{-1,0}=t_{0,-1}=v_{l}$ ($t_{0,1}=t_{1,0}=v_{r}$). The third term introduces the local Hubbard interactions $U$, which is nonzero only for the correlated layer ($z=0$). The last term describes the onsite energies. For the correlated layer $\varepsilon_{z=0}\equiv\varepsilon_{c}$, and for the left and right lead $\varepsilon_{z<0}\equiv\varepsilon_{l}$ and $\varepsilon_{z>0}\equiv\varepsilon_{r}$, respectively.

The nonequilibrium situation is obtained by applying a bias voltage $\Phi$ and/or temperature difference $\Delta T=T_{l}-T_{r}$, between the leads. Here $T_{l}$ and $T_{r}$ is the temperature of the left and right lead, respectively. The usual way to treat such steady state situation (see, e.g.  Okamoto (2007, 2008)) is to start at $t=-\infty$ with three decoupled ($v_{l}=v_{r}=0$) systems consisting of the two leads plus the central region which are separately in equilibrium at different chemical potentials and temperatures. Then one switches on the hybridisation and waits until a steady state is reached. A bias voltage $\Phi$ is obtained by shifting both the onsite energies $\varepsilon_{l/r}$ of the left and right leads as well as their initial chemical potentials $\mu_{l/r}$ by $\pm|e|\Phi/2$, respectively, i.e. $\varepsilon_{l}=\mu_{l}=|e|\Phi/2$, $\varepsilon_{r}=\mu_{r}=-|e|\Phi/2$. Of course, this neglects the long-range Coulomb interaction, which could be added explicitly, for example, by solving self-consistently the Poisson equation, (see, e.g. Freericks and Zlatić (2007); Okamoto and Millis (2004); Chen and Freericks (2007); Hale and Freericks (2012); Bakalov et al. (2016); Charlebois et al. (2013); Lee and MacDonald (2006, 2007)). One major effect of the long-range Coulomb interaction is a voltage drop across the central region. The approximation made here is to assume that the voltage drop takes place only in the central region, which corresponds to the limit of infinite Coulomb screening.

In the following, we use units in which $t_{c}=1,\hbar=1,k_{B}=1,|e|=1$, and $a=1$. The latter is the distance between neighboring sites of the simple cubic lattice. Thus, the current (density) is expressed in units of

[TABLE]

II.2 Method

Here we give a brief overview of the nonequilibrium DMFT approachSchmidt and Monien ; Freericks, Turkowski, and Zlatić (2006); Freericks (2008); Joura, Freericks, and Pruschke (2008); Eckstein, Kollar, and Werner (2009); Okamoto (2007, 2008) and the auxiliary-master-equation approach (AMEA)Arrigoni, Knap, and von der Linden (2013); Dorda et al. (2014); Titvinidze et al. (2015); Dorda et al. (2017), employed to solve the intrinsic impurity problem. For more details we refer to Refs. Arrigoni, Knap, and von der Linden, 2013; Dorda et al., 2014; Titvinidze et al., 2015, 2016; Dorda, Titvinidze, and Arrigoni, 2016  . To describe the behavior of the steady state it is convenient to work in the Keldysh Green’s function formalism Kadanoff and Baym (1962); Schwinger (1961); Keldysh (1965); Haug and Jauho (1998); Rammer and Smith (1986) and introduce $2\times 2$ block Green’s functions

[TABLE]

which we denote by an underscore. $G^{R}$, $G^{K}$ and $G^{A}$ are the retarded, Keldysh, and advanced Green’s functions. $G^{A}$ and $G^{R}$ are related via $G^{A}=(G^{R})^{\dagger}$. While $G^{K}$ can only be obtained from $G^{R}$ in the equilibrium case via the fluctuation dissipation theoremHaug and Jauho (1998), out of equilibrium $G^{K}$ is independent of $G^{R}$ and it needs to be determined separately, which requires the Keldysh Green’s function formalism.

In the steady state, the two-time Green’s functions only depend on the time-difference and it is convenient to switch to the frequency domain $\omega$. Furthermore, due to translational invariance along the $xy$-plane we introduce the momentum variable ${{\boldsymbol{k}}_{\perp}}=(k_{x},k_{y})$. The Green’s function of the correlated layer is given in terms of Dyson’s equation by

[TABLE]

where ${\underline{g}}_{\alpha}(\omega,{{\boldsymbol{k}}_{\perp}})$ denotes the non-interacting, decoupled ($v_{l}=v_{r}=0$) equilibrium Green’s functions for the correlated layer ($\alpha=c$) and for the edge layer of the left ($l$) and the right ($r$) lead.Gre

The retarded part of the equilibrium Green’s functions can be expressed by

[TABLE]

with $\varepsilon({{\boldsymbol{k}}_{\perp}})=-2(\cos k_{x}+\cos k_{y})$. The sign of the square-root for negative argument in Eq. (8) must be chosen such that the Green’s function has the correct $1/\omega$ behavior for $|\omega|\to\infty$. Since the disconnected leads are separately in equilibrium, we can obtain their Keldysh components from the retarded ones via the fluctuation dissipation theoremGre ; Haug and Jauho (1998)

[TABLE]

where $f_{\alpha}(\omega)$ is the Fermi distribution for chemical potential $\mu_{\alpha}$ and temperature $T_{\alpha}$. For the non-interacting isolated central layer the inverse Keldysh Green’s function $[{\underline{g}}_{c}(\omega,{{\boldsymbol{k}}_{\perp}})]^{-1}$ is proportional to $0^{+}$, an infinitesimal imaginary part, and can be neglected in the final steady state in which the leads and the layers are connected to each other.

${\underline{\Sigma}}(\omega)$ in Eq. (6) stands for the self-energy, which in DMFT is approximated by a local and therefore momentum ($\bf k$) independent quantity. It will be determined self-consistently by solving the SIAM problem with the same local parameters $U$ and $\varepsilon_{c}$ as the original model. More specifically we employ the auxiliary master equation approach (AMEA)Arrigoni, Knap, and von der Linden (2013); Dorda et al. (2014); Titvinidze et al. (2015); Dorda et al. (2017) to determine the self-energy $\Sigma(\omega)$ in the equilibrium and nonequilibrium case. The key point of AMEA is to map the infinite SIAM problem to an auxiliary one with a finite number of bath sites $N_{b}$ and two Markovian environments (sink and reservoir), which are crucial to achieve a steady state in a finite system. The resulting open quantum system is described by a Lindblad equation. Its parameters are, however, not taken from a Born-Markov approximation Breuer and Petruccione (2009); Schaller (2014); Nuss et al. (2015) but are used as fit parameters to optimally reproduce the physical hybridization function

[TABLE]

by the hybridization function $\underline{\Delta}_{\mathrm{aux}}(\omega)$, obtained in the auxiliary system. Here $\underline{g}_{0}^{-1}(\omega)$ is the non-interacting Green’s function of the disconnected impurity and the local Green’s function is obtained by

[TABLE]

Clearly, the accuracy of our impurity solver increases exponentially with the number of bath sites $N_{b}$ and becomes exponentially exact in the limit $N_{b}\rightarrow\infty$.Dorda et al. (2017); Schwarz et al. (2016) In practice $N_{b}=4,6$ is sufficient to obtain reliable results for the current. The self consistence cycle proceeds as follows: Starting out with an initial guess for the self energy, we determine the physical hybridization function in Eq. (10). Next the Lindblad parameters of the auxiliary system without interaction are adjusted by minimizing the misfit between the hybridization function of the auxiliary system and that of the physical system. Once the Lindblad parameters are determined, the auxiliary nonequilibrium many-body system is solved by standard numerical techniques, which defines the self energy for the next iteration. This procedure is repeated until convergence is reached. For more details see Refs. Arrigoni, Knap, and von der Linden, 2013; Dorda et al., 2014; Titvinidze et al., 2015, 2016.

II.3 Thermoelectric properties

Once the Green’s function for the central layer is determined, one can calculate the desired thermoelectric properties. Here we shortly overview the corresponding expressions, which we need in the present work and discuss some of their generic features.

The electrical current density ($J_{\text{c}}$) and the energy-current density ($J_{\text{en}}$) through the central layer in case of the linear response regime or zero bias voltage $\Phi=0$ but with finite temperature difference between leads are given byHaug and Jauho (1998); Jauho (2006) (see Appendix A)

[TABLE]

This expression is valid if ${\rm Im\>}g_{l}^{R}(\omega,{{\boldsymbol{k}}_{\perp}})=\lambda{\rm Im\>}g_{r}^{R}(\omega,{{\boldsymbol{k}}_{\perp}})$, where $\lambda$ does not depend on $\omega$ and $\bf k_{\perp}$. Otherwise one has to use more general expression given by the Meir-Wingreen formula (see e.g. Ref.Haug and Jauho (1998); Meir and Wingreen (1992); Jauho (2006); Ren et al. (2012); Titvinidze et al. (2015)).

For the electric current density (energy current density) we have $\zeta_{\text{c}}=1$ ( $\zeta_{\text{en}}=\omega$). The transmission function is given by

[TABLE]

The sum runs over the ${{\boldsymbol{k}}_{\perp}}$-vectors of the first Brillouin zone of the 2D central layer and $N_{\perp}$ is the number of the corresponding vectors. The transmission function depends on the spectral density of the central layer

[TABLE]

and the imaginary part of the surface Green’s functions of the leads via

[TABLE]

Obviously, the ${{\boldsymbol{k}}_{\perp}}$-dependence of both quantities enters through

[TABLE]

with $\gamma_{\alpha}(\omega,{{\boldsymbol{k}}_{\perp}})=\tilde{\gamma}_{\alpha}(\omega,t_{\alpha}\varepsilon({{\boldsymbol{k}}_{\perp}}))$ and $A(\omega,{{\boldsymbol{k}}_{\perp}})=\tilde{A}(\omega,t_{c}\varepsilon({{\boldsymbol{k}}_{\perp}}))$. Therefore, by replacing $\gamma_{\alpha}(\omega,{{\boldsymbol{k}}_{\perp}})$ and $A(\omega,{{\boldsymbol{k}}_{\perp}})$ in Eq. (14) by $\tilde{\gamma}_{\alpha}(\omega,t_{\alpha}\varepsilon)$ and $\tilde{t}A(\omega,t_{c}\varepsilon)$, respectively, we obtain

[TABLE]

where $\rho_{2D}$ is the 2D density of states corresponding to the dispersion relation $\varepsilon({{\boldsymbol{k}}_{\perp}})$, which can be expressed in terms of the complete elliptic integral of the first kind ${\cal K}(x)$Ell as

[TABLE]

As outlined before, the bias voltage enters the chemical potential of the leads and the onsite energies, both are given by $\Phi/2$ ($-\Phi/2$ ) for the left (right) lead. However, the change of the onsite energies enters only in second order and does not contribute to the linear response expressions. Straightforward Taylor expansion then yieldsRen et al. (2012); Kirchner, Zamani, and Muñoz (2013); Benenti et al. (2013)

[TABLE]

The linear response coefficients are the electric conductivity

[TABLE]

the Seebeck coefficient

[TABLE]

and the electronic contribution to the thermal conductance

[TABLE]

These expressions are based on the Lorenz integral

[TABLE]

with $f(\omega)=1/(\exp(\beta\omega)+1)$.

In addition we also calculated the Lorenz ratio (see Eq. (2)), which in the linear response regime can be expressed as

[TABLE]

We will check the validity of the Wiedemann-Franz law

[TABLE]

in the presence of van Hove singularities and interacting electrons. The Wiedemann-Franz law is assumed to be valid for a free electron gas with elastic scattering. By definition, the transmission function ${\cal T}(\omega)$ is non-negative and so is

[TABLE]

The quantities $I_{m}/I_{0}$ can therefore be considered as the $m$-th moment

[TABLE]

of the probability distribution function

[TABLE]

with the proper normalization constant $Z$. Then the Seebeck coefficient is given by the mean

[TABLE]

and the Lorenz ratio

[TABLE]

by the variance of $\rho(\omega)$ and, consequently, the electronic part of the thermal conductance by

[TABLE]

For low temperatures (degenerate electrons gas), the variance of $\rho(\omega)$ is governed by the factor $|f^{\prime}(\omega)|$, which is strongly localized in the interval ${\cal I}_{T}:=[-T,T]$, in which ${\cal T}(\omega)$ is commonly assumed to be slowly varying. Then we obtain

[TABLE]

which corresponds to the Wiedemann-Franz law. We also realize that deviations are to be expected in the case when $T(\omega)$ is strongly varying in the interval $\omega\in{\cal I}_{T}$, which is the case e.g. for the 2D tight-binding density of states, that has a van Hove singularity at $\omega=0$. We see that the electrical conductivity according to Eq. (21) is always positive. The same holds true, based on Eq. (29), for the Lorenz ratio. Then, due to Eq. (2) also the electronic part of the heat conductivity is positive. Eq. (19) then guarantees that, given a bias voltage and no temperature difference, the electrical current flows from the side of higher to the lower potential energy. In the opposite case, of zero bias voltage but finite $\Delta T$, inserting Eq. (19) into yields Eq. (20)

[TABLE]

and, since the bracketed expression is positive, the energy current flows from the side of higher temperature to that of lower temperature. The same is valid for the heat current, which is given by $\kappa_{e}\Delta T$. A sign change can, however, occur in the Seebeck coefficient, which according to Eq. (28) is a measure of the asymmetry of the transmission function. In the present model ${\cal T}(\omega)$ is symmetric about $\omega=0$ for $\varepsilon_{c}=0$. Qualitatively, ${\cal T}(\omega)$ is shifted to positive frequencies by positive values of $\varepsilon_{c}$ and vice versa. Moreover, by definition, $S$ is an anti-symmetric function in $\varepsilon_{c}$.

The key element of the thermoelectric figure of merit is the Seebeck coefficient in Eq. (28), that increases with the asymmetry of the probability distribution function $\rho(\omega)\propto|f^{\prime}(\omega)|{\cal T}(\omega)$. The first factor is a narrow symmetric peak with exponential tails $e^{-\beta|\omega]}$. The asymmetry enters via the transmission function. If the latter is smooth, the mean $\langle\omega\rangle$ is restricted by $|f^{\prime}|$ and the asymmetry is at most $|k_{B}T|$, resulting in $S\lessapprox 1$. In order to have a figure of merit greater than 1, the transmission function has to overcome the peak $|f^{\prime}(\omega)|$. This can be achieved by a transmission function that is zero up to some energy $\omega_{1}>T$ and has a sharp onset at $\omega_{1}$. In the extreme case that ${\cal T}$ is a box function in the interval $(\omega_{1},\omega_{2})$, with $\omega_{2}\gg\omega_{1}\gg T$ then

[TABLE]

In this case $S=1+\beta\omega_{1}$. This is actually the case in the high ZT materials reported in Ouyang et al. (2016).

For non-interacting electrons in the model under consideration, the transmission function is non-zero in the interval ${\cal I}_{0}:=[-2-\gamma(0)+\varepsilon_{c},2+\gamma(0)+\varepsilon_{c}]$ with Lorentzian edges, whose locations are shifted with $\varepsilon_{c}$ (See Fig. 2). Due to the van Hove singularity the transport function ${\cal T}(\omega)$ has sharper peaks compared to the transport function for a system without van Hove singularity, e.g. with a Bethe lattice DOS (compare Figs. 2(a) and 2(b)). It leads to the different behavior of $S$ when the maximum passes the peak of $|f^{\prime}(\omega)|$. Moving the edges of the transmission function via $\varepsilon_{c}$ across the peak of $|f^{\prime}(\omega)|$ will increase $S$ but only until the edge is out of the energy window $\pm T$, then it will decease again and consequently, $S$ has a pronounced peak but the height will be $S\lessapprox 1$ for reasons outlined before.

III Results

The emphases of the present paper is on the influence of strong electronic correlations on the thermoelectric properties, which will be studied in the framework of the device displayed in Fig. 1, consisting of a correlated layer, with local Hubbard interaction $U$ and on-site energy $\varepsilon_{c}$, coupled to two metallic leads. Of particular interest is the impact of the various model parameters. We will first study the linear response behavior of the system, for which only equilibrium calculations are necessary, before we move on to full nonequilibrium studies for strong temperature or potential differences between the leads.

In all our calculations the hopping between neighboring sites inside the correlated layer is considered as energy unit ($t_{c}=1$), while the hopping between neighboring sites of the leads is always chosen as $t_{l}=t_{r}=2$. As far as the other parameters are concerned, we have distinguished between the linear response and the full non-equilibrium calculation. In the first case we have chosen for the hopping into and out of the central region $v_{l}=v_{r}=0.25$ and for the temperature in the leads $T_{l}=T_{r}=0.025$, while in the second case we have used $v_{l}=v_{r}=1$ and $T_{l}=0.9$ and $T_{r}=0.7$, respectively, in order to enhance the effect of the temperature difference.

III.1 Linear response

We start out with the results for the linear response regime. The temperature in the left and right lead is $T=0.025$ and the hopping in and out of the central layer is $v_{l}=v_{r}=0.25$. Due to the particle-hole symmetry for $\varepsilon_{c}=0$, the electrical and thermal conductance are even functions of the onsite energy $\varepsilon_{c}$, while the Seebeck coefficient is an odd function of $\varepsilon_{c}$. Therefore, we restrict our results to $\varepsilon_{c}\leq 0$, which corresponds to fillings $n\geq 1$ in the central layer.

III.1.1 Non-interacting particles

In order to single out the impact of correlations on the thermoelectric properties, we will first study the present model for $U=0$. In this case the self-energy $\Sigma(\omega)$ vanishes and all required quantities for the integral Eq. (24) are given analytically and we can easily determine the integral by numerical means. The results for the electrical and thermal conductance are presented in Fig. 3. We find that both have a maximum at $\varepsilon_{c}=0$ (half-filling), decrease monotonically with increasing $|\varepsilon_{c}|$, and vanish rapidly beyond $|\varepsilon_{c}|=2$ . As pointed out before, both quantities only depend on $|\varepsilon_{c}|$.

Furthermore, we have examined the validity of the Wiedemann-Franz law (see Eq. (25)). In the inset of Fig. 3 we present $L/L_{0}$ as a function of the onsite energy. We observe deviations from the Wiedemann-Franz law only for values $|\varepsilon_{c}|\approx 0.1$ and $|\varepsilon_{c}|\approx 2.0$. The first deviation is due to a pronounced peak in the transmission function ${\cal T}(\omega)$ that stems from the van Hove singularity of the 2D tight-binding density of states (see Eq. (18)) and the second deviation is due to a sharp edge in ${\cal T}(\omega)$ that stems from the finite width of $\gamma_{l/r}(\omega,t_{\alpha}\varepsilon)$. The latter is given by the imaginary part of the Green’s function of the leads. In the wide band limit it is proportional to

[TABLE]

On the other hand, the Seebeck coefficient $S(\varepsilon_{c})$, which is depicted in Fig.  4, is a non-monotonic function and has two maxima: the first maximum is close to half-filling $\varepsilon_{c}\simeq-0.1$, and the second at $\varepsilon_{c}\simeq-2$. As pointed out in section II.3, maxima in $S$ occur, when sharp edges or maxima in the transmission function are moved through the energy window ${\cal I}_{T}$ defined by $|f^{\prime}|$. The situation is illustrated in Fig.  2, where $|f^{\prime}(\omega)|$ and the transmission function ${\cal T}(\omega)$ are shown for different values of $\varepsilon_{c}$. Obviously, the first maximum is due the van Hove singularity, which for $\varepsilon_{c}\simeq-0.1$ crosses the boundary of ${\cal I}_{T}$. The second maximum occurs at $\varepsilon_{c}\simeq-2$, when the edge of the transmission function crosses the boundary. The second maximum of the Seebeck coefficient is higher, because the edge of ${\cal T}$ is more pronounced than the van Hove singularity, resulting in a greater mean $\langle\omega\rangle$.

In order to underpin the impact of the van Hove singularity, we have also performed calculations for the density of states $\rho_{B}(\epsilon)$ of the Bethe lattice, which has no singularity at $\omega=0$. $\rho_{B}(\epsilon)$ is chosen such that it has the same band width as the original density of state. We find that both low-frequency features, the maximum of the Seebeck coefficient and the discrepancy from the Wiedemann-Franz law, disappear.

The features at $\varepsilon_{c}=-2$ survive as they are due to band edge effects.

Finally, still for the non-interacting model, we also present the dependence of the the Seebeck coefficient on the hybridization $v_{l}=v_{r}=v$ in the inset of Fig. 4. Clearly, the Seebeck coefficient increases with decreasing hybridization $v$.

III.1.2 Linear response for finite interaction

Next we discuss the results for the interacting system. In the interacting case, we first have to solve the self-consistent impurity problem and calculate the self-energy $\Sigma(\omega)$, before we can calculate thermoelectric properties introduced in Sec. II.3. In Fig. 5 we show the electrical conductance $\sigma$ and the Lorenz number as a function of $\varepsilon_{c}$ for different values of the local interaction $U$. For small to moderate interactions ($U\lesssim 5.5$) we find, like in the non-interacting case, a decreasing conductance with increasing $\varepsilon_{c}$, which vanishes rapidly beyond $|\varepsilon_{c}|=2+U/2$. As a new aspect, the maximum at $\varepsilon_{c}=0$ decreases with increasing $U$. For strong interactions ($U\gtrsim 5.5$), on the other hand, $\sigma$ behaves differently. A local minimum develops at $\varepsilon_{c}=0$ and a maximum at some finite value $\varepsilon^{*}_{c}$. Interestingly (not depicted in the plot), the filling corresponding to $\varepsilon^{*}_{c}$ is in all cases roughly the same $n\simeq 1.3$. Qualitatively, such non-monotonic behavior can be explained by the fact that there are two competing effects. On the one hand, like in the non-interacting case, with increasing $|\varepsilon_{c}|$ and correspondingly $n$ the conductance should decrease due to the decreasing overlap of the spectral function of the central layer with the imaginary part of the leads Green’s functions. On the other hand, close to half-filling ($\varepsilon_{c}=0$) correlation effects are most dominant and therefore due to backscattering, the conductance should be reduced or even suppressed. The latter is actually the case for very strong interactions ($U\geq 12$), indicating that the isolated correlated interface would be in the Mott insulator phase, but the system always remains metallic when the correlated layer is coupled to the leads. In the limit of strong interaction the conductance is strongly suppressed, but it never becomes equal to zero. The Mott transition is nicely visible in the transmission function, which is shown in Fig.  7, where a gap is found at $U=12$.Gap As far as the Lorenz number is concerned, which is depicted in the lower half of Fig.5, the deviation from the Wiedemann-Franz law, that we found in the non-interacting case at $\varepsilon_{c}\approx 0$, decreases with increasing interaction strength, because the van Hove singularity is smeared out by the interaction.

As a final point on linear response, we turn to the Seebeck coefficient (see Fig. 6). We find that for $U\lesssim 5.5$ the Seebeck coefficient is positive, as in the non-interacting case. With increasing interaction the maxima are shifted to higher $\varepsilon_{c}$-values and are suppressed. Along the same line of reasoning we gave in the non-interacing case, this behavior can easily be explained by the smearing of the transmission function due to interactions. However, a qualitatively new aspect comes into play for strong interactions ($U\gtrsim 5.5$). Now we find that for small onsite energies $\varepsilon_{c}^{\star}\leq\varepsilon_{c}\leq 0$ the Seebeck coefficient becomes negative, before it turns positive again for some value $\varepsilon_{c}^{\star}$. This is in particular the case for $U=8$ and $U=12$ in Fig.  6, while $U=4$ still exhibits a purely positive Seebeck coefficient. To understand the qualitative difference in the Seebeck coefficient corresponding to these $U$ values, we consider the corresponding transmission functions in Fig.  7. We observe that the central peak (van Hove singularity) in the transmission function, which is still predominant for $U=4$, decreases with increasing $U$ and the lower and upper Hubbard band grow, until eventually the Mott-Hubbard gap opens up (See Fig. 7). At the same time, for $U\gtrsim 5.5$, spectral weight is transferred from the lower to the upper Hubbard band, when $\varepsilon_{c}<0$. Consequently, based on Eq. (27), the mean $\langle\omega\rangle$ moves to positive energies, resulting in a negative Seebeck coefficient, absolute value of which increases with increasing $U$. The detailed dependence of the transmission function on $\varepsilon_{c}$ for the case $U=8$ is given in Fig.  8. It is obvious, why the Seebeck coefficient $S=-\beta\langle\omega\rangle$ is increasingly negative for increasing $\varepsilon_{c}$. In addition to the transfer of spectral weight from the lower to the upper Hubbard band, the entire transmission function shifts to the left with increasing $\varepsilon_{c}$, like in the non-interacting case. This shift eventually, for $\varepsilon_{c}<\varepsilon_{c}^{\star}<0$ results in a positive Seebeck coefficient. As we can clearly see in Fig.  8 for the $U=8$ case, the edge close to $\omega=0$ becomes increasingly pronounced with increasing $\varepsilon_{c}$-values. At the same time the Seebeck coefficient increaeses. This is in nice agreement with our general discussion that the Seebeck coefficient is large if a sharp edge or peak is shifted (by changing some parameter) across the energy window, defined by $|f^{\prime}|$. It is particularly interesting to note that the absolute height of the first peak in the Seebeck coefficient is greater for $U=12$ than in the non-interacting case. According to Eq. (1) the thermoelectric figure of merit depends quadratically on $S$, therefore we find that strong correlations would increase the figure of merit, taking into account that the Lorenz number even slightly supports this effect.

The behavior of $ZT$ is completely determined by the behavior of the Seebeck coefficient $S$ and the Lorentz number $L$. Even more, based on the sign of the Seebeck coefficient one can determine whether the conductance is electron like or hole like, which is impossible to do based on the figure of merit. Therefore we prefer to present the results for the Seebeck coefficient $S$ and the Lorentz number $L$ separately.

III.2 Finite temperature difference

Following the discussion of linear response results, we investigate the effect of a finite temperature difference between the leads. We start out with zero bias voltage ($\Phi=0$), where the nonequilibrium situation is driven merely by the temperature difference (${\Delta T>0}$). To emphasize the effect of $\Delta T$ on the the behavior of the system, in our calculations the large temperature difference. Therefore, the temperature values in the leads are $T_{l}=0.9$ and $T_{r}=0.7$, respectively. It should be stressed, that the average temperature is much higher (non-physical) than the one we used in the linear response calculations, which was $T=0.025$. To obtain more smooth results we consider the hopping into and out of the central region is now $v_{l}=v_{r}=1$ instead of $0.25$, which we used in the linear regime. In Fig. 9 we present the current density $J$ as a function of the filling $n$ for different Hubbard interactions $U$. Due to particle-hole symmetry ${J(n)=-J(2-n)}$.

Again, in the non-interacting case, $U=0$, the self-energy is zero and all expressions are analytically available. In this case, the current is positive above the half-filling ($n>1$) and negative in the opposite case. We find that the current has a bell shape $n$-dependence, with a maximum at $n\simeq 1.375$, which corresponds to an onsite energy $\varepsilon_{c}\simeq-1.3$. The direction of the current can be explained by the two factors $\Delta f:=f_{l}(\omega)-f_{r}(\omega)$, which is an antisymmetric function about $\omega=\mu$ (positive above and negative below $\mu$) and ${\cal T}(\omega)$ in Eq. (12). Above half filling ($n>1$), the transmission function ${\cal T}(\omega)$ has negative slope at the chemical potential (see e.g. Fig.  3 for $\varepsilon_{c}<0$) and therefore the integral $\int\Delta f(\omega){\cal T}(\omega)d\omega$ is dominated by the negative part of the integrand. Or in more physical terms: electrons move from the cold right lead to the warm left lead, which can also be described by holes moving from the warm left lead to the cold right lead and $J>0$. Below half-filling, the opposite is the case. The slope of ${\cal T}(\omega)$ at the chemical is positive and, therefore, the current integral is dominated by the positive parts, which can be described by electrons moving from the warm left lead to the cold right lead and $J<0$.

Next we turn to the interacting case. For $U=8$, we have first tested the convergence of our numerical scheme with respect to the number of bath sites in our impurity solver. In Fig.  9 the corresponding results are depicted for $N_{b}=4$ and $N_{b}=6$. We find that the resulting curves can hardly be distinguished. Therefore, we restrict all further results to $N_{b}=4$.

For the weak to intermediate interactions ($U\lesssim 5.5$), as might be expected, the behavior of the current $J$ as a function of the onsite energy $\varepsilon_{c}$ is qualitatively similar to the non-interacting case.

The impact of $U$ is a shift of the current maximum to higher $\varepsilon_{c}$-values, or rather higher fillings, and the value of the maximum decreases slightly. This behavior can qualitatively be explained by our previous findings in the linear response case. Accordingly it the expression for the current, driven by the temperature difference, is $J\propto\sigma S$. Taking into account the much higher average temperature in the present case, the bell shaped structure, the shift of the maximum to higher $\varepsilon_{c}$ values, and the reduction of the maximum, can be inferred from Fig.  5 and Fig.  6. The linear response expression would predict , due to the features of $S$, that for stronger interaction the current is negative for $\varepsilon_{c}$ below a certain threshold $\varepsilon_{c}^{\star}$ and positive above. This is indeed the result, found in the full nonequilibrium many-body calculation, depicted in Fig.  9. We obtain also that for strong coupling, the absolute value of the current increases with increasing interaction strength. In the vicinity of half-filling ($1<n\lesssim 1.05$) the current is larger as compared to the non-interacting system with the same filling. With further increase of the filling $n$, correlation effects become less relevant, therefore the current decreases and changes sign. The change of sign of the current takes place for the higher fillings with increasing the interaction. Finally, for even higher fillings, close to the full filling ($n\gtrsim 1.6$), interaction plays even less important role in the behavior of the system and correspondingly the current is nearly independent from it.

Finally, we consider the case, when in addition to the finite temperature difference mentioned above also a finite bias voltage $\Phi$ is applied, and investigate the current-voltage characteristics. In this case Eq. (12) is not valid anymore, therefore we used the more general Meir-Wingreen formulaHaug and Jauho (1998); Meir and Wingreen (1992); Jauho (2006); Ren et al. (2012); Titvinidze et al. (2015). In the upper panel of Fig.  10 the current density is plotted versus applied voltage for different fillings, which are tuned by $\varepsilon_{c}$. The fillings are chosen such, that the current without bias voltage is negative. We find a linear increase, which agrees with the linear response result, obtained by Eq. (19).

From the current-voltage curve, we can directly compute the power $P=J\cdot\Phi$. Obviously, it is negative for ${0<\Phi<\Phi_{0}(\varepsilon_{c},U)}$ (see Fig. 10[lower panel]), which means that energy conversion from the system takes place. The maximum efficiency of the device is reached for that bias voltage for which the absolute value of the power $P$ reaches its maximum. One can readily see that this is the case for ${\Phi\simeq\frac{1}{2}\Phi_{0}(\varepsilon_{c},U)}$. The corresponding power is $P\simeq\frac{1}{2}J_{0}(\varepsilon_{c},U)\Phi_{0}(\varepsilon_{c},U)$. Our calculations have shown that for $U\gtrsim 5.5$ the maximum value of $|J_{0}(\varepsilon_{c},U)|$ as well as $\Phi_{0}(\varepsilon_{c},U)$ increases with increasing interaction and the optimal filling varies between $1.04<n<1.08$. So, similar to what we have found in the linear response case, the efficiency of the thermoelectric device can be increased in strongly correlated materials close to half-filling.

IV Conclusions

In this work we have investigated linear response and steady state properties of a device, consisting of a correlated mono-layer attached to two metallic leads. In the full many-body nonequilibrium calculation the nonequilibrium properties were driven by a temperature difference in the leads and a bias voltage. The key quantity, as far as the thermoelectric efficiency is concerned, is the thermoelectric figure of merit, which is proportional to the square of the Seebeck coefficient and inverse proportional to the Lorenz number. As the latter is very well described by the Wiedemann-Franz law, the efficiency is governed solely by the Seebeck coefficient. We have given some general arguments that a large Seebeck number requires sharp structures in the transmission function, which indeed occurred for strong interaction strength $U$. We have shown that the efficiency of the thermoelectric device can furthermore be increased by increasing the overall temperature and/or decreasing the coupling strength between device and leads. But more importantly, we have demonstrated, both in the linear response regime and beyond, that the efficiency of the device can be increased by strong correlations close to half-filling.

Appendix A Generalisation of the Meir-Wingreen formula for a layered system

Here we generalize the Meir-Wingreen formula to the case of a layered system and derive equation (12).

We start out from the expressionOkamoto (2007); Not for the current from layer $z$ to layer $z+1$, for which the hopping matrix element in the second term of equation (3) is $t_{z,z+1}$

[TABLE]

We have exploited the fact that the Keldish Green’s function ${\bf G}^{K}$ is antihermitian ($G^{K}_{z+1,z}=-(G^{K}_{z,z+1})^{*}$). Due to the conserved current the results are independent of $z$. Then we find for the current $I_{r/l}$, that flows from the right/ left lead into the central layer, the expression

[TABLE]

Next we will express the off-diagonal elements $G^{K}_{\pm 1,0}$ by diagonal elements of suitable Green’s functions.

Since the first index of the Green’s function ($z=\pm 1$) belongs to sites of the left or right lead where there is no Hubbard interaction, the equation of motion Niu, Lin, and Lin (1999) simply yields

[TABLE]

were we have used that $g^{\gamma}_{\mp 1,\mp 1}=g^{\gamma}_{l/r}$ is the surface Green’s function of the isolated left/ right lead. Based on the general property of Green’s functions $G^{A}=(G^{R})^{\dagger}$ and on definition in equation (16) we find

[TABLE]

where we have introduced the lesser Green’s function $G^{<}$, for which the general relation $G^{<}=G^{K}_{c}-(G_{c}^{R}-G_{c}^{A})$ applies. Here $f_{l/r}$ is Fermi functions for the left and right lead, respectively. Inserting the above expressions for $G^{K}_{\mp 1,0}$ into the Meir-Wingreen formula for the current $I_{l/r}$ (equation (33)) we obtain

[TABLE]

Due to the current conservation we have $I=I_{l}=-I_{r}$ and for arbitrary $x$ we can also express the current as $I=xI_{l}-(1-x)I_{r}$. Following the restricting assumption of JauhoJauho (2006) that $\gamma_{l}=\lambda^{\prime}\gamma_{r}$, we obtain

[TABLE]

Now if we fix the arbitrary parameter by $x=1/(1+\lambda^{\prime})$, so that the term containing $G^{<}$ vanishes, and taking into account the relation $G_{c}^{R}-G_{c}^{A}=-2iA$, we obtain

[TABLE]

In the units, used in this paper, this is nothing else but Eq. (12) along with equation (14), if the $k_{||}$-integration is replaced by the energy integration.

Acknowledgements.

This work was supported by the Austrian Science Fund (FWF): P26508, as well as SfB-ViCoM project F04103, and NaWi Graz. The calculations were partly performed on the D-Cluster Graz and on the VSC-3 cluster Vienna.