CalculateBands

E, C1, … = TightBinding.CalculateBands(HTB, kpath, f1, …) calculates the band structure of a tight-binding object HTB along a path in $k$-space: the eigenvalues of $H(k)$ at every $k$-point of the path and, for every function f given, the character of each eigenstate. G, G1, … = TightBinding.CalculateBands(HTB, kpath, Sigma, f1, …) does the same with a local self-energy $\Sigma(\omega)$. The bands are then no longer sharp, and the result is the band spectral function \begin{equation} A(k,\omega) = -\frac{1}{\pi} \mathrm{Im} \, \mathrm{Tr} \left( \omega + \mathrm{i} 0^{+} - H(k) - \Sigma(\omega) \right)^{-1} \end{equation} at every $k$-point, and its characters. Both results can be plotted with Graphics.PlotBands().

• HTB : Tight-binding object, which can be created using the function NewTightBinding()
• kpath : the path in $k$-space, a table { {Label1, k1}, N1, {Label2, k2}, N2, …, {Labeln, kn} }. Label is a string of which the plot shows the first character (so “Gamma” is drawn as G), k a table of three reals and N the number of $k$-points between two labelled points. The $k$-points are Cartesian, in the inverse length unit of HTB.Cell, such that $k \cdot r$ is a phase; the point X of a simple cubic lattice with lattice constant 1 is {pi,0,0}. The path holds $N_k = N_1 + N_2 + \dots + n$ points. With N equal to 0 the two labelled points follow each other directly, and the plot leaves a gap between them.
• Sigma : (optional) the self-energy, in any of the forms CalculateG() accepts: a response function of the dimension of all orbitals, or a table with one entry per atom (a response function, a table per shell, or 0).
• f1, f2, … : (optional) the characters. Each is either a function that takes an eigenstate $\psi$ (a vector over the orbitals of HTB, indexed psi[1], psi[2], …) and returns a number, usually between 0 and 1, or a matrix $P$ (made with Matrix.ToUserdata) for the character $\mathrm{Re} \sum_{ij} \psi_i^{*} P_{ij} \psi_j$. Without a self-energy 0 to 4 characters can be given, with a self-energy 0, 1 or 3. With a self-energy a function has to be a quadratic form, see below.

A matrix as a character has to be a Matrix userdata. A plain Lua table on the third position is always read as the self-energy.

Without a self-energy:

• E : a matrix of $N_k$ rows and $N_O$ columns, $N_O$ the number of orbitals of HTB. Row $k$ holds the eigenvalues of $H(k)$ in increasing order.
• C1, … : for every character a matrix of the same size, with the character of the eigenstate of the corresponding eigenvalue.

With a self-energy:

• G : a table of $N_k$ response functions (single valued, type ListOfPoles), the band spectral function $A(k,\omega)$ at every $k$-point of the path.
• G1, … : for every character a table of $N_k$ response functions, $A_P(k,\omega) = -\frac{1}{\pi} \mathrm{Im} \, \mathrm{Tr} \, P \left( \omega + \mathrm{i} 0^{+} - H(k) - \Sigma(\omega) \right)^{-1}$, with the same poles as G.

Without a self-energy the band spectral function is a pole of weight one at every eigenvalue, and the pole of G1 has the character as its weight. The result with a self-energy is therefore the same information in a different form: E and C1 become the pole energies and weights of G and G1. A self-energy that is zero gives exactly the eigenvalues with weight one, and the characters.

The character of a spectral function is only defined for a character that is a quadratic form of the eigenstate, $f(\psi) = \sum_{ij} \psi_i^{*} P_{ij} \psi_j$. The weight on atom A, Complex.Re(psi[1]*Conjugate(psi[1])), is one, and so is the bonding character Complex.Re¹⁾/2; a function like Abs(psi[1]) or a threshold is not. With a self-energy the function is evaluated on the unit vectors $e_i$ and on $(e_i + e_j)/\sqrt{2}$ and $(e_i + \mathrm{i} e_j)/\sqrt{2}$, which gives $P$, and then on two more vectors to check that it is such a form. If it is not, CalculateBands stops with an error. The weights of G1 add up to $\mathrm{Tr} \, P$ at every $k$-point, those of G to $N_O$.

For every $k$-point the self-energy is written as a bath, $\Sigma(\omega) = A_0 + C (\omega - H_{\mathrm{bath}})^{-1} C^{\dagger}$, and the matrix with $H(k) + A_0$ in its first block and the bath coupled to it by $C$ is diagonalised, as in CalculateG(). The first $N_O$ components $u$ of every eigenvector give a pole of G with weight $u^{\dagger} u$ and of G1 with weight $u^{\dagger} P u$. The poles are exact, there is no energy grid and no broadening; poles of weight below about $2 \cdot 10^{-13}$ are left out. The size of the matrix is $N_O$ plus the number of bath levels, and the cost per $k$-point grows with its cube: a self-energy with 50 poles on one orbital takes less than a millisecond per $k$-point, one with 400 poles about 0.2 seconds. A self-energy with many poles can be made smaller first with ResponseFunction.ChangeType() to “Tri” with the option “NTriMax”.

The self-energy is added to $H(k)$ in the same basis as in CalculateG, in which the phase of an orbital depends only on its unit cell. The average of G over a regular $k$-mesh is therefore the trace of the Green's function CalculateG returns for that mesh. The eigenvectors are then transformed to the phases CalculateBands uses without a self-energy, in which the phase of an orbital depends on the position of its atom, so that a character means the same with and without a self-energy.

Two p orbitals on a honeycomb lattice, the model of the tight-binding object, with the weight on atom A as character. The band structure is calculated without a self-energy and with a self-energy on atom A only, of two poles at $\pm 2$. At $\Gamma$ each band splits into a quasiparticle band and a satellite. The weights of the four poles add up to 2, their weights on atom A to 1.

Example.Quanty

-- two p orbitals on a honeycomb lattice
dAB = 0.2
tnn = 1.1
HTB = NewTightBinding()
HTB.Name = "dichalcogenide tight binding"
HTB.Cell = {{sqrt(3),0,0},
            {sqrt(3/4),3/2,0},
            {0,0,1}}
HTB.Atoms = { {"A", {0,0,0},       {{"p", {"0"}}}},
              {"B", {sqrt(3),1,0}, {{"p", {"0"}}}}}
HTB.Hopping = {{"A.p","A.p",{         0,   0,0},{{-dAB/2}}},
               {"B.p","B.p",{         0,   0,0},{{ dAB/2}}},
               {"A.p","B.p",{         0,   1,0},{{ tnn  }}},
               {"B.p","A.p",{         0,  -1,0},{{ tnn  }}},
               {"A.p","B.p",{ sqrt(3/4),-1/2,0},{{ tnn  }}},
               {"B.p","A.p",{-sqrt(3/4), 1/2,0},{{ tnn  }}},
               {"A.p","B.p",{-sqrt(3/4),-1/2,0},{{ tnn  }}},
               {"B.p","A.p",{ sqrt(3/4), 1/2,0},{{ tnn  }}}}
 
kpath = { {"M",{0,2*pi/3,0}}, 100, {"G",{0,0,0}}, 100, {"K",{4*pi/(3*sqrt(3)),0,0}}, 100, {"M",{2*pi/(sqrt(3)),0,0}} }
 
-- the character of an eigenstate: its weight on atom A
function OnA(psi)
  return Complex.Re(psi[1]*Conjugate(psi[1]))
end
 
-- without a self-energy: band energies and characters, one row per k-point
E, CA = TightBinding.CalculateBands(HTB, kpath, OnA)
print("without a self-energy, at Gamma (k-point 101):")
print(string.format("  energies %8.4f %8.4f   character A %6.4f %6.4f", E[101][1], E[101][2], CA[101][1], CA[101][2]))
 
-- a self-energy on atom A only, with two poles at -2 and 2
SigmaA = ResponseFunction.New({{0, -2, 2}, {0.5, 0.5}, mu=0, type="ListOfPoles"})
G, GA = TightBinding.CalculateBands(HTB, kpath, {SigmaA, 0}, OnA)
print("with a self-energy on atom A, at Gamma:")
t, tA = ResponseFunction.ToTable(G[101]), ResponseFunction.ToTable(GA[101])
for i = 1, #t[2] do
  print(string.format("  pole %8.4f   weight %6.4f   on A %6.4f", t[1][i+1], t[2][i], tA[2][i]))
end
 
-- the plots, the bands in blue (on A) to red (on B)
Frame = {"Frame",{{"Ymin",-4},{"Ymax",4},{"dYTick",1},{"YLabel","Energy"}}}
pl0 = Graphics.PlotBands(kpath, E, CA, {Frame})
pl1 = Graphics.PlotBands(kpath, G, GA, {Frame, {"GammaL",0.05}})
for name, pl in pairs({BandsWithout=pl0, BandsWithSigma=pl1}) do
  file = io.open(name..".svg", "w")
  file:write(Graphics.ToSVG(pl, {{"RelativeSize",true}}))
  file:close()
end

without a self-energy, at Gamma (k-point 101):
  energies  -3.3012   3.3012   character A 0.5151 0.4849
with a self-energy on atom A, at Gamma:
  pole  -3.5235   weight 0.8875   on A 0.4852
  pole  -1.8712   weight 0.1119   on A 0.0294
  pole   1.8825   weight 0.1089   on A 0.0246
  pole   3.5122   weight 0.8917   on A 0.4608

The two plots are written to BandsWithout.svg and BandsWithSigma.svg.