Verified Commit 030f74c1 authored by Adel Kara Slimane's avatar Adel Kara Slimane
Browse files

Fixes in doc/source/tutorial/energy.rst

parent be30ba86
Pipeline #26686 failed with stages
in 12 minutes and 26 seconds
Energy and heat transport
=========================
The aim of this tutorial is to show how to simulate energy and heat transport with 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>`__.
......@@ -27,7 +27,7 @@ The manybody expectation value of the local energy density operator, in the wave
.. math::
ρ^ε_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 ].
ρ^ε_i = \left \langle \hat ε_i \right \rangle = \sum_{\text{lead } α} \; \sum_{{\text{mode } m_α}} \int \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
......@@ -39,21 +39,21 @@ 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
\frac{\text{d}}{\text{d}t} ρ^ε_i + \sum_j I^ε_{ji} = S^ε_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^ε_{ji}` is the energy current flowing from site :math:`i` to site :math:`j` and :math:`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 :math:`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 :math:`I_{ji}^E(t)` from site :math:`i` to site :math:`j` is the following:
.. math::
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 ]
I_{ji}^ε(t) = \sum_{\text{lead } α} \; \sum_{{\text{mode } m_α}} \int \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}^ε(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^ε(t)` given to particles on site :math:`i` is as follows:
.. math::
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 ]
S_i^ε(t) = \sum_{\text{lead } α} \; \sum_{{\text{mode } m_α}} \int \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 ]
where :math:`V_i(t) = H_{ii}(t) - ε_{ii}(t)` is the time-dependent scalar potential on site :math:`i`.
......@@ -85,8 +85,8 @@ A central site 0 with on-site energy :math:`ε_0(t)=ε_0+Δε\Theta(t-t_0)
Γ = 4 * γc*γc / γ # Energy scaling unit
ε0 = 0.5 * Γ # initial dot energy level
Δε = 2.5 * Γ # change of dot energy level
t0 = 0.5 # time of the heaviside jump
Δε = 2.5 * Γ # change of dot energy level
t0 = 0.5 # time of the heaviside jump
def qdot_potential(site, time, t0, ε0, Δε):
if time > t0:
......@@ -95,13 +95,13 @@ A central site 0 with on-site energy :math:`ε_0(t)=ε_0+Δε\Theta(t-t_0)
return ε0
lat = kwant.lattice.chain(a, norbs = 1)
builder = kwant.Builder()
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
builder[lat.neighbors()] = - γc
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)`.
......@@ -130,19 +130,20 @@ The system is thus the following
plt.rcParams['text.usetex'] = True
plt.rcParams['figure.figsize'] = (11.69,8.27)
plt.rcParams['figure.dpi'] = 40
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
%matplotlib notebook
.. jupyter-execute::
kwant.plot(builder)
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:
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::
......@@ -175,7 +176,7 @@ Now let's declare the other operators. Just like the particle `~kwant.operator.C
We can now initialize the solver
.. jupyter-execute::
# Occupation, for each lead
TL = 1.0 * Γ # Temperature in the left lead
......@@ -194,11 +195,11 @@ 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
particle_current = []
# Have the system evolve forward in time, calculating the operator over the system
particle_current = []
energy_current = []
energy_source = []
energy_density = []
energy_source = []
energy_density = []
heat_current_left = []
heat_current_right = []
......@@ -211,8 +212,9 @@ And run the simulation. Energy operators are evaluated in the same way as the pa
# evolve scattering states in time
solver.evolve(time)
solver.refine_intervals()
energy_current.append(solver.evaluate(energy_current_op))
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))
......@@ -240,21 +242,23 @@ 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) / Γ
particle_current = np.array(particle_current) / Γ**2
heat_current_left_2 = energy_current[:, 0] - energy_source[:, 0] - μL * particle_current[:, 0]
heat_current_right_2 = energy_current[:, 1] - energy_source[:, 1] - μR * particle_current[:, 1]
energy_density = np.array(energy_density) / Γ
# Plot
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 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.plot(times, energy_current[:, 0], 'r-', label = "Left incoming energy current $I^\\varepsilon_{0,-1}$")
plt.plot(times, energy_current[:, 1], 'm-', label = "Right incoming energy current $I^\\varepsilon_{0,1}$")
plt.plot(times, energy_source[:, 0], 'g-', label = "Left source $S_{-1}$")
plt.plot(times, energy_source[:, 1], 'y-', label = "Right source $S_{1}$")
plt.plot(times, heat_current_left, 'c-', label = "Left incoming heat current")
plt.plot(times, heat_current_right, 'b-', label = "Right incoming heat current")
plt.plot(times, heat_current_left_2, 'co', label = "$I^\\varepsilon_{0,-1}-S_{-1}-\mu_L I^N_{0,-1}$", markevery=0.025)
plt.plot(times, heat_current_right_2, 'bo', label = "$I^\\varepsilon_{0,1}-S_{1}-\mu_R I^N_{0,1}$", markevery=0.025)
plt.xlabel('time $[\\hbar / \\Gamma]$')
plt.ylabel('Energy flux $[\\Gamma^2 / \\hbar]$')
......
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