Notebook

HEOM 5a: Fermionic single impurity model¶

Introduction¶

Here we model a single fermion coupled to two electronic leads or reservoirs (e.g., this can describe a single quantum dot, a molecular transistor, etc). Note that in this implementation we primarily follow the definitions used by Christian Schinabeck in his dissertation https://opus4.kobv.de/opus4-fau/files/10984/DissertationChristianSchinabeck.pdf and related publications.

Notation:

$K=L/R$ refers to left or right leads.
$\sigma=\pm$ refers to input/output

We choose a Lorentzian spectral density for the leads, with a peak at the chemical potential. The latter simplifies a little the notation required for the correlation functions, but can be relaxed if neccessary.

$$J(\omega) = \frac{\Gamma W^2}{((\omega-\mu_K)^2 +W^2 )}$$

The Fermi distribution function is:

$$f_F (x) = (\exp(x) + 1)^{-1}$$

Together these allow the correlation functions to be expressed as:

$$C^{\sigma}_K(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} d\omega e^{\sigma i \omega t} \Gamma_K(\omega) f_F[\sigma\beta(\omega - \mu)]$$

As with the bosonic case we can expand these in an exponential series using Matsubara, Pade, or fitting approaches.

The Pade decomposition approximates the Fermi distubition as

$$f_F(x) \approx f_F^{\mathrm{approx}}(x) = \frac{1}{2} - \sum_l^{l_{max}} \frac{2k_l x}{x^2 + \epsilon_l^2}$$

where $k_l$ and $\epsilon_l$ are co-efficients defined in J. Chem Phys 133,10106.

Evaluating the integral for the correlation functions gives,

$$C_K^{\sigma}(t) \approx \sum_{l=0}^{l_{max}} \eta_K^{\sigma_l} e^{-\gamma_{K,\sigma,l}t}$$

where:

$$\eta_{K,0} = \frac{\Gamma_KW_K}{2} f_F^{approx}(i\beta_K W)$$$$\gamma_{K,\sigma,0} = W_K - \sigma i\mu_K$$$$\eta_{K,l\neq 0} = -i\cdot \frac{k_m}{\beta_K} \cdot \frac{\Gamma_K W_K^2}{-\frac{\epsilon^2_m}{\beta_K^2} + W_K^2}$$$$\gamma_{K,\sigma,l\neq 0}= \frac{\epsilon_m}{\beta_K} - \sigma i \mu_K$$

In this notebook we:

compare the Matsubara and Pade approximations and contrast them with the analytical result for the current between the system and the leads.
plot the current through the qubit as a function of the different between the voltages of the leads.

Setup¶

In [1]:

import contextlib
import dataclasses
import time

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

import qutip
from qutip import (
    Options,
    basis,
    destroy,
    expect,
)
from qutip.nonmarkov.heom import (
    HEOMSolver,
    LorentzianBath,
    LorentzianPadeBath,
)

from ipywidgets import IntProgress
from IPython.display import display

%matplotlib inline

Helpers¶

In [2]:

@contextlib.contextmanager
def timer(label):
    """ Simple utility for timing functions:

        with timer("name"):
            ... code to time ...
    """
    start = time.time()
    yield
    end = time.time()
    print(f"{label}: {end - start}")

System and bath definition¶

And let us set up the system Hamiltonian, bath and system measurement operators:

In [3]:

# Define the system Hamiltonian:

# The system is a single fermion with energy level split e1:
d1 = destroy(2)
e1 = 1.0
H = e1 * d1.dag() * d1

In [4]:

# Define parameters for left and right fermionic baths.
# Each bath is a lead (i.e. a wire held at a potential)
# with temperature T and chemical potential mu.

@dataclasses.dataclass
class LorentzianBathParameters:
    lead: str
    Q: object  # coupling operator
    gamma: float = 0.01  # coupling strength
    W: float = 1.0  # cut-off
    T: float = 0.025851991  # temperature
    theta: float = 2.0  # bias

    def __post_init__(self):
        assert self.lead in ("L", "R")
        self.beta = 1 / self.T
        if self.lead == "L":
            self.mu = self.theta / 2.0
        else:
            self.mu = - self.theta / 2.0

    def J(self, w):
        """ Spectral density. """
        return self.gamma * self.W**2 / ((w - self.mu)**2 + self.W**2)

    def fF(self, w, sign=1.0):
        """ Fermi distribution for this bath. """
        x = sign * self.beta * (w - self.mu)
        return fF(x)

    def lamshift(self, w):
        """ Return the lamshift. """
        return 0.5 * (w - self.mu) * self.J(w) / self.W

    def replace(self, **kw):
        return dataclasses.replace(self, **kw)


def fF(x):
    """ Return the Fermi distribution. """
    # in units where kB = 1.0
    return 1 / (np.exp(x) + 1)


bath_L = LorentzianBathParameters(Q=d1, lead="L")
bath_R = LorentzianBathParameters(Q=d1, lead="R")

Spectral density¶

Let's plot the spectral density.

In [5]:

w_list = np.linspace(-2, 2, 100)

fig, ax = plt.subplots(figsize=(12, 7))

spec_L = bath_L.J(w_list)
spec_R = bath_R.J(w_list)

ax.plot(
    w_list, spec_L,
    "b--", linewidth=3,
    label=r"J_L(w)",
)
ax.plot(
    w_list, spec_R,
    "r--", linewidth=3,
    label=r"J_R(w)",
)

ax.set_xlabel("w")
ax.set_ylabel(r"$J(\omega)$")
ax.legend();

Emission and absorption by the leads¶

Next let's plot the emission and absorption by the leads.

In [6]:

w_list = np.linspace(-2, 2, 100)

fig, ax = plt.subplots(figsize=(12, 7))

# Left lead emission and absorption

gam_L_in = bath_L.J(w_list) * bath_L.fF(w_list, sign=1.0)
gam_L_out = bath_L.J(w_list) * bath_L.fF(w_list, sign=-1.0)

ax.plot(
    w_list, gam_L_in,
    "b--", linewidth=3,
    label=r"S_L(w) input (absorption)",
)
ax.plot(
    w_list, gam_L_out,
    "r--", linewidth=3,
    label=r"S_L(w) output (emission)",
)

# Right lead emission and absorption

gam_R_in = bath_R.J(w_list) * bath_R.fF(w_list, sign=1.0)
gam_R_out = bath_R.J(w_list) * bath_R.fF(w_list, sign=-1.0)

ax.plot(
    w_list, gam_R_in,
    "b", linewidth=3,
    label=r"S_R(w) input (absorption)",
)
ax.plot(
    w_list, gam_R_out,
    "r", linewidth=3,
    label=r"S_R(w) output (emission)",
)

ax.set_xlabel("w")
ax.set_ylabel(r"$S(\omega)$")
ax.legend();

Comparing the Matsubara and Pade approximations¶

Let's start by solving for the evolution using a Pade expansion of the correlation function of the Lorentzian spectral density:

In [7]:

# HEOM dynamics using the Pade approximation:

# Solver options, times to solve for and initial system state:
options = Options(nsteps=15000, store_states=True, rtol=1e-14, atol=1e-14)
tlist = np.linspace(0, 100, 1000)
rho0 = basis(2, 0) * basis(2, 0).dag()

Nk = 10  # Number of exponents to retain in the expansion of each bath

bathL = LorentzianPadeBath(
    bath_L.Q, bath_L.gamma, bath_L.W, bath_L.mu, bath_L.T,
    Nk, tag="L",
)
bathR = LorentzianPadeBath(
    bath_R.Q, bath_R.gamma, bath_R.W, bath_R.mu, bath_R.T,
    Nk, tag="R",
)

with timer("RHS construction time"):
    solver_pade = HEOMSolver(H, [bathL, bathR], max_depth=2, options=options)

with timer("ODE solver time"):
    result_pade = solver_pade.run(rho0, tlist, ado_return=True)

with timer("Steady state solver time"):
    rho_ss_pade, ado_ss_pade = solver_pade.steady_state()

RHS construction time: 0.44698143005371094
ODE solver time: 1.6603105068206787
Steady state solver time: 0.11121225357055664

Now let us plot the result which shows the decay of the initially excited impurity. This is not very illuminating, but we will compare it with the Matsubara expansion and analytic solution sortly:

In [8]:

# Plot the Pade results
fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 8))

axes.plot(
    tlist, expect(result_pade.states, rho0),
    'r--', linewidth=2,
    label="P11 (Pade)",
)
axes.axhline(
    expect(rho_ss_pade, rho0),
    color='r', linestyle="dotted", linewidth=1,
    label="P11 (Pade steady state)",
)

axes.set_xlabel('t', fontsize=28)
axes.legend(fontsize=12);

Now let us do the same for the Matsubara expansion:

In [9]:

# HEOM dynamics using the Matsubara approximation:

bathL = LorentzianBath(
    bath_L.Q, bath_L.gamma, bath_L.W, bath_L.mu, bath_L.T,
    Nk, tag="L",
)
bathR = LorentzianBath(
    bath_R.Q, bath_R.gamma, bath_R.W, bath_R.mu, bath_R.T,
    Nk, tag="R",
)

with timer("RHS construction time"):
    solver_mats = HEOMSolver(H, [bathL, bathR], max_depth=2, options=options)

with timer("ODE solver time"):
    result_mats = solver_mats.run(rho0, tlist, ado_return=True)

with timer("Steady state solver time"):
    rho_ss_mats, ado_ss_mats = solver_mats.steady_state()

RHS construction time: 0.4666931629180908
ODE solver time: 1.7071921825408936
Steady state solver time: 0.11576962471008301

We see a marked difference in the Matsubara vs Pade results:

In [10]:

# Plot the Pade results
fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 8))

axes.plot(
    tlist, expect(result_pade.states, rho0),
    'r--', linewidth=2,
    label="P11 (Pade)",
)
axes.axhline(
    expect(rho_ss_pade, rho0),
    color='r', linestyle="dotted", linewidth=1,
    label="P11 (Pade steady state)",
)

axes.plot(
    tlist, expect(result_mats.states, rho0),
    'b--', linewidth=2,
    label="P11 (Mats)",
)
axes.axhline(
    expect(rho_ss_mats, rho0),
    color='b', linestyle="dotted", linewidth=1,
    label="P11 (Mats steady state)",
)

axes.set_xlabel('t', fontsize=28)
axes.legend(fontsize=12);

But which is more correct? The Matsubara or the Pade result?

One advantage of this simple model is that the steady state current to the baths is analytically solvable, so we can check convergence of the result by calculating it analytically (the sum of the currents to and from the system in the steady state must be zero, so the current from one bath is the same as the current to the other).

See the QuTiP-BoFiN paper for a detailed description and references for the analytic result. Below we just perform the required integration numerically.

In [11]:

def analytical_steady_state_current(bath_L, bath_R, e1):
    """ Calculate the analytical steady state current. """

    def integrand(w):
        return (2 / np.pi) * (
            bath_L.J(w) * bath_R.J(w) * (bath_L.fF(w) - bath_R.fF(w)) /
            (
                (bath_L.J(w) + bath_R.J(w))**2 +
                4*(w - e1 - bath_L.lamshift(w) - bath_R.lamshift(w))**2
            )
        )

    def real_part(x):
        return np.real(integrand(x))

    def imag_part(x):
        return np.imag(integrand(x))

    # in principle the bounds for the integral should be rechecked if
    # bath or system parameters are changed substantially:
    bounds = [-10, 10]

    real_integral, _ = quad(real_part, *bounds)
    imag_integral, _ = quad(imag_part, *bounds)

    return real_integral + 1.0j * imag_integral


curr_ss_analytic = analytical_steady_state_current(bath_L, bath_R, e1)

print(f"Analytical steady state current: {curr_ss_analytic}")

Analytical steady state current: (0.0008130726698792024+0j)

To compare the analytical result above with the result from the HEOM, we need to be able to calculate the current from the system to the bath from the HEOM result. In the HEOM description, these currents are captured in the first level auxilliary density operators (ADOs).

In the function state_current(...) below, we extract the first level ADOs for the specified bath and sum the contributions to the current from each:

In [12]:

def state_current(ado_state, bath_tag):
    """ Determine current from the given bath (either "R" or "L") to
        the system in the given ADO state.
    """
    level_1_aux = [
        (ado_state.extract(label), ado_state.exps(label)[0])
        for label in ado_state.filter(level=1, tags=[bath_tag])
    ]

    def exp_sign(exp):
        return 1 if exp.type == exp.types["+"] else -1

    def exp_op(exp):
        return exp.Q if exp.type == exp.types["+"] else exp.Q.dag()

    return -1.0j * sum(
        exp_sign(exp) * (exp_op(exp) * aux).tr()
        for aux, exp in level_1_aux
    )

Now we can calculate the steady state currents from the Pade and Matsubara HEOM results:

In [13]:

curr_ss_pade_L = state_current(ado_ss_pade, "L")
curr_ss_pade_R = state_current(ado_ss_pade, "R")

print(f"Pade steady state current (L): {curr_ss_pade_L}")
print(f"Pade steady state current (R): {curr_ss_pade_R}")

Pade steady state current (L): (-0.0008130302805827138-8.673617379884035e-19j)
Pade steady state current (R): (0.0008130302805827146+3.7947076036992655e-19j)

In [14]:

curr_ss_mats_L = state_current(ado_ss_mats, "L")
curr_ss_mats_R = state_current(ado_ss_mats, "R")

print(f"Matsubara steady state current (L): {curr_ss_mats_L}")
print(f"Matsubara steady state current (R): {curr_ss_mats_R}")

Matsubara steady state current (L): (-0.0011018485316349562-8.673617379884035e-19j)
Matsubara steady state current (R): (0.0011018485316349564-1.0842021724855044e-19j)

Note that the currents from each bath balance as is required by the steady state, but the value of the current is different for the Pade and Matsubara results.

Now let's compare all three:

In [15]:

print(f"Pade current (R): {curr_ss_pade_R}")
print(f"Matsubara current (R): {curr_ss_mats_R}")
print(f"Analytical curernt: {curr_ss_analytic}")

Pade current (R): (0.0008130302805827146+3.7947076036992655e-19j)
Matsubara current (R): (0.0011018485316349564-1.0842021724855044e-19j)
Analytical curernt: (0.0008130726698792024+0j)

In this case we observe that the Pade approximation has converged more closely to the analytical current than the Matsubara.

The Matsubara result could be improved by increasing the number of terms retained in the Matsubara expansion (i.e. increasing Nk).

Current as a function of bias voltage¶

Now lets plot the current as a function of bias voltage (the bias voltage is the parameter theta for the two baths).

We will calculate the steady state current for each theta both analytically and using the HEOM with the Pade correlation expansion approximation.

In [16]:

# Theta (bias voltages)

thetas = np.linspace(-4, 4, 100)

# Setup a progress bar:

progress = IntProgress(min=0, max=2 * len(thetas))
display(progress)

# Calculate the current for the list of thetas


def current_analytic_for_theta(e1, bath_L, bath_R, theta):
    """ Return the analytic current for a given theta. """
    current = analytical_steady_state_current(
        bath_L.replace(theta=theta),
        bath_R.replace(theta=theta),
        e1,
    )
    progress.value += 1
    return np.real(current)


def current_pade_for_theta(H, bath_L, bath_R, theta, Nk):
    """ Return the steady state current using the Pade approximation. """
    bath_L = bath_L.replace(theta=theta)
    bath_R = bath_R.replace(theta=theta)

    bathL = LorentzianPadeBath(
        bath_L.Q, bath_L.gamma, bath_L.W, bath_L.mu, bath_L.T,
        Nk, tag="L",
    )
    bathR = LorentzianPadeBath(
        bath_R.Q, bath_R.gamma, bath_R.W, bath_R.mu, bath_R.T,
        Nk, tag="R",
    )

    solver_pade = HEOMSolver(H, [bathL, bathR], max_depth=2, options=options)
    rho_ss_pade, ado_ss_pade = solver_pade.steady_state()
    current = state_current(ado_ss_pade, bath_tag="R")

    progress.value += 1
    return np.real(current)


curr_ss_analytic_thetas = [
    current_analytic_for_theta(e1, bath_L, bath_R, theta)
    for theta in thetas
]

# The number of expansion terms has been dropped to Nk=6 to speed
# up notebook execution. Increase to Nk=10 for more accurate results.
curr_ss_pade_theta = [
    current_pade_for_theta(H, bath_L, bath_R, theta, Nk=6)
    for theta in thetas
]

IntProgress(value=0, max=200)

Below we plot the results and see that even with Nk=6, the HEOM Pade approximation gives good results for the steady state current. Increasing Nk to 10 gives very accurate results.

In [17]:

fig, ax = plt.subplots(figsize=(12, 7))

ax.plot(
    thetas, 2.434e-4 * 1e6 * np.array(curr_ss_analytic_thetas),
    color="black", linewidth=3,
    label=r"Analytical",
)
ax.plot(
    thetas, 2.434e-4 * 1e6 * np.array(curr_ss_pade_theta),
    'r--', linewidth=3,
    label=r"HEOM Pade $N_k=10$, $n_{\mathrm{max}}=2$",
)


ax.locator_params(axis='y', nbins=4)
ax.locator_params(axis='x', nbins=4)

ax.set_xticks([-2.5, 0, 2.5])
ax.set_xticklabels([-2.5, 0, 2.5])
ax.set_xlabel(r"Bias voltage $\Delta \mu$ ($V$)", fontsize=28)
ax.set_ylabel(r"Current ($\mu A$)", fontsize=28)
ax.legend(fontsize=25);

About¶

In [18]:

qutip.about()

QuTiP: Quantum Toolbox in Python
================================
Copyright (c) QuTiP team 2011 and later.
Current admin team: Alexander Pitchford, Nathan Shammah, Shahnawaz Ahmed, Neill Lambert, Eric Giguère, Boxi Li, Jake Lishman, Simon Cross and Asier Galicia.
Board members: Daniel Burgarth, Robert Johansson, Anton F. Kockum, Franco Nori and Will Zeng.
Original developers: R. J. Johansson & P. D. Nation.
Previous lead developers: Chris Granade & A. Grimsmo.
Currently developed through wide collaboration. See https://github.com/qutip for details.

QuTiP Version:      4.7.1
Numpy Version:      1.22.4
Scipy Version:      1.8.1
Cython Version:     0.29.33
Matplotlib Version: 3.5.2
Python Version:     3.10.4
Number of CPUs:     2
BLAS Info:          Generic
OPENMP Installed:   False
INTEL MKL Ext:      False
Platform Info:      Linux (x86_64)
Installation path:  /home/runner/work/qutip-tutorials/qutip-tutorials/qutip/qutip
================================================================================
Please cite QuTiP in your publication.
================================================================================
For your convenience a bibtex reference can be easily generated using `qutip.cite()`

Testing¶

This section can include some tests to verify that the expected outputs are generated within the notebook. We put this section at the end of the notebook, so it's not interfering with the user experience. Please, define the tests using assert, so that the cell execution fails if a wrong output is generated.

In [19]:

assert np.allclose(curr_ss_pade_L + curr_ss_pade_R, 0)
assert np.allclose(curr_ss_mats_L + curr_ss_mats_R, 0)
assert np.allclose(curr_ss_pade_R, curr_ss_analytic, rtol=1e-4)