On the delta function broadening in the Kubo-Greenwood equation

TL;DR

This paper presents a new nonlinear extrapolation method for calculating DC electrical conductivity from first principles, reducing finite-size effects and ambiguities in the Kubo-Greenwood approach.

Contribution

A novel nonlinear extrapolation scheme for the Kubo-Greenwood equation that improves accuracy in finite systems and aligns well with Landauer-based methods.

Findings

Method predicts conductivity close to the thermodynamic limit for small systems.

Overcomes common ambiguities in finite-system conductivity calculations.

Results are consistent with Landauer equation-based approaches.

Abstract

Understanding DC electrical conductivity is crucial for the study of materials. Macroscopic DC conductivity can be calculated from first principles using the Kubo-Greenwood equation. The procedure involves finding the thermodynamic limit of the current response to an electric field that is slowly switched on, and then taking the limit of the switching rate to zero. We introduce a nonlinear extrapolation procedure executed in systems with periodic boundary conditions, which predicts conductivity close to the thermodynamic limit even for very small systems. The scheme also overcomes a large part of the usual ambiguities of the DC conductivity definition for finite systems. We numerically compare our method to the Landauer equation-based approach and find both techniques to be consistent with each other.

Tables

Table 1. Table 1: Predicted DC conductivity as a function of the broadening function δ η ​ ( ω ) subscript 𝛿 𝜂 𝜔 \delta_{\eta}(\omega) assuming Drude equation for conductivity ( 12 ).

Function name	${\delta }_{\eta }\left(\omega \right)\eta$	${\delta }_{\eta }\left(t\right)$	Predicted ${\sigma }_{DC}\left(\eta \right)$
Sinc	$\frac{\sin \left(2\omega /\eta \right)}{2\pi \omega /\eta }$	$\theta \left(2{\eta }^{-1}-\left|t\right|\right)$	$A\left(1-e^{2/B\eta }\right)$
Lorentzian	$\frac{1}{\pi \left({\omega }^2/{\eta }^2+1\right)}$	$e^{-\eta \left|t\right|}$	$\frac{A}{1+B\eta }$
Gaussian	$\frac{1}{\sqrt{\pi }}{e}^{-{\omega }^2/{\eta }^2}$	$e^{-\frac{1}{4}{x}^2{\eta }^2}$	$A\frac{\sqrt{\pi }{e}^{\frac{1}{B^2{\eta }^2}}\text{erfc}\frac{1}{B\eta }}{B\eta }$

Equations

$$\mathrm{Re}[\sigma (\omega )]=\sum _{mn}{\Gamma }_{mn}\delta (E_m-E_n-\omega ),$$

$$j(t)={\int }_{-\infty }^t\sigma (t-t^{\prime })E(t^{\prime })d{t}^{\prime }+O(E^2),$$

$${\sigma }_{DC}\equiv \underset{\eta \to 0}{\lim }\left(\underset{L\to \infty }{\lim }\text{Re}[\sigma (0){]}_{\eta }\right),$$

$$\text{Re}[\sigma (\omega )]=A/(1+B^2{\omega }^2),$$

$$j(t)={\int }_{-\infty }^t\sigma (t-t^{\prime })E(t^{\prime })d{t}^{\prime },$$

$$j(0)=2\pi {\int }_{-\infty }^{\infty }\text{Re}[\sigma (\omega )]E(\omega )d\omega ,$$

$$\text{Re}\left[\sigma (\omega )\right]=-\frac{\pi }{Z}\sum _{nm}\frac{e^{-\beta E_n}-e^{-\beta E_n}}{E_n-E_m}{j}_{nm}{j}_{mn}\delta (E_n-E_m-\omega ),$$

$$\text{Re}\left[\sigma (\omega )\right]=-\pi \sum _{nm}\frac{f_n-f_m}{{\epsilon }_n-{\epsilon }_m}{j}_{nm}{j}_{mn}\delta ({\epsilon }_n-{\epsilon }_m-\omega ),$$

$${\Gamma }_{mn}=-\pi \frac{f_n-f_m}{{\epsilon }_m-{\epsilon }_n}{j}_{nm}{j}_{mn}.$$

$$j(0)=\sum _{mn}{\Gamma }_{mn}2\pi E({\epsilon }_n-{\epsilon }_m-\omega ){|}_{\omega =0}.$$

$$\underset{L\to \infty }{\lim }\text{Re}[\sigma (0){]}_{\eta }\approx \text{Re}[\sigma (0){]}_{\text{Num}\eta }$$

$$\text{Re}[\sigma (\omega ){]}_{\text{Drude}}=A/(1+B^2{\omega }^2).$$

$$\mathrm{Re}[\sigma (0){]}_{\text{Drude}\eta }=\int \frac{A}{1+B^2{\omega }^2}{\delta }_{\eta }(\omega )d\omega .$$

$$t_{\text{crit}}=\frac{R}{2v_{\max }},$$

$$⟨\sigma ⟩(\omega )\equiv \int P(\mathbf{\Gamma },\mathbf{\Delta })\text{Re}[\sigma (\omega )]d\mathbf{\Gamma }d\mathbf{\Delta }$$

$$⟨\sigma ⟩(\omega )=\int P(\mathbf{\Gamma },\omega )\mathbf{\Gamma }d\mathbf{\Gamma }.$$

$$⟨{\sigma }_{\eta }⟩(\omega )=(⟨\sigma ⟩\ast {\delta }_{\eta })(\omega ),$$

$${\delta }_{\eta }(\Delta E-\omega )=(\delta \ast {\delta }_{\eta })(\Delta E-\omega )$$

$$⟨{\sigma }_{\eta \beta }⟩(\mu )=(⟨{\sigma }_{\eta \infty }⟩\ast f_{\beta }^{\prime })(\mu ).$$

$$⟨{\sigma }_{\eta \beta }⟩(\mu ,\omega )=(⟨{\sigma }_{\infty }⟩\ast (f_{\beta }^{\prime }{\delta }_{\eta }))(\mu ,\omega ).$$

$$R=R_{\text{bulk}}+R_{\text{contacts}}$$

$${\sigma }_L\equiv S{\left(\frac{d{R}_{\text{bulk}}}{dL}\right)}^{-1}=S{\left(\frac{dR}{dL}\right)}^{-1}.$$

$$⟨\hat{X}⟩(t)=⟨{\underrightarrow{e}}^{i\epsilon {\int }_{-\infty }^t{\hat{V}}_I(t_1)f(t_1)d{t}_1}\hat{X}{\underleftarrow{e}}^{-i\epsilon {\int }_{-\infty }^t{\hat{V}}_I(t_1)f(t_1)d{t}_1}{⟩}_0.$$

$${\underleftarrow{e}}^{{\int }_0^t\hat{V}(t_1)d{t}_1}\equiv \underset{\delta \to 0}{\lim }(1+\delta \hat{V}(t-\delta ))...(1+\delta \hat{V}(\delta ))(1+\delta \hat{V}(0)).$$

$$⟨\hat{X}⟩=⟨\hat{X}{⟩}_0-i\epsilon {\int }_{-\infty }^0⟨[\hat{X},{\hat{V}}_I(t_1)]{⟩}_{0}f(t_1)d{t}_1$$

$$⟨\hat{X}⟩=⟨\hat{X}{⟩}_0+\frac{1}{Z}\sum _{mn}\frac{e^{-\beta E_n}-e^{-\beta E_m}}{E_n-E_m-\omega +i\eta }{X}_{nm}{V}_{mn}$$

$$t(\vec{A})=t{e}^{i(\vec{A},\vec{a})}$$

$$\hat{j}=-\frac{1}{V}{\partial }_A\hat{H},$$

$$\hat{H}={\hat{H}}_0-AV{\hat{j}}_0,$$

$$\hat{j}={\hat{j}}_0-A\hat{M}$$

$$⟨\hat{j}⟩=-⟨\hat{M}{⟩}_{0}A(t)+iV{\int }_{-\infty }^0⟨[{\hat{j}}_0,{\hat{j}}_{0I}(t_1)]{⟩}_{0}A(t_1)d{t}_1.$$

Full text

††journal: Computer Physics Communication

1 Introduction

This work was motivated by the issue of ambiguity of the definition of DC conductivity in finite-sized quantum systems, calculated by Kubo’s linear response theory [1]. The Kubo-Greenwood formula [2, 3] expresses the real part of AC electrical conductivity $\sigma(\omega)$ as a sum of delta functions:

[TABLE]

where $E_{n}$ are eigenstate energies, $\omega$ is the frequency of the external electric field, and $\Gamma_{mn}$ depends on the Hamiltonian and temperature of the system, but not on $\omega$. Multiple equivalent forms of Eq. (1) can be found in the literature [1, 3, 4, 5, 6, 7, 8, 9, 10, 11]. Conductivity $\sigma(\omega)$ is defined as a Fourier transform of the time dependent $\sigma(t)$, which a part of the linear response equation for current density $j(t)$:

[TABLE]

where $E(t)$ is a time dependent electric field.

Due to the delta function in Eq. (1), DC conductivity is always zero for closed and finite systems. A non-zero DC conductivity can be defined rigorously as a result of taking two consecutive limits, one of which includes taking the size of the system to infinity:

First, one defines a pseudo-conductivity $\text{Re}[\sigma(\omega)]_{\eta}$, where each delta function in Eq. (1) is exchanged for a smooth representation $\delta(\omega)\rightarrow\delta_{\eta}(\omega)$, where $\eta$ is an effective width of the representation. For example $\delta_{\eta}(\omega)=\eta^{-1}\pi^{-1/2}e^{-\omega^{2}/\eta^{2}}$ is a valid smearing function. 2. 2.

A finite temperature DC conductivity is defined as a result of two consecutive limits:

[TABLE]

where $L$ is the linear size of system.

In experiments, DC conductivity of finite systems is not zero, which pushes us to modify Eq. (1). Theoretically, finite-sized DC conductivity can be defined by employing additional constructs like infinite electric leads or a thermostat [12]. Although these definitions allow us to talk about DC conductivity of a finite system, each of them involves some freedom and can yield different results, i.e., there is no unambiguous definition of DC conductivity for finite systems.

In this paper, we are primarily concerned with the definitions based on the artificial broadening of the delta function in Eq. (1) [1, 3, 4, 5, 6, 7, 8, 9, 13]. The methods are usually reasoned empirically: a physically meaningful DC conductivity can be obtained from the Kubo-Greenwood equation only if the delta functions are broadened. If the broadening is small, $\eta\ll\Delta E$ ($\Delta E$ is the distance between energy levels), conductivity $\text{Re}[\sigma(\omega)]_{\eta}$ will experience strong oscillations due to the discrete energy spectrum. If the broadening is large, the features of conductivity, like resonant peaks, may get washed out. Therefore, one hopes to have a region of the broadening $\eta$, where the oscillations from the discrete spectrum are suppressed, but macroscopic features are not washed out. DC conductivity then can be defined as a value $\text{Re}[\sigma(0)]_{\eta}$ obtained from this region. Unfortunately, as Calderin et al. demonstrated in [1], such an approach may yield significantly different values of DC conductivity for small systems depending on the form of the broadening function $\delta_{\eta}(\omega)$ and the chosen value of $\eta$. The problem is more pronounced for systems with a sharp conductivity peak at $\omega=0$. The peak gets quickly suppressed if the delta function is broadened but has a width comparable to the distance between energy levels. Obtaining DC conductivity for such systems may be difficult and requires consideration of very large systems.

As an illustration of the scope of the issue, we plotted DC conductivity calculated by equation (1) with different types of smearing for a tight-binding nearest-neighbor hopping 3D system with a small onsite disorder. Conductivity as a function of the broadening parameter $\eta$ is presented on Fig. 1a. From the picture, apart from the arguable absence of the plateau for each curve, there is another problem: the peak heights for different types of broadening differ up to 30%, and do not converge as the system size increases (see Fig. 1b). Considering larger systems quickly becomes computationally intractable without radical optimisations (see Kernel Polynomial Method [14]), since the eigenvalue problem scales as $O(L^{9})$. Extracting DC conductivity by the basic method from small systems seems impossible.

In this paper we provide a solution to the aforementioned problem which allows extraction of macroscopic DC conductivity even from tiny systems. Instead of a peak height, we propose to define a finite-sized DC conductivity based on a specific non-linear extrapolation of $\text{Re}[\sigma(0)]_{\eta}$ from finite $\eta$ to zero. The scheme is based on the equivalence between delta function broadening and current response to variable electric field: $\text{Re}[\sigma(0)]_{\eta}$ is a current response to electric field $E(\omega)=\delta_{\eta}(\omega)$ at $t=0$. If $\text{Re}[\sigma(0)]_{\eta}$ has reached the thermodynamic limit and macroscopic conductivity follows the Drude equation

[TABLE]

where $A$ and $B$ are constants, we can predict $\text{Re}[\sigma(0)]_{\eta}$ by calculating response to the corresponding electric field and then obtain DC conductivity $\sigma_{\text{DC}}=A$ by extrapolating from finite $\eta$ to zero. We call the procedure a Drude-based extrapolation because of the assumption of macroscopic conductivity to follow Drude equation.

It turns out if the temperature is high enough, the result of the extrapolation is practically independent of the choice $\delta_{\eta}(\omega)$, and the value of finite-size DC conductivity converges rapidly to the thermodynamic limit, i.e., we get close to the limit even for very small systems. Conductivity obtained by such extrapolation also matches well with values obtained from the Laundauer-based approach [12], where the system is attached to infinite electric leads, and DC conductivity is extracted from the change of resistance, as the size of the system increases. In summary, for finite system sizes our approach overcomes a large part of the usual ambiguities and provides a reasonable approximation of the macroscopic DC conductivity.

The paper is structured as follows: We start with a quick review of the Kubo-Greenwood equation and various associated definitions in Section 2. We then proceed to the description of the method and numerical tests in Section 3. Then we apply the method to two other toy models: one is a 3D system with point defects, and the other is a 1D disordered chain.

2 Summary of the Kubo-Greenwood equation

In this section, we quickly summarise the general notions of linear response and the Kubo-Greenwood equation.

2.1 General results from linear response

Linear response of the current of the system to an electric field can be expressed as follows:

[TABLE]

where $E(t)$ is external electric field and $\sigma(t)$ is conductivity. Generally $\sigma(t)$ is a tensor, but we consider only an isotropic case for simplicity. For the same reason, we also assume electric field and current to be directed along the $x$ axis. Without loss of generality we can also choose $t=0$ and assume $E(t^{\prime})=E(-t^{\prime})$, since due to causality, fields at $t>0$ should not influence current at $t=0$. Under the above condition, it is possible to express current at $t=0$ in Eq. (5) using only real parts of $\sigma$ and $E$ in Fourier space

[TABLE]

where $\sigma(\omega)=\frac{1}{2\pi}\int\sigma(t)e^{-i\omega t}dt$. Notice how Eq. (6) provides a direct way of expressing current using only real part of conductivity.

2.2 Kubo-Greenwood equation and associate definitions summary

In all equations we assume $\hbar=e=1$, where $\hbar$ is Plank constant and $e$ is electron charge.

Kubo-Greenwood equation predicts $\text{Re}[\sigma(\omega)]$ for a quantum system with finite number of energy levels (see Appendix A for derivation):

[TABLE]

where $E_{n}$ are energies of the system, $Z\equiv\sum_{n}e^{-\beta E_{n}}$ is the partition function, $j_{nm}\equiv\langle n|\hat{j_{x}}|m\rangle$ are matrix elements of the current density operator $\hat{j}_{x}=-V^{-1}\partial{\hat{H}}/\partial{A_{x}}$, where $A_{x}$ is a vector potential along the x-axis, $\hat{H}$ is the Hamiltonian and $V$ is the volume of the system. The dimensionality of the system is arbitrary.

If the system consists of non-interacting particles, equation (7) reduces to the following form

[TABLE]

where $\varepsilon_{n}$ are one-electron energies, and $f_{n}=(1+e^{\beta(\varepsilon_{n}-\mu)})^{-1}$ is the Fermi distribution.

Connecting this equation to the Eq. (1) in the introduction, we obtain

[TABLE]

Substituting equation for conductivity (8) into (6) yields the prediction for the current at $t=0$

[TABLE]

The form of this equation coincides with the equation for DC conductivity (1), where each delta function is substituted by $\delta(\omega)\rightarrow 2\pi E(\omega)$. Smearing delta function by a narrow function $\delta(\omega)_{\eta}$ corresponds to a wide function $E(t)$ in time domain, or, in other words, slowly switched electric field.

In thermodynamic limit, one usually expects a constant current in response to the constant electric field while in finite systems the current starts dropping after the time $t\gtrsim 1/\Delta E$, where $\Delta E$ is an energy difference between nearby energy levels. This can be seen from the form of equation (10). A long-lasting electric field corresponds to a narrow delta function. Due to the repulsion of the energy levels, probability density that a delta function is located at $\omega=0$ is zero, so the DC conductivity will tend to zero, as the width of delta function diminishes.

3 Formulation of the method

This section explains the motivation behind our method and demonstrates its numerical advantage over the classical approach.

3.1 Drude based extrapolation

Our prime interest was macroscopic DC conductivity obtained from the Eq. (3). Suppose we found a numerical estimate for thermodynamic limit

[TABLE]

for some values of $\eta$. In order to find the second limit $\eta\to 0$, let us assume that the “true” conductivity of a macroscopic system follows the Drude equation for small $\omega$:

[TABLE]

Remembering the duality between variable electric field and broadening of the delta function, and using Eq. (6), the Drude equation (12) predicts conductivity for arbitrary $\delta_{\eta}(\omega)$:

[TABLE]

We can now compare the prediction to our numerical results by first finding the best-fit parameters for $A$ and $B$, and then analyzing how well the Drude model describes the data and whether other theories can describe the data better. We call the fitting procedure a Drude-based extrapolation because the estimate of the DC conductivity $\sigma_{\text{DC}}=A$ is the result of the extrapolation of ${\text{Re}}[\sigma(0)]_{\eta}$ from finite $\eta$ to zero.

We performed the extrapolation for the example described on Fig. 1. Predicted DC conductivity as a function of the fitting region for several types of broadening $\delta_{\eta}(\omega)$ (see Table 1) is presented in Fig. 2(a). The plot contains a region of plateau, where predicted DC conductivity is almost independent of the type of broadening and the extrapolation region. In the calculation, we assumed that the thermodynamic limit for ${\text{Re}}[\sigma(0)]_{\eta}$ is already reached for each $\eta$. This condition seems to be violated for small $\eta$, where predicted conductivity suddenly drops. Certain physical arguments (see Appendix B) suggest that the thermodynamic limit of ${\text{Re}}[\sigma(0)]_{\eta}$ is practically reached if the corresponding excitation time of the system is below

[TABLE]

where $v_{\max}$ is the maximum group velocity of the band in the direction of the electric field, $R=(L-\lambda_{\beta})/2$ is the effective radius of the system. Here $L$ is the linear size of the system and $\lambda_{\beta}\approx\sqrt{|t|\beta/2}$ is a diffusion length, defined from the evolution of a one-particle equation $\dot{\psi}=-\hat{H}\psi$ at a time $t=\beta/2$ (see Appendix B for details). The critical time explains why for the Gaussian and Lorentzian broadening the plateaus are somewhat less flat compared to using $\mathrm{sinc}$. The former two functions always have tails in the time domain $t>t_{\text{crit}}$ (see Table 1), while the sinc function is exactly zero for $|t|>2\eta^{-1}$. Note that the relevant values of $\eta$ in our method are much larger than $\Delta E$. Also, the values of the DC conductivity that our method predicts are almost an order of magnitude larger compared to the classical approach, and, as will become clear soon, much closer to the thermodynamic limit.

The consistency of the Drude-based extrapolation can be confirmed further by studying the convergence of the predicted conductivity with the increase of the size of the system. Fig. 2(b) shows DC conductivity as a function of the chemical potential for different sizes of the system. For the considered example, the Drude-based extrapolation converges to the thermodynamic limit within 5$\%$ for a cube as small as $6\times 6\times 6$. The strong oscillations of conductivity for the system of the size $4\times 4\times 4$ can be suppressed by increasing the temperature.

Let us emphasize that the quality of the extrapolation depends strongly on the chosen model. Using a Gaussian model for conductivity, or simply extrapolating linearly does not produce a plateau and results in predictions highly dependent on the extrapolation region (see Fig. 2(c)).

3.2 Averaging over disorder

A complementary perspective on the Drude-based extrapolation can be gained through a general connection of conductivity averaged over disorder to the pseudo-conductivity $\text{Re}[\sigma(\omega)]_{\eta}$.

In a mathematical sense, conductivity of a single finite system according to the Kubo-Greenwood Eq. (1) is not a function, but a functional. Conductivity can be defined as a proper function if averaged over disorder. Indeed, suppose two energy states $m$ and $n$ have a probability distribution $P(\mathbf{\Gamma},\mathbf{\Delta})$, where $\mathbf{\Gamma}\equiv\Gamma_{mn}$ and $\mathbf{\Delta}\equiv\Delta_{mn}$ from Eq. (1). Then we can define averaged conductivity as

[TABLE]

Using Eq. (1) and Performing integration over $\mathbf{\Delta}$ results in

[TABLE]

So, if $P(\mathbf{\Gamma},\omega)$ is a smooth function of $\omega$, then $\langle\text{Re}[\sigma(\omega)]\rangle$ is a smooth function as well.

Suppose our goal is to find $\langle\sigma\rangle(\omega)$ from sample systems. One of the ways to do so is to find $\langle\sigma_{\eta}\rangle(\omega)\equiv\langle\text{Re}[\sigma_{\eta}(\omega)]\rangle$ first, and then extrapolate it to the limit $\eta\to 0$. Such extrapolation may benefit from the fact that $\langle\sigma_{\eta}\rangle(\omega)$ is a convolution $\langle\sigma\rangle(\omega)$ with $\delta_{\eta}(\omega)$:

[TABLE]

This can be easily seen from Eq. (1), and relation

[TABLE]

Similarly, averaged conductivity at finite temperature as a function of chemical potential $\mu$ can be expressed as convolution of averaged conductivity at zero temperature with the derivative of Fermi distribution with respect to $\beta$:

[TABLE]

Here $f_{\beta}^{\prime}(\mu)\equiv\beta\cosh(\beta\mu/2)^{-2}$ is the derivative of the Fermi distribution over $\beta$. The relation can be obtained by manipulating equation (8). Joining Eq. (17) and (19) results in the expression for conductivity with finite broadening and temperature expressed through convolution of averaged conductivity at zero temperature and broadening with the function $f_{\beta}^{\prime}(\mu)\delta_{\eta}(\omega)$:

[TABLE]

Notice, that Eq. (20) is valid for any system.

Convolution in the real space corresponds to multiplication in the Fourier space, and therefore convolution suppresses high-frequency harmonics. When we assume a model for $\langle\sigma_{\infty}\rangle(\mu,\omega)$ and then fit the parameters of the model, we effectively fit the low-frequency region in the Fourier space of the model to the low-frequency region of the numerical data. The fitting region has an effective size $\eta^{-1}\times\beta^{-1}$ with the error between the model and the data weighted proportionally to the Fourier transform of $f_{\beta}^{\prime}(\mu)\delta_{\eta}(\omega)$.

4 Method comparison

We compare our method to two other possible approaches:

The first approach follows the logic described in the introduction: we first find $\text{Re}[\sigma(\omega)]_{\eta}$ for several values of $\omega$ at $\eta=W/N$, where $W$ is the band width, and $N$ is the number of energy levels. We then fit the results with a Drude equation (12) to obtain DC conductivity $\sigma_{\text{DC}}=A$. We already know that conductivity obtained by this method is sensitive to the representation of the delta function and considered $\eta$, but we wanted to explore how much the method underestimates conductivity for different disorders. 2. 2.

The second approach is based on the Landauer equation [12]. Here one finds DC conductance of the system placed between infinite electric leads. Empirically, the resistance of the system can be thought of as a sum of resistances of the bulk and contacts:

[TABLE]

If the size of the system is much larger than the correlation length, contact resistance can be assumed constant, while bulk resistance will increase linearly with the length of the system. The conductivity of the system can then be defined as

[TABLE]

We made a comparison of conductivity as a function of chemical potential at three different values of the onsite disorder obtained by different methods. The results are presented in Fig. 3. The values of conductivity obtained by our method and the Landauer-based approach are within the error-bars, while the classical approach underestimates conductivity by a factor of 5 for the weak disorder, and becomes closer to other methods for high disorders.

5 Applications

We applied our method to two test systems. The first is a 3D system with a low concentration of impurities, and the second is a 1D disordered chain. The essential features of the models are the appearance of new length scales: the typical distance between impurities and localization length.

5.1 3D system with impurities

The model of the uncorrelated weak onsite disorder, which we used in the previous examples, is not very realistic for real materials. Another common approach is to model disorder by a small concentration of strong impurities in an otherwise perfect system. We were interested whether the addition of the new length scale would influence the convergence rate of the method. Intuitively it seems impossible that the conductivity of a system smaller than an average distance between impurities can be close to the thermodynamic limit. On the other hand, the intuition about the convergence of $\text{Re}[\sigma(\omega)]_{\eta}$ described in Appendix B did not rely on any specific properties of the disorder distribution.

We considered a tight binding 3D system of the size $L\times L\times L$ and nearest neighbor hopping similar to the system described in Fig. 1, but with the difference that onsite potential could be either $V(1-p)$ with a small probability $p$ or $-Vp$ with probability $1-p$. With such choice the average potential is zero, while the difference between ”impurities” and background potential is $V$. We considered two examples with $p=0.003$ and $p=0.001$. Calculated conductivity as a function of chemical potential for different sizes of the system is presented on Fig. 4. Comparing the results to Fig. 2(b) it is evident that convergence rate is not influenced by the introduction of a new length scale.

5.2 1D system with disorder

Often an intuition about a 3D or 2D system can be obtained by looking at a 1D system with a similar structure. Unfortunately, 1D systems with disorder have zero DC conductivity in thermodynamic limit due to the localization of the wave-functions [15]. On the other hand, short chains have an empirically well defined conductivity, since their resistance increases linearly with the size of the system. This result is demonstrated on Fig. 5(a)), where we plotted Landauer resistance vs length of a 1D tight binding chain with nearest neighbor hopping $t=-1$, onsite disorder $\Delta\varepsilon=0.1$, chemical potential $\mu=-1$ and temperature corresponding to $\beta=30$.

We were interested to what extend the Drude-based extrapolation would reproduce these results. It appears conductivity obtained by our method predicts resistance of short chains very well (see blue dots on Fig. 5(a)), but underestimates conductivity for long chains. This is expected, since long chains are insulators and conductivity cannot follow the Drude equation for long excitation times. This becomes also evident from the plot of conductivity vs fitting region (Fig. 5(b)): while in 3D systems conductivity was almost independent of the fitting region, in 1D conductivity predicted from long and short excitation times is different.

The important result is that we can predict static response of systems connected to the leads from dynamic response of the system with periodic boundary conditions. This is not the case with the classical approach, which noticeably underestimates the DC response.

6 Summary and outlook

In the current paper, we proposed a procedure for estimating the thermodynamic limit of DC conductivity from a dynamic current response of a microscopic quantum system. The following steps roughly describe the procedure:

First, conductivity $\sigma_{\mathrm{DC}\eta}$ of the microscopic system is found using Kubo-Greenwood equation (7) with smeared delta-function $\delta(\omega)\rightarrow\delta_{\eta}(\omega)$ for several values of the broadening parameter $\eta$. The conductivity $\sigma_{\mathrm{DC}\eta}$ is then averaged over many realizations of the disorder. 2. 2.

A particular model for conductivity $\sigma(\omega)$ is adopted (Drude model in our case), and parameters of the model (and thus DC conductivity) are found by fitting $\sigma_{\mathrm{DC}\eta}$ calculated by Eq. (17) to numerical results.

Consistency of the procedure can be verified by considering the convergence of the obtained DC conductivity as the size of the system increases similarly as in the classical approach. Importantly, the convergence is accelerated greatly if the broadening parameter $\eta$ and the temperature are above a certain threshold. Physical considerations (see B) suggest that the critical value of $\eta$ is such that the Fourier transform of $\delta_{\eta}(\omega)$ – $\delta_{\eta}(t)$, has a width of $L/v_{\mathrm{max}}$, where $L$ is the size of the system and $v_{\mathrm{max}}$ - maximum group velocity in the corresponding direction. The critical temperature is $T_{\mathrm{crit}}=2D/L^{2}$, where $D$ is the effective diffusion coefficient of the Schrodinger equation in the imaginary time.

Another major factor influencing the convergence rate is the chosen model for conductivity and the type of the broadening function. In the considered example of a 3D tight binding model, the best performers were the Drude model for conductivity and the sinc broadening. In real systems, the situation may change due to the presence of optical resonances. Gaussian broadening, because of being more localized, may perform better, making the final result less sensitive to a wrong choice of conductivity model. The question, though, requires further investigation.

So far, we considered Drude-based extrapolation in the context of the thermodynamic limit. Another way of looking at the extrapolation is through a lens of data compression: Drude model predicts an electric current response of a finite system to a family of slowly varying electric fields acting over the times $t<t_{\text{crit}}=L/2v_{\text{max}}$, thus compressing the information about the responsiveness of the system to two numbers. Notably, the value of conductivity obtained by our method also matches the Landauer-based approach, which implies a static response of a system connected to electrical leads can be predicted from the dynamic response of the system detached from the leads. It is not the case with a classical approach, which predicts a much lower value for DC conductivity than Landauer.

While the Drude-based extrapolation works well in practice, a rigorous understanding is still lacking. The arguments presented in B provide a strong motivation for the convergence conditions of our method to thermodynamic limit, but we do not have mathematical proof. Looking at the plateau in our plots, one may wonder if there is a way to define DC conductivity unambiguously, even for finite systems. Unfortunately, we do not see any physical basis for this claim, since the Kubo-Greenwood equation is exact. Experimental definition of finite-sized conductivity contains similar ambiguities to the theoretical approaches: one needs to exclude contact resistance and interaction of the system with the thermostat. The empirically well-defined plateau only means that $\sigma(\omega)_{\eta}$ is indeed very close to the thermodynamic limit and that our assumptions about the shape of conductivity as a function of frequency are close to the actual shape.

We believe that an important application of our method is its extension to the Kernel Polynomial based Methods (KPM) [14] and other versions of the Kubo-Greenwood formula [16]. The KPM methods allow calculation of conductivity of very large systems at low temperatures by removing the necessity of calculating the whole spectrum and scaling linearly with the system size.

The application of our method to the existing approaches in DFT packages [1] should be straightforward. Calculations of conductivity at different broadenings of the delta function are performed routinely as a consistency check of the classical approach. Adding the extrapolation step on the top is a computationally efficient way to greatly improve convergence to the thermodynamic limit.

7 Acknowledgements

The author thanks Jorge Sofo for proposing this problem and numerous helpful discussions. We thank Lucas Hackl for the help with structuring and preparing this manuscript. We thank Barbara Kennedy for the help with the language and outlook. We thank Hyokun Yun, William Faugno, Jacob Robins, Brett Green, Likun Zhang, Dmytro Oliinychenko and Xiaoxiao Li for productive discussions, which helped to clarify various aspects of the work.