Theory
Info
Here, we briefly introduce the main concepts that underlie all spectral methods currently implemented in KITE. Focus is placed on the role of the simulation parameters that can be adjusted by the enduser in order to suit its specific purposes.
The central object to any calculation done in KITE is the lattice singleparticle Hamiltonian (SPH) — \(H\) — which is always a sparse \(D\!\times\!D\) hermitian matrix (\(D\) being the total number of orbitals in the lattice). This fully embodies the simulated system, from its underlying bravais lattice structure and local orbital basis, to all the nonperiodic terms that realize disordered potentials, specific boundary conditions, complex structural defects and external magnetic fields.
Any finitedimensional SPH has a bounded realvalued spectrum which must be shifted and rescaled to suitably fit within \([1,1]\), the convergence interval of the method. This conversion is performed internally by KITE, which transforms \(H\to\mathcal{H}=(H\!\!\varepsilon_{0})/\delta\varepsilon\) and rescales all energy variables by \(\delta\varepsilon\).
This step requires an early (over)estimation of \(H\)'s spectral bandwidth, which can be set manually (using spectrum_range
= \([\,\varepsilon_{0}\!\!\delta\varepsilon/2\),\(\,\varepsilon_{0}\!+\!\delta\varepsilon/2]\) in kite.Configuration()
) or be automatically done by KITE upon generation of the hdf5 configuration file. At this stage, the user must also specify:
 The dimensions of the simulated lattice,
length=[lx,ly,(lz)]
;  The number of subdivisions for parallelization of the matrixvector operation,
divisions=[nx,ny,(nz)]
, where \(n_{x}n_{y}n_{z}\) is the available CPUcores;  The type of boundary conditions;
 The type of data to be handled in internal arithmetic operations (for further information see Settings).
Generally, a target function \(\mathcal{Q}\) evaluated by KITEx
fits one of the following forms:
where \(\mathcal{O}_{j}\) are sparse lattice operators (identities, velocity operators or spin operators) and \(F_{j}\) are functions of \(\mathcal{H}\), as well as other scalar parameters \(\left\{ \lambda_{i}^{j}\right\} _{i}\), such as the energy or a frequency. Furthermore, the \(\left\Psi\right\rangle\) in the above equation is a specific state vector of the system that depends on the observable that is being computed. Hence, there are two district categories of observables that can be computed with KITEx
, namely
 Traces Over The Entire Hilbert Space

This includes global observables such as the Density of States (DoS) and DCconductivity, as well as the \(1^{\text{st}}\) and \(2^{\text{nd}}\)order optical conductivity (see João et al.^{1} for further details). For all these cases, the trace is evaluated stochastically as an average of expectation values for \(R\) normalized random vectors^{2}, i.e.,
 Within the user interface, the number of independent random vectors is specified by the parameter
num_random
, which must be large enough to ensure a wellestimated trace. The associated error scales as \(1/\sqrt{R\,D}\), and thus requires very few random vectors if the simulated system is very large^{2}. On top of this averaging, if \(\mathcal{H}\) has a random component (by hosting disorder or featuring randomly twisted boundaries), it is often the case that the results are to be averaged over an ensemble of random Hamiltonians. Such averaging is also done insideKITEx
and the number of random configurations is specified by user with the parameternum_disorder
.  Diagonal Matrix Elements

This class of target functions includes local observables such as the local density of states (LDoS) and the \(\mathbf{k}\)space spectral function (for ARPES's response), as well as the timeevolution of Gaussian wavepackets. Note that
num_random
is no longer a relevant parameter for these target functions. 
In both classes, the core of the method is to expand the functions \(F_{j}\) as a truncated Chebyshev series of \(\mathcal{H}\), which allows one to write
 where \(T_{n}(x)\) are Chebyshev polynomials of the \(1^{\text{st}}\)kind and \(G_{j}\) are the expansion coefficients of \(F_{j}\). The advantage gained by using the Chebyshev expansion is that both \(\text{Tr}\left[\cdots\right]\) and \(\left\langle\Psi\right\cdots \left\Psi\right\rangle\) can be evaluated recursively using only matrixvector operations^{1}. For all observables implemented in KITE, the functions \(F_{j}\) are of three types:

 SingleParameter Dirac\(\delta\) Functions.—\(F_{j}(\lambda_{1}^{i};\mathcal{H})\to\delta\left(\lambda\mathcal{H}\right)\)
 Broadened SingleParticle Green's Functions.—\(F_{j}(\lambda_{1}^{i},\lambda_{2}^{i};\mathcal{H})\to\left[\lambda+i\eta\mathcal{H}\right]^{1}\)
 Quantum TimeEvolution Operators.—\(F_{j}\left(\lambda_{1}^{i},\mathcal{H}\right)\to\exp\left(\frac{i\,t\,\mathcal{H}}{\hbar}\right)\)

For these functions, analytical forms of the Chebyshev expansion coefficients are known^{2}^{3}^{4}^{5}^{6}^{7} and used in KITE. In the user interface, the truncation order \(M\) is specified by the parameter
num_moments
, and always impacts the validity of the expanded results. Nevertheless, its precise effect depends crucially on the specific case, as shown in the Figure below. In particular, one has the following cases:  Dirac\(\delta\) Function

An order\(M\) expansion (regularized by the Jackson kernel) produces a Gaussian approximation of \(\delta(\lambda\!\!\mathcal{H})\) endowed by a width \(\sigma_{\lambda}\!\approx\!\delta \varepsilon\,\pi/M\) in \(\lambda\)^{2}^{7}. The choice of \(M\) then fixes the effective spectral broadening, \(\sigma_{\lambda}\), which must be sufficiently narrow to accurately describe all relevant features of the calculated curve. However, if it becomes too narrow (\(M\) too high), the discrete eigenvalues of the SPH are wellresolved and the obtained data start suffering from large (finitesize) fluctuations. For information on other available kernels see Weisse et al.^{2}.
Rule of Thumb
If \(\Delta\varepsilon\) is the meanlevel spacing of the simulated system (that depends on the system size), then \(M\) must be kept smaller than \(\frac{\pi\,\delta\varepsilon}{\Delta\varepsilon}\) in order to avoid resolving individual energy levels. Simultaneously, for obtaining highresolution results, the arificial broadening much remain much smaller than the total bandwidth, i.e., \(M\gg \pi\,\delta\varepsilon\).
 SingleParticle Green's Functions

Usually, no kernel is required here^{3}^{4} as singularities of the Green's function are naturally broadened by a finite \(\eta\). Provided \(\eta\) exceeds the spacing between eigenvalues of the SPH, the truncation order may be arbitrarily increased and convergence is achieved when the data ceases to depend on M. Smaller values of \(\eta\) typically require higher values of \(M\) to converge.
Rule of Thumb
Since the energy resolution is fixed by \(\eta\), the number of polynomials must be larger enough to resolve such a broadening. As shown below, in Figure (b), an apt rule of thumb is to have \(M\gtrsim 100/\eta\).
 Quantum TimeEvolution Operators

The truncation error introduced here translates into a limitation of the available timeinterval of accurate unitary timeevolution. In the panels of Figure (c) it is demonstrated that the timeevolution operator is converged as long as \(t\!\lesssim\!\hbar\,M/\delta\varepsilon\).
Rule of Thumb
A safe empirical rule of thumb (used in Santos Pires et al.^{6}) is to have \(M\!\gtrsim\!8\,\delta\varepsilon\,\hbar^{1}t_{\text{max}}\!\!\), where \(t\!<\!t_{\text{max}}\) is the length of the timeinterval intended for the evolution.

For some target functions, the output of KITE will be a function of energy (e.g., the DoS) with maybe additional spacial coordinates (as it happens with the LDoS or the ARPES response). In contrast, most response functions are not properties of the Fermi level (the zerotemperature longitudinal conductivity is the only exception here) and therefore require that the raw output of KITE is numerically integrated in energy (typically weighted by the FermiDirac distribution at a given finite temperature). This step is always done at the postprocessing level by
KITEtools
.PostProcessing Integration
For target functions that require energy integrations, a key parameter is the number of energy points for which the integrand is evaluated. In order to adjust it at the postprocessing level, the user can use the
E
flag of theKITEtools
executable.

KITE: highperformance accurate modelling of electronic structure and response functions of large molecules, disordered crystals and heterostructures, S. M. João, M. Anđelković, L. Covaci, T. G. Rappoport, João M. Viana Parente Lopes, and A. Ferreira, R. Soc. open sci. 7, 191809 (2020). ↩↩

Kernel polynomial method, Alexander Weiße, Gerhard Wellein, Andreas Alvermann, and Holger Fehske. Rev. Mod. Phys. 78, 275 (2006). ↩↩↩↩↩

Critical delocalization of chiral zero energy modes in graphene, A. Ferreira and E. Mucciolo, Phys. Rev. Lett. 115, 106601 (2015). ↩↩

Numerical evaluation of Green's functions based on the Chebyshev expansion, A. Braun and P. Schmitteckert, Phys. Rev. B 90, 165112 (2014). ↩↩

An accurate and efficient scheme for propagating the time dependent Schrödinger equation, H. TalEzer and R. Kosloff, J. Chem. Phys. 81, 39673971 (1984). ↩

Landauer transport as a quasisteady state on finite chains under unitary quantum dynamics, J. P. Santos Pires, B. Amorim, and J. M. Viana Parente Lopes, Phys. Rev. B 101, 104203 (2020). ↩↩

Spectral functions of onedimensional systems with correlated disorder, N. A. Khan, J. M. Viana Parente Lopes, J. P. Santos Pires, J. M. B. Lopes dos Santos, J. Phys.: Condens. Matt. 31, 175501 (2019). ↩↩