# Comparision between bias chemical potential and bias electrical potential¶

We consider a circular shaped central scattering region with two leads attached on the left and on the lower site as shown below. The left lead has an additional bias potential $$V_b$$. In second quantization, the Hamiltonian is

$\hat{H}(t) = \sum_{i,j} \gamma_{ij} \, \hat{c}^\dagger_i \hat{c}_j + \sum_{i < i_p} V_b \, \hat{c}^\dagger_i \hat{c}_i$

where $$\gamma_{ii} = 4$$ and for nearest-neighbors $$\gamma_{ij} = -1$$. We are interested in the time-dependent current $$I(t)$$ through the coupling element (highlighted in red in the system plot) from the left lead to the central scattering region:

$I(t) = - 2 \text{Im} \sum_{\alpha} \int \frac{dE}{2 \pi} f_\alpha(E) \psi_{\alpha E}^*(t,i_p) \mathbf{H}_{i_p, i_p + 1} \psi_{\alpha E}(t,i_p + 1)$

Two different results are obtained whether the potential $$V_b$$ is explicitly time dependent or not. We refer to Ref. 1 for more information.

Static voltage

In the static case, the potential has a constant value $$V_b \equiv V = 0.5$$. One can calculate the constant value of the current directly from the solution of above equations at the initial time $$t = 0$$.

Dynamic voltage pulse

In the dynamic case the potential is zero at the initial time, but is switched on smoothly to reach the similar value as in the static case. We parametrize the potential as

$\begin{split}V_b(t) = \begin{cases} 0, & \text{for } t < 0\\ \frac{V}{2} \left ( 1 - \cos\left (\frac{\pi t}{\tau} \right) \right) , & \text{for } 0 \leq t \leq \tau \\ V , & \text{for } t > \tau \end{cases}\end{split}$

such that the phase is

$\begin{split}\phi(t) = e/\hbar \int_{-\infty}^t d t' V_b(t') = \begin{cases} 0, & \text{for } t < 0\\ \frac{e V}{2 \hbar} \left ( t - \frac{\tau}{\pi} \sin\left (\frac{\pi t}{\tau} \right) \right), & \text{for } 0 \leq t \leq \tau \\ e V / \hbar (t - \tau / 2) , & \text{for } t > \tau \end{cases}\end{split}$

Performing a gauge transform, the time-dependent lead potential can be absorbed in the time-dependent system-lead coupling element between $$i_b$$ and $$i_b + 1$$:

$\hat{H}(t) = \sum_{i,j} \gamma_{ij} \, \hat{c}^\dagger_i \hat{c}_j - [e^{- i \phi(t)} - 1 ] \, \hat{c}^\dagger_{i_b + 1} \hat{c}_{i_b} + \text{h.c.}$

The result obtained from the static and from the dynamic case are shown below.

tkwant features highlighted

• Use of tkwant.manybody.State to solve the manybody Schrödinger equation.

• Use of refine_intervals for adaptive refinement the manybody integral.

• Use of tkwant.leads.add_voltage to add time-dependence to leads.

import tkwant
import kwant
from math import sin, pi
import matplotlib.pyplot as plt

def am_master():
"""Return true for the MPI master rank"""
return tkwant.mpi.get_communicator().rank == 0

def plot_currents(times, dynamic_current, static_current):
plt.plot(times, dynamic_current, lw=3, label='dynamic')
plt.plot([times[0], times[-1]], [static_current] * 2, lw=2, label='static')
plt.legend(loc=4)
plt.xlabel(r'time $t$')
plt.ylabel(r'current $I(t)$')
plt.show()

def circle(pos):
(x, y) = pos
return x ** 2 + y ** 2 < 5 ** 2

def make_system():
gamma = 1
lat = kwant.lattice.square(norbs=1)
syst = kwant.Builder()
syst[lat.shape(circle, (0, 0))] = 4 * gamma
syst[lat.neighbors()] = -1 * gamma

def voltage(site, V_static):
return 2 * gamma + V_static

leadr[lat(0, 0)] = 2 * gamma

return lat, syst

omega = 0.1
t_upper = pi / omega
if 0 <= time < t_upper:
return V_dynamic * (time - sin(omega * time) / omega) / 2
return V_dynamic * (time - t_upper / 2)

def tkwant_calculation(syst, operator, tmax, chemical_potential, V_static, V_dynamic):

params = {'V_static': V_static, 'V_dynamic': V_dynamic}

# occupation -- for each lead
occupations = [occup_left, occup_lower]

# Initialize the time-dependent manybody state. Here we use a lower
# accuracy for the initial adaptive refinement (at time=0) to speed up the calculation.
# For a realistic simulation, set rtol and atol to smaller values.
state = tkwant.manybody.State(syst, tmax, occupations, params=params, refine=False)
state.refine_intervals(rtol=0.05, atol=0.05)

# loop over time, calculating the current as we go
expectation_values = []
times = range(tmax)
for time in times:
state.evolve(time)
# For a realistic simulation, one must set rtol and atol to smaller values,
# or remove both arguments completely, to converge the manybody integral.
state.refine_intervals(rtol=0.05, atol=0.05)
result = state.evaluate(operator)
expectation_values.append(result)
return times, expectation_values

def main():

lat, syst = make_system()

# add the "dynamic" part of the voltage
for a in syst.neighbors(b) if a not in added_sites]
if am_master():
hop_lw=lambda a, b: 0.3 if (a, b) in left_lead_interface else 0.1,
hop_color=lambda a, b: 'red' if (a, b) in left_lead_interface else 'k')

syst = syst.finalized()

# observables

chemical_potential = 1
tmax = 100

# all the voltage applied via the time-dependenece
V_static = 0
V_dynamic = 0.5
times, dynamic_current = tkwant_calculation(syst, current_operator, tmax,
chemical_potential,
V_static, V_dynamic)

# all the voltage applied statically
# even though we are doing a "tkwant calculation", we only calculate
# the results at t=0.
V_static = 0.5
V_dynamic = 0
_, (static_current, *_) = tkwant_calculation(syst, current_operator, 1,
chemical_potential,
V_static, V_dynamic)

if am_master():
plot_currents(times, dynamic_current, static_current)

if __name__ == '__main__':
main()


The complete source code of this example can be found in voltage_raise.py.