HEOM 5b: Discrete boson coupled to an impurity and fermionic leads¶

import contextlib
import dataclasses
import time

import matplotlib.pyplot as plt
import numpy as np

import qutip
from qutip import (
    Options,
    destroy,
    qeye,
    tensor,
)
from qutip.nonmarkov.heom import (
    HEOMSolver,
    LorentzianPadeBath,
)

from ipywidgets import IntProgress
from IPython.display import display

%matplotlib inline

@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}")

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
    )

# Define the system Hamiltonian:

@dataclasses.dataclass
class SystemParameters:
    e1: float = 0.3  # fermion mode energy splitting
    Omega: float = 0.2  # bosonic mode energy splitting
    Lambda: float = 0.12  # coupling between fermion and boson
    Nbos: int = 2

    def __post_init__(self):
        d = tensor(destroy(2), qeye(self.Nbos))
        a = tensor(qeye(2), destroy(self.Nbos))
        self.H = (
            self.e1 * d.dag() * d +
            self.Omega * a.dag() * a +
            self.Lambda * (a + a.dag()) * d.dag() * d
        )
        self.Q = d

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


sys_p = SystemParameters()

# 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
    gamma: float = 0.01  # coupling strength
    W: float = 1.0  # cut-off
    T: float = 0.025851991  # temperature (in eV)
    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)


# We set W = 1e4 to investigate the wide-band limit:

bath_L = LorentzianBathParameters(W=10**4, lead="L")
bath_R = LorentzianBathParameters(W=10**4, lead="R")

def steady_state_pade_for_theta(sys_p, bath_L, bath_R, theta, Nk, Nc, Nbos):
    """ Return the steady state current using the Pade approximation. """
    options = Options(nsteps=15000, store_states=True, rtol=1e-14, atol=1e-14)

    sys_p = sys_p.replace(Nbos=Nbos)
    bath_L = bath_L.replace(theta=theta)
    bath_R = bath_R.replace(theta=theta)

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

    solver_pade = HEOMSolver(
        sys_p.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")

    return np.real(2.434e-4 * 1e6 * current)

assert 1 == 1