Kwant: a software package for quantum transport

Kwant is a Python package for numerical quantum transport calculations. It aims to be an user-friendly, universal, and high-performance toolbox for the simulation of physical systems of any dimensionality and geometry that can be described by a tight-binding model. Kwant has been designed such that the natural concepts of the theory of quantum transport (lattices, symmetries, electrodes, orbital/spin/electron-hole degrees of freedom) are exposed in a simple and transparent way: Defining a new simulation setup is very close to describing the corresponding mathematical model. Kwant offers direct support for calculations of transport properties (conductance, noise, scattering matrix), dispersion relations, modes, wave functions, various Green's functions, and out-of-equilibrium local quantities. Other computations involving tight-binding Hamiltonians can be implemented easily thanks to its extensible and modular nature. Kwant is free software available at http://kwant-project.org/.


INTRODUCTION
Solving the scattering problem is one of the most common and general tasks in condensed matter physics. Instead of describing states in a closed geometry, one considers the scattering of particles in a finite system coupled (possibly strongly) to infinite leads. Its solution by itself directly yields the conductance and various other transport properties, but it can also be used as a building block for the calculation of more complicated physical phenomena, such as supercurrent, nonequilibrium density of states at a high voltage bias, or the evaluation of the topological properties of a topological insulator.
The history of numerical simulation of the scattering problem goes back to the early days of mesoscopic physics [1][2][3] when the first algorithms were developed. The most popular one of these is the recursive Green's function algorithm (RGF). Various groups created their own implementations of it, which quickly became an invaluable tool to verify, extend, or even replace the analytical approach even despite being restricted to quasi-one-dimensional geometries and to a particular type of tight-binding Hamiltonian. Besides quantum transport, the scattering problem naturally emerges in other contexts and many packages with a different focus (such as density functional theory, or transistor simulations) were developed [4][5][6][7][8][9][10]. Nevertheless, up to now no package existed whose main emphasis is efficiently solving with comparatively little effort the scattering problem for arbitrary single-particle tight-binding Hamiltonians.
Here we introduce Kwant, a publicly available package that is designed to • solve the scattering problem in a robust and highly efficient way, • exhibit a high degree of inter-operability with other packages and algorithms from any part of the code, including both defining and solving the scattering problems, • support an easy and expressive way to define a broad range of tight-binding systems as required for exploratory research.
Kwant uses highly efficient and robust algorithms that allow one to (i) significantly outperform the most commonly used recursive Green's function method and (ii) avoid usual instabilities occurring with many commonly used algorithms (for instance in dealing with the evanescent modes of complex electrodes). Inter-operability removes the need from specialized packages to reimplement the solution of scattering problem, while benefiting from the advanced and efficient algorithms used in Kwant. Finally, expressiveness is an especially important feature for mesoscopic physics, since it allows one to define a broad range of physical systems using the associated physics concepts directly. In short, the way one writes down a Hamiltonian in Kwant is very close to what one would write on a blackboard. The definition of a physical system amounts to writing a simple Python program that operates with physical concepts such as lattices, shapes, symmetries, and potentials. We hope that the free availability of a user-friendly, generic and high-performance code for quantum transport calculations will help to advance the field by allowing researchers to concentrate more on the physics and to perform computations that were considered out-of-reach due to their complexity. An example of a device that was simulated with Kwant is shown in Fig. 1: a cylindrical semiconducting wire with spin-orbit interaction, partially covered by a superconductor, used to create Majorana fermions [11][12][13].  [14]. The balls represent lattice sites. The red region is a tunnel barrier, used to measure tunneling conductance, the blue region is a superconductor.
to which a few semi-infinite periodic electrodes are connected. Within the Landauer-Büttiker formalism these leads act as wave guides leading plane waves into and out of the scattering region and correspond to the contacts of a quantum transport experiment. The Hamiltonian for such a system takes the form where c † i (c j ) are the usual fermionic creation (destruction) operators, i and j label the different degrees of freedom of the system, and H i j are the elements of an infinite Hermitian matrix. Alternatively, the same Hamiltonian can be written in first quantization asĤ = ∑ i, j H i j |i j| . (2) The degrees of freedom usually take the form |i = |r r rα , where r r r corresponds to the lattice coordinates of a site, and α labels its internal degrees of freedom. These can include any combination of spin, atomic orbital, and Nambu electron-hole degree of freedom for superconductivity. We may expressĤ as a sum over all the Hamiltonian fragmentŝ H r r rr r r = ∑ αα H r r rαr r r α |r r rα r r r α , that couple sites r r r and r r r . Hamiltonians of this form can arise directly from an approximate atomic description of a physical system, in which case sites correspond to atoms or molecules. Alternatively, a finite-difference discretization of a continuum Hamiltonian also results in a tight-binding Hamiltonian.
Kwant represents such Hamiltonians as annotated infinite graphs like the one shown in Fig. 2. Each node of the graph corresponds to a site r r r and is annotated with the typically small Hermitian matrix H r r rr r r that is a representation ofĤ r r rr r r . Each edge between sites r r r and r r r corresponds to a non-zero H r r rr r r = H † r r r r r r . The periodicity of the leads allows a finite representation of these infinite objects.
In summary, defining a scattering geometry amounts to defining a graph and the matrices H r r rr r r associated with it.

Scattering theory
We focus on the wave function formulation of the scattering problem due to its simpler structure compared to nonequilibrium Green's functions, and since Kwant's default solver is based on the wave function approach. The nonequilibrium Green's function formalism is mathematically equivalent to the wave function approach due to the Fisher-Lee relation, however it is less stable [15].
Without loss of generality, we can consider the case with a single lead. This is possible because several leads can always be considered as a single effective lead with disjoint sections. In the basis in which the sites are ordered according to the reverse distance to the scattering region (scattering region S last, first unit cell of the lead before that, second unit cell before the first one, etc.), the Hamiltonian of such a system has the tridiagonal block form where H S is the (typically large) Hamiltonian matrix of the scattering region S. H L is the (typically much smaller) Hamiltonian of one unit cell of the lead, while the block submatrix is the Hamiltonian V L connecting one unit cell of the lead to the next. Finally, V LS is the hopping from the system to the leads. We define the wave function of an infinite system as (. . . , ψ L (2), ψ L (1), ψ S ), where ψ S is the wave function in the scattering region, and ψ L (i) the wave function in the i-th unit cell away from the scattering region in the lead. Due to the translational invariance of the leads, the general form of the wave function in them is a superposition of plane waves. The eigenstates of the translation operator in the lead take the form such that they obey the Schrödinger equation with χ n the n-th eigenvector, and λ n the n-th eigenvalue. The normalizability requirement on the wave function reads |λ n | ≤ 1. The modes with |λ n | < 1 are evanescent, the rest are propagating and λ n = e ik n can be expressed in term of the longitudinal momentum k n of mode (or channel) n. The propagating modes are normalized according to the expectation value of the particle current, such that The modes are further sorted into incoming ones φ in n ( I = +1), outgoing ones φ out n ( I = −1) and evanescent ones φ ev n ( I = 0). With these notations, the scattering states in the leads take the form and the scattering wave function inside the system The scattering matrix S nm and the wave function inside the scattering region φ S n are the main raw outputs of Kwant. Their calculation can be done by matching the wave function in the leads with the one in the scattering region which amounts to inserting the above form of the wave function into the tightbinding equations Hψ n = εψ n , with H given by Eq. (4).
Examples of transport properties that can be obtained from the scattering matrix include conductance, shot noise, spin currents, Peltier and Seebeck coefficients and many other quantities. For instance, the differential conductance G ab = dI a /dV b (where a and b label two electrodes) is given by the Landauer formula The internal properties of the system such as local density of states or current density can be obtained from the φ S n using the general relation where f n (E) = 1/[1 + e (E−µ n )/kT n ] is the Fermi function of the lead to which channel n is associated.

THE DESIGN OF KWANT
As explained in the introduction, we have designed Kwant for performance and interoperability on the one hand, and flexibility and ease of use on the other. Combining all of these requirements is a nontrivial task since flexibility and ease of use can be best achieved using features of a high-level language (Python in our case), while performance and interoperability require a universal low-level language interface as well as simple data structures.
This apparent contradiction is resolved once we notice that while the time necessary for defining the tight-binding Hamiltonian will scale linearly with system size (if implemented efficiently), solving the scattering problem or applying any other relatively complicated numerical algorithm will scale less well. This allows us to separate the work into two phases: During the first the tight-binding Hamiltonian is constructed in Python using concepts that are natural in physics such as lattices, shapes, and functions of coordinates (e.g. an electrostatic potential V = tanh(αx) or a hopping in presence of magnetic field, t = t 0 exp[iB z (x 1 −x 2 )(y 1 +y 2 )]). Subsequently, the Hamiltonian is prepared for the second phase by transforming it into a low-level representation. During the following second phase, this optimized representation is used as input for high-performance numerical calculations.
As we show later in Sec. 7, this two-phase approach indeed does not cause any significant drop in performance. On the contrary, using the nested dissection algorithm [16] implemented in sparse linear algebra libraries, such as MUMPS [17,18], allows Kwant to significantly outperform a reference implementation of the RGF algorithm written in pure C. An additional advantage of the separation of defining tight-binding systems and solving the scattering problem is that the default implementations of either of the two phases can be substituted by other ones that are more adapted to a specific problem.

Defining tight-binding systems
The most natural way to think about a finite tight-binding system is to consider it as a mapping from the vertices and edges of a graph to the corresponding values of the Hamiltonian for the sites and hoppings. In Kwant such a mapping is represented by a Builder object. Its implementation as a hash table (Python dictionary) allows to efficiently add or remove sites and hoppings that are present in the system, thereby changing the geometry of the system incrementally.
Sites can be often classified by type of atom or the lattice to which they belong. We represent this in Kwant by defining a site to be a combination of a family and a tag. The site family identifies the class of the site, while the site tag identifies a unique site within this family. An example of a tight-binding system with only a single site family is shown in Fig. 3(a). The Builder object then maps every site (identified by tag α) to a value, in this case the onsite Hamiltonian term. Hoppings are represented as pairs of sites and mapped to the corresponding hopping Hamiltonian. Note that in this example, the site tag could be any object, such as an integer or a string. This way of defining a tight-binding system is completely universal. However, in the common case when many sites belong to the same crystal structure the notion of site family becomes more useful. In that case we can define a site family to represent each of the sublattices and the site tags to be tuples of integer coefficients (n 1 , . . . , n d ) that describe the site position coefficients in the basis of the Bravais lattice vectors a a a 1 , . . . , a a a d , with d the dimensionality of the lattice. This enables us to also use operations in real space, since the site tags are directly related to the real-space position x x x = ∑ i n i a a a i . For example we may operate with all sites of a lattice that lie within a certain geometrical shape, or make the Hamiltonian values depend on the real space position of the sites. This approach applied to a honeycomb lattice with two site families corresponding to its two sublattices is shown in Fig. 3

(b).
A tight-binding model on a regular lattice typically contains only a few 'kinds' of hoppings, with hoppings of a single kind being related to each other by lattice translations. In order to make use of that, Kwant allows to manipulate all hoppings of a single kind in a single operation.
While it would be sufficiently general to only allow complex constants as elements of the Hamiltonian, Kwant allows two important generalizations motivated by common usage. First, it is natural to identify several degrees of freedom that occupy the same position in space, such as orbitals, with a single site. In this case the values of the onsite Hamiltonians and the hoppings become matrices, as defined in Eq. (3). Second, the Hamiltonian matrix elements can be often seen as functions of position and other parameters on which the calculation may depend. To accommodate this, Kwant allows the values of a Hamiltonian to be given by program subroutines that are evaluated at a later stage.
The final feature of Builder that we are going to discuss is symmetry support. In order to reduce the amount of bookkeeping, Kwant enforces the Hermiticity of tight-binding sys-tems: H r r rr r r = H † r r r r r r . In addition, tight-binding systems in Kwant are allowed to have a certain real space symmetry that is automatically enforced. Builders with a translational symmetry may be attached to other builders, thereby forming the leads in a scattering geometry, as described in Sec. 2.2.
The details of the Builder interface are described in Appendix A.

Low-level representation of tight-binding systems
The mapping format used by Builder objects is very useful for defining tight-binding systems, but it is a rather poor choice for performing computations. First of all, much of it is specific to Python, and hence hard to interface with code written in other languages. It is also not memory-efficient, and does not allow to easily calculate certain crucial properties of the tightbinding system, such as the Hilbert space dimension. In order to avoid these problems, we define a low-level representation of a tight-binding system suitable for efficient calculations. In essence this representation is a sparse graph of numbered sites and an array of values of sites and hoppings of this graph. The aspects that distinguish such a system from a regular sparse matrix are that the values may be functions, and that semiinfinite leads may be attached to it.
A low-level system may be created from a Builder by finalizing it. In the case where Kwant is interfaced with an external package, which performs e.g. a DFT calculation, a low-level system may also be defined directly by the other package, fully avoiding the usage of Kwant or even Python for defining the system. Once a low-level system has been created, it can be used to perform numerical calculations. Kwant uses low-level systems as an input to quantum transport solvers that calculate various quantities of interest for tight-binding systems with multiple attached semi-infinite leads. Kwant contains efficient implementations of algorithms for computing the scattering matrix, the retarded Green's function, the scattering wave functions, modes in the leads, and the local density of states.
Since low-level system is an universal format, any computation with tight-binding Hamiltonians can be trivially adapted to use Kwant low-level systems as input. If, for example, one is interested in the eigenstates of a quantum system, the full Hamiltonian matrix can be requested from the low-level system for a specific set of parameters and passed on to a standard eigenstate calculation routine as provided for example by ARPACK [19] (bundled for Python by SciPy [20]). On the other hand, the fact that the Hamiltonian values are implemented as functions allows to naturally integrate Kwant with a Poisson solver, and implement Hatree-Fock mean field calculations.
The solving phase and low-level systems are described in more detail in Appendix B.

COMPARISON WITH OTHER QUANTUM TRANSPORT PACKAGES
The quantum transport problem appears in many physical settings, and hence it is not surprising that it is addressed by various software packages from different domains. In particular, there are quite a few packages, including commercial ones, for computing transport in molecular junctions. Examples of these packages are TranSiesta/Atomistix Toolkit [4], SMEAGOL [5], OpenMX [6] or nanodcal/nanodsim [7]. These combine density functional theory (DFT) with the non-equilibrium Green's function technique. Another group of codes that deal with the transport problem is mainly geared towards the simulation of transistors on the nano-scale. It includes NEMO5 [8], nextnano [9], or NanoTCAD Vides [10], and TB_Sim [21]. These packages focus on more complicated physical effects, and often go beyond the scope of Kwant by including phonon effects, Boltzmann transport, or self-consistent electrostatic potential calculation. However, all these packages are specialized to a certain class of tight-binding systems, and extending them to include an additional physical effect, such as superconductivity is often impossible or requires a lot of work.
In contrast, Kwant focuses on generality in order to allow simulating the broad variety of complex models and geometries that are encountered in mesoscopic physics. There exist several private codes for solving the scattering problem (for example Refs. [22][23][24]) in addition to the publicly available KNIT package [25] in which one of us was involved. To the best of our knowledge Kwant significantly outperforms these packages since it uses the nested dissection algorithm instead of RGF (see Sec. 7 for details). Another advantage of Kwant, as described in the Sec. 3, is that it provides both an advanced set of tools for defining tight-binding models, and a universal interface that allows to interact with other codes.

ILLUSTRATION OF KWANT USAGE: UNIVERSAL CONDUCTANCE FLUCTUATIONS IN A QUANTUM BILLIARD
In order to illustrate how various aspects of Kwant work together when applied to a scattering problem we turn to the classic example of deterministic chaos in a stadium billiard. Despite their regular shape, stadium billiards are not integrable, showing an irregular density of states inside the scattering region, and universal conductance fluctuations. Thanks to Kwant's expressiveness, the Python program that defines the billiard system, performs numerical calculations, and creates two figures is less than 30 lines long. In the following the complete program is presented together with explanations.
The first step is to make Kwant's functionality available within Python, import kwant As described in Sec. 3.1, we need to create the builder object that will contain the information about the system being constructed. We also create the lattice that is used. We continue with defining a function that specifies the shape of the the scattering region. Given a point (x, y) this function returns True if the point is inside the shape and False otherwise.
def stadium(position): x, y = position x = max(abs(x) -7 , ) return x**2 + y**2 < 1 **2 We proceed to add these sites to the scattering region and set the corresponding values of the onsite potential to −4t (we use t = −1).
sys[sqlat.shape(stadium, ( , ))] = 4 The expression sqlat.shape(stadium, ( , )) represents all the sites of the lattice that belong to the stadium (provided they can be reached from the central point (0, 0)). We then set the hopping Hamiltonian matrix elements between nearest neighbors to t = −1: The scattering region is now fully defined.
To complete the scattering problem, as explained in Sec. 2.2, we need to define the leads. The procedure for a lead is very similar to that for the scattering region. The only difference between the two is that upon creation of each lead its symmetry is specified. All further operations with the lead will automatically respect this symmetry: For example, adding a single site to the lead will also add all the image sites under the symmetry as well. Thus, in order to provide enough information, specifying the structure of a single unit cell is sufficient. We construct two leads by defining the sites which belong to a unit cell of each lead, and attach them to the scattering region. This finishes the definition of the scattering problem. The next required step, as explained in Sec. 3.2, is to transform the system into a form suitable for efficient numerical calculations: We can now use the Kwant solvers to obtain various physical observables. We compute the conductance, given by Eq. The function kwant.smatrix return the scattering matrix of the system at a given energy. This scattering matrix is then used to calculate the transmission from one lead to another. The function kwant.ldos returns the local density of states in the scattering region. We finish the program by plotting the calculated data as shown in Fig. 4 Here, in order to make the leads visible, we have plotted the first ten lead unit cells.
The resulting plots (see Fig. 4) show the universal conductance fluctuations of the stadium billiard and its irregular, chaotic density of states. We see that defining the scattering problem and calculating the relevant physical observables is transparent, logical, and involves a minimal number of steps.

A FULL-SCALE APPLICATION: HANLE EFFECT IN A GRAPHENE-BASED NON-LOCAL SPIN VALVE
We continue with a simulation of a graphene-based nonlocal spin valve as sketched in the inset of Fig. 5. The device consists of a graphene nanoribon where the two sides serve as contacts 0 and 3. Two additional magnetic metallic leads (1, 2) are deposited on top of the nanoribon. This setup allows to measure the non-local resistance R 01,23 = V /I of the device: a current I is passed through contacts 0 and 1 and the voltage V is measured between 2 and 3. Such a device allows to measure a "pure" spin signal since no electrical current is flowing through the electrodes 2 and 3 where the voltage drop is measured, and was studied experimentally recently [26,27]. We model the system using a tight-binding Hamiltonian, where the sum over i j is restricted over nearest neighbours and the corresponding hopping amplitude t ss a takes different values inside the graphene layer (a = G), the magnetic electrodes (a = F) and at the graphene-ferromagnet interface (a = GF). s and s denote spin indices. The metallic magnetic leads are modeled using a cubic lattice; it is attached to the honeycomb lattice of graphene by adding a hopping t GF from each lattice point in the last slice of the cubic lattice to the nearest atom in the graphene lattice.
We take t F = t G = 1 while the spin-filtering due to the presence of the magnetic electrodes is included in the interfacial hopping, where −1 < β < 1 characterizes the spin polarization of the graphene-ferromagnet interface. Lastly, the graphene onsite potential contains static disorder plus an in-plane magnetic field H x perpendicular to the magnetization of the electrodes where the v i are random numbers uniformly distributed inside [−0.5, 0.5], W characterizes the strength of the disorder potential and σ x , σ z are the Pauli matrices. This is the minimal model that allows to simulate the Hanle effect in a lateral spin valve. Additional ingredients, such as magnetic disorder, could be added to introduce a finite spindiffusion length in the system.
Even though we present a minimal model, it is far more complex than the two-terminal geometries typically studied in numerical simulations. In particular, it is necessary to study a four-terminal geometry involving different Bravais lattices, and to include the spin degree of freedom. Kwant's design allows a natural implementation of this system: The different Bravais lattices are represented by different site families, and the spin degree of freedom by matrix Hamiltonians associated with sites and hoppings in Builder. In addition, the solver for computing the scattering matrix is general enough to allow for arbitrary geometries with an arbitrary number of contacts.
Apart from native Kwant features, the present example also benefits from the fact that Kwant is a Python package and as such can be easily embedded and extended by custom code. and must be custom-coded by the user. This is easily achieved by combining coordinate information provided by Kwant with a kd-tree (a data structure designed to find the nearest lattice point) provided by the SciPy package [20]. Also, the calculation of the non-local resistance requires solving a 3 × 3 linear equation, which is done in one line of code using a call to the NumPy package [28]. Finally, a few additional lines of code allow to parallelize the script for a multi-core workstation or cluster. Fig. 5 shows the result of a numerical simulation using Kwant: the non-local resistance as a function of H x for the two configurations where the magnetizations are parallel P and anti-parallel AP (the latter is obtained by setting β → −β in one of the magnetic electrodes). We find a typical Hanle signal: when the magnetic field is increased, the spin precesses around the x axis resulting in a change of the sign of the spindependent signal ∆R = R P − R AP . Since different trajectories have different lengths, the precession angle is spread over a finite window. Eventually, this makes ∆R vanish at large magnetic fields.
Hanle precession is typically described within a semiclassical diffusion description, and the current study is to our knowledge the first to take into account quantum coherence. In fact, the full quantum model presented here allows to go beyond a semi-classical description and study the effect of phase coherence and/or ballistic propagation (W = 0).

BENCHMARK
We now show a comparison of the performance of Kwant with a C implementation of the RGF algorithm applied to a prototypical quantum transport problem: the calculation of the conductance of a square tight-binding system. The system consists of L × L single-orbital sites that belong to a square lattice. On two opposite sides leads of width L are attached. The energy was kept fixed, so that the number of propagating modes in the leads was proportional to L. The measurements were performed on a computer with an Intel i5-2520M 64-bit processor and 8 GiB of main memory running a variant of GNU/Linux with single-threaded OpenBLAS [29]. Figure 6 shows the dependence of the running time on L and compares Kwant with an alternative code written entirely in C (using the same BLAS and LAPACK) that implements the RGF method. One can see that for large system sizes, Kwant with the MUMPS-based solver is up to ten times faster than the C RGF code. For small systems the situation is inverted but less dramatic, the cross-over occurs around L = 50. It is evident that for small systems Kwant's construction step takes up a considerable fraction of the total time. The much simpler RGF code only supports quasi-1-d geometries and therefore requires no system construction in the sense of Kwant. Construction is typically performed considerably less often (one time per geometry) than solving (once per geometry and set of parameters). Figure 7 compares the memory footprint of Kwant's MUMPS-based solver with that of an RGF solver implemented in C. Evidently, increased memory usage is the price for the superior speed of Kwant's solver. Note, however, that with main memory sizes commonly available today even on lowend computers (a few GiB) the MUMPS-based solver is able to tackle systems of more than 10 6 sites. To be able to handle even larger systems RGF has been implemented as well.  small by today's standards but whose inclusion would have obstructed the power-law character of the curves. It was measured as the increase of the maximum "resident set size" taken up by a given computation over that reported for an empty run of the respective software.

CONCLUSION
The Kwant project has a double objective. First, to gather high-performance algorithms that are useful in the field of quantum transport. Second, to provide a simple and clear but powerful user interface for defining and working with tight-binding models. We refrained from designing such an interface from scratch (an approach that usually does not age well) but rather chose to extend an existing language (Python) with new capabilities. In this approach, the usual input files present in most scientific software disappear and are replaced by small programs. This results in more flexibility, since tight integration with other packages becomes possible, as well as pre-and post-processing of the data. We believe that such a modern approach to scientific programming has a strong impact on the usefulness of a code. This paper describes version 1.0 of Kwant, the first version released to the public. Kwant is a work in progress and there are many ideas for improvements collected in a to-do-list that is available with the source code. For instance, more solvers could be added and Kwant's low-level system format could be modified to allow general symmetries.
Kwant is free (open source) software [30] (distributed under the liberal "simplified BSD license"), and we hope that the project's website http://kwant-project.org/ and mailing list http://kwant-project.org/community will develop into a hub for a community of users and developers. Contribu-tions to Kwant as well as the sharing of related code modules are welcome.
Site objects are Kwant's abstraction for the sites of tightbinding systems. Each site belongs to exactly one site family, typically a crystal lattice (monatomic crystals) or a crystal sublattice (polyatomic crystals). Within a site family, individual sites are distinguished by a tag. In the common case where the site family is a regular lattice, the site tag is simply given by its integer lattice coordinates.
Any crystal lattice with any basis can be created easily by specifying the primitive vectors of its Bravais lattice and the coordinates of its sites inside the lattice unit cell. For example, the honeycomb lattice of graphene can be created using Even though most examples in this article are 2-dimensional, Kwant is not limited to any specific dimensionality and uses a dimensionality-independent graph representation of tightbinding systems. Site families are in principle more general then Bravais lattices. For example, one could create a site family for an amorphous material or even label sites with names consisting of characters. In practice, however, Bravais lattices are sufficient to construct most systems.

A.2. Tight-binding systems as Python mappings
Having introduced an important prerequisite, sites, we now proceed to discuss the structure of builders. The main idea here is to represent a tight-binding system as a mapping from sites and hoppings (=pairs of sites) to the corresponding submatrices of the Hamiltonian. Similar to any other mapping in Python, builder items are set using the syntax sys[key] = value. The following example shows how to create a simple system by setting its sites and hoppings one by one.
First, a square lattice and an empty builder are initialized. The resulting system is shown in Fig. 8.
Builders are fully-fledged Python mappings. This means that in addition to setting values of sites and hoppings like in the above example, it is possible to query values, e.g. print sys[lat(1, )]. It is also possible to delete items: del sys[lat( , 1)].
Builder objects ensure their consistency during manipulation: They automatically remain Hermitian, and hoppings may only exist between sites present in the builder. Hence, the value of each hopping (i, j) is invariably bound to be the Hermitian conjugate of the hopping ( j, i), and when a site is deleted, all of the hoppings to it and from it are deleted as well.

A.3. Values of sites and hoppings
In the previous subsection the onsite Hamiltonians and hopping values were just scalar constants. While this allows to define any tight-binding system by introducing a sufficient number of sublattices, in many cases it is useful to also allow these values to be matrices. Thanks to the flexibility of Python, this works in the most natural way in Kwant. If instead of a scalar value of the Hamiltonian one needs to use a matrix value, one may just write Here we set the Hamiltonian of two sites to be equal to the Pauli matrices σ y and σ z . We used two different ways to define a matrix-like object: as a NumPy array, and a tinyarray. NumPy [28] is the standard Python array library, and tinyarray is a library developed for Kwant, specifically optimized to be fast with small arrays, as discussed in Appendix C. In this application, the only difference between these two is performance, which is better for tinyarray arrays. The number of orbitals may be different for different sites as long as the sizes of all the value matrices are consistent: the hopping integral H i j from site j to i must have the size n i × n j , with n i and n j the sizes of onsite Hamiltonian of these two sites. (In practice all the Hamiltonian element matrices are very often square and of the same size, e.g. 2 × 2 for a model with spin.) In order to simplify the definition of more complicated Hamiltonians, Kwant also supports values that are functions. If a function is assigned as a value of a builder for some site or hopping, its evaluation is postponed until the last possible moment, that is when the Hamiltonian is actually used in a calculation. (Performing these calculations is the topic of Appendix B.) Setting functions as values with a builder works exactly like setting constant values: sys[key] = value, with the only difference that the value is now a function, not a number or a matrix. The value function is called with the corresponding site or hopping as arguments, and optionally some additional parameters.
Value functions help to elegantly define numerical values that depend on position. Furthermore, they allow to use different numerical values with a single finalized system (varying parameters such as magnetic field or electrostatic potential).

A.3.1. Example: value functions for hoppings and sites
As an example let us use value functions to model a system with constant perpendicular magnetic field and disorder. The following code snippets are part of the quantum Hall effect example shown fully in Appendix D.5 and further used in Appendix B.2.
First, we consider the hoppings. We choose the vector potential in the Landau gauge with B B B = Bẑ, so that A A A = −Byx. Including a vector potential in a tight-binding system is done using Peierls substitution [31,32]. The hopping integral t i j (Φ) from the point (x j , y j ) to the point (x i , y i ) is then given by This function receives four arguments: The first two are the sites that are connected by the hopping for which the value is requested. The remaining two are user-specified additional parameters. phi corresponds to Φ. salt is not used here, but it will be used in the function onsite below. Declaring that second parameter is necessary because each value function of a given system receives the same user-specified parameters.
Having specified the hoppings, we define an uncorrelated Gaussian disorder potential for the sites. The usual approach for defining disorder is to use a random number generator and evaluate the disordered potential on each site in sequence. With Kwant, however, solvers are allowed to evaluate the value of a site or hopping more than once and expect that it does not change during a single invocation of a solver. Additionally, since the order in which the values of sites and hoppings are evaluated is undefined and depends on internal details of the 1 It is also possible to add a magnetic field to a system with non-parallel leads [33]: For each lead, an adopted Landau gauge is chosen. For the scattering region, a Landau gauge is also adopted with local gauge transformations applied to the interface sites of each lead that mediate between the different Landau gauges. A Kwant module that automates this procedure is planned to be made available.
builder and the used solver. One solution to this problem would be to generate a large Since this is a value function for sites, the first argument has to be the site on which it is going to be evaluated. phi is the user-specified parameter that was used with hopping above and is ignored here. Passing different string values to the salt parameter results in different realizations of the disorder, so salt plays a role similar to a random seed.

A.3.2. Example: value functions that return matrices
The following second example, quite different from the preceding one, demonstrates the flexibility of value functions. We consider a section of the Majorana fermion script of Appendix D. 6 and some Kronecker-products of them tau_z = tinyarray.array(numpy.kron(s_z, s_ )) tau_x = tinyarray.array(numpy.kron(s_x, s_ )) sigma_z = tinyarray.array(numpy.kron(s_ , s_z)) tau_zsigma_x = tinyarray.array(numpy.kron(s_z, s_x)) The calls to tinyarray.array are optional. Their purpose is to improve the performance of the value functions. These constants are used within the value functions that define the site Hamiltonians Note the following features of these value functions: • Their return values are matrices. A value function may return anything that would be valid as a constant value for a builder.
• They do not depend on their site parameters. Still, Kwant will call them for each site and hopping individually such that the same values will be re-calculated many times. Even though this may seem wasteful, as discussed in Appendix C such inefficiencies do not play a role in practice.
• A single extra parameter p is passed, as opposed to the many parameters that define the Hamiltonian. p is a namespace object (see Appendix D.6 for its definition) that contains all the actual parameters as its attributes. This technique is useful to avoid mistakes when passing many parameters to value functions.

A.4. Acting on multiple sites/hoppings at once
The features of builders that have been introduced so far are in principle sufficient to define (site-by-site and hoppingby-hopping) systems with an arbitrarily complex geometry. In order to facilitate a higher level of abstraction, in addition to the simple keys (single sites and single hoppings) builders can also be indexed by composite keys that combine multiple simple keys. This feature has been inspired by the "slicing" of some Python containers and by NumPy's "fancy indexing".
The most basic kind of these advanced keys is a list of simple keys. With such a list it is possible to assign a value to multiple sites or hoppings in one go. We now create a list of sites of the lattice lat belonging to a disk. for y in range(-1 , 11): if x**2 + y**2 < 13**2: sites.append(lat(x, y)) This list of sites can be added to the system in one step: Any non-tuple object that supports iteration can be used just like a list. For example, the following code snippet adds the same sites using a generator expression. The syntax is quite self-explanatory. We refer to the Python tutorial [37] for further information.
sys [(lat(x, y) for x in range(-1 , 11) for y in range(-1 , 11) if x**2 + y**2 < 13**2)] = .5 Finally, the most high-level way to add many sites at once is to use a method of the lattice named shape. This method finds all the sites of a lattice that fit inside a certain region and can be reached from a given starting point. Having defined a function that specifies the circular region def disk(pos): x, y = pos return x**2 + y**2 < 13**2 the code required for adding the same sites as before to the system is very compact: sys[lat.shape(disk, ( , ))] = .5 Note that we do not have to explicitly loop over the possible candidate sites anymore. The advantage of using shape becomes greater for non-square lattices and for more complicated shapes. We will see an example of such usage further below. With a Kwant builder, one will often first add all the sites to a system, and then all the hoppings that connect them. If we regard all hoppings that only differ by a lattice translation as a single kind, that second step typically adds only a few kinds of hoppings. These are represented by objects of the type HoppingKind from the module kwant.builder that can be used directly as keys: HoppingKind(displacement, lattice2, lattice1) generates all the hoppings inside the system that start at lattice1, end at lattice2, and whose lattice coordinates differ by displacement. Using HoppingKind to set all the nearest neighbor hoppings of a system on a square lattice can be done as To go even further, Kwant can calculate a list of all the n-th nearest neighbor hopping kinds on any lattice. That list, returned by the method neighbors can be used directly as a key: sys[lat.neighbors(1)] = 1 (The argument 1 to lat.neighbors means that nearest neighbors are to be found and could have been omitted since it is the default value.) This high-level key allows to access the n-th nearest neighbors of any regular lattice in a single line of code. The complete yet very compact construction script for the circular quantum dot is listed in Appendix D.2 and its output is shown in Fig. 9.
For further information about keys we refer to the documentation of Kwant, specifically the documentation of the method expand of Builder. Here, we will just show another example that demonstrates the flexibility of the approach. Consider Fig. 10 that has been generated by the script shown in Appendix D.3. This program is identical to the one for the disk except that it features a different type of lattice lat = kwant.lattice.honeycomb() and a different region-defining function: def bean(pos): x, y = pos rr = x**2 + y**2 return rr**2 < 15 * y * rr + x**2 * y**2 Note that the function defining the shape (bean in the example) receives points in real space. In Kwant, each site "lives" in two coordinate systems: the coordinates system of the lattice (with  integer coefficients) and the real-space "world" coordinates. These two coordinate systems are only identical in the simplest square lattice (that is actually used in several examples of this article), but not in general, like in the example above.

A.5. Symmetries
So far all examples in this section involved systems with a finite number of sites. Kwant supports infinite periodic systems as well: Upon creation of a builder a spatial symmetry can be specified that will be automatically enforced for that system. Infinite systems with a translation symmetry can be attached to finite systems as leads.
Kwant's abstraction of symmetries closely matches the mathematical concept of spatial symmetry: A symmetry in Kwant is specified by a set of spatial transformations of sites and hoppings that generate the symmetry group and a set of sites and hoppings called fundamental domain that contains a single representative from each set of mutually symmetrical sites/hoppings 3 . In order to completely describe an infinite periodic tight-binding system it is thus sufficient to know the symmetry and those sites and hoppings of the system that lie within its fundamental domain.
In practice, users of Kwant never specify symmetries in this low-level way, they use a predefined symmetry class instead. (Translational symmetry is by far the most important for defining scattering problems, and this is why it is the only class predefined in Kwant.) A kwant.TranslationalSymmetry object is initialized by one or several real-space vectors, each representing a translation under which the system is to remain invariant. Currently, the user has no control over the fundamental domain -it is deduced from the translation symmetry periods using a simple algorithm.
Sites and hoppings can be added and manipulated in the same way as in finite systems, with the difference that each site/hopping is now treated as a representative of all the sites/hoppings symmetrical to it. The following example first creates a builder with a symmetry. Then, all the sites symmetrical to lat(3, ) are added. Finally, the sites that have been just added are deleted using a site as key that differs from lat(3, ) but is symmetrical to it. To achieve this behavior, a builder with a symmetry internally maps all sites and hoppings to the fundamental domain. Each serves as a unique representative for itself and all of its images under the symmetry. This mapping is often invisible to the user, but can lead to confusion when the sites or hoppings of a builder with a symmetry are printed and seem to differ from the ones that were set.
Due to limitations of the current low-level system format, in Kwant 1.0 it is only possible to finalize systems with 1-d translational symmetry.

A.6. Leads
Infinite periodic systems can be attached to finite systems as leads. Kwant uses the convention that the period of the lead symmetry must point away from the scattering region. Attaching a lead works by specifying a set of sites of the finite system to which the lead is to be attached: the lead interface. The sites of the lead interface must belong to an image of a single unit cell of the lead under the lead symmetry, and the hopping between the lead and the lead interface is assumed to be equal to the hopping between the neighboring lead unit cells. This can always be achieved, requiring sometimes to add an extra unit cell of the lead to the scattering region.
To make attaching of leads easier, builders provide the method attach_lead. This method automatically adds sites and hoppings to the scattering region to make it compatible with the lead shape and the requirements imposed on the lead interface. The added sites and hoppings are assigned the values of the corresponding sites and hoppings in the lead. The lead interface is then calculated automatically.
The following example (with complete code shown in Appendix D.4) creates a ring-shaped scattering region. It then defines an infinite system with period (−2, 1) and adds to it a horizontal line of 12 sites with coordinates ranging from (−6, 0) to (5, 0). Thanks to the symmetry mechanism, the infinite system contains also all images of that horizontal line under all multiples of the period. This infinite system is attached twice as two different leads to the scattering region: once outside and once inside. The argument lat( , ) to the second call of attach_lead specifies the starting point for the algorithm: the lead is attached at the nearest possible place in the direction opposite to the lead direction. As they are attached, leads are labeled with consecutive integers starting from zero. A plot of the resulting system is shown in Fig. 11. The scattering region is represented by the black dots. The lead is represented by its first two unit cells (colored dots). Note that the lead unit cells are not perpendicular to the lead direction, as their shape is determined by the fundamental domain and the latter is chosen implicitly. This does not affect the physics of the system, it does, however, affect which sites have to be added to the system by attach_lead. Since tight-binding Hamiltonians are sparse, it is natural to store them as annotated graphs. That representation of the Hamiltonian matrix is illustrated in Fig. 12: A tight-binding degree of freedom i corresponds to a node of a graph, every non-zero hopping matrix element H i j to an edge connecting nodes i and j. The graph thus represents the structure of the non-zero entries of the Hamiltonian. Together with the values H i j the complete tight-binding Hamiltonian is defined.
Such annotated graphs are called low-level systems within Kwant and are implemented by a hierarchy of classes of objects in the module kwant.system. As explained in Sec. 2, low-level systems form the common interface between the two phases of construction and solving. Within a low-level system, the graph structure itself is stored in a compressed sparse row format [38]. The values of the onsite Hamiltonians and hoppings H i j are provided upon request by the method hamiltonian given i, j and possibly other parameters that are required to evaluate the Hamiltonian matrix elements (see Appendix A.3). Periodic infinite systems that can be attached as leads to a finite system are represented in a similar way. Since no Python-specific features are used, low-level systems are compatible with C, Fortran and similar languages. This makes Kwant independent of a single programming language: even though the current version of Kwant is mostly written in Python, systems and tools for working with systems can be implemented in other programming languages as well.
Most often, a low-level system will be obtained from a Builder instance by calling the method finalized. In this case, the mapping from system sites to the integers numbering the graph is not controlled by the user.
It is also perfectly possible to implement a low-level system directly, deriving from FiniteSystem: revealing the matrix structure of the example shown in Fig. 12. Here we have used the convenience method hamiltonian_submatrix that can return any submatrix of the Hamiltonian as a NumPy array, or as a SciPy sparse matrix, and that returns the full Hamiltonian matrix by default. The numbering of the graph and the numbering of the indices of the entries in the full Hamiltonian matrix coincide when the hopping matrix elements H i j are all scalar. In general, however, the hopping matrix elements can be n i × n j matrices H i j . The Hamiltonian matrix should then be interpreted as a block matrix, with the graph indices as block indices.

B.2. Quantum transport
Given the Hamiltonian matrix of the system and of the leads attached to it, various methods may be used to compute transport properties. The default algorithm used in Kwant is the wave function approach that has been introduced in Sec. 2.2 and amounts to setting-up and solving of a sparse system of linear equations. The full algorithm, including the case of noninvertible hopping matrices, a numerically stabilized version of the algorithm and a discussion of the connection to the Green's function formalism will be published elsewhere [15].
The sparse matrix of the full linear system to be solved contains the full Hamiltonian of the scattering region and additional rows and columns for outgoing and evanescent modes of the leads. At first glance, a direct solution might seem inefficient and even infeasible for large systems. However, there is a large variety of established libraries for solving general sparse systems of linear equations that allow nevertheless for a very efficient and stable solution. The efficiency of these sparse linear solvers depends crucially on choosing a good ordering of the coefficient matrix. For systems arising from regular grids, nested dissection orderings have been shown to be optimal. For such orderings, the time needed for solving a two-dimensional tight-binding system with the geometry of a rectangle of length L and width W scales as O(LW 2 ) [16]. 4 This is a more favorable scaling compared to the RGF algorithm whose execution time scales as O(LW 3 ) [1 -3]. Indeed, in practice we find that the transport algorithm in Kwant is considerably faster than RGF for systems of intermediate and large size, as shown in Sec. 7. Solving the linear system of equations yields the scattering matrix. Within the framework of Kwant, the function kwant.smatrix performs this calculation given a low-level system, the Fermi energy and optional user-defined parameters of the system. It returns a scattering matrix object that can be queried for e.g. the values of transmission or noise between leads. The following code, for example, calculates conductance as a function of a range of magnetic fluxes. Together with the value functions of Appendix A.3.1 and a few lines of code that construct the system that are shown in the full listing in Appendix D.5, the quantum Hall effect conductance plateaus of Fig. 13 are generated. Note that the list passed as the args parameter contains the values of the two user-defined parameters of the value functions. The arguments to transmission specify that the conductance from lead 0 to lead 1 is to be returned. In fact, kwant.smatrix is a shortcut for the function smatrix of kwant.solvers.default, the default solver module. This module offers additional functionality as well, for example the calculation of the scattering wave function (9) of each scattering state (8). When kwant.smatrix is used, these wave functions are discarded since keeping them may require an excessive amount of memory for large systems. Instead, when it is needed, the wave function can be obtained using kwant.wave_function. This routine prepares another function that must be supplied the lead number as the only argument to finally get an array of the scattering wave functions corresponding to the incoming modes of that lead. The following code uses this mechanism to calculate the local density of all the states incoming from a given lead. It returns an array that contains a density value for each site of the system. Such an array can be color-plotted with the function map of Kwant's plotter module:

B.3. Exact diagonalization of a finite system Hamiltonian
One of the design objectives of Kwant was to make it easy to perform arbitrary user-specified computations with complex tight-binding systems. When the included solvers do not provide the desired calculation it can be often implemented in just a few lines of Python. (In fact, this is nothing else than writing a solver.) As an example, we apply the ARPACK library [19] to solve a large scale eigenvalue problem and calculate the lowest eigenenergies of a finite system. ARPACK is easily accessible from Python thanks to SciPy. Kwant's hamiltonian_submatrix can return the Hamiltonian matrix in a sparse format understood by SciPy, so there remains almost nothing to be done: The output of this code is shown in Fig. 15.
We have demonstrated that the use of the high-level dynamic language Python for the interface of a quantum transport library can offer considerable advantages in terms of usability. Since the most straightforward way to create a library with a Python interface is to write it in Python and also due to the expressiveness of this language, we have striven to not only provide a Python interface, but to also use Python as much as possible for the implementation of the library itself. That, however, brings in the risk of drastically decreased performance compared to the compiled languages used traditionally for similar tasks -a carelessly written Python program can be about 100 times slower than its C counterpart. By employing a number of measures that are the focus of this section we have ensured that the performance of Kwant remains competitive with other quantum transport codes. In fact, due to a novel approach to solving the scattering problem (see Appendix B and Ref. [15]) Kwant can be more than an order of magnitude faster than traditional quantum transport codes for large systems, as shown in Sec. 7.

C.1. Resource usage of quantum transport calculations
The individual steps of a quantum transport calculation have running times that scale differently with n, the number of sites in the system. The asymptotically most expensive step is the solving: the execution time scales as O(n 3−2/d ) for a d-dimensional system in the case of the RGF algorithm, for instance. Most other parts of the calculation such as initialization, definition of the system, preparation of solving, and (typical) post-processing of the results, have an asymptotic execution time between O(1) and O(n). Similar considerations are valid for memory usage with the memory cost of the solving step dominating as well.
The parts of a quantum transport code with asymptotically highest cost are typically small in terms of code size compared to the rest. This has the important consequence that one may allow oneself to implement most of the code with less attention to performance without a significant penalty. This is not merely an excuse to work sloppily: Freeing oneself from the constraint to strive for the highest possible computer efficiency allows to focus on other aspects like human-time efficiency both of the users and the authors of the code. This insight is the guiding principle behind the technical choices that were made when implementing Kwant.

C.2. Programming languages used
Kwant is written primarily in Python. While some parts of it are implemented in other languages, all these pieces are held together by Python code. Care has been taken to optimize all the asymptotically most expensive components. For the less critical parts pragmatic choices were made: If some component was measured to be too inefficient for typical work loads, it was optimized gradually until a satisfactory performance was obtained. Overall, all of the following approaches are present in Kwant: 1. Usage of pure Python (often delegating low-level operations to libraries like NumPy and SciPy).
For the quantum transport solvers, a combination of approaches 1 and 2 is used. There are currently two solvers, both based on solving the sparse linear system defined by Eqs. (4,8): One utilizes the libraries UMFPACK [40,41] or SuperLU [42,43] that are wrapped by SciPy. The other uses the more efficient MUMPS library [17,18] that, however, is not included in SciPy and needs to be installed separately. The MUMPS-based solver makes use of the nested dissection orderings provided by Metis [44] and Scotch [45].
The part of Kwant that calculates modes and self-energies of leads employs the same approach: Pure Python code uses the services of the LAPACK library [46] that was specially wrapped with Cython as NumPy and SciPy do not provide all the necessary LAPACK functions.
Because highly optimized libraries perform the "heavy lifting" in these asymptotically most expensive parts of Kwant, the fraction of total running time spent in these libraries slowly approaches 100% for all of Kwant as system size grows. In practice more than half of all time is spent in such libraries for all but small systems which means that even if all of Kwant would have been implemented in highly optimized C or Fortran, the additional speed-up one could hope for would be less than a factor of two. Given that this would be much more work, and that it would prevent one from profiting from the dynamic features of Python, we believe that such optimization would not be reasonable.
The routine hamiltonian_submatrix that creates a representation of a (sub)matrix of the Hamiltonian of a Kwant system is an example of a component that has been optimized using Cython even though its execution time is only linear in terms of the number of sites (for sparse matrices). The fact that this routine is called each time a system is solved for a given set of parameters makes it more performance-relevant than the construction and finalization a system that need to be performed only once for a given geometry.
The most low-level optimization of Kwant is embodied by the module tinyarray that has been implemented in pure C++ and is available as a stand-alone library for Python. Kwant uses many instances of small vectors and matrices that cannot be collected into larger arrays, the most frequently used example being the tags and coordinates of lattice sites. These small vectors and matrices could be represented by NumPy arrays or by Python tuples but both solutions are unsatisfactory: NumPy arrays were not designed for this application and are hence too resource-hungry when used in great numbers. Additionally, they cannot be used as keys for Python dictionaries. Python tuples, on the other hand, do not provide mathematical operations and have a high memory-overhead as well, since each element of the tuple exists as an individual Python object. The tinyarray library offers a solution by providing arrays that unlike those of NumPy are optimized for "tiny" sizes. These arrays can be used as dictionary keys as they are hashable and immutable.  This example shows (see Fig. 15), as a function of magnetic field strength, the lowest eigenenergies of a system that supports Majorana fermions. Parts of it are discussed in Appendix A.3.2 and Appendix B.3. from matplotlib import pyplot import kwant import numpy import tinyarray import scipy.sparse.linalg # Python >= 3.3 provides SimpleNamespace in the # standard library so we can simply import it: # >>> from types import SimpleNamespace # (Kwant does not yet support Python 3.) class SimpleNamespace(object): """A simple container for parameters.""" def __init__(self, **kwargs): [2] D. J. Thouless and S. Kirkpatrick, J. Phys. C 14, 235 (1981).