2.6. Solving the many-body problem

Note

The examples in this section take several minutes on a single core desktop computer. To speed up the computation the example scrips can be run in parallel, see section Parallelization with MPI.

We like to study manybody problem for an infinite one-dimensional chain. In second quantization the Hamiltonian reads

\[\hat{H}(t) = \sum_{i,j} \gamma_{ij} \, \hat{c}^\dagger_i \hat{c}_j + \sum_{i} w(t) \theta(i_b - i) \, \hat{c}^\dagger_i \hat{c}_i\]

where \(\gamma_{ii} = 1\) and \(\gamma_{ij} = -1\) and \(w(t)\) is taken similar to section Solving one-body problems as

\[w(t) = \theta(t) v_p e^{- 2 (t / \tau)^2}.\]

We are interested in the time evolution of the electron density

\[n_i(t) = \langle \hat c^\dagger_i \hat c_i \rangle (t)\]

The time-dependent pulse \(w(t)\) which act on all lattice sites to the left of \(i_b\) is taken into account by a gauge transform similar to Time-dependent potentials and pulses and translates to a time-dependent coupling between \(i_b\) and \(i_b+1\). The Kwant code to build the system is

import kwant
import cmath
from scipy.special import erf

def gaussian(time, t0=40, A=1.57, sigma=24):
    return A * (1 + erf((time - t0) / sigma))

# time dependent coupling with gaussian pulse
def coupling_nn(site1, site2, time):
    return - cmath.exp(- 1j * gaussian(time))

def make_system(L=400):

    # system building
    lat = kwant.lattice.square(a=1, norbs=1)
    syst = kwant.Builder()

    # central scattering region
    syst[(lat(x, 0) for x in range(L))] = 1
    syst[lat.neighbors()] = -1
    # time dependent coupling between two sites in the center
    syst[lat(L // 2, 0), lat(L // 2 - 1, 0)] = coupling_nn

    # add leads
    sym = kwant.TranslationalSymmetry((-1, 0))
    lead_left = kwant.Builder(sym)
    lead_left[lat(0, 0)] = 1
    lead_left[lat.neighbors()] = -1
    syst.attach_lead(lead_left)
    syst.attach_lead(lead_left.reversed())

    return syst

Note that the code to build the system is basically the same as in the previous examples of the onebody problem which were treated in first quantization. The system looks similar to

../_images/manybody_2_1.png

For representation purpose, the central scattering system has been shrinked to only 20 sites in the plot and the time-dependent coupling is highlighed in red.

Two approaches are possible to obtain the density exectation value: Either a high-level approach using manybody.State where the preprocessing is done automatically and which provides additional functionality. Alternatively a low-level approach using manybody.WaveFunction, where the different preprocessing steps must be handled manually. Both ways are shown below.

2.6.1. High-level automatic approach

The high-level approach comprises all preprocessing steps. The entire code is:

import tkwant
import kwant
import matplotlib.pyplot as plt

syst = make_system().finalized()
sites = [site.pos[0] for site in syst.sites]
times = [40, 80, 120, 160]

density_operator = kwant.operator.Density(syst)

state = tkwant.manybody.State(syst, max(times))

density0 = state.evaluate(density_operator)

for time in times:
    state.evolve(time=time)
    if time == 40:
        state.refine_intervals()
    error = state.estimate_error()
    print('time={}, error={:10.4e}'.format(time, error))
    density = state.evaluate(density_operator)
    plt.plot(sites, density - density0, label='time={}'.format(time))

plt.legend()
plt.xlabel(r'site position $i$')
plt.ylabel(r'charge density $n$')
plt.show()

time=40, error=5.1438e-07

time=80, error=5.0540e-04

time=120, error=3.6453e-03

time=160, error=1.1463e-02
../_images/manybody_3_4.png

Note that this approach is much simpler and provides additional methods to fascilitate the numerical procedure without the need to fine-tune the quadrature by hand. While the high-level approach is less flexible, it can still be adapted in various ways. In the following we show how to change the lead occupation.

See also

The complete example script including MPI directives for parallel execution can be found in 1d_wire_high_level.py.

Chemical potential and temperature of the leads

By default, the chemical potential and the temperature in all leads are identical and equal zero. To set them in all leads to the same, non-zero value, is possible via

occupations = tkwant.manybody.lead_occupation(chemical_potential=0.5, temperature=0.1)
state = tkwant.manybody.State(syst, max(times), occupations)

One can also set different values in each lead as

occup_left = tkwant.manybody.lead_occupation(chemical_potential=0.5, temperature=0.1)
occup_right = tkwant.manybody.lead_occupation(chemical_potential=0.2)
occupations = [occup_left, occup_right]

state = tkwant.manybody.State(syst, max(times), occupations)

Instead of a scalar value one can also provide a sequence of values for the chemical_potential, one value for each band of the lead spectrum:

occupations = tkwant.manybody.lead_occupation(chemical_potential=[-0.5, 0.5])

In this case the number of values must match the number of bands in that lead. The ordering of the bands in the same obtained with kwantSpectrum.

Adaptive refinement and error estimate

The class manybody.State provides methods to estimate the quadrature error of the manybody integral and to adaptively refine the approximation to a given accuracy.

The error ist estimated via

error = state.estimate_error()
print('estimated integration error= {:10.4e}'.format(error))

estimated integration error= 1.3478e-08

By default, the error is estimated on the density expectation value. One can obtain the error also for other expectation values, as e.g. the current:

current_operator = kwant.operator.Current(syst)
error = state.estimate_error(error_op=current_operator)
print('estimated integration error= {:10.4e}'.format(error))

estimated integration error= 6.7801e-10

The quadrature intervals can be refined via

state.refine_intervals();

By default, the refinement is done up to a certain accuracy of the density expectation value. Again, the behavior can be changed

current_operator = kwant.operator.Current(syst)
state.refine_intervals(rtol=1E-3, atol=1E-3, error_op=current_operator);

Note

Adaptive refinement is computationally expensive. Exploring initially at low precision is often a good idea.

2.6.2. Low-level manual approach

The low-level approach is close to the algorithm to solve the manybody problem which described in the Tkwant paper. The code is:

from tkwant import leads, manybody
import kwant
import kwantspectrum

import functools
import numpy as np
import matplotlib.pyplot as plt

syst = make_system().finalized()
sites = [site.pos[0] for site in syst.sites]
times = [40, 80, 120, 160]

density_operator = kwant.operator.Density(syst)

# calculate the spectrum E(k) for all leads
spectra = kwantspectrum.spectra(syst.leads)

# estimate the cutoff energy Ecut from T, \mu and f(E)
# All states are effectively empty above E_cut
occupations = manybody.lead_occupation(chemical_potential=0, temperature=0)
emin, emax = manybody.calc_energy_cutoffs(occupations)

# define boundary conditions
bdr = leads.automatic_boundary(spectra, tmax=max(times), emin=emin, emax=emax)

# calculate the k intervals for the quadrature
interval_type = functools.partial(manybody.Interval, order=20,
                                  quadrature='gausslegendre')
intervals = manybody.calc_intervals(spectra, occupations, interval_type)
intervals = manybody.split_intervals(intervals, number_subintervals=10)

# calculate all onebody scattering states at t = 0
tasks = manybody.calc_tasks(intervals, spectra, occupations)
psi_init = manybody.calc_initial_state(syst, tasks, bdr)

# set up the manybody wave function
wave_function = manybody.WaveFunction(psi_init, tasks)

density0 = wave_function.evaluate(density_operator)

for time in times:
     wave_function.evolve(time)
     density = wave_function.evaluate(density_operator)
     plt.plot(sites, density - density0, label='time={}'.format(time))

plt.legend()
plt.xlabel(r'site position $i$')
plt.ylabel(r'charge density $n$')
plt.show()
../_images/manybody_12_0.png

The role of each function can be deduced from the Tkwant paper and the function documentation. While most lines of above code are generic, a few lines are responsible for the numerical accuracy of the result and must be fine tuned for each problem in question.

The numerical accuracy is controled by the integration order (given by the variable order) of a quadrature interval and by the number of sub intervals (by the variable number_subintervals), in which each initial quadrature interval is divided. The actual value of the variable order, is less crucial and typically ranges in between 10 and 20. The value of number_subintervals is highly system dependent and must be tuned.

Note

The numerical precision of the manybody expectation value is mainly determined by the integer variable number_subintervals in above example. Larger values lead to a more precise result on the cost of longer compute time. The actual value is highly system dependent. It is a good practice to start with a low value and to gradually increase it until the result converges.

To better understand the logic between these two parameters, let us state that for Gaussian quadrature rules as Gauss-Legendre or Gauss-Kronrod, the sampling points are not distributed equidistantly over the quadrature interval. The purpose of the function manybody.split_intervals() is to split a quadrature interval with a given order equidistantly into number_subintervals with similar order. From this follows, that order=2 and number_subintervals=10 and order=10 and number_subintervals=2 will both lead to the same number of sampling points to approximate the integral, but with a very different distribution of the points.

See also

The complete example script including MPI directives for parallel execution can be found in 1d_wire_low_level.py

2.6.3. Summary

To summarize, we like to highlight the similarity between the onebody and the manybody approach. The first one is the definition of the system using Kwant, which is the same, whether the Hamiltonian is written in first quantization (onebody) or in second quantization (manybody). The second similarity is the API of the solvers for the onebody and the manybody Schrödinger equation. We will show this on the example of the two classes onebody.ScatteringStates() and manybody.State(). After defining an observable, as e.g.

density_operator = kwant.operator.Density(syst)

both states can be evolved forward in time and evaluate expectation values similarly

Onebody

mode = 0
psi = tkwant.onebody.ScatteringStates(syst, energy=1, lead=0, tmax=10)[mode]
psi.evolve(time=5)
density = psi.evaluate(density_operator)

Manybody

state = tkwant.manybody.State(syst, tmax=10)
state.evolve(time=5)
density = state.evaluate(density_operator)

See also

More customization options for the high- and the low-level approach be found in the section Advanced manybody settings. Additional examples for solving the manybody Schrödinger equation can be found in the section More examples.