Non-equilibrium density with dephasing in a linear chain

In this example we calculate the non-equilibrium density for a linear chain adding an electron-phonon elastic dephasing model.

import numpy
import pynegf
import scipy
import matplotlib.pyplot as plt
from scipy import sparse

We define the hamiltonian for a linear chain similar to other example.

def orthogonal_linear_chain(nsites=100, contact_size=2, coupling=1.0):
    mat = numpy.zeros(shape=(nsites, nsites), dtype='complex128')

    for i in range(nsites - contact_size):
        mat[i - 1, i] = coupling
    for i in range(nsites - contact_size, nsites):
        mat[i - 1, i] = coupling
    mat[0, nsites - contact_size] = coupling

    mat_csr = sparse.csr_matrix(mat)
    mat_csr = mat_csr + mat_csr.getH()
    mat_csr.sort_indices()

    return mat_csr

ideal_hamiltonian = orthogonal_linear_chain()

ham = ideal_hamiltonian.copy()

We define the system structure in the same way as in the previous examples.

The non-equilibrium calculation is very similar to the equilibrium case. We need to add a real-axis integration which covers energy point between the minimum and maximum electrode chemical potentials. The energy point density should be smaller than the thermal broadening, otherwise electronic structure features which are sharp in energy like resonant levels might be missed.

Note that we are setting a small chemical potential and this example is a basic proof of concept. In a realistic calculation we need to properly ensure consistency between the imposed boundary conditions and the chemical potential, and possibly solve self-consistently with a poisson solver.

We will solve the density for different coupling with the dephasing bath and look at the carrier density along the chain.

for coupling in [0.0, 0.05, 0.1]:
    negf = pynegf.PyNegf()
    negf.set_hamiltonian(ham)
    negf.set_identity_overlap(100)
    negf.init_structure(
        ncont=2,
        contend=numpy.array([97, 99]),
        surfend=numpy.array([95, 97]))
    negf.params.ec = -3.5
    negf.params.kbt_dm = (.01, .01)
    negf.params.g_spin = 2.0
    negf.params.mu[0] = -0.05
    negf.params.mu[1] = 0.05
    # Only the first value is used.
    #negf.params.np_real[0] = 200
    negf.params.verbose = 100
    negf.set_params()
    negf.set_diagonal_elph_dephasing(numpy.array([coupling]*100))
    negf.solve_density()
    density_matrix = negf.density_matrix()
    diagonal = numpy.real(density_matrix.diagonal())
    plt.plot(diagonal[:96], label='Coupling {}'.format(coupling))
    plt.legend()
    plt.grid()
plt.show()

INFO:root:Running libnegf on 1 processes
INFO:root:Running libnegf on 1 processes
INFO:root:Running libnegf on 1 processes

The electron-phonon dephasing model introduces back-scattering in the system and the non-equilibrium component of the density decreases linearly along the chain.