The aim of this tutorial is to introduce gauge invariant energy related quantities (following Ref `[1] <#references>`__) and how to calculate them using tKwant.

Theory

------

The need of gauge invariance comes from the existence of time dependent electrodynamic fields. Indeed, these are accounted for by the introduction of the electromagnetic scalar potential :math:`V(\vec r, t)` and vector potential :math:`\vec A(\vec r, t)` in the hamiltonian :math:`\hat H`. However, this makes the hamiltonian's expectation value gauge dependent since an electromagnetic gauge transformation doesn't leave its expectation value invariant (refer to Ref `[1] <#references>`__ for more information).

The gauge invariant energy operator :math:`\hat ε` is defined as being the kinetic energy operator plus any eventual stationary and physical scalar potential. The local energy density operator :math:`\hat ε_i` on site :math:`i` writes:

where its matrix elements :math:`ε_{ij}` are derived from the hamiltonian's matrix elements :math:`H_{ij}`:

* Values on hoppings coincide with the hamiltonian's :math:`\forall i \! \neq \! j \; ε_{ij}(t) = H_{ij}(t)`

* Values on sites are the hamiltonian's at a time :math:`t_0` where it is still stationary: :math:`\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:

Where :math:`I^E_{ji}` is the energy current flowing from site :math:`i` to site :math:`j` and :math:`S^E_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 :math:`I^E_{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 :math:`I_{ji}^E(t)` from site :math:`i` to site :math:`j` is the following:

It verifies :math:`I_{ji}^E(t) = - I_{ij}^E(t)`, necessary for interpreting the quantity as the net current flowing from site :math:`i` to site :math:`j` ; and :math:`H_{ij} = 0 \implies I_{ji}^E(t) = 0`, which means that no energy current flows between disconnected sites.

The expression of the energy source :math:`S_i^E(t)` given to particles on site :math:`i` is as follows:

The heat current leaving a lead :math:`α` to the central system is defined as:

.. math::

I_α^H(t) = I^E_α - S^E_α - μ I^N_α

Calculating energy transport with tKwant

-----------------------------------------

Each one of the previously defined energy quantities have been implemented as classes in tKwant under the `~tkwant.operator.energy` module. Namely `~tkwant.operator.energy.Density`, `~tkwant.operator.energy.Current`, `~tkwant.operator.energy.Source` and `~tkwant.operator.energy.LeadHeatCurrent`. An additional class `~tkwant.operator.energy.CurrentDivergence` exists to calculate easily the outgoing energy flux off a set of sites. These classes enable using :math:`\hat ε`, :math:`\hat H` or a user defined operator (with specified values on sites) as energy operators, more information is available in the class' documentation.

Let's illustrate the use of these classes in a simple system, the quantum dot. Modeled by a 1D chain with a heaviside potential applied on one site, connected with the leads with a coupling constant :math:`γ_c`:

.. jupyter-execute::

from mpi4py import MPI

import kwant

import tkwant

γ = 1.0

γc = 0.2

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 = 0.0 # 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(1) is the dot, the rest belongs already formally to the leads

builder[(lat(0))] = 0

builder[(lat(1))] = qdot_potential

builder[(lat(2))] = 0

builder[lat.neighbors()] = - γc

We want to calculate the time dependent energy density on the site ``lat(1)`` (the site on which the time dependent potential has been applied), the heat and energy current flowing through ``(lat(0), lat(1))`` (from ``lat(1)`` to ``lat(0)``) and ``(lat(2), lat(1))``. For that matter, two leads (left and right) need to be added to the system: they are taken with a zero onsite potential and with a hopping constant :math:`γ`. One site from each lead is added to the system to be able to calculate properly the energy currents since nearest neighbors of a chosen hopping are needed.

.. jupyter-execute::

# Define and attach the leads

lead = kwant.Builder(kwant.TranslationalSymmetry((-a,)))

# Append lat(0) to the list, to calculate the heat current between lat(2) and lat(1)

added_sites_right.append(lat(2))

The system is thus the following

.. jupyter-execute::

:hide-code:

import matplotlib.pyplot as plt

plt.rcParams['text.usetex'] = True

plt.rcParams['figure.figsize'] = (11.69,8.27)

plt.rcParams['axes.formatter.useoffset'] = True

plt.rcParams['axes.formatter.limits'] = (-2, 2)

plt.rcParams['axes.grid'] = True

plt.rcParams['font.size'] = 20

%matplotlib inline

.. jupyter-execute::

kwant.plot(builder)

``added_sites_left`` and ``added_sites_right`` are a list of sites we will use when declaring `~tkwant.operator.energy.LeadHeatCurrent` operator instances by giving them to the ``added_lead_sites`` parameter:

.. jupyter-execute::

# Create finalized system

syst = builder.finalized()

μL = 0.5 * Γ # Chemical potential in the left lead

μR = -0.5 * Γ # Chemical potential in the right lead

The ``added_lead_sites`` parameter is used by the `~tkwant.operator.energy.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``, excluding hoppings that connect to lead sites. The calculated hoppings are stored in its `~tkwant.operator.energy.LeadHeatCurrent.hoppings` member. We can check that these correspond to ``(lat(1), lat(0))`` for the left lead and ``(lat(1), lat(2))`` for the right lead:

Now let's declare the other operators. Just like the particle `~kwant.operator.Current` and `~kwant.operator.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 operators, feel free to visit its respective documentation.