The aim of this tutorial is to show how to simulate energy and heat transport with tKwant. The corresponding theoretical framework can be found in Ref [1].
Readers who are not familiar with t-Kwant are invited to visit its documentation webpage first.
In the following, we take \(e=\hbar=1\).
In the presence of time-dependent electromagnetic fields (described by an electromagnetic scalar potential \(V(\vec r, t)\) and an electromagnetic vector potential \(\vec A(\vec r, t)\)), the Hamiltonian’s expectation value is (in general) gauge dependent i.e. an electromagnetic gauge transformation does not leave its expectation value invariant.
The gauge invariant energy operator \(\hat ε\) is defined as being the kinetic energy operator plus any stationary scalar potential, that has a physical origin and is present from the remote past to the remote future. The local energy density operator \(\hat ε_i\) on site \(i\) reads:
where its matrix elements \(ε_{ij}\) are derived from the Hamiltonian’s matrix elements \(H_{ij}\):
Values on hoppings coincide with the Hamiltonian’s \(\forall i \! \neq \! j \; ε_{ij}(t) = H_{ij}(t)\)
Values on sites are the Hamiltonian’s before the application of the time-dependent electromagnetic field (i.e. for times \(t \leq t_0\)) \(\forall i \, \forall t \; ~ ε_{ii}(t) = H_{ii}(t_0)\). This captures the initial scalar potential with the assumption that it physically exists all along the experiment, the time dependent control stacking on top of it.
The manybody expectation value of the local energy density operator, in the wavefunction formalism using scattering states, reads:
The equation of motion of this quantity
can be interpreted as an energy conservation equation with a source term:
where \(I^ε_{ji}\) is the energy current flowing from site \(i\) to site \(j\) and \(S^ε_i\) is the power given to the charged particles by the time dependent electromagnetic fields that are embedded in the Hamiltonian. Note that the expression taken for the energy current \(I^ε_{ji}\) is not unique and only its divergence has a physical meaning: the energy flux out of a closed surface. A similar situation exists for the Poynting vector in classical electrodynamics.
The taken expression for the energy current \(I_{ji}^ε(t)\) from site \(i\) to site \(j\) is the following:
It verifies \(I_{ji}^ε(t) = - I_{ij}^ε(t)\), necessary for interpreting the quantity as the net current flowing from site \(i\) to site \(j\) ; and \(H_{ij} = 0 \implies I_{ji}^ε(t) = 0\), which means that no energy current flows between disconnected sites.
The expression of the energy source \(S_i^ε(t)\) given to particles on site \(i\) is as follows:
where \(V_i(t) = H_{ii}(t) - ε_{ii}(t)\) is the time-dependent scalar potential on site \(i\).
The heat current leaving a lead \(α\) to the central system \(C\) is defined as:
where \(I^ε_α=\sum_{i \in α}\sum_{j \in C}I^ε_{ji}\), \(I^N_α=\sum_{i \in α}\sum_{j \in C}I^N_{ji}\), \(S^ε_α=\sum_{i \in α}S^ε_i\), and \(μ_α\) is chemical potential of lead \(α\).
Each one of the previously defined energy quantities have been implemented tkwantoperator module and are to be used with tkwant. The classes are: EnergyDensity, EnergyCurrent, EnergySource and LeadHeatCurrent. An additional class EnergyCurrentDivergence exists to calculate easily the outgoing energy flux off a set of sites. These classes enable using \(\hat ε\), \(\hat H\) or a user defined operator (with specified values on sites) as energy operators, more information is available in the class’ documentation.
A 1D chain, whose sites have zero on-site energy and are connected with each other through nearest neighbor hoppings \(γ\), has a central site \(0\) with an on-site energy \(ε_0(t)=ε_0+Δε\Theta(t-t_0)\). The central site is connected to its two nearest neighbors through a different hopping term \(γ_c\). We want to calculate the time dependent energy density on the central site, the energy and heat currents flowing to its right and left side.
In kwant, such a system can be created as follows:
from mpi4py import MPI
import kwant
import tkwant
import tkwantoperator
γ = 1.0 # nearest-neighbor hopping term in the leads
γc = 0.2 # nearest-neighbor hopping term between the central site and the leads
a = 1.0 # lattice constant
Γ = 4 * γc*γc / γ # Energy scaling unit
ε0 = 0.5 * Γ # initial dot energy level
Δε = 2.5 * Γ # change of dot energy level
t0 = 1.5 # time of the heaviside jump
def qdot_potential(site, time, t0, ε0, Δε):
if time > t0:
return ε0 + Δε
else:
return ε0
lat = kwant.lattice.chain(a, norbs = 1)
builder = kwant.Builder()
# lat(0) is the dot, the rest belongs already formally to the leads
builder[(lat(-1))] = 0
builder[(lat(0))] = qdot_potential
builder[(lat(1))] = 0
builder[lat.neighbors()] = - γc
The left currents will be computed on the hopping (lat(0), lat(-1)) (from lat(-1) to lat(0)) and the right currents on (lat(0), lat(1)) (from lat(1) to lat(0)). We will create and attach each lead to the system while adding its edge site to the system, by using the add_cells parameter in the attach_lead method. Given that wave functions are not accessible in the leads in kwant and tkwant, adding one site from each lead to the system needs to be done to be able to calculate properly the energy and heat currents, since the wave functions also need to be accessed on lat(-2) (respectively lat(-2)) to compute \(I^ε_{0,-1}(t)\) (respectively \(I^ε_{0,1}(t)\)).
# Define and attach the leads
lead = kwant.Builder(kwant.TranslationalSymmetry((-a,)))
lead[lat(0)] = 0
lead[lat.neighbors()] = - γ
# Attach the left lead and add the first cell of the lead
# (i.e. the site lat(-2)) to the central system. The added site
# is returned by `attach_lead()` and stored in `added_sites_left`
added_sites_left = builder.attach_lead(lead, add_cells=1)
# Append lat(-1) to the list, to calculate the heat current between lat(-1) and lat(0)
added_sites_left.append(lat(-1))
# Attach the right lead and add the first cell of the lead
# (i.e. the site lat(2)) to the central system. The added site
# is returned by `attach_lead()` and stored in `added_sites_right`
added_sites_right = builder.attach_lead(lead.reversed(), add_cells=1)
# Append lat(1) to the list, to calculate the heat current between lat(1) and lat(0)
added_sites_right.append(lat(1))
The system is thus the following
kwant.plot(builder);
The added sites from each lead are stored in the added_sites_left and added_sites_right variables. These variables are then used to instance LeadHeatCurrent by giving them to the added_lead_sites parameter:
# Create finalized system
syst = builder.finalized()
μL = 0.5 * Γ # Chemical potential in the left lead
μR = -0.5 * Γ # Chemical potential in the right lead
lead_heat_current_right_op = tkwantoperator.LeadHeatCurrent(syst, chemical_potential=μR, added_lead_sites=added_sites_right)
lead_heat_current_left_op = tkwantoperator.LeadHeatCurrent(syst, chemical_potential=μL, added_lead_sites=added_sites_left)
The added_lead_sites parameter is used by the LeadHeatCurrent class to find the hopping(s) on which it will calculate the heat current: they are the outgoing hoppings from the set of sites added_lead_sites to the system, excluding the hoppings that connect to lead sites. The calculated hoppings are stored in its hoppings member. This approach (attach_lead with add_cells=1 then giving the added cells to LeadHeatCurrent) has the benefit to make it easier for the user to calculate lead currents when there are many hoppings connecting the central system to the leads.
We can check, in our specific case, that each LeadHeatCurrent instance computed the right hoppings: they should be (lat(0), lat(-1)) for the left lead and (lat(0), lat(1)) for the right lead
print(lead_heat_current_left_op.hoppings == [(lat(0), lat(-1))])
print(lead_heat_current_right_op.hoppings == [(lat(0), lat(1))])
True
True
Now let’s declare the other operators. Just like the particle Current and Density, energy operators receive a where argument upon initialization (but can’t be left blank or set to None because of lead interface issues). For more information on any energy operator, feel free to visit its respective documentation.
particle_current_op = kwant.operator.Current(syst, where=[(lat(0), lat(-1)), (lat(0), lat(1))])
energy_current_op = tkwantoperator.EnergyCurrent(syst, where=[(lat(0), lat(-1)), (lat(0), lat(1))])
energy_source_op = tkwantoperator.EnergySource(syst, where=[lat(-1),lat(1)])
energy_density_op = tkwantoperator.EnergyDensity(syst, where=[lat(0)])
We can now initialize the solver
# Occupation, for each lead
TL = 1.0 * Γ # Temperature in the left lead
TR = 0.0 * Γ # Temperature in the right lead
occupation = [None] * len(syst.leads)
occupation[0] = tkwant.manybody.lead_occupation(chemical_potential=μL, temperature=TL)
occupation[1] = tkwant.manybody.lead_occupation(chemical_potential=μR, temperature=TR)
# Maximum simulation time
tmax = 6. / Γ
# Initialize the solver (to solve t-dep SEQ)
solver = tkwant.manybody.State(syst, tmax, occupation, params={'t0': t0, 'ε0': ε0, 'Δε': Δε}, error_op=energy_density_op)
And run the simulation. Energy operators are evaluated in the same way as the particle operators.
# Have the system evolve forward in time, calculating the operator over the system
particle_current = []
energy_current = []
energy_source = []
energy_density = []
heat_current_left = []
heat_current_right = []
dt = 0.01 / Γ # Interval between two measurements of the operators
import numpy as np
times = np.arange(0, tmax, dt)
for time in times:
# evolve scattering states in time
solver.evolve(time)
solver.refine_intervals()
particle_current.append(solver.evaluate(particle_current_op))
energy_current.append(solver.evaluate(energy_current_op))
energy_source.append(solver.evaluate(energy_source_op))
energy_density.append(solver.evaluate(energy_density_op))
heat_current_left.append(solver.evaluate(lead_heat_current_left_op))
heat_current_right.append(solver.evaluate(lead_heat_current_right_op))
Which gives the following plots:
A python file containing all the above lines of code is available: quantum_dot.py. It can be ran in parallel to shorten the simulation time, as explained in Parallel calculations
[1] A. Kara Slimane, P. Reck, G. Fleury, Simulating time-dependent thermoelectric transport in quantum systems [ Phys. Rev. B | arXiv ]