A Python package for exactly solving low-dimensional fermionic Hamiltonians. It was developed with the intention of constructing and diagonalizing matrix representations of low-dimensional, interacting fermionic Hamiltonians, in particular effective (surrogate) models that describe the low-energy (sub-gap) physics of the superconducting Anderson impurity model (see the arXiv preprint).
FermPy is based on NumPy and SciPy. Small systems use dense NumPy arrays as operator representations while large systems employ sparse SciPy arrays.
Functionality of this package includes:
- Construction of matrix representations of operators (
$H, S_z$ , etc.) from a basis (chosen from symmetry considerations) - Diagonalization (numpy.linalg.eigh, scipy.sparse.linalg.eigh) of operators
- Calculation of expectation values and thermal averages of operators
FermPy is found on PyPI and easily installed by running
pip install fermpy
from the command-line.
As an example of the intended usage of FermPy, we calculate and plot the excitation spectrum and average spin on a quantum dot (QD), coupled to a superconducting (SC) lead which is described by an effective three-level discretization (see the arXiv preprint for more details).
Import the necessary packages.
import numpy as np
import matplotlib.pyplot as plt
from fermpy.operator_sum import OpSum
from fermpy.hamiltonian_constructors import build_ham_op_sum_S_D
from fermpy.operator_representation import OpRep, operator_exp_val
from fermpy.bases import BasisS2Construct the operator sum for the Hamiltonian gamma and level positions xi are collected in Table 1 of the Supplemental Material of arXiv preprint. Everything is in units of the SC gap Delta.
Gamma = np.linspace(0, 10, 201)
h_os = build_ham_op_sum_S_D(epsilon_d=-7.5,
U=15,
Gamma=Gamma,
Delta=1,
gamma=[0.5080, 1.8454, 1.8454],
xi=[0, 2.7546, -2.7546])Define the basis. Here, we choose the eigenstates of total spin
basis = BasisS2(1 + 3) # One quantum dot level + three lead levelsBuild the matrix representation of the Hamiltonian from the operator sum h_os and basis. We only consider the singlet, doublet, and triplet states (
h_op = OpRep(h_os, basis, qns_to_consider=(0, 1, 2)) # quantum numbers: 2*sDiagonalize energies and eigenvectors are dictionaries with quantum numbers as keys (
energies: dict[int, np.ndarray]
eigenvectors: dict[int, np.ndarray]
energies, eigenvectors = h_op.diagonalize_all_qns(return_eigenvectors=True)Extract the ground state quantum number and energy.
s_gs = np.argmin([E[:, 0] for E in energies.values()], axis=0)
E_gs = np.array([energies[s][idx, 0] for idx, s in enumerate(s_gs)])Calculate the z-projection of spin on the quantum dot in the ground state. The syntax of OpSum is similar to ITensor.
sz_op = OpRep(opsum=OpSum([('sz', 0)]), basis=basis)
sz_exp = operator_exp_val(op_rep=sz_op, state_rep=eigenvectors, axis=1)
sz_exp_gs = [sz_exp[s][idx, 0] for idx, s in enumerate(s_gs)]Plot the excitation spectrum and average dot spin in the ground state.
fig, (ax_top, ax_bot) = plt.subplots(2, 1,
figsize=(5.5, 5.5),
sharex=True,
constrained_layout=True
)
ax_top.set(ylabel=r'$(E - E_\mathrm{GS}) / \Delta$',
xlim=(0, 10),
ylim=(-0.05, 2)
)
ax_bot.set(xlabel=r'$\Gamma / \Delta$',
ylabel=r'$\langle S_{z, \mathrm{dot}} \rangle$',
xlim=(0, 10),
ylim=(-0.05, 0.55)
)
ax_top.hlines(1, 0, 10, linestyle='dashed', color='k', linewidth=0.5)
ax_bot.hlines(0, 0, 10, linestyle='dashed', color='k', linewidth=0.5)
colors = {0: 'tab:blue', 1: 'tab:orange', 2: 'tab:green'}
lines_to_label = []
for s, energies_s in energies.items():
lines_to_label.append(ax_top.plot(Gamma,
energies_s - E_gs[:, None],
color=colors[s])[0]
)
ax_top.legend(handles=lines_to_label,
labels=[0, '1/2', 1],
loc='upper right',
title='s'
)
ax_bot.plot(Gamma, sz_exp_gs)Here's the output figure.