class tkwant.manybody.State(syst, tmax=None, occupations=None, params=None, spectra=None, boundaries=None, intervals=<class 'tkwant.manybody.Interval'>, refine=True, combine=False, error_op=None, scattering_state_type=<class 'tkwant.onebody.onebody.ScatteringStates'>, manybody_wavefunction_type=<class 'tkwant.manybody.WaveFunction'>, mpi_distribute=<function round_robin>, comm=None)[source]

Solve the time-dependent many-particle Schrödinger equation.

  • syst (kwant.builder.FiniteSystem or tkwant.system.ExtendedSystem) – The low level system for which the wave functions are to be calculated.

  • tmax (float, optional) – The maximum time up to which to simulate. Sets the boundary conditions such that they are accurate up to tmax. Must be set if boundaries are not provided. Mutually exclusive with boundaries.

  • occupations (tkwant.manybody.Occupation or sequence thereof, optional) –

    Lead occupation. By default (or if occupations is set to False), all leads are taken into account and are considered as equally occupied with: chemical potential \(\mu = 0\), temperature \(T = 0\) and the non-interacting Fermi-Dirac distribution as distribution function \(f(E)\). To change the default values (\(\mu, T, f(E)\)), a tkwant.manybody.Occupation instance is precalculated with tkwant.manybody.lead_occupation and passed as occupations argument. If occupations is only one element, respectively a sequence with only one element, (\(\mu, T, f(E)\)) will be identical in each lead. In the most general case, if (\(\mu, T, f(E)\)) is different for each lead, occupations must be a sequence of tkwant.manybody.Occupation instances with an ordering similar to syst.leads. In that case, occupations must have the same length as syst.leads, respectively spectra. A lead is not considered, if the corresponding occupations element is set to False. Otherwise, for a lead to be considered, an element of the occupations sequence must have at least the following attributes:

    • energy_range : energy integration range

    • bands : int or list of int, bands (n) to be considered, all bands considered if None

    • distribution : callable, distribution function. Calling signature: (energy).

  • params (dict, optional) – Extra arguments to pass to the Hamiltonian of syst, excluding time.

  • spectra (sequence of spectrum, optional) – Energy dispersion \(E_n(k)\) for the leads. Must have the same length as syst.leads. If needed but not present, it will be calculated on the fly from syst.leads.

  • boundaries (sequence of BoundaryBase, optional) – The boundary conditions for each lead attached to syst. Must have the same length as syst.leads. Mutually exclusive with tmax.

  • intervals (tkwant.manybody.Interval sequence or class, optional) –

    Momentum intervals and quadrature rules on these intervals. If intervals is a sequence, it represents the momentum intervals. In that case, initial integration intervals are not calculated from occupations but intervals is used instead. Each element of the intervals sequence must have at least the following attributes:

    • lead : int, lead index

    • band : int, band index (n)

    • kmin : float, lower momentum bound

    • kmax : float, upper momentum bound, must be larger kmin

    • integration_variable : string, energy vs. momentum integration

    • order : int, quadrature order

    • quadrature : string, quadrature rule to use. See tkwant.integration.calc_abscissas_and_weights

    If intervals is a (data) class, it is passed to the tkwant.manybody.calc_intervals routine as Interval argument. By default, intervals are calculated from tkwant.manybody.calc_intervals and tkwant.manybody.Interval.

  • refine (bool, optional) – If True, intervals are refined at the initial time.

  • combine (bool, optional) – If True, intervals are grouped by lead indices.

  • error_op (callable or kwant.operator, optional) – Observable used for the quadrature error estimate. Must have the calling signature of kwant.operator. Default: Error estimate with density expectation value.

  • scattering_state_type (tkwant.onebody.ScatteringStates, optional) – Class to calculate time-dependent onebody wavefunctions starting in an equilibrium scattering state. Name of the time argument and initial time are taken from this class. If this is not possible, default values are used as a fallback.

  • manybody_wavefunction_type (tkwant.manybody.WaveFunction, optional) – Class to evolve a many-particle wavefunction in time.

  • mpi_distribute (callable, optional) – Function to distribute the tasks dict keys over all MPI ranks. By default, keys must be integer and are distributed round-robin like.

  • comm (mpi4py.MPI.Intracomm, optional) – The MPI communicator over which to parallelize the computation.


The name of the time argument (time_name) and the initial time of the evolution (time_start) are taken from the default values of the scattering_state_type.__init__ method. Changing the default values by partial prebind (e.g. via functools) is possible.


add_boundstate(psi, energy, weight=1)[source]

Add a boundstate to the manybody solver

  • psi (ndarray) – Projection of the boundstate wavefunction to the central scattering region

  • energy (float) – Energy of the boundstate

  • weight (float, optional) – Weight of the boundstate in the discrete sum. Set to one by default.

calc_integrand(operator, interval, root=0)[source]

Return the integrand and the corresponding abscissa values on an interval.

This routine is mainly for debugging purpose.

  • operator (callable or kwant.operator) – Observable to be calculated

  • interval (tkwant.manybody.Interval) – Interval over which the manybody sum should be calculated. Must be present in the solver.

  • root (int or None, optional) – root receive the result, other rank receive None. If root is None all ranks receive the result.


  • abscissa (ndarray) – The integration points (abcissa) of the quadrature. Corresponds to energy values if integration is performed over energy (interval.integration_variable == ‘energy’) or momentum (if interval.integration_variable == ‘momentum’).

  • integrand (ndarray) – The expectation value of operator on all abcissa points.

estimate_error(intervals=None, error_op=None, estimate=None, full_output=False)[source]

Estimate the numerical error of the integration quadrature.

  • intervals (tkwant.manybody.Interval or sequence thereof, optional) – If present, the error is estimated on given intervals \(I_n\), otherwise the total error over all integrals is evaluated.

  • error_op (callable or kwant.operator, optional) – Observable used for the quadrature error estimate. Must have the calling signature of kwant.operator. Default: error_op from initialization.

  • estimate (callable, optional) – Function to estimate an error on an interval \(I_n\). Default: _error_estimate_quadpack, which estimates abserr of QUADPACK (see refine_intervals())

  • full_output (bool, optional) – If the expectation value of error_op is an array and if full_output is true, then the error is estimated on each point of the array. By default, we only return the maximum value of the array error.


error – Error estimate of the momentum (energy) quadrature.

  • Evaluating this function without an argument for interval is a measure of the total error and evaluates error = \(max(\sum_n^N I_n)\), where the sum runs over all N interval errors \(I_n\) stored in the solver.

  • Evaluating this function on a sequence of intervals passed as argument via interval, error is a sequence and the elements are \(max(I_n)\). (first index for intervals) The individual interval errors do NOT sum up to the total error mentioned above. This is due to the triangle inequality \(max(\sum_n^N I_n) \leq \sum_n^N max(I_n)\).

  • If full_output is true, no max() element is taken but the full array-like error of the observable is returned.

  • By default, the error \(I_n\) is identical to abserr of QUADPACK, see refine_intervals().

Return type

ndarray of floats


An interval that is passed as argument via intervals must already be present in the solver, othewise a KeyError is raised.

evaluate(observable, root=0)[source]

Evaluate the expectation value of an operator at the current time.

  • observable (callable or kwant.operator) – An operator to evaluate the expectation value. Must have the calling signature of kwant.operator.

  • root (int or None, optional) – MPI return rank on which to return the result. If root is an integer, it must be in the range of valid MPI ranks 0 <= root < self.comm.size. In that case, the calculated result is returned only on that specific MPI rank where rank == root, whereas the result is None on all other MPI ranks with rank != root. Alternatively, if root is None, the calculated result is returned on all MPI ranks. Setting root to None is generally slower. By default, the result is returned on MPI rank zero only.


result – The expectation value of observable, integrated over all occupied bands. For Kronrod-like quadratures, the expectation value corresponding to the higher order rule is returned by default. Set the instance variable return_element to None or call the evaluate method of the base class directly, in order to get both, the higher and the lower order result. The result might not be returned on all MPI ranks; note the explanation above for input parameter root.

Return type



Returning the result on all MPI ranks (by setting root=None), might be slower, as an additional broadcast step is needed.


Evolve all wavefunctions up to time.


time (int or float) – time argument up to which the solver should be evolved


Return a list of all intervals stored in the solver.


intervals – List of all momentum intervals stored in the solver.

Return type


refine_intervals(atol=1e-05, rtol=1e-05, limit=2000, error_op=None, intervals=None)[source]

Refine intervals until the quadrature error is below tolerance.

  • atol (float, optional) – Absolute accuracy requested.

  • rtol (float, optional) – Relative accuracy requested.

  • limit (integer, optional) – Maximum number of intervals stored in the solver. A warning is raised and the refinement stops if limit is reached.

  • error_op (callable or kwant.operator, optional) – Observable used for the quadrature error estimate. Must have the calling signature of kwant.operator. Default: error_op from initialization.

  • intervals (sequence of tkwant.manybody.Interval, optional) – Apply the refinement process only to the (Gauss-Kronrod) intervals given in the sequence. Note that in this case, all intervals must be present in the solver already. By default, the refinement is done on all Gauss-Kronrod intervals stored in the solver.


  • abserr (float) – Estimate of the modulus of the absolute error, which should equal or exceed abs(i-result), where i is the exact integral value. If error_op has an array-like output, we report the maximal value of the error. (sum errors over all intervals and take the maximum element).

  • intervals (list of tkwant.manybody.Interval) – All subintervals J taken into accound. Intervals are ordered according to errors.

  • errors (ndarray of floats) – Error estimates E(J) on the intervals in descending order. If error_op has an array-like output, the error is returned on all array points. The shape of errors is like error_op (its expectation value) with an additional first dimension for the interval index.


This routine implements a globally adaptive strategy based on quadpacks QAG algorithm 1. It attemps to reduce the absolute error abserr by subdividing (bisecting) the interval with the largest error estimate.

result corresponds to the quadrature estimate of the integral \(\int_a^b f(x) dx\):

\(result = \sum_{J \in \mathcal{P}[a, b]} Q_f(J)\).

Here J is a subinterval of the [a, b] interval and \(Q_f(J)\) is a quadrature rule applied on function f(x) on interval J. In this algorithm, a (2*n + 1) Kronrod estimate is used to calculate Q.

abserr corresponds to the sum of errors:

\(abserr = \sum_{J \in \mathcal{P}[a, b]} E_f(J)\).

We use the quadpack error estimate, described in method _error_estimate_quadpack, for \(E_f(J)\). If the refinement procedure is successful, the following inequality holds:

\(abserr \leq \max \left\{ atol, rtol \cdot result \right\}\)

In the physical sense, the integral we estimate is the manybody integral

\(\langle \hat{A}_{ij}(t) \rangle = \sum_{\alpha} \int_{- \pi}^{\pi} \frac{dk}{2 \pi} v_{\alpha}(k) \theta(v_{\alpha}(k)) f_\alpha(k) [\psi_{\alpha, k}(t)]_i \hat{A} [\psi^\dagger_{\alpha, k}(t)]_j\)

where \(\hat{A}\) corresponds to the error_op and the error is estimated for the expectation value \(\langle \hat{A}_{ij}(t) \rangle\). Note that above inequality condition must be fulfilled on each site i and j individually. This is the case if error_op generates an array-like output. However, the inequality condition must be fulfilled only at the current time t of the solver.

Moreover, note that only intervals stored in the solver are refined. One-body states, that are not part of an interval, are not altered by this method.


Piessens, R., de Doncker-Kapenga, E., Ueberhuber, C. W., and Kahaner, D. K., QUADPACK A Subroutine Package for Automatic Integration, Springer-Verlag, Berlin, (1983).