A general method for calculating lattice Green functions on the branch cut

TL;DR

This paper introduces a versatile Chebyshev expansion method for accurately computing lattice Green functions at any real frequency, applicable to various lattice types, and improves convergence by modeling van Hove singularities.

Contribution

The authors develop a general Chebyshev expansion approach for calculating lattice Green functions, including a technique to accelerate convergence by subtracting modeled singularities.

Findings

Achieves 6-9 significant figure accuracy with 1000 series terms.

Applicable to multiple lattice types including square, cubic, and diamond.

Method effectively captures van Hove singularities in Green functions.

Abstract

We present a method for calculating the complex Green function $G_{ij} (\omega)$ at any real frequency $\omega$ between any two sites $i$ and $j$ on a lattice. Starting from numbers of walks on square, cubic, honeycomb, triangular, bcc, fcc, and diamond lattices, we derive Chebyshev expansion coefficients for $G_{ij} (\omega)$. The convergence of the Chebyshev series can be accelerated by constructing functions $f(\omega)$ that mimic the van Hove singularities in $G_{ij} (\omega)$ and subtracting their Chebyshev coefficients from the original coefficients. We demonstrate this explicitly for the square lattice and bcc lattice. Our algorithm achieves typical accuracies of 6--9 significant figures using 1000 series terms.

Tables

Table 1. Table 3: Numbers of walks of length n 𝑛 n that return to the origin on various lattices. These correspond to existing entries in the Online Encyclopedia of Integer Sequences (OEIS).

$n$	$W_n^{\text{chain}}$	$W_n^{\text{sq}}$	$W_n^{\text{bcc}}$	$W_n^{\text{hon}}$	$W_n^{\text{diam}}$	$W_n^{\text{cubic}}$	$W_n^{\text{hcub}}$	$W_n^{\text{tri}}$	$W_n^{\text{fcc}}$
0	1	1	1	1	1	1	1	1	1
1	0	0	0	0	0	0	0	0	0
2	2	4	8	3	4	6	8	6	12
3	0	0	0	0	0	0	0	12	48
4	6	36	216	15	28	90	168	90	540
5	0	0	0	0	0	0	0	360	4320
6	20	400	8000	93	256	1860	5120	2040	42240
7	0	0	0	0	0	0	0	10080	403200
8	70	4900	343000	639	2716	44730	190120	54810	4038300
9	0	0	0	0	0	0	0	290640	40958400
10	252	63504	16003008	4653	31504	1172556	7939008	1588356	423550512
OEIS#	A000984	A002894	A002897	A002893	A002895	A002896	A039699	A002898	A002899

Table 2. Table 4: Values of local Chebyshev moments g n = ⟨ 𝟎 | T n ​ ( H ^ ) | 𝟎 ⟩ = ∑ k = 0 n a n ​ k ​ z − k ​ W k subscript 𝑔 𝑛 quantum-operator-product 0 subscript 𝑇 𝑛 ^ 𝐻 0 superscript subscript 𝑘 0 𝑛 subscript 𝑎 𝑛 𝑘 superscript 𝑧 𝑘 subscript 𝑊 𝑘 g_{n}=\left<\mathbf{0}\right|T_{n}(\hat{H})\left|\mathbf{0}\right>=\sum_{k=0}^{n}a_{nk}z^{-k}W_{k} on various lattices. We quote the integer-valued quantities z n ​ g n superscript 𝑧 𝑛 subscript 𝑔 𝑛 z^{n}g_{n} , where z 𝑧 z is the coordination number of each lattice. The last eight series have been submitted to the Online Encyclopedia of Integer Sequences as OEIS A288454, A288455, A288456, A288457, A288458, A288459, A288460, and A288461.

Table 3. Table 5: Coefficients of powers in Chebyshev polynomials, a n ​ k subscript 𝑎 𝑛 𝑘 a_{nk} , such that T n ​ ( x ) = ∑ k = 0 n a n ​ k ​ x k subscript 𝑇 𝑛 𝑥 superscript subscript 𝑘 0 𝑛 subscript 𝑎 𝑛 𝑘 superscript 𝑥 𝑘 T_{n}(x)=\sum_{k=0}^{n}a_{nk}x^{k} . Blank elements are zero.

Table 4. Table 6: Chebyshev transform pairs, i.e., functions f ​ ( ω ) 𝑓 𝜔 f(\omega) for ω ∈ [ − 1 , 1 ] 𝜔 1 1 \omega\in[-1,1] and Chebyshev coefficients f n subscript 𝑓 𝑛 f_{n} for n = 0 , 1 , 2 , … 𝑛 0 1 2 … n=0,1,2,\dotsc such that f ​ ( ω ) = 1 π ​ 1 − ω 2 ​ ∑ n = 0 ∞ ( 2 − δ n ) ​ T n ​ ( ω ) ​ f n 𝑓 𝜔 1 𝜋 1 superscript 𝜔 2 superscript subscript 𝑛 0 2 subscript 𝛿 𝑛 subscript 𝑇 𝑛 𝜔 subscript 𝑓 𝑛 f(\omega)=\frac{1}{\pi\sqrt{1-\omega^{2}}}\sum_{n=0}^{\infty}(2-\delta_{n})T_{n}(\omega)f_{n} and f n = ∫ − 1 1 𝑑 ω ​ T n ​ ( ω ) ​ f ​ ( ω ) subscript 𝑓 𝑛 superscript subscript 1 1 differential-d 𝜔 subscript 𝑇 𝑛 𝜔 𝑓 𝜔 f_{n}=\int_{-1}^{1}d\omega~{}T_{n}(\omega)f(\omega) . The harmonic numbers are defined as H n = ∑ k = 1 n 1 / k subscript 𝐻 𝑛 superscript subscript 𝑘 1 𝑛 1 𝑘 H_{n}=\sum_{k=1}^{n}1/k . Generally, stronger divergences in f ​ ( ω ) 𝑓 𝜔 f(\omega) correspond to more slowly decaying tails in f n subscript 𝑓 𝑛 f_{n} .

Equations

$$\hat{H}$$

$$G_{ij}(\omega )$$

$$G(\omega ,\mathbf{r})$$

$$W_{xyzn}^{\text{bcc}}$$

$$\frac{n!}{n_x!m_x!n_y!m_y!n_z!m_z!}.$$

$$W_{xyzn}^{\text{cubic}}$$

$$W_{xyzn}^{\text{hon}}$$

$$W_{xyzn}^{\text{tri}}$$

$$=⟨\mathbf{r}|\sum _{j=0}^n\left(\genfrac{}{}{0ex}{}{n}{j}\right)(3{\hat{H}}^{\text{hon}}{)}^{2j}(-3{)}^{n-j}|0⟩$$

$$=\sum _{j=0}^n\left(\genfrac{}{}{0ex}{}{n}{j}\right)(-3{)}^{n-j}{W}_{x,y,z,2j}^{\text{hon}}.$$

$$W_{uvwsn}^{\text{fcc}}$$

$$x$$

$$y$$

$$z$$

$$G(\omega )=\sum _{n=0}^{\infty }\frac{⟨i|{\hat{H}}^n|j⟩}{{\omega }^{n+1}}.$$

$$T_n(\cos \theta )$$

$$U_{n-1}(\cos \theta )$$

$${\int }_{-1}^{1}dx{\varphi }_n(x){\varphi }_{n^{\prime }}(x)$$

$$\sum _{n=0}^{\infty }{\varphi }_n(x){\varphi }_n(x^{\prime })$$

$$T_n(x)$$

$$a_{00}$$

$$a_{nk}$$

$${\int }_{-1}^{1}d\nu \frac{T_n(\nu )}{(\nu -\omega )\sqrt{1-{\nu }^2}}=\pi U_{n-1}(\omega ).$$

$$g_n=⟨i|T_n(\hat{H})|j⟩.$$

$$g_n=\sum _{k=0}^n{a}_{nk}{z}^{-k}{W}_k$$

$$g(\omega )=⟨i|\delta (\omega -\hat{H})|j⟩.$$

$$g(\omega )=\frac{1}{\pi \sqrt{1-{\omega }^2}}\sum _{n=0}^{\infty }(2-{\delta }_n)T_n(\omega )g_n.$$

$$G(\omega )$$

$$z^n{g}_n=\sum _{k=0}^n{a}_{nk}{z}^{n-k}{W}_k$$

$$G(\omega )$$

$$f_n^{\text{win}}$$

$$\displaystyle g(\omega)=\int_{\text{BZ}}\frac{d^{d}k~{}}{(2\pi)^{d}}\delta\big{(}\omega-\varepsilon(\mathbf{k})\big{)}$$

$$g(\omega )=\int \frac{d^{d-1}{k}_{\perp }}{(2\pi {)}^d}\frac{1}{\left|\nabla \epsilon (\mathbf{k})\right|}.$$

Full text

Consider the quantum mechanical tight-binding Hamiltonian

[TABLE]

where $t$ is the hopping amplitude between nearest-neighbor sites $i$ and $j$ on a lattice. Define the “Greenian” operator $\hat{G}(\omega)=(\omega-\hat{H})^{-1}$ where $\omega-\hat{H}\equiv\omega\hat{1}-\hat{H}$ and $\hat{1}$ is the identity operator. The matrix elements of the Greenian

[TABLE]

are called the lattice Green function (LGF). In this paper we consider $i$ and $j$ to be fixed, and we will often omit these indices for brevity.

Lattice Green functions are not limited to quantum mechanics, but arise frequently in many other areas of physics. The ability to compute LGFs can be useful, for example, for simulations of Hubbard modelsbloch2008review and for non-perturbative renormalization group studies of scalar boson models.caillol2012 ; caillol2013 For $i=j$, $G_{ij}(\omega)$ can be expressed as closed forms in terms of named special functions (mainly elliptic integrals or generalized hypergeometric functions) for square, bcc, honeycomb, diamond, cubic, hypercubic, triangular, and fcc lattices.guttmann2010 ; maradudin1960 ; morita1971 ; katsura1971 ; joyce1972 ; joyce1994 ; joyce2002 ; joyce2001 ; joyce2003 ; delves2001 ; ray2014arxiv ; schwalm1988 ; schwalm1992 For $i\neq j$, $G_{ij}(\omega)$ can be expressed as closed forms for square, bcc,ray2014arxiv honeycomb,ray2014arxiv triangular,ray2014arxiv kagome,ray2014arxiv diced,ray2014arxiv and cubicjoyce2002 lattices; spatial recurrence relations exist but are often numerically unstable.morita1971recurrence ; morita1975 ; berciu2009 In this paper we develop a general numerical method applicable to LGFs for which no closed form is known.

In the complex $\omega$ plane, the Green function has a branch cut running from $\omega_{\text{min}}$ to $\omega_{\text{max}}$ along the real axis, where $\omega_{\text{min}}$ and $\omega_{\text{max}}$ are the lowest and highest eigenvalues of $\hat{H}$. Wherever $\omega$ occurs as an argument of a Green function, it is to be interpreted as including an infinitesimal imaginary shift, i.e., as $\omega+i0^{+}$. We will always scale $\hat{H}$ such that all its eigenvalues lie within the interval $[-1,1]$ (see Fig. 1). For the lattices treated in this paper, we ensure this by choosing $t=1/z$ where $z$ is the coordination number (the number of neighbors of each site). 111For a physical tight-binding model, the hopping amplitudes are negative, so it would be more appropriate to choose $t=-1/z$.

It is well known that the lattice Green function $G_{ij}(\omega)$ can be written as an inverse power series about $\omega=\infty$, where the coefficient of the $\omega^{-n-1}$ term is related to the number of paths of length $n$ from site $i$ to site $j$. maassarani2000 ; mamedov2008 This series is useful for numerically evaluating the Green function outside the unit disk ($\left\lvert{\omega}\right\rvert>1$), and sometimes to evaluate the quantities $G(\pm 1)$, which are known as Watson integrals. (See Ref. guttmann2010, for a review.) However, for certain applications, one requires values of the Green function “on the cut” (i.e., on either side of the branch cut), where the power series diverges. This is a more difficult problem.

The most direct approach is to write the Green function as a $d$-dimensional integral over the Brillouin zone, typically of the form

[TABLE]

This expression does not lend itself well to numerical implementation because it involves high-dimensional integration of singular integrands.berciu2010 In previous papers loh2011hlgf ; loh2013hlgf we presented an approach for calculating Green functions accurately ($>12$ s.f.) based on time-frequency Fourier transformation and contour integration techniques. Unfortunately, that approach relies on a factorization property of hypercubic lattice Green functions in the time domain, and we have not been able to generalize it to non-hypercubic lattices.

Some of the literature focuses on finding ordinary differential equations that are satisfied by each lattice Green function.guttmann2010 This does not appear to help directly with evaluating the LGFs on the cut.

The recursion method and the continued-fraction methodberciu2010 ; moller2012 can be used to evaluate LGFs on the cut, although accuracy appears to be limited to 2–6 decimal places.

In this paper we evaluate LGFs by converting their power series into Chebyshev series. This amounts to analytic continuation from the region $\left\lvert{\omega}\right\rvert>1$ to the region $\left\lvert{\omega}\right\rvert\leq 1$, as depicted in Fig. 1.

I Counting walks on lattices

Tables 1 and 3 show combinatorial formulas and explicit values for the numbers of polygons (closed walks) on various lattices. Most of these results are well known in the literature.guttmann2010 ; domb1960 ; joyce1972 ; joyce1994 ; joyce2001 ; joyce2002 ; joyce2003 ; bailey2008 We reproduce them here for the reader’s convenience, organized to reveal that the formulas fall into various families: the bcc family (1D chain, 2D square, 3D body-centered cubic); honeycomb family (2D honeycomb, 3D diamond); cubic family (1D chain, 2D square, 3D cubic, 4D hypercubic); and triangular family (2D triangular, 3D face-centered cubic).

Table 1 also gives formulas for the numbers of open walks beginning at the origin and ending at position $\mathbf{r}=(x,y,z)$ on the lattice. We outline the derivations below.

BCC family: Consider a walk on a bcc lattice. Let $n_{x}$, $m_{x}$, $n_{y}$, $m_{y}$, $n_{z}$, and $m_{z}$ be the numbers of steps in the $\pm x$, $\pm y,$ and $\pm z$ directions. Let $n$ be the total number of steps. During each step, the walker moves simultaneously in the $\pm x$, $\pm y,$ and $\pm z$ directions by either $+1$ or $-1$ units. Thus $n=n_{x}+m_{x}=n_{y}+m_{y}=n_{z}+m_{z}$. Let the net displacements in each direction be $x=n_{x}-m_{x}$, $y=n_{y}-m_{y}$, and $z=n_{z}-m_{z}$. Then $n_{x}=\frac{n+x}{2}$, $n_{y}=\frac{n+y}{2}$, and $n_{z}=\frac{n+z}{2}$. Thus the number of walks of length $n$ with total displacement $(x,y,z)$ is

[TABLE]

This derivation can easily be generalized to a $d$-dimensional bcc lattice.

Cubic family: Consider a walk on a cubic lattice. Suppose the numbers of steps in the $\pm x$, $\pm y,$ and $\pm z$ directions are given by the six integers $n_{x}$, $m_{x}$, $n_{y}$, $m_{y}$, $n_{z}$, and $m_{z}$. The number of such walks is given by the multinomial coefficient

[TABLE]

The total number of steps is $n=n_{x}+m_{x}+n_{y}+m_{y}+n_{z}+m_{z}$ and the net displacements in each direction are $x=n_{x}-m_{x}$, $y=n_{y}-m_{y}$, and $z=n_{z}-m_{z}$. Let $j=m_{x}$, $k=m_{y}$, and $l=m_{z}$. Then $n_{x}=j+x$, $n_{y}=k+y$, and $n_{z}=l+z$. Furthermore, $n=2(j+k+l)+x+y+z$. Thus the total number of walks of length $n$ with total displacement $(x,y,z)$ is

[TABLE]

where $s\equiv\frac{n-x-y-z}{2}$ and $l\equiv s-j-k$. This derivation is easily generalized to a $d$-dimensional hypercubic lattice.

Honeycomb family: The honeycomb lattice can be viewed as the projection of a puckered subset of a cubic lattice, as shown in Fig. 2. Consider a walk starting at the origin. Suppose the numbers of steps in the $\pm x$, $\pm y,$ and $\pm z$ directions are given by the six integers $n_{x}$, $m_{x}$, $n_{y}$, $m_{y}$, $n_{z}$, and $m_{z}$. The total number of steps is $n=n_{x}+m_{x}+n_{y}+m_{y}+n_{z}+m_{z}$ and the net displacements in each direction are $x=n_{x}-m_{x}$, $y=n_{y}-m_{y}$, and $z=n_{z}-m_{z}$. Let $j=m_{x}$, $k=m_{y}$, and $l=m_{z}$. Then $n_{x}=j+x$, $n_{y}=k+y$, and $n_{z}=l+z$.

On the odd-numbered steps of the walk, the walker can only travel in the $+x$, $+y$, or $+z$ directions. Likewise, on even-numbered steps, the walker can only travel in the $-x$, $-y$, or $-z$ directions.

If $n$ is even, then there must be exactly $n/2$ steps along “positive” directions, and $n/2$ steps along “negative” directions. So $n/2=n_{x}+n_{y}+n_{z}=m_{x}+m_{y}+m_{z}$. Thus the number of walks of length $n$ is the number of permutations of step displacement vectors along odd-numbered steps, times the number of ways of permutations for even-numbered steps:

[TABLE]

where $l\equiv n/2-j-k$, and it is assumed that $x+y+z=0$.

If $n$ is odd, then there are $(n+1)/2$ positive steps and $(n-1)/2$ negative steps, so

[TABLE]

where $l\equiv(n-1)/2-j-k$.

Walks on a diamond lattice can be counted in a similar fashion.

Triangular family: If one starts at the origin of the honeycomb lattice and performs two successive hops, there are three paths that return to the origin, and one path to each of the 6 A sites surrounding the origin. Thus $6\hat{H}^{\text{tri}}=(3\hat{H}^{\text{hon}})^{2}-3(\hat{1})$. Therefore the number of paths on a triangular lattice of length $n$ with displacement $(x,y,z)$ (where $x+y+z=0$), using the same coordinate scheme, is

[TABLE]

For closed walks $(x=y=z=0)$ this reduces to the formula for $W^{\text{tri}}_{n}$ shown in Table 1, derived in Ref. guttmann2010, using a different approach.

Similarly, using $12\hat{H}^{\text{fcc}}=(4\hat{H}^{\text{diam}})^{2}-4(\hat{1})$, one obtains

[TABLE]

in terms of the $(u,v,w,s)$ coordinates for diamond and fcc lattices embedded in a 4D grid. The “physical” 3D Cartesian coordinates $(x,y,z)$ are given by the projection

[TABLE]

In the above discussion we have derived the $W_{n}$ using combinatorics. There are other methods to obtain $W_{n}$, such as contour integration techniques.ray2014arxiv

II Basic approach for calculating Green functions

Power moments: Starting from the definition of $G(\omega)$, Eq. (2), and expanding in powers of $\hat{H}$ shows that the Green function can be written in an inverse power series

[TABLE]

From Eq. (1) it is easy to see that $\left<i\right|\hat{H}^{n}\left|j\right>=t^{n}W_{n}=z^{-n}W_{n}$, where $W_{n}$ is the number of walks of length $n$ that begin at site $i$ and end at site $j$. As discussed earlier, Eq. (14) converges only for $\left\lvert{\omega}\right\rvert>1$.

Chebyshev polynomials: The Chebyshev polynomials of the first and second kinds are defined as

[TABLE]

They satisfy orthogonality and completeness relations

[TABLE]

where $\phi_{n}(x)=\sqrt{\frac{2-\delta_{n}}{\pi\sqrt{1-x^{2}}}}T_{n}(x)$. Here $\delta_{nn^{\prime}}$ is the Kronecker delta function and $\delta_{n}\equiv\delta_{n,0}$. The Chebyshev polynomials can be written in terms of monomials as

[TABLE]

where $n=0,1,2,\dotsc,\infty$; $k=0,1,2,\dotsc,n$; and $n-k$ is even. The first few coefficients $a_{nk}$ are shown in Table 5. In Mathematica we have found it fastest to to evaluate $a_{nk}$ using the recursion $a_{nk}=2a_{n-1,k-1}-a_{n-2,k}$ for $n\geq 2$.

The functions $T_{n}(x)/\sqrt{1-x^{2}}$ and $\pi U_{n-1}(x)$ are Hilbert transforms of each other; that is, they obey Kramers-Kronig relations:

[TABLE]

Chebyshev moments: Define the Chebyshev moment $g_{n}$ to be the matrix element of the Chebyshev polynomial of the Hamiltonian operator,

[TABLE]

Using Eq. (19) we obtain

[TABLE]

where the coefficients $a_{nk}$ are given by Eq. (20). Tables 2 and 4 give formulas and values for Chebyshev moments on various lattices.

Spectral function: Define the spectral function $g(\omega)$ as the matrix elements of a Dirac delta function,

[TABLE]

This is a generalization of the local density of states. Expanding the Dirac delta function using the Chebyshev polynomial completeness relation, Eq. (18), we see that the spectral function is the sum of an infinite series of Chebyshev polynomials weighted by the Chebyshev moments,

[TABLE]

Green function: The spectral function and the Green function are related by $g(\omega)=-\frac{1}{\pi}\mathop{\mathrm{Im}}G(\omega+i0^{+})$ and $G(\omega)=\int_{-1}^{1}d\nu~{}\tfrac{1}{\omega-\nu+i0^{+}}g(\nu)$. Thus

[TABLE]

From Eq. (25) and Eq. (21) we then obtain

[TABLE]

Basic algorithm: Our approach for numerical computation of Green function may now be outlined as follows:

Calculate $\{W_{n}\}$ using combinatorial formulas. 2. 2.

Calculate $\{g_{n}\}$ using Eq. (23). 3. 3.

Calculate $G(\omega)$ using Eq. (27).

Illustration: Figure 3 shows the bcc lattice Green function $G_{\mathbf{i}\mathbf{j}}(\omega)\equiv G_{xyz}(\omega)$ for various displacements ${\mathbf{j}-\mathbf{i}}=(x,y,z)$, computed using this approach.

Practical considerations: The sum in Eq. (23) involves large cancellations between terms, so it must be calculated using exact arithmetic. Rewriting the sum as

[TABLE]

allows us to perform the computation exactly using integer arithmetic.

Computing $g_{n}$ for $n=0,1,2,3,\dotsc,N$ directly from values of $W_{k}$ would require at least $O(N^{2})$ time. By using expressions for $W_{k}$ we can derive expressions for $g_{n}$ that may be quicker to evaluate. See Table 2.

The sum has an infinite number of terms. If we truncate the series after $N$ terms, we necessarily introduce some error. The Chebyshev series on $\omega\in[-1,1]$ is equivalent to a Fourier series on $\cos\omega\in[0,\pi]$, and so truncation error takes the form of Gibbs oscillations. These can be mitigated by multiplying $g_{n}$ by an appropriate window function such as a Kaiser-Bessel windowloh2014 ; karki2016 before inserting it into Eq. (27):

[TABLE]

where $\beta$ is a tuning parameter. Nevertheless, accuracy is still limited.

III Van Hove singularities

As illustrated in Fig. 3, lattice Green functions contain van Hove singularities, which are difficult to expand in terms of Chebyshev polynomials. They cause the Chebyshev coefficients $g_{n}$ to decay slowly as some inverse power of $n$. We can accelerate the convergence of the series by identifying the functional forms of the singularities and subtracting the corresponding tails from the sequence of Chebyshev coefficients.

First let us review the theory of van Hove singularities. The density of states on a lattice is

[TABLE]

where $\varepsilon(\mathbf{k})$ is the dispersion relation and $\delta$ is the $d$-dimensional Dirac delta function. Roughly speaking, this may be written as

[TABLE]

When the wavevector $\mathbf{k}$ approaches a critical point $\mathbf{k}_{0}$ where $\nabla\varepsilon=\mathbf{0}$, the integrand diverges, producing a van Hove singularity in $g(\omega)$ at $\omega=\varepsilon_{0}=\varepsilon(\mathbf{k}_{0})$. The dispersion relation may be expanded as a power series about the critical point. Typically

[TABLE]

where $(\lambda_{1},\dotsc,\lambda_{d})$ are the eigenvalues of the Hessian matrix $\nabla\nabla\varepsilon$ and $(q_{1},\dotsc,q_{d})$ are coordinates along orthogonal eigenvectors of the Hessian. Suppose there are $m$ positive and $n$ negative eigenvalues. The numbers $m$ and $n$ determine the shape of the constant-energy surfaces of the dispersion relation; for example, for $d=3$, these surfaces may be ellipsoids, one-sheeted hyperboloids, or two-sheeted hyperboloids. Rescaling coordinates to $(p_{1},\dotsc,p_{n},P_{1},\dotsc,P_{m})$ and $\xi=\omega-\varepsilon_{0}$, and using $\left\lvert{\lambda_{1}\dots\lambda_{d}}\right\rvert=\left\lvert{\det\nabla\nabla\varepsilon(\mathbf{k}_{0})}\right\rvert$, one obtains

[TABLE]

To obtain meaningful results as $\xi\rightarrow 0$, the integratiosn over $P_{1},\dotsc,P_{n}$ should be extended to infinity, but the integrals over $p_{1},\dotsc,p_{m}$ must be cut off at a wavenumber $L$ on the scale of the Brillouin zone.

For $\xi>0$ the leading singular terms in $I_{mn}$, excluding constant backgrounds, are

[TABLE]

For $\xi<0$ we may use $I_{mn}(\xi)=I_{nm}(-\xi)$.

Equation (32) assumed that the Hessian, $\nabla\nabla\varepsilon$, is finite at $\mathbf{k}_{0}$. For tight-binding models that exhibit flat bands where the Hessian is zero, Dirac cones where the Hessian is infinite, Weyl cones, or other unusual features in the band structure, the critical points should be treated on a case-by-case basis.

If there are two or more critical points $\mathbf{k}_{1},\mathbf{k}_{2},\dotsc$ contributing to a van Hove singularity at the same energy $\varepsilon_{0}$, the singular forms simply combine additively:

[TABLE]

Singularity subtraction for the DoS: We may exploit knowledge of the van Hove singularities as follows:

Choose an approximate DoS, $f(\omega)$, whose singularities mimic the van Hove singularities in the density of states, $g(\omega)$. 2. 2.

Calculate its Chebyshev coefficients $f_{n}$ analytically. 3. 3.

Subtract the approximation in the Chebyshev domain to obtain the residual Chebyshev coefficients, $h_{n}=g_{n}-f_{n}$. 4. 4.

Construct the residual function $h(\omega)$ from the coefficients $h_{n}$, for certain values of $\omega$. 5. 5.

Construct $g(\omega)=f(\omega)+h(\omega)$ for these $\omega$ values.

Since $g(\omega)$ and $f(\omega)$ have similar singularities, $g_{n}$ and $f_{n}$ will have the same tails at large $n$, and $h_{n}$ will decay faster. Thus the Chebyshev series for $h(\omega)$ will converge faster.

Illustration for square lattice DoS: The square lattice DoS has band-edge step discontinuities and a mid-band logarithmic divergence,

[TABLE]

Construct an approximation using the functions listed in Table 6, and write down its Chebyshev coefficients:

[TABLE]

The exact Chebyshev coefficients, from Table 2, are

[TABLE]

Let $h_{n}=g_{n}-f_{n}$ and calculate

[TABLE]

The functions $f_{n}$, $h_{n}$, $g_{n}$, $f(\omega)$, $h(\omega)$, and $g(\omega)$ are shown in Figs. 4. By dealing with the leading-order van Hove singularities analytically, we have decreased the truncation error of the 1000-term series from $10^{-3}$ to $10^{-9}$.

Illustration for bcc lattice DoS: The van Hove singularities for the bcc lattice density of states are difficult to derive, because the Hessian is zero at the critical points and one has to expand the dispersion relation to third order, and because it is not clear how to cut off the logarithmic divergence correctly. Here we “cheat” by using the leading terms in the series expansion of the closed form involving elliptic integrals:katsura1971 ; guttmann2010

[TABLE]

We now construct an approximation using the functions listed in Table 6, and write down its Chebyshev coefficients:

[TABLE]

where $H_{n}=\sum_{k=1}^{n}1/k$ are the harmonic numbers. Since $H_{n}\sim\ln n$ for large $n$, we see that $f_{n}\sim(\ln n)/n$. The exact Chebyshev coefficients, from Table 2, are

[TABLE]

As before, we compute $h_{n}=g_{n}-f_{n}$ and $h(\omega)$. The results are shown in Figs. 4. Again, the singularity subtraction has decreased the truncation error of the 1000-term series from $10^{-3}$ to $10^{-9}$.

The results agree with closed forms for $g^{\text{sq}}(\omega)$ and $g^{\text{bcc}}(\omega)$ in terms of complete elliptic integrals.guttmann2010

Singularity subtraction for nonlocal spectral functions: We now consider the spectral function, Eq. (24), for injecting a particle at the origin and removing it at position $\mathbf{r}$. We have

[TABLE]

Van Hove singularities in $g_{\mathbf{r}}(\omega)$ arise from regions $\mathbf{k}~{}\approx~{}\mathbf{k}_{\alpha}$ where $\left\lvert{\nabla\varepsilon}\right\rvert\approx 0$. We can generalize Eq. (35) such that the contribution from each critical point $\mathbf{k}_{\alpha}$ is weighted by a different phase factor:

[TABLE]

Square lattice spectrum: For the square lattice, the dispersion relation $\varepsilon(\mathbf{k})$ has two saddle points at $(\pi,0)$ and $(0,\pi)$. Thus the logarithmic singularity at $\omega=0$ is modified to

[TABLE]

We have verified that subtracting this singularity reduces the truncation error in $g_{xy}(\omega)$ to about $10^{-9}$, similar to the case of $g_{00}(\omega)$.

BCC lattice spectrum: For the bcc lattice, $\varepsilon(\mathbf{k})$ has eight third-order saddle points at $(\pm\frac{\pi}{2},\pm\frac{\pi}{2},\pm\frac{\pi}{2})$. The nature of the dominant $(\ln^{2}\left\lvert{\omega}\right\rvert)$ singularities at $\omega=0$ in Fig. 3 can be explained in terms of Eq. (45). Unfortunately, we are unable to predict the coefficient of the subdominant $(\ln\left\lvert{\omega}\right\rvert)$ singularity, because it depends delicately on the cutoff of a logarithmic integral. Thus we proceed as follows. For simplicity we focus on the case where $x$, $y$, and $z$ are multiples of 4. We know the dominant singularity to be

[TABLE]

Calculate the dominant residual coefficients $g^{A}_{n}=g_{n}-f^{A}_{n}$. Assume the subdominant singularity is of the form

[TABLE]

Determine the coefficient $c$ by performing a least-squares fit, over a suitable range, such that $f^{B}_{n}\approx g^{A}_{n}$. Calculate the secondary residual coefficients ${g^{B}_{n}=g^{A}_{n}-f^{B}_{n}}$. Perform a Chebyshev transform to obtain $g^{B}(\omega)$. Finally, reconstruct the spectral function as

[TABLE]

We have verified, in a few cases, that this singularity subtraction approach reduces the truncation error in $g(\omega)$ to about $10^{-9}$. Figure 6 shows the case of $g^{\text{bcc}}_{400}(\omega)$. (According to Ref. ray2014arxiv, , $g^{\text{bcc}}_{xyz}(\omega)$ is expressible in terms of hypergeometric functions.)

Real part of Green function: In this section we have dealt with the DoS $g(\omega)$, which is proportional to $\mathop{\mathrm{Im}}G(\omega)$. To obtain the complex function $G(\omega)$ we may write $G(\omega)=F(\omega)+H(\omega)$ where

[TABLE]

Unfortunately, for many of the singular forms for $f(\omega)$ tabulated in Table 6, $F(\omega)$ cannot be calculated in closed form. Thus Eq. (53) may need to be evaluated numerically.

IV Discussion

In this paper we calculate the Chebyshev moments combinatorially. Another way to obtain the moments is to calculate $\left|\psi_{n}\right>\equiv T_{n}(\hat{H})\left|j\right>$ in the site basis by repeated application of the Chebyshev recursion formula $\left|\psi_{n+1}\right>=2\hat{H}\left|\psi_{n}\right>-\left|\psi_{n-1}\right>$. The Chebyshev moments are obtained by taking the inner product with $\left<i\right|$. This is referred to as the “spectral method,” “equation-of-motion method,” or “kernel polynomial method” for calculating the density of states.wangChebyshev1994 ; silver1994 ; silver1996 ; loh2014 ; karki2016 . With this approach, computing the $L$th Chebyshev moment requires storing values of three wavefunctions on $O(L^{d})$ sites, which may be memory-intensive. The “recursion method”berciu2010 is similar. The continued fraction methodberciu2010 may be more efficient but still requires $O(L^{d-1})$ storage.

In place of Chebyshev polynomials, one can use Legendre polynomials or any other family of orthogonal polynomials. Chebyshev polynomials have the advantage that Eq. (25) can be implemented using fast Fourier transform methods.

We have attempted to accelerate the convergence of the Chebyshev series using techniques such as Borel summation or Wynn’s epsilon rule. For a fixed value of $\omega$, we computed the partial sums of the Chebyshev series and applied convergence acceleration transformations. We only achieved limited success. In our opinion, it is more effective to fit the Chebyshev coefficients as in Eq. (50), which allows us to “accelerate” the convergence of the series “globally” for all values of $\omega$ simultaneously.

We have considered other methods of analytically continuing the inverse power series into the unit disk (Fig. 1). For example, one can use the inverse power series to calculate the derivatives of the Green function, $G^{(k)}(\Omega)$ ($k=0,1,2,\dotsc,L$), at a “pivot” point $\Omega=\omega+i\sqrt{1-\omega^{2}}$ in the complex plane. One can then calculate the Green function using the power series about the pivot point, $G(\omega)\approx\sum_{k=0}^{L}\frac{G^{(k)}(\Omega)}{k!}(\omega-\Omega)^{k}$. We have found that $G^{(k)}(\Omega)$ must be evaluated extremely accurately, up to very large values of $k$, in order to obtain $G(\omega)$ with modest accuracy. The “Chebyshev analytic continuation” method in this paper is preferable, as it naturally lends itself to exact integer arithmetic.

In many cases $G(\omega)$ is related to Gauss, Appell, or Lauricella hypergeometric functions.guttmann2010 ; ray2014arxiv In those cases our method may be viewed as a method (albeit an indirect one) for analytic continuation of hypergeometric functions.

V Conclusions

We have developed and demonstrated a general and efficient method for calculating lattice Green functions. The method relies on combinatorial formulas for the numbers of walks on the lattice, which are available for bcc-like, cubic-like, and honeycomb-like lattices. The method can be used to calculate imaginary parts (spectra) as well as real parts of the Green functions. The basic algorithm (Sec. II) gives Green functions to about 3 decimal places by summing Chebyshev series to 1000 terms. Singularity subtraction (Sec. III) increases the accuracy to about 6–9 decimals with little extra computational effort. Arbitrary-precision integer arithmetic is required. Fast cosine transforms and least-squares fitting routines may be useful when implementing the algorithms.