Commit be30ba86 authored by Genevieve Fleury's avatar Genevieve Fleury
Browse files

Modification of the tutorial on energy transport

parent 0aa60c09
Pipeline #26680 failed with stages
in 12 minutes and 23 seconds
Energy and heat transport
=========================
The aim of this tutorial is to introduce gauge invariant energy related quantities (following Ref `[1] <#references>`__) and how to calculate them using tKwant.
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] <#references>`__.
Theory
------
In the following, we take :math:`e=\hbar=1`.
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).
In the presence of time-dependent electromagnetic fields (described by an electromagnetic scalar potential :math:`V(\vec r, t)` and an electromagnetic vector potential :math:`\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 :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:
The gauge invariant energy operator :math:`\hat ε` is defined as being the kinetic energy operator plus any eventual stationary and physical scalar potential that is present from the remote past. The local energy density operator :math:`\hat ε_i` on site :math:`i` reads:
.. math::
\hat ε_i = \frac{1}{2} \sum_j ε_{ij} \hat c^\dagger_i \hat c_j + ε_{ji} \hat c^\dagger_j \hat c_i
where its matrix elements :math:`ε_{ij}` are derived from the hamiltonian's matrix elements :math:`H_{ij}`:
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.
* 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 before the application of the time-dependent electromagnetic field (i.e. for times :math:`t \leq t_0`) :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:
.. math::
ρ^E_i = \left \langle \hat ε_i \right \rangle = \sum_{\text{lead } α} \; \sum_{{\text{mode } m_α}} \frac{{\text{d}ε}}{2π} f_α(E) \sum_{j} \text{Re} \left [ \left [ ψ_k^{m_α ε} \right ]^\dagger ε_{ij} ψ_j^{m_α ε} \right ]
ρ^ε_i = \left \langle \hat ε_i \right \rangle = \sum_{\text{lead } α} \; \sum_{{\text{mode } m_α}} \frac{{\text{d}E}}{2π} f_α(E) \sum_{j} \text{Re} \left [ \left [ ψ_i^{m_α E} \right ]^\dagger ε_{ij} ψ_j^{m_α E} \right ].
The equation of motion of this quantity
......@@ -31,38 +35,41 @@ The equation of motion of this quantity
\frac{\text{d}}{\text{d}t} \langle \hat ε_i \rangle = \text{i} \left \langle \left [ \hat H, \hat ε_i \right ] \right \rangle + \langle ∂_t \hat ε_i \rangle
Can be interpreted as an energy conservation equation with a source term:
can be interpreted as an energy conservation equation with a source term:
.. math::
\frac{\text{d}}{\text{d}t} ρ^E_i + \sum_j I^E_{ji} = S^E_i
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.
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:
.. math::
I_{ji}^E(t) = \sum_{\text{lead } α} \; \sum_{{\text{mode } m_α}} \frac{{\text{d}ε}}{h} f_α(ε) \sum_{k} - \text{Im} \left [ \left [ ψ_k^{m_α ε} \right ]^\dagger ε_{ki} ε_{ij} ψ_j^{m_α ε} - \left [ ψ_k^{m_α ε} \right ]^\dagger ε_{kj} ε_{ji} ψ_i^{m_α ε} \right ]
I_{ji}^ε(t) = \sum_{\text{lead } α} \; \sum_{{\text{mode } m_α}} \frac{{\text{d}E}}{2\pi} f_α(E) \sum_{k} - \text{Im} \left [ \left [ ψ_k^{m_α E} \right ]^\dagger ε_{ki} ε_{ij} ψ_j^{m_α E} - \left [ ψ_k^{m_α E} \right ]^\dagger ε_{kj} ε_{ji} ψ_i^{m_α E} \right ]
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.
It verifies :math:`I_{ji}^ε(t) = - I_{ij}^ε(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}^ε(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 expression of the energy source :math:`S_i^ε(t)` given to particles on site :math:`i` is as follows:
.. math::
S_i^E(t) = \sum_{\text{lead } α} \; \sum_{{\text{mode } m_α}} \frac{{\text{d}ε}}{h} f_α(E) \sum_{k} - \text{Im} \left [ \left [ ψ_i^{m_α ε} \right ]^\dagger V'_i ε_{ik} ψ_k^{m_α ε} + \left [ ψ_k^{m_α ε} \right ]^\dagger V'_k ε_{ki} ψ_i^{m_α ε} \right ] + \text{Re} \left [ \left [ ψ_i^{m_α ε} \right ]^\dagger ∂_t ε_{ik} ψ_k^{m_α ε} \right ]
S_i^ε(t) = \sum_{\text{lead } α} \; \sum_{{\text{mode } m_α}} \frac{{\text{d}E}}{2\pi} f_α(E) \sum_{k} - \text{Im} \left [ \left [ ψ_i^{m_α E} \right ]^\dagger V'_i ε_{ik} ψ_k^{m_α E} + \left [ ψ_k^{m_α E} \right ]^\dagger V'_k ε_{ki} ψ_i^{m_α E} \right ] + \text{Re} \left [ \left [ ψ_i^{m_α E} \right ]^\dagger ∂_t ε_{ik} ψ_k^{m_α E} \right ]
The heat current leaving a lead :math:`α` to the central system is defined as:
where :math:`V_i(t) = H_{ii}(t) - ε_{ii}(t)` is the time-dependent scalar potential on site :math:`i`.
.. math::
I_α^H(t) = I^E_α - S^E_α - μ I^N_α
The heat current leaving a lead :math:`α` to the central system :math:`C` is defined as:
.. math::
I_α^H(t) = I^ε_α - S^ε_α - μ_α I^N_α
Calculating energy transport with tKwant
-----------------------------------------
where :math:`I^ε_α=\sum_{i \in α}\sum_{j \in C}I^ε_{ji}`, :math:`I^N_α=\sum_{i \in α}\sum_{j \in C}I^N_{ji}`, :math:`S^ε_α=\sum_{i \in α}S^ε_i`, and :math:`μ_α` is chemical potential of lead :math:`α`.
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`:
Example : Quantum dot embedded in a 1D chain
--------------------------------------------
A central site 0 with on-site energy :math:`ε_0(t)=ε_0+Δε\Theta(t-t_0)` is connected through nearest-neighbor hopping term :math:`γ_c` to two left and right semi-infinite 1D chains. We start by defining the central part of the system as follows
.. jupyter-execute::
......@@ -70,8 +77,8 @@ Let's illustrate the use of these classes in a simple system, the quantum dot. M
import kwant
import tkwant
γ = 1.0
γc = 0.2
γ = 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
......@@ -79,7 +86,7 @@ Let's illustrate the use of these classes in a simple system, the quantum dot. M
ε0 = 0.5 * Γ # initial dot energy level
Δε = 2.5 * Γ # change of dot energy level
t0 = 0.0 # time of the heaviside jump
t0 = 0.5 # time of the heaviside jump
def qdot_potential(site, time, t0, ε0, Δε):
if time > t0:
......@@ -90,13 +97,13 @@ Let's illustrate the use of these classes in a simple system, the quantum dot. M
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
# 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
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.
We want to calculate the time dependent energy density on the site ``lat(0)`` (the site on which the time dependent potential has been applied), the left heat and energy currents flowing through the hopping ``(lat(-1), lat(0))`` (from ``lat(-1)`` to ``lat(0)``) and the right heat and energy currents flowing through ``(lat(1), lat(0))`` (from ``lat(1)`` to ``lat(0)``). For that matter, two (left and right) leads 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 and heat currents since e.g. all four sites ``lat(-2)``, ``lat(-1)``, ``lat(0)``, and ``lat(1)``are needed to calculate the energy current :math:`I^ε_{0,-1}(t)`.
.. jupyter-execute::
......@@ -106,12 +113,12 @@ We want to calculate the time dependent energy density on the site ``lat(1)`` (t
lead[lat.neighbors()] = - γ
added_sites_left = builder.attach_lead(lead, add_cells=1)
# Append lat(0) to the list, to calculate the heat current between lat(0) and lat(1)
added_sites_left.append(lat(0))
# Append lat(-1) to the list, to calculate the heat current between lat(-1) and lat(0)
added_sites_left.append(lat(-1))
added_sites_right = builder.attach_lead(lead.reversed(), add_cells=1)
# Append lat(0) to the list, to calculate the heat current between lat(2) and lat(1)
added_sites_right.append(lat(2))
# 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
......@@ -135,7 +142,7 @@ The system is thus the following
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:
The lists of sites ``added_sites_left`` and ``added_sites_right`` are used when declaring `~tkwant.operator.energy.LeadHeatCurrent` operator instances by giving them to the ``added_lead_sites`` parameter:
.. jupyter-execute::
......@@ -149,21 +156,21 @@ The system is thus the following
lead_heat_current_left_op = tkwant.operator.energy.LeadHeatCurrent(syst, chemical_potential=μL, added_lead_sites=added_sites_left)
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:
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(0), lat(-1))`` for the left lead and ``(lat(0), lat(1))`` for the right lead:
.. jupyter-execute::
print(lead_heat_current_left_op.hoppings == [(lat(1), lat(0))])
print(lead_heat_current_right_op.hoppings == [(lat(1), lat(2))])
print(lead_heat_current_left_op.hoppings == [(lat(0), lat(-1))])
print(lead_heat_current_right_op.hoppings == [(lat(0), lat(1))])
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.
.. jupyter-execute::
energy_current_op = tkwant.operator.energy.Current(syst, where=[(lat(0), lat(1)), (lat(2), lat(1))])
energy_source_op = tkwant.operator.energy.Source(syst, where=[lat(1)])
energy_density_op = tkwant.operator.energy.Density(syst, where=[lat(1)])
energy_current_divergence_op = tkwant.operator.energy.CurrentDivergence(syst, where=[lat(1)])
particle_current_op = kwant.operator.Current(syst, where=[(lat(0), lat(-1)), (lat(0), lat(1))])
energy_current_op = tkwant.operator.energy.Current(syst, where=[(lat(0), lat(-1)), (lat(0), lat(1))])
energy_source_op = tkwant.operator.energy.Source(syst, where=[lat(-1),lat(1)])
energy_density_op = tkwant.operator.energy.Density(syst, where=[lat(0)])
We can now initialize the solver
......@@ -187,13 +194,13 @@ And run the simulation. Energy operators are evaluated in the same way as the pa
.. jupyter-execute::
# Have the system evolve forward in time, calculating the operator over the system
# 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 = []
energy_current_divergence = []
dt = 0.01 / Γ # Interval between two measurements of the operators
......@@ -210,7 +217,6 @@ And run the simulation. Energy operators are evaluated in the same way as the pa
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))
energy_current_divergence.append(solver.evaluate(energy_current_divergence_op))
Which gives the following plots:
......@@ -234,24 +240,28 @@ Which gives the following plots:
energy_source = np.array(energy_source) / Γ**2
heat_current_left = np.array(heat_current_left) / Γ**2
heat_current_right = np.array(heat_current_right) / Γ**2
heat_current_left_2 = energy_current[:, 0] - energy_source[:, 0] - μL * particle_current[:, 0] / Γ**2
heat_current_right_2 = energy_current[:, 1] - energy_source[:, 1] - μR * particle_current[:, 1] / Γ**2
energy_density = np.array(energy_density) / Γ
energy_current_divergence = np.array(energy_current_divergence) / Γ**2
# Plot
plt.plot(times, energy_current[:, 0], label = "Left outgoing current $I^E_{0,1}$")
plt.plot(times, energy_current[:, 1], label = "Right outgoing current $I^E_{2,1}$")
plt.plot(times, energy_source, label = "Source $S_1$")
plt.plot(times, energy_current[:, 0], label = "Left incoming energy current $I^E_{0,-1}$", 'r-')
plt.plot(times, energy_current[:, 1], label = "Right incoming energy current $I^E_{0,1}$", 'm-')
plt.plot(times, energy_source[:, 0], label = "Left source $S_{-1}$", 'g-')
plt.plot(times, energy_source[:, 1], label = "Right source $S_{1}$", 'y-')
plt.plot(times, heat_current_left, label = "Left lead outgoing heat current")
plt.plot(times, heat_current_right, label = "Right lead outgoing heat current")
plt.plot(times, heat_current_left, label = "Left incoming heat current", 'c-')
plt.plot(times, heat_current_right, label = "Right incoming heat current", 'b-')
plt.plot(times, heat_current_left2, label = "$I^E_{0,-1}-S_{-1}-μL*I^N_{0,-1}$", 'co')
plt.plot(times, heat_current_right2, label = "$I^E_{0,1}-S_{1}-μL*I^N_{0,1}$", 'bo')
plt.xlabel('time $[\\hbar / \\Gamma]$')
plt.ylabel('Energy flux $[\\Gamma^2 / \\hbar]$')
plt.legend()
plt.show()
plt.plot(times, energy_density, label = "Energy density $\\varepsilon_1$")
plt.plot(times, energy_density, label = "Energy density $\\varepsilon_0$")
plt.xlabel('time $[\\hbar / \\Gamma]$')
plt.ylabel('Energy density $[\\Gamma]$')
plt.legend()
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment