Theory

Linear scaling method

The linear scaling method we use is based on the expansion of the quantities of interest we want to study in terms of Chebyshev polynomials. These polynomials have been routinely employed in chemistry and physics and exhibit greater accuracy compared to other sets of orthogonal polynomials. In what follows, we give a brief introduction to the topic. For more details, see this reference.

The expansion of a function \(f: [-1, 1] \rightarrow {\rm R}\) in terms of Chebyshev polynomials of the first kind reads \[ \begin{aligned} f(x) & = \left[ \alpha_0 + 2 \sum\limits_{n=1}^{\infty} \alpha_n T_n(x)\right] \\ \alpha_n & = \int_{-1}^1 dx \frac{f(x) T_n(x)}{\pi \sqrt{1 - x^2}}, \end{aligned} \] where \(T_n(x)\) is the n-th Chebyshev polynomial of the first kind and \(\alpha_n\) is its corresponding moment. The Chebyshev polynomials benefit from the following advantageous recursion relation

\[ T_0(x) = 1, \quad T_1(x) = x, \quad T_{n+1}(x) = 2 x T_n(x) - T_{n-1}(x). \]

Although the above-mentioned expansion is correct and efficient, it can be problematic when evaluating the moments due to the form of the denominator. A considerable simplification can be introduced by working with the modified orthogonal functions and the corresponding dot product:

\[ \begin{aligned} & \psi_n(x) = \frac{T_n(x)}{\pi \sqrt{1 - x^2}} \\ & \langle f\rvert\psi_n\rangle \equiv \int_{-1}^1 dx \sqrt{1 - x^2} f(x) \psi_n(x) = \int_{-1}^1 dx f(x) T_n(x). \end{aligned} \] Rewriting the expansion in this modified framework yields

\[ \begin{aligned} f(x) & = \frac{1}{\pi \sqrt{1 - x^2}} \left[ \mu_0 + 2 \sum\limits_{n=1}^{\infty} \mu_n T_n(x)\right]\\ \mu_n & = \int_{-1}^1 dx f(x) T_n(x), \end{aligned} \] which, as we will see later, allows for an efficient recursive calculation of the moments.

For our purposes, the quantities we aim to represent using the Chebyshev polynomials basis are functions of the Hamiltonian \(\hat{H}\) rather than of a simple independent variable. Before expressing these functions as a Chebyshev expansion, it is important to ensure that the spectrum of the Hamiltonian lies within the interval \([-1, 1]\), as the expansion is only defined there. If needed, the Hamiltonian should be rescaled by

\[ \tilde{H} = \frac{(H - E_C)}{E_W},\]

where \(E_C\) is the center of the original spectrum and \(E_W\) is its full width. For simplicity, in what follows, we omit the ~ sign on the Hamiltonian; however, we will always refer to the scaled Hamiltonian. The results obtained will need to be appropriately scaled back to fit the original spectrum, hence the Hamiltonian ¹. The quantities of interest, which are functions of the (rescaled) Hamiltonian, can be expressed as a Chebyshev series by observing

\[ \begin{aligned} f(\hat{H}) & = \sum\limits_n f(\hat{H}) \lvert n\rangle \langle n\rvert = \sum\limits_n f(E_n) \lvert n\rangle \langle n\rvert = \sum\limits_n \sum\limits_m \mu_m T_m(E_n)\lvert n\rangle \langle n\rvert \\ & = \sum\limits_n \sum\limits_m \mu_m T_m(\hat{H})\lvert n\rangle \langle n\rvert = \sum\limits_m \mu_m T_m(\hat{H}) \end{aligned} \]

where the set \({\lvert n\rangle }\) consists of the single-particle eigenfunctions of the Hamiltonian with eigenvalues \({E_n}\). As an example, let us explicitly write the expansion of the density of states, i.e.

\[ \rho(E) = {\rm Tr}\left[ \delta(E-\hat{H}) \right] = \sum\limits_n \langle n\rvert \delta(E-\hat{H}) \lvert n\rangle . \] The operator to be traced is, as we just pointed out, a function of the Hamiltonian. By inserting the resolution of the identity and using the manipulations shown above, we obtain

\[ \int dE' \delta(E' - \hat{H}) T_n(E') = \sum\limits_i \int dE' \delta(E' - \hat{H}) T_n(E') \lvert i\rangle \langle i\rvert = \sum\limits_i T_n(E_i) \lvert i\rangle \langle i\rvert = T_n(\hat{H}). \] The expansion of the density of states in terms of Chebyshev moments is then

\[ \begin{aligned} \rho(E) &= \sum\limits_n \mu_n T_n(E_n) \\ \mu_n &= \langle n\rvert T_n(\hat{H}) \lvert n\rangle . \end{aligned} \] This expression, though exact, is not particularly convenient, as it requires knowledge of all single-particle eigenfunctions \(\lvert n\rangle \). This would defeat the purpose of any other calculation, since knowledge of these eigenstates and their corresponding eigenvalues defines any property of the system ². Nonetheless, the form presented above is used in our code to evaluate specific matrix elements of functions of the Hamiltonian - e.g., Green’s function matrix elements in the site basis, for which only a single evaluation is necessary. Fortunately, the calculation of the trace can be approximated in an efficient and accurate way through a stochastic evaluation ^{3 4}. In this approach, the trace of a general single-particle operator \(\hat{A}\) is approximated as

\[ {\rm Tr}[ \hat{A} ] \simeq \sum\limits_{r=1}^R \langle r\rvert \hat{A} \lvert r\rangle ,\]

where the set of \({\lvert r\rangle }\) consists of random states constructed from independent and identically distributed random coefficients with zero mean. The vectors are normalized to 1, which ensures that they can indeed be considered states for the purpose of our calculations. This means, however, that when evaluating quantities such as the density of states, the average over the random states will be normalized to 1 rather than to the number of particles. In terms of accuracy, it has been shown ⁵ that the error of the approximation is proportional to \(1/\sqrt{N R}\), where N is the system size and R is the number of random vectors employed in the stochastic evaluation. This works in our favor, as the accuracy increases with the system size - so much so that, for a system with a billion sites, averaging over only a few \((<5)\) states is sufficient for most purposes ⁶: in this case \(1/\sqrt{N R} \lesssim 10^{-9}\).

Kubo formula for the conductivity

For the analysis of linear response, we use a method based on the Kubo formula. In this context, we focus on the case of responses to electrical signals. According to the Kubo formula, the response of an operator \(\hat{A}\) to a zero-frequency (DC) perturbation is givem by

\[ \langle \hat{A} \rangle = \Omega \lim\limits_{\tau, t \rightarrow \infty} \int^t_0 dt' e^{-t'/\tau }\int^\beta_0 d\beta' {\rm Tr}\left[ \hat{\rho}_{\rm eq} \hat{A}(0) \hat{\mathbf{J}}(-t' - i \hbar \beta') \right] \cdot \mathbf{E_0}, \]

where \(\Omega\) is the volume of the system, \( \tau \) is the intrinsic lifetime, \( \hat{\rho}_{\rm eq} \) is the equilibrium density, \( \mathbf{E_0} \) is the DC perturbation, and \( \hat{J} \) is the current operator. If the general operator \( \hat{A}=\hat{\mathbf{J}}_\alpha \), we have

\[ \label{eq_kubo_mb} \langle \hat{J}_\alpha \rangle = E_{0, \alpha}\left\{ \lim\limits_{\tau, t \rightarrow\infty} \int^t_0 dt' e^{-t'/\tau}\int^\beta_0 d\beta' {\rm Tr}\left[ \hat{\rho}_{\rm eq} \hat{J}_\alpha(t' + i \hbar \beta') \hat{\mathbf{J}}_\alpha(0) \right] \right\}, \]

where the volume has been ‘absorbed’ into the current, and we are considering the spatially averaged quantity. According to Ohm’s law, the quantity in the brackets in the equation above represents the DC conductivity \(\sigma\). This formula simplifies considerably in the non-interacting case. Following the steps detailed in ⁷, we can rewrite the above equation as

\[ \langle \hat{A}(\mu) \rangle = \lim\limits_{\eta \rightarrow 0} - i \hbar \Omega \int dE' f(E' - \mu) {\rm Tr}\left[ \delta(E' - \hat{H}) \hat{A} \frac{1}{(\hat{H} - E') (\hat{H} -E' - i\eta)}(\hat{\mathbf{J}} \cdot \mathbf{E_0}) - h.c. \right], \] where the Fermi occupation function \(f(E)\) and the denominator of the fraction come from the density matrix. Additionally, the following propertiy has been used

\[ \begin{aligned} \int dE' F(E')\delta(E' - E_n) = F(E_n)\ \sum\limits_n F(E_n) \lvert n\rangle \lvert n \rangle \langle n\rvert = F(\hat{H}), \end{aligned} \]

where \(\lvert n\rangle \) are the single particle eigenfunction of the Hamiltonian with energy \(E_n\). Rewriting this equation in terms of Green’s functions, one obtains

\[ \langle \hat{A}(\mu) \rangle = i \hbar \Omega \int dE' f(E' - \mu) {\rm Tr}\left[ \delta(E' - \hat{H}) \hat{A} \frac{d G^+(E')}{dE'}(\hat{\mathbf{J}} \cdot \mathbf{E_0}) - h.c. \right], \]

where \(G^\pm(E')\) are the retarded and advanced Green’s functions. This formula is referred to as the Kubo-Bastin formula and is equivalent to the Kubo formula for non-interacting systems.

If we focus on the conductivity, i.e. \(\hat{A} = \hat{J} = q \hat{V}/\Omega\), where q is the carrier charge and \(\hat{V}\) is the velocity operator, then with some simple manipulations, we arrive at

\[ \sigma(\mu, T) = - \frac{2 \hbar e^2}{\Omega} \int dE' f(E' - \mu) {\rm Tr}\left[ \delta(E' - \hat{H}) \hat{V} \frac{d \operatorname{Im}\left[G^+(E')\right]}{dE'} \hat{V}\right], \] where we make use of \(G^+(E') - G^-(E') = -2i\operatorname{Im}\left[G^+(E')\right]\). Integrating the last equation by parts leads to

\[ \sigma(\mu, T) = - \frac{\hbar e^2}{\Omega} \int dE' \left[- \frac{\partial f(E' - \mu)} {\partial E'}\right] {\rm Tr}\left[ \delta(E' - \hat{H}) \hat{V} \operatorname{Im} \left[G^+(E')\right]\hat{V}\right], \] where the derivative of the Fermi function acts as a window, selecting the energies around the chemical potential that contribute to the dissipative conductivity, as expected. The zero-temperature limit of the above equation is the so-called Kubo-Greenwood formula

\[ \sigma(E) = \frac{\pi \hbar e^2}{\Omega} {\rm Tr}\left[ \delta(E' - \hat{H}) \hat{V} \delta(E - \hat{H}) \hat{V}\right], \] where we make use of the identity \(\operatorname{Im}\left[G^+(E)\right] = - \pi \delta(E - \hat{H})\). For simplicity, the direction of the velocity operators has been omitted, but for calculating any element of the conductivity tensor, one simply needs to use the velocity operators in the corresponding directions, i.e.

\[ \sigma_{xy} \propto {\rm Tr}[ \delta(E' - \hat{H}) \hat{V_x} \delta(E - \hat{H}) \hat{V_y}]. \] This formula can then be expanded in terms of Chebyshev polynomials by transforming the delta functions back into Green’s functions and evaluating the corresponding moments for the stochastic samples, as shown in ⁷, with the only difference being the normalization of the random states and, consequently, the volume prefactor ⁸.

Velocity operator

As shown in the previous sections, to compute the conductivity of the system, we need the velocity operator. In this section, we demonstrate how this is handled when the system is initialized via the Lattice Validator input. The velocity operator can be computed using Heisenberg’s equation of motion, \(\hat{{V}} = i/\hbar [\hat{H}, \hat{R}]\), where \(\hat{R}\) is the position operator. While this is, in principle, a simple operation requiring two matrix-matrix multiplications and an addition, it becomes ill-defined for generic boundary conditions other than hard-wall. In fact, the position operator cannot be defined for periodic systems, and the above-mentioned formula fails. This problem, however, can be circumvented by using the 'distance' information provided by the various bonds that make up the Hamiltonian. Specifically, the translation vector associated with each bond will give us information about the real-space displacement between the two unit cells it connects, i.e. all the spatial information we need to compute the commutator \([H, R]\). The periodic boundary conditions are then automatically handled when implementing the recipe given by the translation associated with each bond to construct the full operator, much like the construction of the Hamiltonian. Note that it is not important where spatial origin is set, as it would simply introduce a contribution proportional to the identity in the position operator, which does not contribute to the value of the commutator. As an example, let us consider the case of an unevenly connected chain, as shown in the following figure:

This system has the following Hamiltonian (on the lattice)

\[ H = \begin{pmatrix} 0 & t_0 & 0 & 0\\ t_0 & 0 & t_1 & 0\\ 0 & t_1 & 0 & t_0\\ 0 & 0 & t_0 & 0\\ \end{pmatrix}. \]

Considering a unit cell made of two different sites, this Hamiltonian can equivalently be represented by the following bond information

\[ \begin{aligned} & (0, 0, 0)\, , \, \bar{t}_0 = \begin{pmatrix} 0 & t_0\\ 0 & 0 \end{pmatrix} \\ & (1, 0, 0) \, , \, \bar{t}_1 = \begin{pmatrix} 0 & 0\\ t_1 & 0 \end{pmatrix}, \end{aligned} \] where the vector represents the translation in terms of unit cells associated with the bond (hence \((0, 0, 0)\) is the bond within the unit cell), and \(\bar{t}_{0, 1}\) are matrices spanning the different sites of the unit cell, indicating the connections within the different sites and their corresponding strengths. Since the Hamiltonian is the sum of the contributions from the two bonds, we can sum the velocity contributions associated with each bond. In this system only the x-component of the velocity operator will be non-zero. For the trivial translation associated with the bond \(t_0\), the commutator \([\bar{t}_0, {r}]\) can be directly evaluated, where \(\bar{t}_0\) is the first bond Hamiltonian from the previous equation and the matrix \(r\) is simply the position operator within the unit cell. For the second bond, a similar approach can be used, but care must be taken with the position operator. Let us now consider the subspace made of the original unit cell and the cell reached by the translation \((1, 0, 0)\), i.e. the displacement associated with the bond

\[ \bar{t}_1 = \begin{pmatrix} 0 & 0\\ t_1 & 0 \end{pmatrix}. \] The position operator in the x-direction, in the site basis within this subspace ⁹, reads:

\[ R_x = \begin{pmatrix} p^0_x & 0 & 0 & 0 \\ 0 & p^1_x & 0 & 0 \\ 0 & 0 & p^0_x + d_x & 0 \\ 0 & 0 & 0 & p^1_x + d_x \end{pmatrix} = \begin{pmatrix} \bar{r} & \bar{0} \\ \bar{0} & \bar{r}_d \end{pmatrix}, \] where \(\bar{0}\) is the \(2 \times 2\) null matrix, \(p^i_x\) is the position in the x-direction of the site \(i\) and \(d_x\) is the displacement in the x-direction corresponding to the lattice vector. The matrices \(\bar{r}\) and \(\bar{r}_d\) are the position operators in the original unit cell and the displaced one, respectively. Note that, in general, the displacement vector is

\[ (l_1, l_2, l_3) \cdot (\mathbf{v}_1, \mathbf{v}_2, \mathbf{v}_3), \] where \(l_i\) is the translation associated with the bond for the i-th lattice vector \(\mathbf{v}_i\). Our example is the trivial case of a translation by a single lattice vector of magnitude one. In the same subspace, the Hamiltonian contribution of the (directional) bond is

\[ H^{t_1} = \begin{pmatrix} \bar{0} & \bar{t}_1 \\ \bar{0} & \bar{0} \end{pmatrix}. \] The position of non-zero block is determined by the directionality of the bond definition: if the translation were \((-1, 0, 0)\), i.e. hopping to the left, the non-zero block would be the lower left one. In terms of the block representation, it is straightforward to evaluate the commutator

\[ V^{t_1}_x = \frac{i}{\hbar} [H^{t_1}, R_x] = \frac{i}{\hbar} \left[ \begin{pmatrix} \bar{0} & \bar{t}_1 \cdot \bar{r}_d \\ \bar{0} & \bar{0} \end{pmatrix} - \begin{pmatrix} \bar{0} & \bar{r} \cdot \bar{t}_1\\ \bar{0} & \bar{0} \end{pmatrix}\right], \] where the superscript \({t_1}\) denotes that the contribution to the corresponding quantities comes from the bond \(t_1\). Since the block structure of the resulting contribution to the velocity operator is the same as the bond, one can construct the velocity associated with each bond as \(v^{t_i}_x = i/ \hbar (\bar{t}_i \cdot \bar{r}_{d^i} - \bar{r} \cdot \bar{t}_i)\), and distribute these across the whole system in the same way that the Hamiltonian is constructed. The storage strategy will be the same as well, simply keeping track of the (intra-)inter-cell velocity contributions and the associated translations. In this sense, the boundary conditions do not represent a problem from this perspective, because the correct behavior of the velocity operator is imposed by means of the distribution functions, i.e. how the sub-blocks are repeated in the matrix corresponding to the whole system.

Transport properties

HQS Qolossal offers the possibility of calculating carrier mobility and diffusion coefficient. These are quantities derived from the DC conductivity and the density of states. Carrier mobility $$\sigma =e{n}_{e,h}{\mu }_{e,h}$$ with e electron charge, $n_{e,h}$ the electron (hole) density, $\sigma$ the DC conductivity and, finally, ${\mu }_{e,h}$ the mobility.

Diffusion coefficient $$\sigma (\mu )=e^2\rho (\mu )D(\mu )$$ with $\sigma$ the DC conductivity, $\rho$ the density of states and $D$ the diffusion coefficient. All quantities refer to their values at the specified chemical potential $\mu$.

Optical properties

The optical electrical conductivity is computed via the Chebyshev expansion of the well-known Kubo-Bastin formula¹⁰

$${\sigma }_{\alpha \beta }(\mu ,T)=\frac{i{e}^2\hslash }{\Omega }{\int }_{-\infty }^{\infty }dϵf(ϵ)\mathrm{Tr}\left[{\hat{V}}_{\alpha }\delta (ϵ-\hat{H}){\hat{V}}_{\beta }\frac{d{G}^{+}(ϵ)}{dϵ}-{\hat{V}}_{\alpha }\frac{d{G}^{-}(ϵ)}{dϵ}{\hat{V}}_{\beta }\delta (ϵ-\hat{H})\right]$$ with $f(ϵ)$ the fermi function, $\hat{V}$ the velocity operator and $G^{\pm }(ϵ)$ the $\pm$ components of the Greens function. This is expanded into Chebyshev polynomials the same way as for the DC conductivity. All other optical properties are derived from the optical conductivity. Below all definitions.

Permittivity $${\epsilon }_{\alpha \beta }(\omega )={\delta }_{\alpha \beta }+\frac{i\sigma (\omega )}{{\epsilon }_0\omega }$$

Reflectivity $$\begin{gathered} \frac{(1-n{)}^2+k^2}{(1+n{)}^2+k^2}, \\ N(\omega )=n(\omega )+ik(\omega )=\sqrt{\epsilon (\omega )} \end{gathered}$$ Absorption coefficient $$\alpha =\frac{2\omega k(\omega )}{c}$$

Bibliography

1. This means that any prefactor coming from the rescaling of the Hamiltonian has to be included. For example, the expansion of the density of states will be presented later in the text. The delta function contained in its definition will yield \(\delta(E - H) = E_W^{-1} \delta(\tilde{E} - \tilde{H})\). The factor \(E_W^{-1}\) will have to be included in the final rescaling to obtain the correct result. ↩

2. Aside from representing a computational effort which scales like the cubic power (or square, if the Hamiltonian is sparse) of the system's size, since it is the diagonalization of the whole Hamiltonian. ↩

3. Drabold, David A., and Otto F. Sankey. 1993. “Maximum Entropy Approach for Linear Scaling in the Electronic Structure Problem.” Phys. Rev. Lett. 70 (June): 3631–34. https://doi.org/10.1103/PhysRevLett.70.3631. ↩

4. Silver, R. N., and H. Röder. 1994. “DENSITIES OF STATES OF MEGA-DIMENSIONAL HAMILTONIAN MATRICES.” International Journal of Modern Physics C 05 (04): 735–53. https://doi.org/10.1142/S0129183194000842. ↩

5. Weiße, Alexander, Gerhard Wellein, Andreas Alvermann, and Holger Fehske. 2006. “The Kernel Polynomial Method.” Rev. Mod. Phys. 78 (March): 275–306. https://doi.org/10.1103/RevModPhys.78.275. ↩

6. The appropriate number of random states to average over also depends on the order of the expansion: the higher the number of moments, the lower the frequency resolution of the expansion, hence the larger the sampling pool must become. ↩

7. Fan, Zheyong, José H. Garcia, Aron W. Cummings, Jose Eduardo Barrios-Vargas, Michel Panhans, Ari Harju, Frank Ortmann, and Stephan Roche. 2021. “Linear Scaling Quantum Transport Methodologies.” Physics Reports 903: 1–69. https://doi.org/https://doi.org/10.1016/j.physrep.2020.12.001. ↩ ↩2

8. The corresponding volume will be the volume of the unit cell instead of the total volume. ↩

9. Note that, in this case, the subspace mentioned is equivalent to the whole system, but in any case with more than two unit cells, this will not be the case. ↩

10. Jose H. Garcı́a, Lucian Covaci, and Tatiana G. Rappoport. 2015. "Real-space calculation of the conductivity tensor for disordered topological matter". Phys. Rev. Lett., 114, p. 116602. https://doi.org/10.1103/PhysRevLett.114.116602. ↩