4.3. kwant.kpm – Kernel Polynomial Method

class kwant.kpm.SpectralDensity(hamiltonian, params=None, operator=None, num_vectors=10, num_moments=None, energy_resolution=None, vector_factory=None, bounds=None, eps=0.05, rng=None)[source]

Bases: object

Calculate the spectral density of an operator.

This class makes use of the kernel polynomial method (KPM), presented in [R12], to obtain the spectral density \(ρ_A(e)\), as a function of the energy \(e\), of some operator \(A\) that acts on a kwant system or a Hamiltonian. In general

\[ρ_A(e) = ρ(e) A(e),\]

where \(ρ(e) = \sum_{k=0}^{D-1} δ(E-E_k)\) is the density of states, and \(A(e)\) is the expectation value of \(A\) for all the eigenstates with energy \(e\).


hamiltonian : FiniteSystem or matrix Hamiltonian

If a system is passed, it should contain no leads.

params : dict, optional

Additional parameters to pass to the Hamiltonian and operator.

operator : operator, dense matrix, or sparse matrix, optional

Operator for which the spectral density will be evaluated. If it is callable, the densities at each energy will have the dimension of the result of operator(bra, ket). If it has a dot method, such as numpy.ndarray and scipy.sparse.matrices, the densities will be scalars. If no operator is provided, the density of states is calculated with a faster algorithm.

num_vectors : positive int, default: 10

Number of random vectors for the KPM method.

num_moments : positive int, default: 100

Number of moments, order of the KPM expansion. Mutually exclusive with ‘energy_resolution’.

energy_resolution : positive float, optional

The resolution in energy of the KPM approximation to the spectral density. Mutually exclusive with ‘num_moments’.

vector_factory : function, optional

The user defined function f(n) generates random vectors of length n that will be used in the algorithm. If not provided, random phase vectors are used. The default random vectors are optimal for most cases, see the discussions in [R12] and [R22].

bounds : pair of floats, optional

Lower and upper bounds for the eigenvalue spectrum of the system. If not provided, they are computed.

eps : positive float, default: 0.05

Parameter to ensure that the rescaled spectrum lies in the interval (-1, 1); required for stability.

rng : seed, or random number generator, optional

Random number generator used by vector_factory. If not provided, numpy’s rng will be used; if it is an Integer, it will be used to seed numpy’s rng, and if it is a random number generator, this is the one used.


When passing an operator defined in operator, the result of operator(bra, ket) depends on the attribute sum of such operator. If sum=True, densities will be scalars, that is, total densities of the system. If sum=False the densities will be arrays of the length of the system, that is, local densities.

[R12](1, 2) Rev. Mod. Phys., Vol. 78, No. 1 (2006).
[R22]Phys. Rev. E 69, 057701 (2004)


In the following example, we will obtain the density of states of a graphene sheet, defined as a honeycomb lattice with first nearest neighbors coupling.

We start by importing kwant and defining a FiniteSystem,

>>> import kwant
>>> def circle(pos):
...     x, y = pos
...     return x**2 + y**2 < 100
>>> lat = kwant.lattice.honeycomb()
>>> syst = kwant.Builder()
>>> syst[lat.shape(circle, (0, 0))] = 0
>>> syst[lat.neighbors()] = 1

and after finalizing the system, create an instance of SpectralDensity

>>> fsyst = syst.finalized()
>>> rho = kwant.kpm.SpectralDensity(fsyst)

The energies and densities can be accessed with

>>> energies, densities = rho()


>>> energies, densities = rho.energies, rho.densities


energies (array of floats) Array of sampling points with length 2 * num_moments in the range of the spectrum.
densities (array of floats) Spectral density of the operator evaluated at the energies.

Return the spectral density evaluated at energy.


energy : float or sequence of float, optional


float, if energy is float, array of float if energy

is a sequence, a tuple (energies, densities) if

energy was not provided.


If energy is not provided, then the densities are obtained by Fast Fourier Transform of the Chebyshev moments.

add_moments(num_moments=None, *, energy_resolution=None)[source]

Increase the number of Chebyshev moments.


num_moments: positive int

The number of Chebyshev moments to add. Mutually exclusive with ‘energy_resolution’.

energy_resolution: positive float, optional

Features wider than this resolution are visible in the spectral density. Mutually exclusive with ‘num_moments’.


Increase the number of random vectors.


num_vectors: positive int

The number of random vectors to add.


Returns the total spectral density.

Returns the integral over the whole spectrum with an optional distribution function. distribution_function should be able to take arrays as input. Defined using Gauss-Chebyshev integration.

Previous topic

4.2. kwant.rmt – Random matrix theory Hamiltonians

Next topic

4.4. kwant.continuum – Tools for continuum systems

This Page