.. _tutorial: Tutorial: Energy and heat transport =================================== 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>`__. Readers who are not familiar with t-Kwant are invited to visit its `documentation webpage `_ first. Theory ------ In the following, we take :math:`e=\hbar=1`. 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 stationary scalar potential, that has a physical origin and is present from the remote past to the remote future. 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}`: * 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:: ρ^ε_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 .. math:: \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: .. math:: \frac{\text{d}}{\text{d}t} ρ^ε_i + \sum_j I^ε_{ji} = S^ε_i 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}^ε(t)` from site :math:`i` to site :math:`j` is the following: .. math:: 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_α}} \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`. The heat current leaving a lead :math:`α` to the central system :math:`C` is defined as: .. math:: I_α^H(t) = I^ε_α - S^ε_α - μ_α I^N_α 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 `~tkwantoperator` module and are to be used with `tkwant `_. The classes are: `~tkwantoperator.EnergyDensity`, `~tkwantoperator.EnergyCurrent`, `~tkwantoperator.EnergySource` and `~tkwantoperator.LeadHeatCurrent`. An additional class `~tkwantoperator.EnergyCurrentDivergence` 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. Example : 1D Quantum dot ------------------------ A 1D chain, whose sites have zero on-site energy and are connected with each other through nearest neighbor hoppings :math:`γ`, has a central site :math:`0` with an on-site energy :math:`ε_0(t)=ε_0+Δε\Theta(t-t_0)`. The central site is connected to its two nearest neighbors through a different hopping term :math:`γ_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: .. jupyter-execute:: 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 `~kwant.builder.Builder.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 :math:`I^ε_{0,-1}(t)` (respectively :math:`I^ε_{0,1}(t)`). .. jupyter-execute:: # 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 .. jupyter-execute:: :hide-code: import matplotlib.pyplot as plt %matplotlib inline plt.rcParams['text.usetex'] = True plt.rcParams['figure.figsize'] = (11.69,8.27) plt.rcParams['figure.dpi'] = 150 plt.rcParams['axes.formatter.useoffset'] = True plt.rcParams['axes.formatter.limits'] = (-2, 2) plt.rcParams['axes.grid'] = True plt.rcParams['font.size'] = 20 .. jupyter-execute:: 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 `~tkwantoperator.LeadHeatCurrent` 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 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 `~tkwantoperator.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 `~tkwantoperator.LeadHeatCurrent.hoppings` member. This approach (``attach_lead`` with ``add_cells=1`` then giving the added cells to `~tkwantoperator.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 `~tkwantoperator.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 .. jupyter-execute:: 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 operator, feel free to visit its respective documentation. .. jupyter-execute:: 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 .. jupyter-execute:: # 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. .. jupyter-execute:: # 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: .. jupyter-execute:: :hide-code: # Rescale results times = np.array(times) * Γ energy_current = np.array(energy_current) / Γ**2 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 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], '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]$') plt.legend() plt.show() plt.plot(times, energy_density, label = "Energy density $\\varepsilon_0$") plt.xlabel('time $[\\hbar / \\Gamma]$') plt.ylabel('Energy density $[\\Gamma]$') plt.legend() plt.show() A python file containing all the above lines of code is available: :download:`quantum_dot.py `. It can be ran in parallel to shorten the simulation time, as explained in :ref:`parallel-calculations` Reference --------- [1] A. Kara Slimane, P. Reck, G. Fleury, Simulating time-dependent thermoelectric transport in quantum systems [ `Phys. Rev. B `_ | `arXiv `_ ]