Pymablock: an algorithm and a package for quasi-degenerate perturbation theory

Isidora Araya Day1,2$\star$, Sebastian Miles1, Hugo K. Kerstens2, Daniel
Varjas3,4, Anton R. Akhmerov2$\dagger$

1 QuTech, Delft University of Technology, Delft 2600 GA, The Netherlands

2 Kavli Institute of Nanoscience, Delft University of Technology, 2600 GA
Delft, The Netherlands

3 Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Strasse
38, 01187 Dresden, Germany

4 Institute for Theoretical Solid State Physics, IFW Dresden and Würzburg-
Dresden Cluster of Excellence ct.qmat, Helmholtzstr. 20, 01069 Dresden,
Germany

⋆<EMAIL_ADDRESS><EMAIL_ADDRESS>
April 4, 2024

## Abstract

A common technique in the study of complex quantum-mechanical systems is to
reduce the number of degrees of freedom in the Hamiltonian by using quasi-
degenerate perturbation theory. While the Schrieffer–Wolff transformation
achieves this and constructs an effective Hamiltonian, its scaling is
suboptimal, and implementing it efficiently is both challenging and error-
prone. We introduce an algorithm for constructing an equivalent effective
Hamiltonian as well as a Python package, Pymablock, that implements it. Our
algorithm combines an optimal asymptotic scaling with a range of other
improvements. The package supports numerical and analytical calculations of
any order and it is designed to be interoperable with any other packages for
specifying the Hamiltonian. We demonstrate how the package handles
constructing a k.p model, analyses a superconducting qubit, and computes the
low-energy spectrum of a large tight-binding model. We also compare its
performance with reference calculations and demonstrate its efficiency.

###### Contents

  1. 1 Introduction
  2. 2 Constructing an effective model
    1. 2.1 k.p model of bilayer graphene
    2. 2.2 Dispersive shift of a transmon qubit coupled to a resonator
    3. 2.3 Induced gap in a double quantum dot
  3. 3 Perturbative block-diagonalization algorithm
    1. 3.1 Problem statement
    2. 3.2 Existing solutions
    3. 3.3 Pymablock’s algorithm
    4. 3.4 Equivalence to Schrieffer–Wolff transformation
    5. 3.5 Extra optimization: common subexpression elimination
  4. 4 Implementation
    1. 4.1 The data structure for block operator series
    2. 4.2 The implicit method for large sparse Hamiltonians
  5. 5 Benchmark
  6. 6 Conclusion

##  1 Introduction

Effective models enable the study of complex quantum systems by reducing the
dimensionality of the Hilbert space. Their construction separates the low and
high-energy subspaces by block-diagonalizing a perturbed Hamiltonian

$\mathcal{H}=\begin{pmatrix}H_{0}^{AA}&0\\\
0&H_{0}^{BB}\end{pmatrix}+\mathcal{H}^{\prime},$ (1)

where $H_{0}^{AA}$ and $H_{0}^{BB}$ are separated by an energy gap, and
$\mathcal{H}^{\prime}$ is a series in a perturbative parameter. This procedure
requires finding a series of the basis transformation $\mathcal{U}$ that is
unitary and that also cancels the off-diagonal block of the transformed
Hamiltonian order by order, as shown in Fig. 1. The low-energy effective
Hamiltonian $\tilde{\mathcal{H}}^{AA}$ is then a series in the perturbative
parameter, whose eigenvalues and eigenvectors are approximate solutions of the
complete Hamiltonian. As a consequence, the effective model is sufficient to
describe the low-energy properties of the original system while also being
simpler and easier to handle.

A common approach to constructing an effective Hamiltonian is the
Schrieffer–Wolff transformation [1, 2], also known as Löwdin partitioning [3],
or quasi-degenerate perturbation theory. This method parameterizes the unitary
transformation $\mathcal{U}=e^{-\mathcal{S}}$ and finds the series
$\mathcal{S}$ that decouples the $A$ and $B$ subspaces of
$\tilde{\mathcal{H}}=e^{\mathcal{S}}\mathcal{H}e^{-\mathcal{S}}$. This idea
enabled advances in multiple fields of quantum physics. As an example, all the
k.p models are a result of treating crystalline momentum as a perturbation
that only weakly mixes atomic orbitals separated in energy [4]. More broadly,
this method serves as a go-to tool in the study of superconducting circuits
and quantum dots, where couplings between circuit elements and drives are
treated as perturbations to reproduce the dynamics of the system [5, 6].

Figure 1:  Block-diagonalization of a Hamiltonian with a first order
perturbation.

Constructing effective Hamiltonians is, however, both algorithmically complex
and computationally expensive. This is a consequence of the recursive
equations that define the unitary transformation, which require an
exponentially growing number of matrix products in each order. In particular,
already a 4-th order perturbative expansion that is necessary for many
applications may require hundreds of terms. While the computational complexity
is only a nuisance when analysing model systems, it becomes a bottleneck
whenever the Hilbert space is high-dimensional. Several alternative approaches
improve the performance of the Schrieffer–Wolff algorithm by either using
different parametrizations of the unitary transformation [7, 8, 3, 9, 10],
adjusting the problem setting to density matrix perturbation theory [11, 12],
or a finding a similarity transform instead of a unitary [13]. A more recent
line of research even applies the ideas of Schrieffer–Wolff transformation to
quantum algorithms for the study of many-body systems [14, 15]. Despite these
advances, neither of the approaches combines an optimal scaling with the
ability to construct effective Hamiltonians.

We introduce an algorithm to construct effective models with optimal scaling,
thus making it possible to find high order effective models for systems with
millions of degrees of freedom. This algorithm exploits the efficiency of
recursive evaluations of series satisfying polynomial constraints and obtains
the same effective Hamiltonian as the Schrieffer–Wolff transformation. We make
the algorithm available via the open source package Pymablock 111The
documentation and tutorials are available in
https://pymablock.readthedocs.io/(PYthon MAtrix BLOCK-diagonalization), a
versatile tool for the study of numerical and symbolic models.

##  2 Constructing an effective model

We illustrate the construction of effective models by considering several
representative examples. The simplest application of effective models is the
reduction of finite symbolic Hamiltonians, which appear in the derivation of
low-energy dispersions of materials. Starting from a tight-binding model, one
performs Taylor expansions of the Hamiltonian near a $k$-point, and then
eliminates several high-energy states [4, 16]. In the study of superconducting
qubits, for example, the Hamiltonian contains several bosonic operators, so
its Hilbert space is infinite-dimensional, and the coupling between bosons
makes the Hamiltonian impossible to diagonalize. The effective qubit model
describes the analytical dependence of qubit frequencies and couplings on the
circuit parameters [17, 18, 5, 19, 20, 21]. This allows to design circuits
that realize a desired qubit Hamiltonian, as well as ways to understand and
predict qubit dynamics, for which computational tools are being actively
developed [22, 23, 24]. Finally, mesoscopic quantum devices are described by a
single particle tight-binding model with short range hoppings. This produces a
numerical Hamiltonian that is both big and sparse, which allows to compute a
few of its states but not the full spectrum [25]. Because only the low-energy
states contribute to observable properties, deriving how they couple enables a
more efficient simulation of the system’s behavior.

Pymablock treats all the problems, including the ones above, using a unified
approach that only requires three steps:

  * •

Define a Hamiltonian

  * •

Call `pymablock.block_diagonalize`

  * •

Request the desired order of the effective Hamiltonian

The following code snippet shows how to use Pymablock to compute the fourth
order correction to an effective Hamiltonian $\tilde{\mathcal{H}}$:

⬇

# Define perturbation theory

H_tilde, *_ = block_diagonalize([H_0, H_1], subspace_eigenvectors=[vecs_A,
vecs_B])

\par# Request 4th order correction to the effective Hamiltonian

H_AA_4 = H_tilde[0, 0, 4]

The function `block_diagonalize` interprets the Hamiltonian $H_{0}+H_{1}$ as a
series with two terms, zeroth and first order and calls the block
diagonalization routine. The subspaces to decouple are spanned by the
eigenvectors `vecs_A` and `vecs_B` of $H_{0}$. This is the main function of
Pymablock, and it is the only one that the user ever needs to call. Its first
output is a multivariate series whose terms are different blocks and orders of
the transformed Hamiltonian. Calling `block_diagonalize` only defines the
computational problem, whereas querying the elements of `H_tilde` does the
actual calculation of the desired order. This interface treats arbitrary
formats of Hamiltonians and system descriptions on the same footing and
supports both numerical and symbolic computations.

###  2.1 k.p model of bilayer graphene

To illustrate how to use Pymablock with analytic models, we consider two
layers of graphene stacked on top of each other, as shown in Fig. 2. Our goal
is to find the low-energy model near the $\mathbf{K}$ point [16]. To do this,
we first construct the tight-binding model Hamiltonian of bilayer graphene.

Figure 2: Crystal structure and hoppings of AB-stacked bilayer graphene.

The main features of the model are its 4-atom unit cell spanned by vectors
$\mathbf{a}_{1}=(1/2,\sqrt{3}/2)$ and $\mathbf{a}_{2}=(-1/2,\sqrt{3}/2)$, and
with wave functions $\phi_{A,1},\phi_{B,1},\phi_{A,2},\phi_{B,2}$, where $A$
and $B$ indices are the two sublattices, and $1,2$ are the layers. The model
has hoppings $t_{1}$ and $t_{2}$ within and between the layers, respectively,
as shown in Fig. 2. We also include a layer-dependent onsite potential $\pm
m$.

We define the Bloch Hamiltonian using the Sympy package for symbolic Python
[26].

⬇

t_1, t_2, m = sympy.symbols(”t_1 t_2 m”, real=True)

alpha = sympy.symbols(r”\alpha”)

\parH = Matrix([

[m, t_1 * alpha, 0, 0],

[t_1 * alpha.conjugate(), m, t_2, 0],

[0, t_2, -m, t_1 * alpha],

[0, 0, t_1 * alpha.conjugate(), -m]]

)

$H=\begin{pmatrix}m&t_{1}\alpha&0&0\\\ t_{1}\alpha^{*}&m&t_{2}&0\\\
0&t_{2}&-m&t_{1}\alpha\\\ 0&0&t_{1}\alpha^{*}&-m\end{pmatrix}$

where
$\alpha(\mathbf{k})=1+e^{i\mathbf{k}\cdot\mathbf{a}_{1}}+e^{\mathbf{k}\cdot\mathbf{a}_{2}}$,
with $k$ the wave vector. We consider $\mathbf{K}=(4\pi/3,0)$ the reference
point point in $\mathbf{k}$-space: $\mathbf{k}=(4\pi/3+k_{x},k_{y})$ because
$\alpha(\mathbf{K})=0$, making $k_{x}$ and $k_{y}$ small perturbations.
Additionally, we consider $m\ll t_{2}$ a perturbative parameter.

To call `block_diagonalize`, we need to define the subspaces for the block
diagonalization, so we compute the eigenvectors of the unperturbed Hamiltonian
at the $\mathbf{K}$ point, $H(\alpha(\mathbf{K})=m=0)$. Then, we substitute
$\alpha(\mathbf{k})$ into the Hamiltonian, and call the block diagonalization
routine using that $k_{x}$, $k_{y}$, and $m$ are perturbative parameters via
the `symbols` argument.

⬇

vecs = H.subs({alpha: 0, m: 0}).diagonalize(normalize=True)[0]

\parH_tilde, U, U_adjoint = block_diagonalize(

H.subs({alpha: alpha_k}),

symbols=(k_x, k_y, m),

subspace_eigenvectors=[vecs[:, :2], vecs[:, 2:]] # AA, BB

)

The order of the variables in the perturbative series will be that of

    symbols}.
    For example, requesting the term $\propto k_x^

i k_y^j m^l$fromtheeffectivemodelisdonebycalling\verb|H_tilde[0, 0, i, j,
l]|,wherethefirsttwoindicesaretheblockindices(AA).Theseriesoftheunitarytransformation$U$and$U^†$arealsodefined,andwemayusethemtotransformotheroperators.\par
Wecollectcorrectionsuptothirdorderinmomentumtocomputethestandardquadraticdispersionofbilayergrapheneandtrigonalwarping.Wequerythesetermsfrom\verb|H_tilde|andthoseproportionaltomasstoobtaintheeffectiveHamiltonian(shownasproducedbythecode):{\small\begin{gathered}\tilde{H}_{\textrm{eff}}=\begin{bmatrix}m&\frac{3t_{1}^{2}}{4t_{2}}(-k_{x}^{2}-2ik_{x}k_{y}+k_{y}^{2})\\\
\frac{3t_{1}^{2}}{4t_{2}}(-k_{x}^{2}+2ik_{x}k_{y}+k_{y}^{2})&-m\end{bmatrix}+\\\
\begin{bmatrix}\frac{3mt_{1}^{2}}{2t_{2}^{2}}(-k_{x}^{2}-k_{y}^{2})&\frac{\sqrt{3}t_{1}^{2}}{8t_{2}}(k_{x}^{3}-5ik_{x}^{2}k_{y}+9k_{x}k_{y}^{2}+3ik_{y}^{3})\\\
\frac{\sqrt{3}t_{1}^{2}}{8t_{2}}(k_{x}^{3}+5ik_{x}^{2}k_{y}+9k_{x}k_{y}^{2}-3ik_{y}^{3})&\frac{3mt_{1}^{2}}{2t_{2}^{2}}(k_{x}^{2}+k_{y}^{2})\end{bmatrix}\end{gathered}}Thefirsttermisthestandardquadraticdispersionofgappedbilayergraphene.Thesecondtermcontainstrigonalwarpingandthecouplingbetweenthegapandmomentum.Allthetermstakelessthantwosecondsinapersonalcomputertocompute.\par$222The
full code is available at
https://pymablock.readthedocs.io/en/latest/tutorial/bilayer˙graphene.html.

###  2.2 Dispersive shift of a transmon qubit coupled to a resonator

The need for analytical effective Hamiltonians often arises in circuit quantum
electrodynamics (cQED) problems, which we illustrate by studying a transmon
qubit coupled to a resonator [5]. Specifically, we choose the standard problem
of finding the frequency shift of the resonator due to its coupling to the
qubit, a phenomenon used to measure the qubit’s state [17]. The Hamiltonian of
the system is given by

$\mathcal{H}=-\omega_{t}(a^{\dagger}_{t}a_{t}-\frac{1}{2})+\frac{\alpha}{2}a^{\dagger}_{t}a^{\dagger}_{t}a_{t}a_{t}+\omega_{r}(a^{\dagger}_{r}a_{r}+\frac{1}{2})-g(a^{\dagger}_{t}-a_{t})(a^{\dagger}_{r}-a_{r}),$
(2)

where $a_{t}$ and $a_{r}$ are bosonic annihilation operators of the transmon
and resonator, respectively, and $\omega_{t}$ and $\omega_{r}$ are their
frequencies. The transmon has an anharmonicity $\alpha$, so that its energy
levels are not equally spaced. In presence of both the coupling $g$ between
the transmon and the resonator and the anharmonicity, this Hamiltonian admits
no analytical solution. We therefore treat $g$ as a perturbative parameter. To
deal with the infinite dimensional Hilbert space, we observe that the
perturbation only changes the occupation numbers of the transmon and the
resonator by $\pm 1$. Therefore computing $n$-th order corrections to the
$n_{0}$-th state allows to disregard states with any occupation numbers larger
than $n_{0}+n/2$. We want to compute the second order correction to the levels
with occupation numbers of either the transmon or the resonator being $0$ and
$1$. We accordingly truncate the Hilbert space to the lowest 3 levels of the
transmon and the resonator. The resulting Hamiltonian is a $9\times 9$ matrix
that we construct using Sympy [26]. Finally, to compute the energy corrections
of the lowest levels, we call `block_diagonalize` for each state separately,
replicating a regular perturbation theory calculation for single
wavefunctions. To do this, we observe that $H_{0}$ is diagonal, and use
`subspace_indices` to assign the elements of its eigenbasis to the $A$ (0) or
$B$ (1) subspace. For example, to find the qubit-dependent frequency shift of
the resonator, $\chi$, we start by computing the second order correction to
$\ket{0_{t}0_{r}}$:

⬇

indices = [0, 1, 1, 1, 1, 1, 1, 1, 1] # 00 is the first state in the basis

H_tilde, *_ = block_diagonalize(H, subspace_indices=indices, symbols=[g])

H_tilde[0, 0, 2][0, 0] # 2nd order correction

$E^{(2)}_{00}=\frac{g^{2}}{-\omega_{r}+\omega_{t}}.$ (3)

Repeating this process for the other levels requires changing
`subspace_indices` according to the basis of $H$, and yields the desired
resonator frequency shift:

$\displaystyle\chi$
$\displaystyle=(E^{(2)}_{11}-E^{(2)}_{10})-(E^{(2)}_{01}-E^{(2)}_{00})$ (4)
$\displaystyle=-\frac{2g^{2}}{\alpha+\omega_{r}-\omega_{t}}+\frac{2g^{2}}{-\alpha+\omega_{r}+\omega_{t}}-\frac{2g^{2}}{\omega_{r}+\omega_{t}}+\frac{2g^{2}}{\omega_{r}-\omega_{t}}$
$\displaystyle=-\frac{4\alpha
g^{2}\left(\alpha\omega_{t}-\omega_{r}^{2}-\omega_{t}^{2}\right)}{\left(\omega_{r}-\omega_{t}\right)\left(\omega_{r}+\omega_{t}\right)\left(-\alpha+\omega_{r}+\omega_{t}\right)\left(\alpha+\omega_{r}-\omega_{t}\right)}.$

In this example, we have not used the rotating wave approximation, including
the frequently omitted counter-rotating terms $\sim a_{r}a_{t}$ to illustrate
the extensibility of Pymablock. Computing higher order corrections to the
qubit frequency only requires increasing the size of the truncated Hilbert
space and requesting `H_tilde[0, 0, n]` to the desired order $n$.

###  2.3 Induced gap in a double quantum dot

Large systems pose an additional challenge due to the cubic scaling of linear
algebra routines with matrix size. To overcome this, Pymablock is equipped
with an implicit method, which utilizes the sparsity of the input and avoids
the construction of the full transformed Hamiltonian. We illustrate the
efficiency of this method by applying it to a system of two quantum dots
coupled to a superconductor between them, shown in Fig. 3, and described by
the Bogoliubov-de Gennes Hamiltonian:

$H_{BdG}=\begin{cases}(\mathbf{k}^{2}/2m-\mu_{sc})\sigma_{z}+\Delta\sigma_{x}\quad&\text{for
}L/3\leq x\leq 2L/3,\\\
(\mathbf{k}^{2}/2m-\mu_{n})\sigma_{z}\quad&\text{otherwise},\end{cases}$ (5)

where the Pauli matrices $\sigma_{z}$ and $\sigma_{x}$ act in the electron-
hole space, $\mathbf{k}$ is the 2D wave vector, $m$ is the effective mass, and
$\Delta$ the superconducting gap. We use the Kwant package [27] to build the
Hamiltonian of the system 333The full code is available at
https://pymablock.readthedocs.io/en/latest/tutorial/induced˙gap.html., which
we define over a square lattice of $L\times W=200\times 40$ sites. On top of
this, we consider two perturbations: the barrier strength between the quantum
dots and the superconductor, $t_{b}$, and an asymmetry of the dots’
potentials, $\delta\mu$. The system is large: it is a sparse array of size
$63042\times 63042$, with $333680$ non-zero elements, so even storing all the
eigenvectors would take $60$ GB of memory. The perturbations are also sparse,
with $632$, and $126084$ non-zero elements for the barrier strength and the
potential asymmetry, respectively. The sparsity structure of the Hamiltonian
and the perturbations is shown in the left panel of Fig. 3, where we use a
smaller system of $L\times W=8\times 2$ for visualization. Therefore, we use
sparse diagonalization [28] and compute only four eigenvectors of the
unperturbed Hamiltonian closest to zero energy, which are the Andreev states
of the quantum dots.

We now call the block diagonalization routine and provide the computed
eigenvectors.

Because we only provide the low-energy subspace, Pymablock uses the implicit
method. Calling `block_diagonalize` is now the most time-consuming step
because it requires pre-computing several decompositions of the full
Hamiltonian. It is, however, manageable and it only produces a constant
overhead of less than three seconds. To compute the spectrum, we collect the
lowest three orders in each parameter in an appropriately sized tensor.

This takes two more seconds to run, and we can now compute the low-energy
spectrum after rescaling the perturbative corrections by the magnitude of each
perturbation.

Finally, we plot the spectrum of the $2$ Andreev states in Fig. 3.

Figure 3:  Hamiltonian (left) and Andreev levels (right) of two quantum dots
coupled to a superconductor (inset). The barrier $t_{b}$ between the dots and
the superconductor, $H_{10}$, and the asymmetry $\delta\mu$ between the dots’
potential, $H_{01}$, are perturbations.

As expected, the crossing at $E=0$ due to the dot asymmetry is lifted when the
dots are coupled to the superconductor. In addition, we observe how the
proximity gap of the dots increases with the coupling strength. Computing the
spectrum of the system for $3$ points in parameter space would require the
same time as the total runtime of Pymablock in this example. This demonstrates
the speed of the implicit method and the efficiency of Pymablock’s algorithm.

##  3 Perturbative block-diagonalization algorithm

###  3.1 Problem statement

Pymablock finds a series of the unitary transformation $\mathcal{U}$ (we use
calligraphic letters to denote series) that block-diagonalizes the Hamiltonian

$\mathcal{H}=H_{0}+\mathcal{H}^{\prime},\quad
H_{0}=\begin{pmatrix}H_{0}^{AA}&0\\\ 0&H_{0}^{BB}\end{pmatrix},$ (6)

with $\mathcal{H}^{\prime}=\mathcal{H}^{\prime}_{D}+\mathcal{H}^{\prime}_{O}$
containing an arbitrary number and orders of perturbations with block-diagonal
and block-offdiagonal components, respectively. The series here may be
multivariate, and they represent sums of the form

$\mathcal{A}=\sum_{n_{1}=0}^{\infty}\sum_{n_{2}=0}^{\infty}\cdots\sum_{n_{k}=0}^{\infty}\lambda_{1}^{n_{1}}\lambda_{2}^{n_{2}}\cdots\lambda_{k}^{n_{k}}A_{n_{1},n_{2},\ldots,n_{k}},$
(7)

where $\lambda_{i}$ are the perturbation parameters and
$A_{n_{1},n_{2},\ldots,n_{k}}$ are linear operators. The problem statement,
therefore, is finding $\mathcal{U}$ and $\tilde{\mathcal{H}}$ such that

$\tilde{\mathcal{H}}=\mathcal{U}^{\dagger}\mathcal{H}\mathcal{U},\quad\tilde{\mathcal{H}}^{AB}=0,\quad\mathcal{U}^{\dagger}\mathcal{U}=1,$
(8)

which is schematically shown in Fig. 1. Series multiply according to the
Cauchy product:

$\mathcal{C}=\mathcal{A}\mathcal{B}\Leftrightarrow
C_{\mathbf{n}}=\sum_{\mathbf{m}+\mathbf{p}=\mathbf{n}}A_{\mathbf{m}}B_{\mathbf{p}}.$

The Cauchy product is the most expensive operation in perturbation theory,
because it involves a large number of multiplications between potentially
large matrices. For example, evaluating $\mathbf{n}$-th order of $\mathcal{C}$
requires $\sim\prod_{i}n_{i}\equiv N$ multiplications of the series
elements.444If both $\mathcal{A}$ and $\mathcal{B}$ are known in advance, fast
Fourier transform-based algorithms can reduce this cost to $\sim N\log N$. In
our problem, however, the series are constructed recursively and therefore
this optimization is impossible. A direct computation of all the possible
index combinations in a product between three series
$\mathcal{A}\mathcal{B}\mathcal{C}$ would have a higher cost $\sim N^{2}$,
however, if we use associativity of the product and compute this as
$(\mathcal{A}\mathcal{B})\mathcal{C}$, then the scaling of the cost stays
$\sim N$. There are many ways to solve the problem (8) that give identical
expressions for $\mathcal{U}$ and $\tilde{\mathcal{H}}$. We are searching for
a procedure that satisfies two additional constraints:

  * •

It has the same complexity scaling as a Cauchy product, and therefore $\sim N$
multiplications per additional order.

  * •

It does not require multiplications by $H_{0}$.

  * •

It requires only one Cauchy product by $\mathcal{H}_{D}$, the block-diagonal
part of $\mathcal{H}$.

The first requirement is that the algorithm scaling is optimal: the desired
expression at least contains a Cauchy product of $\mathcal{U}$ and
$\mathcal{H}$. Therefore the complexity scaling of the complete algorithm may
not become lower than the complexity of a Cauchy product and we aim to reach
this lower bound. The second requirement is because in perturbation theory,
$n$-th order corrections to $\tilde{\mathcal{H}}$ carry $n$ energy
denominators $1/(E_{i}-E_{j})$, where $E_{i}$ and $E_{j}$ are the eigenvalues
of $H_{0}$ belonging to different subspaces. Therefore, any additional
multiplications by $H_{0}$ must cancel with additional energy denominators.
Multiplying by $H_{0}$ is therefore unnecessary work, and it gives longer
intermediate expressions. The third requirement we impose by considering a
case in which $\mathcal{H}_{O}=0$, where $\mathcal{H}_{D}$ must at least enter
$\tilde{\mathcal{H}}$ as an added term, without any products. Moreover,
because $\mathcal{U}$ depends on the entire Hamiltonian, there must be at
least one Cauchy product by $\mathcal{H}^{\prime}_{D}$. The goal of our
algorithm is thus to be efficient and to produce compact results that do not
require further simplifications.

###  3.2 Existing solutions

A common approach to constructing effective Hamiltonians is to use the
Schrieffer–Wolff transformation [1]:

$\tilde{\mathcal{H}}=e^{\mathcal{S}}\mathcal{H}e^{-\mathcal{S}},\\\
e^{\mathcal{S}}=1+\mathcal{S}+\frac{1}{2!}\mathcal{S}\mathcal{S}+\frac{1}{3!}\mathcal{S}\mathcal{S}\mathcal{S}+\cdots,$
(9)

where $\mathcal{S}=\sum_{n}S_{n}$ is an antihermitian polynomial series in the
perturbative parameter, making $e^{\mathcal{S}}$ a unitary transformation.
Requiring that $\tilde{\mathcal{H}}^{AB}=0$ gives a recursive equation for
$S_{n}$, whose terms are nested commutators between the series of
$\mathcal{S}$ and $\mathcal{H}$. Similarly, the transformed Hamiltonian is
given by a series of nested commutators

$\tilde{\mathcal{H}}=\sum_{j=0}^{\infty}\frac{1}{j!}\Big{[}\mathcal{H},\sum_{n=0}^{\infty}S_{n}\Big{]}^{(j)},$
(10)

where the superscript $(j)$ denotes the $j$-th nested commutator
$[A,B]^{(j)}=[[A,B]^{(j-1)},B]$, with $[A,B]^{(0)}=A$ and $[A,B]^{(1)}=AB-BA$.
Regardless of the specific implementation, this expression does not meet
either of our two requirements:

  * •

The direct computation of the series elements requires $\sim\exp N$
multiplications, and even an optimized one has a $\sim N^{2}$ scaling.

  * •

Evaluating Eq. (10) contains multiplications by $H_{0}$.

Alternative parametrizations of the unitary transformation $\mathcal{U}$
require solving unitarity and block diagonalization conditions too, but give
rise to a different recursive procedure for the series elements. For example,
using hyperbolic functions

$\displaystyle\mathcal{U}=\cosh{\mathcal{G}}+\sinh{\mathcal{G}},\quad\mathcal{G}=\sum_{i=0}^{\infty}G_{i},$
(11)

leads to different recursive expressions for $G_{i}$ [8], but does not change
the algorithm’s complexity. On the other hand, using a polynomial series
directly

$\mathcal{U}=\sum_{i=0}^{\infty}U_{i},$ (12)

gives rise to another recursive equation for $U_{i}$ [7, 3, 9, 10]. Still,
this choice results in an expression for $\tilde{\mathcal{H}}$ whose terms
include products by $H_{0}$, and therefore requires additional
simplifications. The following three algorithms satisfy both of our
requirements while solving a related problem. First, density matrix
perturbation theory [11, 29, 12] constructs the density matrix
$\mathcal{\rho}$ of a perturbed system as a power series with respect to a
perturbative parameter:

$\mathcal{\rho}=\sum_{i=0}^{\infty}\rho_{i}.$ (13)

The elements of the series are found by solving two recursive conditions,
$\mathcal{\rho}^{2}=\mathcal{\rho}$ and $[\mathcal{H},\mathcal{\rho}]=0$,
which avoid multiplications by $H_{0}$ and require a single Cauchy product
each. This approach, however, deals with the entire Hilbert space, rather than
the low-energy subspace, and does not provide an effective Hamiltonian.
Second, the perturbative similarity transform by C. Bloch [13, 2] constructs
the effective Hamiltonian in a non-orthogonal basis, which preserves the
Hamiltonian spectrum while breaking its hermiticity. Finally, the recursive
Schrieffer–Wolff algorithm [24] applies the Schrieffer–Wolff transformation to
the output of lower-order iterations, and calculates the effective Hamiltonian
at a fixed perturbation strength, rather than as a series. We thus identify
the following open question: can we construct an effective Hamiltonian with a
linear scaling algorithm that produces compact expressions?

###  3.3 Pymablock’s algorithm

The first idea that Pymablock exploits is the recursive evaluation of the
operator series, which we illustrate by considering the unitarity condition.
Let us separate the transformation $\mathcal{U}$ into an identity and
$\mathcal{U}^{\prime}=\mathcal{W}+\mathcal{V}$:

$\mathcal{U}=1+\mathcal{U}^{\prime}=1+\mathcal{W}+\mathcal{V},\quad\mathcal{W}^{\dagger}=\mathcal{W},\quad\mathcal{V}^{\dagger}=-\mathcal{V}.$
(14)

We use the unitarity condition $\mathcal{U}^{\dagger}\mathcal{U}=1$ by
substituting $\mathcal{U}^{\prime}$ into it:

$\displaystyle
1=(1+\mathcal{U}^{\prime\dagger})(1+\mathcal{U}^{\prime})=1+\mathcal{U}^{\prime\dagger}+\mathcal{U}^{\prime}+\mathcal{U}^{\prime\dagger}\mathcal{U}^{\prime}.$
(15)

This immediately yields

$\displaystyle\mathcal{W}=\frac{1}{2}(\mathcal{U}^{\prime\dagger}+\mathcal{U}^{\prime})=-\frac{1}{2}\mathcal{U}^{\prime\dagger}\mathcal{U}^{\prime}.$
(16)

Because $\mathcal{U}^{\prime}$ has no $0$-th order term,
$(\mathcal{U}^{\prime\dagger}\mathcal{U}^{\prime})_{\mathbf{n}}$ does not
depend on the $\mathbf{n}$-th order of $\mathcal{U}^{\prime}$ nor
$\mathcal{W}$, and therefore Eq. (16) allows to compute $\mathcal{W}$ using
the already available lower orders of $\mathcal{U}^{\prime}$. Alternatively,
using Eq. (14) we could define $\mathcal{W}$ as a Taylor series in
$\mathcal{V}$:

$\mathcal{W}=\sqrt{1+\mathcal{V}^{2}}-1\equiv
f(\mathcal{V})\equiv\sum_{n}a_{n}\mathcal{V}^{2n}.$

A direct computation of all possible products of terms in this expression
requires $\sim\exp N$ multiplications. A more efficient approach for
evaluating this expression introduces each term in the sum as a new series
$\mathcal{A}^{n+1}=\mathcal{A}\mathcal{A}^{n}$ and reuses the previously
computed results. This optimization brings the exponential cost down to $\sim
N^{2}$. However, we see that the Taylor expansion approach is both more
complicated and more computationally expensive than the recurrent definition
in Eq. (16). Therefore, we use Eq. (16) to efficiently compute $\mathcal{W}$.
More generally, a Cauchy product $\mathcal{A}\mathcal{B}$ where $\mathcal{A}$
and $\mathcal{B}$ have no $0$-th order terms depends on
$\mathcal{A}_{1},\ldots,\mathcal{A}_{n-1}$ and
$\mathcal{B}_{1},\ldots,\mathcal{B}_{n-1}$. This makes it possible to use
$\mathcal{AB}$ in a recurrence relation, a property that we exploit throughout
the algorithm. To compute $\mathcal{U}^{\prime}$ we also need to find
$\mathcal{V}$, which is defined by the requirement
$\tilde{\mathcal{H}}^{AB}=0$. Additionally, we constrain $\mathcal{V}$ to be
block off-diagonal: $\mathcal{V}^{AA}=\mathcal{V}^{BB}=0$, a choice we make to
ensure that the resulting unitary transformation is equivalent to the
Schrieffer–Wolff transformation (see section 3.4). In turn, this makes
$\mathcal{W}$ block-diagonal and makes the norm of $\mathcal{U}^{\prime}$
minimal. The remaining condition for finding a recurrent relation for
$\mathcal{U}^{\prime}$ is that the transformed Hamiltonian

$\tilde{\mathcal{H}}=\mathcal{U}^{\dagger}\mathcal{H}\mathcal{U}=\mathcal{H}_{D}+\mathcal{U}^{\prime\dagger}\mathcal{H}_{D}+\mathcal{H}_{D}\mathcal{U}^{\prime}+\mathcal{U}^{\prime\dagger}\mathcal{H}_{D}\mathcal{U}^{\prime}+\mathcal{U}^{\dagger}\mathcal{H}^{\prime}_{O}\mathcal{U},$
(17)

is block-diagonal, a condition that determines $\mathcal{V}$. Here we used
$\mathcal{U}=1+\mathcal{U}^{\prime}$ and
$\mathcal{H}=\mathcal{H}_{D}+\mathcal{H}^{\prime}_{O}$, since $H_{0}$ is
block-diagonal by definition. Because we want to avoid products by
$\mathcal{H}_{D}$, we need to get rid of the terms that contain it by
replacing them with an alternative expression. Our strategy is to define an
auxiliary operator $\mathcal{X}$ that we can compute without ever multiplying
by $\mathcal{H}_{D}$. Like $\mathcal{U}^{\prime}$, $\mathcal{X}$ needs to be
defined via a recurrence relation, which we determine later. Because Eq. (17)
contains $\mathcal{H}_{D}$ multiplied by $\mathcal{U}^{\prime}$ from the left
and from the right, eliminating $\mathcal{H}_{D}$ requires moving it to the
same side. To achieve this, we choose $\mathcal{X}=\mathcal{Y}+\mathcal{Z}$ to
be the commutator between $\mathcal{U}^{\prime}$ and $\mathcal{H}_{D}$:

$\mathcal{X}\equiv[\mathcal{U}^{\prime},\mathcal{H}_{D}]=\mathcal{Y}+\mathcal{Z},\quad\mathcal{Y}\equiv[\mathcal{V},\mathcal{H}_{D}]=\mathcal{Y}^{\dagger},\quad\mathcal{Z}\equiv[\mathcal{W},\mathcal{H}_{D}]=-\mathcal{Z}^{\dagger},$
(18)

where $\mathcal{Y}$ is therefore block off-diagonal and $\mathcal{Z}$, block
diagonal. We use
$\mathcal{H}_{D}\mathcal{U}^{\prime}=\mathcal{U}^{\prime}\mathcal{H}_{D}-\mathcal{X}$
to move $\mathcal{H}_{D}$ through to the right and find

$\displaystyle\tilde{\mathcal{H}}$
$\displaystyle=\mathcal{H}_{D}+\mathcal{U}^{\prime\dagger}\mathcal{H}_{D}+(\mathcal{H}_{D}\mathcal{U}^{\prime})+\mathcal{U}^{\prime\dagger}\mathcal{H}_{D}\mathcal{U}^{\prime}+\mathcal{U}^{\dagger}(\mathcal{H}^{\prime}_{O}\mathcal{U})$
(19)
$\displaystyle=\mathcal{H}_{D}+\mathcal{U}^{\prime\dagger}\mathcal{H}_{D}+\mathcal{U}^{\prime}\mathcal{H}_{D}-\mathcal{X}+\mathcal{U}^{\prime\dagger}(\mathcal{U}^{\prime}\mathcal{H}_{D}-\mathcal{X})+\mathcal{U}^{\dagger}\mathcal{H}^{\prime}_{O}\mathcal{U}$
$\displaystyle=\mathcal{H}_{D}+(\mathcal{U}^{\prime\dagger}+\mathcal{U}^{\prime}+\mathcal{U}^{\prime\dagger}\mathcal{U}^{\prime})\mathcal{H}_{D}-\mathcal{X}-\mathcal{U}^{\prime\dagger}\mathcal{X}+\mathcal{U}^{\dagger}\mathcal{H}^{\prime}_{O}\mathcal{U}$
$\displaystyle=\mathcal{H}_{D}-\mathcal{X}-\mathcal{U}^{\prime\dagger}\mathcal{X}+\mathcal{U}^{\dagger}\mathcal{H}^{\prime}_{O}\mathcal{U},$

where the terms multiplied by $\mathcal{H}_{D}$ cancel according to Eq. (15).
The transformed Hamiltonian does not contain multiplications by
$\mathcal{H}_{D}$ anymore, but it does depend on $\mathcal{X}$, an auxiliary
operator whose recurrent definition we do not know yet. To find it, we first
focus on its anti-Hermitian part, $\mathcal{Z}$. Since recurrence relations
are expressions whose right-hand side contains Cauchy products between series,
we need to find a way to make a product appear. We do so by using the
unitarity condition
$\mathcal{U}^{\prime\dagger}+\mathcal{U}^{\prime}=-\mathcal{U}^{\prime\dagger}\mathcal{U}^{\prime}$
to obtain the recursive definition of $\mathcal{Z}$:

$\displaystyle\mathcal{Z}$
$\displaystyle=\frac{1}{2}(\mathcal{X}-\mathcal{X}^{\dagger})$ (20)
$\displaystyle=\frac{1}{2}\Big{[}(\mathcal{U}^{\prime}+\mathcal{U}^{\prime\dagger})\mathcal{H}_{D}-\mathcal{H}_{D}(\mathcal{U}^{\prime}+\mathcal{U}^{\prime\dagger})\Big{]}$
$\displaystyle=\frac{1}{2}\Big{[}-\mathcal{U}^{\prime\dagger}(\mathcal{U}^{\prime}\mathcal{H}_{D}-\mathcal{H}_{D}\mathcal{U}^{\prime})+(\mathcal{U}^{\prime}\mathcal{H}_{D}-\mathcal{H}_{D}\mathcal{U}^{\prime})^{\dagger}\mathcal{U}^{\prime}\Big{]}$
$\displaystyle=\frac{1}{2}(-\mathcal{U}^{\prime\dagger}\mathcal{X}+\mathcal{X}^{\dagger}\mathcal{U}^{\prime}).$

Similar to computing $W_{\mathbf{n}}$, computing $Z_{\mathbf{n}}$ requires
lower-orders of $\mathcal{X}$ and $\mathcal{U}^{\prime}$. Then, we compute the
Hermitian part of $\mathcal{X}$ by requiring that $\tilde{\mathcal{H}}^{AB}=0$
in the Eq. (19) and find

$\mathcal{X}^{AB}=(\mathcal{U}^{\dagger}\mathcal{H}^{\prime}_{O}\mathcal{U}-\mathcal{U}^{\prime\dagger}\mathcal{X})^{AB}.$
(21)

Once again, despite $\mathcal{X}$ enters the right hand side, because all the
terms lack 0 order, this defines a recursive relation for $\mathcal{X}^{AB}$,
and therefore $\mathcal{Y}$. The final part is straightforward: using
$\mathcal{H}_{D}=H_{0}+\mathcal{H}^{\prime}_{D}$ and the definition of
$\mathcal{Y}$ in Eq. (18) fixes $\mathcal{V}$ as a solution of:

$\mathcal{V}^{AB}H_{0}^{BB}-H_{0}^{AA}\mathcal{V}^{AB}=\mathcal{Y}^{AB}-[\mathcal{V},\mathcal{H}^{\prime}_{D}]^{AB},$
(22)

a Sylvester’s equation, which we only need to solve once for every new order.
In the eigenbasis of $H_{0}$, the solution of Sylvester’s equation is
$V^{AB}_{\mathbf{n},ij}=(\mathcal{Y}-[\mathcal{V},\mathcal{H}^{\prime}_{D}])^{AB}_{\mathbf{n},ij}/(E_{i}-E_{j})$,
where $E_{i}$ are the eigenvalues of $H_{0}$. We now have the complete
algorithm:

  1. 1.

Define series $\mathcal{U}^{\prime}$ and $\mathcal{X}$ and make use of their
block structure and Hermiticity.

  2. 2.

To define the diagonal blocks of $\mathcal{U}^{\prime}$, use
$\mathcal{W}=-\mathcal{U}^{\prime\dagger}\mathcal{U}^{\prime}/2$.

  3. 3.

To find the off-diagonal blocks of $\mathcal{U}^{\prime}$, solve Sylvester’s
equation
$\mathcal{V}^{AB}H_{0}^{BB}-H_{0}^{AA}\mathcal{V}^{AB}=\mathcal{Y}^{AB}-[\mathcal{V},\mathcal{H}^{\prime}_{D}]^{AB}$.
This requires $\mathcal{X}$.

  4. 4.

To find the diagonal blocks of $\mathcal{X}$, define
$\mathcal{Z}=(-\mathcal{U}^{\prime\dagger}\mathcal{X}+\mathcal{X}^{\dagger}\mathcal{U}^{\prime})/2$.

  5. 5.

For the off-diagonal blocks of $\mathcal{X}$, use
$\mathcal{Y}^{AB}=(-\mathcal{U}^{\prime\dagger}\mathcal{X}+\mathcal{U}^{\dagger}\mathcal{H}^{\prime}\mathcal{U})^{AB}$.

  6. 6.

Compute the effective Hamiltonian as
$\tilde{\mathcal{H}}=\mathcal{H}_{D}-\mathcal{X}-\mathcal{U}^{\prime\dagger}\mathcal{X}+\mathcal{U}^{\dagger}\mathcal{H}^{\prime}_{O}\mathcal{U}$.

###  3.4 Equivalence to Schrieffer–Wolff transformation

Pymablock’s algorithm and the Schrieffer–Wolff transformation both find a
unitary transformation $\mathcal{U}$ such that $\tilde{\mathcal{H}}^{AB}=0$.
They are therefore equivalent up to a gauge choice in each subspace, $A$ and
$B$. We establish the equivalence between the two by demonstrating that this
gauge choice is the same for both algorithms. The Schrieffer–Wolff
transformation uses $\mathcal{U}=\exp\mathcal{S}$, where
$\mathcal{S}=-\mathcal{S}^{\dagger}$ and
$\mathcal{S}^{AA}=\mathcal{S}^{BB}=0$, this restriction makes the result
unique [2]. On the other hand, our algorithm produces the unique block-
diagonalizing transformation with a block structure
$\mathcal{U}_{AA}=\mathcal{U}_{AA}^{{\dagger}}$,
$\mathcal{U}_{BB}=\mathcal{U}_{BB}^{{\dagger}}$ and
$\mathcal{U}_{AB}=-\mathcal{U}_{BA}^{{\dagger}}$. The uniqueness is a
consequence of the construction of the algorithm, where calculating every
order gives a unique solution satisfying these conditions. To see that the two
solutions are identical, we expand the Taylor series of $\exp\mathcal{S}$.
Every even order gives a Hermitian, block-diagonal matrix, while every odd
order gives an anti-Hermitian block off-diagonal matrix, showing that
$\exp\mathcal{S}$ has the same structure as $\mathcal{U}$ above. The reverse
statement about the structure of $\log\mathcal{U}$ can be seen similarly,
using the Taylor series of the logarithm around $1$. Using a series expansion
is justified by the perturbative nature of the result, meaning that
$\mathcal{S}$ is close to $0$ and $\mathcal{U}$ is close to $1$. Because of
the uniqueness of both results, we find that $\exp\mathcal{S}$ from
conventional Schrieffer–Wolff transformation is identical to $\mathcal{U}$
found by our algorithm, which remains true if both power series are truncated
at a finite order.

###  3.5 Extra optimization: common subexpression elimination

While the algorithm of Sec. 3.3 satisfies our requirements, we improve it
further by reusing products that are needed in several places, such that the
total number of matrix multiplications is reduced. Firstly, we rewrite the
expressions for $\mathcal{Z}$ in Eq. (20) and $\tilde{\mathcal{H}}$ in Eq.
(19) by utilizing the Hermitian conjugate of
$\mathcal{U}^{\prime\dagger}\mathcal{X}$ without recomputing it:

$\displaystyle\mathcal{Z}=\frac{1}{2}\left[(-\mathcal{U}^{\prime\dagger}\mathcal{X})-\textrm{h.c.}\right],$
$\displaystyle\tilde{\mathcal{H}}=\mathcal{H}_{D}+\mathcal{U}^{\dagger}\mathcal{H}^{\prime}_{O}\mathcal{U}-(\mathcal{U}^{\prime\dagger}\mathcal{X}+\textrm{h.c.})/2,$

where h.c. is the Hermitian conjugate, and $\mathcal{X}$ drops out from the
diagonal blocks of $\tilde{\mathcal{H}}$ because the diagonal blocks of
$\mathcal{X}$ are anti-Hermitian. Additionally, we reuse the repeated
$\mathcal{A}\equiv\mathcal{H}^{\prime}_{O}\mathcal{U}^{\prime}$ in

$\mathcal{U}^{\dagger}\mathcal{H}^{\prime}_{O}\mathcal{U}=\mathcal{H}^{\prime}_{O}+\mathcal{A}+\mathcal{A}^{\dagger}+\mathcal{U}^{\prime\dagger}\mathcal{A}.$
(23)

Next, we observe that some products from the
$\mathcal{U}^{\dagger}\mathcal{H}_{O}\mathcal{U}$ term appear both in
$\mathcal{X}$ in Eq. (21) and in $\tilde{\mathcal{H}}$ (19). To avoid
recomputing these products, we introduce
$\mathcal{B}=\mathcal{X}-\mathcal{H}^{\prime}_{O}-\mathcal{A}$ and define the
recursive algorithm using $\mathcal{B}$ instead of $\mathcal{X}$. With this
definition, we compute the off-diagonal blocks of $\mathcal{B}$ as:

$\displaystyle\mathcal{B}^{AB,BA}$
$\displaystyle=\left[\mathcal{X}-\mathcal{H}^{\prime}_{O}-\mathcal{A}\right]^{AB,BA}$
(24)
$\displaystyle=\left[\mathcal{A}^{\dagger}+\mathcal{U}^{\prime\dagger}\mathcal{A}-\mathcal{U}^{\prime\dagger}\mathcal{X}\right]^{AB,BA}$
$\displaystyle=\left[\mathcal{U}^{\prime\dagger}\mathcal{H}^{\prime}_{O}+\mathcal{U}^{\prime\dagger}\mathcal{A}-\mathcal{U}^{\prime\dagger}\mathcal{X}\right]^{AB,BA}$
$\displaystyle=-(\mathcal{U^{\prime}}^{\dagger}\mathcal{B})^{AB,BA},$

where we also used Eq. (21) and the definition of $\mathcal{A}$. The diagonal
blocks of $\mathcal{B}$, on the other hand, are given by

$\displaystyle\mathcal{B}^{AA,BB}$
$\displaystyle=\left[\mathcal{X}-\mathcal{H}^{\prime}_{O}-\mathcal{A}\right]^{AA,BB}$
(25)
$\displaystyle=\left[\frac{1}{2}[(-\mathcal{U}^{\prime\dagger}\mathcal{X})-\textrm{h.c.}]-\mathcal{A}\right]^{AA,BB}$
$\displaystyle=\left[\frac{1}{2}[(-\mathcal{U}^{\prime\dagger}[\mathcal{X}-\mathcal{H}^{\prime}_{O}-\mathcal{A}])-\textrm{h.c.}]-\frac{1}{2}[\mathcal{A}^{\dagger}+\mathcal{A}]+{\frac{1}{2}[(-\mathcal{U}^{\prime\dagger}\mathcal{A})-\textrm{h.c.}]}\right]^{AA,BB},$
$\displaystyle=\left[\frac{1}{2}[(-\mathcal{U}^{\prime\dagger}\mathcal{B})-\textrm{h.c.}]-\frac{1}{2}[\mathcal{A}^{\dagger}+\textrm{h.c.}]\right]^{AA,BB},$

where we used Eq. (20) and that $\mathcal{U}^{\prime\dagger}\mathcal{A}$ is
Hermitian. Using $\mathcal{B}$ changes the relation for $\mathcal{V}^{AB}$ in
Eq. (22) to

$\mathcal{V}^{AB}H_{0}^{BB}-H_{0}^{AA}\mathcal{V}^{AB}=\left(\mathcal{B}-\mathcal{H}^{\prime}-\mathcal{A}-[\mathcal{V},\mathcal{H}^{\prime}_{D}]\right)^{AB}.$
(26)

Finally, we combine Eq. (19), Eq. (23), Eq. (25) and Eq. (24) to obtain the
final expression for the effective Hamiltonian:

$\tilde{\mathcal{H}}_{D}=\mathcal{H}_{D}+(\mathcal{A}+\textrm{h.c.})/2-(\mathcal{U}^{\dagger}\mathcal{B}+\textrm{h.c.})/2.$
(27)

Together with the series $\mathcal{U}^{\prime}$ in Eqs. (16,26),
$\mathcal{A}=\mathcal{H}^{\prime}_{O}\mathcal{U}^{\prime}$, and $\mathcal{B}$
in Eqs. (25,24), this equation defines the optimized algorithm.

##  4 Implementation

###  4.1 The data structure for block operator series

The optimized algorithm from the previous section requires constructing $14$
operator series, whose elements are computed using a collection of recurrence
relations. This warrants defining a specialized data structure suitable for
this task that represents a multidimensional series of operators. Because the
recurrent relations are block-wise, the data structure needs to keep track of
separate blocks. In order to support varied use cases, the actual
representation of the operators needs to be flexible: the block may be dense
arrays, sparse matrices, symbolic expressions, or more generally any object
that defines addition and multiplication. Finally, the series needs to be
queryable by order and block, so that it supports a block-wise multivariate
Cauchy product—the main operation in the algorithm. The most straightforward
way to implement a perturbation theory calculation is to write a function that
has the desired order as an argument, computes the series up to that order,
and returns the result. This makes it hard to reuse already computed terms for
a new computation, and becomes complicated to implement in the
multidimensional case when different orders in different perturbations are
needed. We find that a recursive approach addresses these issues: within this
paradigm, each series needs to define how its entries depend on lower-order
terms. To address these requirements, we define a `BlockSeries` Python class
and use it to represent the series of $\mathcal{U}$, $\mathcal{H}$, and
$\tilde{\mathcal{H}}$, as well as the intermediate series used to define the
algorithm. The objects of this class are equipped with a function to compute
their elements and it stores the already computed results in a dictionary.
Storing the results for reuse is necessary to optimize the evaluation of
higher order terms and it allows to request additional orders without
restarting the computation. For example, the definition of the `BlockSeries`
for $\tilde{\mathcal{H}}$ has the following form:

⬇

H_tilde = BlockSeries(

shape=(2, 2), # 2x2 block matrix

n_infinite=n, # number of perturbative parameters

eval=compute_H_tilde, # function to compute the elements

name=”H_tilde”,

dimension_names=(”lambda”, …), # parameter names

)

Here `compute_H_tilde` is a function implementing Eq. (27) by querying other
series objects. Calling `H_tilde[0, 0, 2]`, the second order perturbation
$\sim\lambda^{2}$ of the $AA$ block, then does the following:

  1. 1.

Evaluates `compute_H_tilde(0, 0, 2)` if it is not already computed.

  2. 2.

Stores the evaluation result in a dictionary.

  3. 3.

Returns the result.

To conveniently access multiple orders at once, we implement NumPy array
indexing so that `H_tilde[0, 0, :3]` returns a NumPy masked array array with
the orders $\sim\lambda^{0}$ , $\sim\lambda^{1}$, and $\sim\lambda^{2}$ of the
$AA$ block. The masking allows to support a common use case where some orders
of a series are zero, so that they are omitted from the computations. We
expect that the `BlockSeries` data structure is suitable to represent a broad
class of perturbative calculations, and we plan to extend it to support more
advanced features in the future. We utilize `BlockSeries` to implement
multiple other optimizations. For example, we exploit Hermiticity when
computing the Cauchy product of $U^{\prime\dagger}U^{\prime}$ in Eq. (16), by
only evaluating half of the matrix products, and then complex conjugate the
result to obtain the rest. Similarly, for Hermitian and anti-Hermitian series,
like the off-diagonal blocks of $\mathcal{U}^{\prime}$, we only compute the
$AB$ blocks, and use the conjugate transpose to obtain the $BA$ blocks. This
approach should also allow us to implement efficient handling of symmetry-
constrained Hamiltonians, where some blocks either vanish or are equal to
other blocks due to a symmetry. Moreover, using `BlockSeries` with custom
objects yields additional information about the algorithm and accommodates its
further development. Specifically, we have used a custom object with a counter
to measure the algorithm complexity (see also Sec. 5) and to determine which
results are only used once so that they can be immediately discarded from
storage.

###  4.2 The implicit method for large sparse Hamiltonians

A distinguishing feature of Pymablock is its ability to handle large sparse
Hamiltonians, that are too costly to diagonalize, as illustrated in Sec. 2.3.
Specifically, we consider the situations when the size $N_{A}$ of the $A$
subspace is small compared to the entire Hilbert space, so that obtaining the
basis $\Psi_{A}$ of the $A$ subspace is feasible using sparse diagonalization.
The projector on the $A$ subspace $P_{A}=\Psi_{A}^{\dagger}\Psi_{A}$ is then a
low-rank matrix, a property that we exploit to avoid constructing the $B$
subspace explicitly. Furthermore, the solution of Sylvester’s equation in Eq.
22 amounts to multiplying $N_{A}$ large vectors, rows of
$Y_{\mathbf{n}}^{AB}$, by the energy denominators $E_{i}-E_{j}$, where $E_{i}$
are the $N_{A}$ eigenvalues of the $A$ subspace provided by sparse
diagonalization. The key tool to solve this problem is the projector approach
introduced in Ref. [30], which introduces an equivalent extended Hamiltonian
using the projector $P_{B}=1-P_{A}$ onto the $B$ subspace:

$\displaystyle\bar{\mathcal{H}}=\begin{pmatrix}\Psi_{A}^{\dagger}\mathcal{H}\Psi_{A}&\Psi_{A}^{\dagger}\mathcal{H}P_{B}\\\
P_{B}\mathcal{H}\Psi_{A}&P_{B}\mathcal{H}P_{B}\end{pmatrix}.$ (28)

In other words, the subspace $\bar{A}$ is written in the basis of $\Psi_{A}$,
while the basis of the $\bar{B}$ subspace is the same as the original complete
basis of $\mathcal{H}$ to preserve its sparsity. We also project out the
$A$-degrees of freedom from the $\bar{B}$ subspace to avoid duplicate
solutions in $\bar{\mathcal{H}}$, which introduces $N_{A}$ eigenvectors with
zero eigenvalues. Introducing $\bar{\mathcal{H}}$ allows to multiply by
operators of a form $P_{B}H_{\mathbf{n}}P_{B}$ efficiently by using the low-
rank structure of $P_{A}$. In the code we represent the $\bar{B}\bar{B}$
operators as `LinearOperator` objects from the SciPy package [28], enabled by
the ability of the `BlockSeries` to store arbitrary objects. Storing the
$\bar{A}\bar{A}$ and $\bar{A}\bar{B}$ blocks as dense matrices—efficient
because these are small and dense—finishes the implementation of the
Hamiltonian. To solve the Sylvester’s equation we write it for every row of
$V_{\mathbf{n}}^{\bar{A}\bar{B}}$ separately:

$\displaystyle
V_{\mathbf{n},ij}^{\bar{A}\bar{B}}(E_{i}-H_{0})=Y^{\bar{A}\bar{B}}_{\mathbf{n},j}$
(29)

This equation has a solution despite $E_{i}-H_{0}$ not being invertible
because $Y^{\bar{A}\bar{B}}_{\mathbf{n}}P_{A}=0$. We solve this equation using
the MUMPS sparse solver [31, 32], which prepares an efficient sparse LU-
decomposition of $E_{i}-H_{0}$, or the KPM approximation of the Green’s
function [33]. Both methods work on sparse Hamiltonians with millions of
degrees of freedom.

##  5 Benchmark

To the best of our knowledge, there are no other packages implementing
arbitrary order quasi-degenerate perturbation theory. Literature references
provide explicit expressions for the effective Hamiltonian up to fourth order,
together with the procedure for obtaining higher order expressions [34].
Because the full reference expressions are lengthy, we do not provide them,
but for example at $4$-th order the effective Hamiltonian is a sum of several
expressions of the form

$\sum_{m^{{}^{\prime\prime}}m^{{}^{\prime\prime\prime}}l}\frac{H^{\prime}_{mm^{{}^{\prime\prime}}}H^{\prime}_{m^{{}^{\prime\prime}}m^{{}^{\prime\prime\prime}}}H^{\prime}_{m^{{}^{\prime\prime\prime}}l}H^{\prime}_{lm^{{}^{\prime}}}}{(E_{m^{{}^{\prime\prime}}}-E_{l})(E_{m^{{}^{\prime\prime\prime}}}-E_{l})(E_{m}-E_{l})},$
(30)

where the $m$-indices label states from the $A$-subspace and $l$-indices label
the states from the $B$-subspace. More generally, at $n$-th order each term is
a product of $n$ matrix elements of the Hamiltonian and $n-1$ energy
denominators. Directly carrying out the summation over all the states requires
$\mathcal{O}(N_{A}^{2}N_{B}^{n-1})$ operations, where $N_{A}$ and $N_{B}$ are
the number of states in the two subspaces. In other words, the direct
computation scales worse than a matrix product with the problem size.
Formulating Eq. (30) as $n-1$ matrix products combined with $n-1$ solutions of
Sylvester’s equation, brings this complexity down to $\mathcal{O}((n-1)\times
N_{A}N_{B}^{2})$. This optimization, together with the hermiticity of the sum,
allows us to evaluate the reference expressions for the effective Hamiltonian
for $2$-nd, $3$-rd, and $4$-th order using $1$, $4$, and $27$ matrix products,
respectively. Pymablock’s algorithm utilizes $1$, $3$, and $14$, matrix
products to obtain the same orders of the effective Hamiltonian. Its advantage
becomes even more pronounced at higher orders due to the exponential growth of
the number of terms in the reference expressions. While finding the optimized
implementation from the reference expressions is possible for the $3$-rd
order, we expect it to be extremely challenging for the $4$-th order, and
essentially impossible to do manually for higher orders. Moreover, because the
`BlockSeries` class tracks absent terms, in practice the number of matrix
products depends on the sparsity of the block structure of the perturbation,
as shown in Fig. 4.

Figure 4:  Matrix products required to compute $\tilde{H}^{AA}_{n}$ for a
dense and block off-diagonal first-order perturbation (left) and a dense and
block off-diagonal perturbative series with terms of all orders present
(right).

The efficiency of Pymablock becomes especially apparent when applied to sparse
numerical problems, similar to Sec. 2.3. We demonstrate the performance of the
implicit method by using it to compute the low-energy spectrum of a large
tight-binding model, and comparing Pymablock’s time cost to that of sparse
diagonalization. We define a 2D square lattice of $52\times 52$ sites with
nearest-neighbor hopping and a random onsite potential $\mu(\mathbf{r})$. The
perturbation $\delta\mu(\mathbf{r})$ interpolates between two different
disorder realizations. For the sake of an illustration, we choose the system’s
parameters such that the dispersion of the lowest few levels with $\delta\mu$
features avoided crossings and an overall nonlinear shape, whose details are
not relevant. Similar to Sec. 2.3, constructing the effective Hamiltonian
involves three steps. First, we compute the $10$ lowest states of the
unperturbed Hamiltonian using sparse diagonalization. Second,
`block_diagonalize` computes a sparse LU decomposition of the Hamiltonian at
each of the $10$ eigenenergies. Third, we compute corrections $\tilde{H}_{1}$,
$\tilde{H}_{2}$, and $\tilde{H}_{3}$ to the effective Hamiltonian, each being
a $10\times 10$ matrix. Each of these steps is a one-time cost, see Fig. 5.
Finally, to compare the perturbative calculation to sparse diagonalization, we
construct the effective Hamiltonian
$\tilde{H}=H_{0}+\delta\mu\tilde{H}_{1}+\delta\mu^{2}\tilde{H}_{2}+\delta\mu^{3}\tilde{H}_{3}$
and diagonalize it to obtain the low-energy spectrum for a range of
$\delta\mu$. This has a negligible cost compared to constructing the series.
The comparison is shown in Fig. 5. We observe that while the second order
results are already very close to the exact spectrum, the third order
corrections fully reproduce the sparse diagonalization. At the same time, the
entire cost of computing the perturbative band structure for a range of
$\delta\mu$ is lower than computing a single additional sparse
diagonalization.

Figure 5:  Top panels: band structure of the perturbative effective
Hamiltonian (black) of a tight-binding model compared to exact sparse
diagonalization (gray). Bottom panel: a comparison of the Pymablock’s time
cost with sparse diagonalization. Most of the time is spent in the LU
decomposition of the Hamiltonian (red). The entire cost of the implicit method
is lower than a single sparse diagonalization (gray). The operations of
negligible cost are not shown. The bars length corresponds to the average time
cost over $40$ runs, and the error bars show the standard deviation.

##  6 Conclusion

We developed an algorithm for constructing an effective Hamiltonian that
combines advantages of different perturbative expansions. The main building
block of our approach is a set of recurrence relations that define several
series that depend on each other and combine into the effective Hamiltonian.
Our algorithm constructs the same effective Hamiltonians as the
Schrieffer–Wolff transformation [1], while keeping the linear scaling per
extra order similar to the density matrix perturbation theory [11, 12] or the
non-orthogonal perturbation theory [13]. Its expressions minimize the number
of matrix multiplications per order, making it appealing both for symbolic and
numerical computations. We provide a Python implementation of the algorithm in
the Pymablock package [35]. The package is thoroughly tested (94% test
coverage as of version 2.0), becoming a reliable tool for constructing
effective Hamiltonians that combine multiple perturbations to high orders. The
core of the Pymablock interface is the `BlockSeries` class that handles
arbitrary objects as long as they support algebraic operations. This enables
Pymablock’s construction of effective models for large tight-binding models
using its implicit method. It also allows Pymablock to solve both symbolic and
numerical problems in diverse physical settings, and potentially to
incorporate it into existing packages, such as scqubits [22], QuTiP [36, 37],
or dft2kp [38]. Beyond the Schrieffer–Wolff transformation, the Pymablock
package provides a foundation for defining other perturbative expansions. We
anticipate extending it to time-dependent problems, where the different
regimes of the time-dependent drive modify the recurrence relations that need
to be solved. Applying the same framework to problems with weak position
dependence would allow to construct a nonlinear response theory of quantum
materials. Finally, we expect that in the many-particle context the same
framework supports implementing different flavors of diagrammatic expansions.

## Acknowledgements

We thank Valla Fatemi and Antonio Manesco for feedback on the manuscript.

## Data availability

The code used to produce the reported results is available on Zenodo [35].

#### Author contributions

A. R. A. had the initial idea and oversaw the project. All authors developed
the algorithm. I. A. D., S. M., H. K. K, and A. R. A. wrote the package. I. A.
D. and A. R. A. wrote the paper.

#### Funding information

This research was supported by the Netherlands Organization for Scientific
Research (NWO/OCW) as part of the Frontiers of Nanoscience program, a NWO VIDI
grant 016.Vidi.189.180, and OCENW.GROOT.2019.004. D.V. acknowledges funding
from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)
under Germany’s Excellence Strategy through the Würzburg-Dresden Cluster of
Excellence on Complexity and Topology in Quantum Matter – ct.qmat (EXC 2147,
project-ids 390858490 and 392019).

## Glossary

## References

  * [1] J. R. Schrieffer and P. A. Wolff,  _Relation between the anderson and kondo hamiltonians_ ,  Phys. Rev. 149, 491 (1966),  10.1103/PhysRev.149.491.
  * [2] S. Bravyi, D. P. DiVincenzo and D. Loss,  _Schrieffer–wolff transformation for quantum many-body systems_ ,  Annals of Physics 326(10), 2793 (2011),  10.1016/j.aop.2011.06.004.
  * [3] P. Löwdin,  _Studies in Perturbation Theory. IV. Solution of Eigenvalue Problem by Projection Operator Formalism_ ,  Journal of Mathematical Physics 3(5), 969 (1962),  10.1063/1.1724312.
  * [4] J. M. Luttinger and W. Kohn,  _Motion of electrons and holes in perturbed periodic fields_ ,  Phys. Rev. 97, 869 (1955),  10.1103/PhysRev.97.869.
  * [5] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson and W. D. Oliver,  _A quantum engineer’s guide to superconducting qubits_ ,  Applied physics reviews 6(2) (2019),  10.1063/1.5089550.
  * [6] J. Romhányi, G. Burkard and A. Pályi,  _Subharmonic transitions and bloch-siegert shift in electrically driven spin resonance_ ,  Phys. Rev. B 92, 054422 (2015),  10.1103/PhysRevB.92.054422.
  * [7] J. H. Van Vleck,  _On $\sigma$-type doubling and electron spin in the spectra of diatomic molecules_,  Phys. Rev. 33, 467 (1929),  10.1103/PhysRev.33.467.
  * [8] I. Shavitt and L. T. Redmon,  _Quasidegenerate perturbation theories. A canonical van Vleck formalism and its relationship to other approaches_ ,  The Journal of Chemical Physics 73(11), 5711 (1980),  10.1063/1.440050.
  * [9] D. J. Klein,  _Degenerate perturbation theory_ ,  The Journal of Chemical Physics 61(3), 786 (1974),  10.1063/1.1682018.
  * [10] K. Suzuki and R. Okamoto,  _Degenerate Perturbation Theory in Quantum Mechanics_ ,  Progress of Theoretical Physics 70(2), 439 (1983),  10.1143/PTP.70.439.
  * [11] R. McWeeny,  _Perturbation theory for the fock-dirac density matrix_ ,  Phys. Rev. 126, 1028 (1962),  10.1103/PhysRev.126.1028.
  * [12] L. A. Truflandier, R. M. Dianzinga and D. R. Bowler,  _Notes on density matrix perturbation theory_ ,  The Journal of Chemical Physics 153(16), 164105 (2020),  10.1063/5.0022244.
  * [13] C. Bloch,  _Sur la théorie des perturbations des états liés_ ,  Nuclear Physics 6, 329 (1958),  10.1016/0029-5582(58)90116-0.
  * [14] J. Wurtz, P. W. Claeys and A. Polkovnikov,  _Variational schrieffer-wolff transformations for quantum many-body dynamics_ ,  Phys. Rev. B 101, 014302 (2020),  10.1103/PhysRevB.101.014302.
  * [15] Z. Zhang, Y. Yang, X. Xu and Y. Li,  _Quantum algorithms for schrieffer-wolff transformation_ ,  Phys. Rev. Res. 4, 043023 (2022),  10.1103/PhysRevResearch.4.043023.
  * [16] E. McCann and M. Koshino,  _The electronic properties of bilayer graphene_ ,  Reports on Progress in Physics 76(5), 056503 (2013),  10.1088/0034-4885/76/5/056503.
  * [17] A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin and R. J. Schoelkopf,  _Cavity quantum electrodynamics for superconducting electrical circuits: An architecture for quantum computation_ ,  Phys. Rev. A 69, 062320 (2004),  10.1103/PhysRevA.69.062320.
  * [18] G. Zhu, D. G. Ferguson, V. E. Manucharyan and J. Koch,  _Circuit qed with fluxonium qubits: Theory of the dispersive regime_ ,  Phys. Rev. B 87, 024510 (2013),  10.1103/PhysRevB.87.024510.
  * [19] X. Li, T. Cai, H. Yan, Z. Wang, X. Pan, Y. Ma, W. Cai, J. Han, Z. Hua, X. Han, Y. Wu, H. Zhang _et al._ ,  _Tunable coupler for realizing a controlled-phase gate with dynamically decoupled regime in a superconducting circuit_ ,  Phys. Rev. Appl. 14, 024070 (2020),  10.1103/PhysRevApplied.14.024070.
  * [20] A. Blais, A. L. Grimsmo, S. M. Girvin and A. Wallraff,  _Circuit quantum electrodynamics_ ,  Rev. Mod. Phys. 93, 025005 (2021),  10.1103/RevModPhys.93.025005.
  * [21] E. A. Sete, A. Q. Chen, R. Manenti, S. Kulshreshtha and S. Poletto,  _Floating tunable coupler for scalable quantum computing architectures_ ,  Phys. Rev. Appl. 15, 064063 (2021),  10.1103/PhysRevApplied.15.064063.
  * [22] P. Groszkowski and J. Koch,  _Scqubits: a Python package for superconducting qubits_ ,  Quantum 5, 583 (2021),  10.22331/q-2021-11-17-583.
  * [23] S. P. Chitta, T. Zhao, Z. Huang, I. Mondragon-Shem and J. Koch,  _Computer-aided quantization and numerical analysis of superconducting circuits_ ,  New Journal of Physics 24(10), 103020 (2022),  10.1088/1367-2630/ac94f2.
  * [24] B. Li, T. Calarco and F. Motzoi,  _Nonperturbative analytical diagonalization of hamiltonians with application to circuit qed_ ,  PRX Quantum 3, 030313 (2022),  10.1103/PRXQuantum.3.030313.
  * [25] A. Melo, T. Tanev and A. R. Akhmerov,  _Greedy optimization of the geometry of Majorana Josephson junctions_ ,  SciPost Phys. 14, 047 (2023),  10.21468/SciPostPhys.14.3.047.
  * [26] A. Meurer, C. P. Smith, M. Paprocki, O. Čertík, S. B. Kirpichev, M. Rocklin, A. Kumar, S. Ivanov, J. K. Moore, S. Singh, T. Rathnayake, S. Vig _et al._ ,  _Sympy: symbolic computing in python_ ,  PeerJ Computer Science 3, e103 (2017),  10.7717/peerj-cs.103.
  * [27] C. W. Groth, M. Wimmer, A. R. Akhmerov and X. Waintal,  _Kwant: a software package for quantum transport_ ,  New Journal of Physics 16(6), 063065 (2014),  10.1088/1367-2630/16/6/063065.
  * [28] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright _et al._ ,  _Scipy 1.0: fundamental algorithms for scientific computing in python_ ,  Nature methods 17(3), 261 (2020),  10.1038/s41592-019-0686-2.
  * [29] R. McWeeny,  _Self-consistent perturbation thèory_ ,  Chemical Physics Letters 1(12), 567 (1968),  10.1016/0009-2614(68)85047-X.
  * [30] M. Irfan, S. R. Kuppuswamy, D. Varjas, P. M. Perez-Piskunow, R. Skolasinski, M. Wimmer and A. R. Akhmerov,  _Hybrid kernel polynomial method_ ,  10.48550/arXiv.1909.09649 (2019).
  * [31] P. R. Amestoy, I. S. Duff, J.-Y. L’Excellent and J. Koster,  _A fully asynchronous multifrontal solver using distributed dynamic scheduling_ ,  SIAM Journal on Matrix Analysis and Applications 23(1), 15 (2001),  10.1137/S0895479899358194.
  * [32] P. R. Amestoy, A. Guermouche, J.-Y. L’Excellent and S. Pralet,  _Hybrid scheduling for the parallel solution of linear systems_ ,  Parallel Computing 32(2), 136 (2006),  10.1016/j.parco.2005.07.004,  Parallel Matrix Algorithms and Applications (PMAA’04).
  * [33] A. Weiße, G. Wellein, A. Alvermann and H. Fehske,  _The kernel polynomial method_ ,  Rev. Mod. Phys. 78, 275 (2006),  10.1103/RevModPhys.78.275.
  * [34] W. Roland,  _Spin-orbit coupling effects in two-dimensional electron and hole systems_ ,  Springer Tracts in Modern Physiscs: Springer, Berlin, Heidelberg 191 (2003),  10.1007/b13586.
  * [35] I. Araya Day, S. Miles, H. K. Kerstens, D. Varjas and A. R. Akhmerov,  _Pymablock_ ,  10.5281/zenodo.10629243 (2024).
  * [36] J. Johansson, P. Nation and F. Nori,  _Qutip: An open-source python framework for the dynamics of open quantum systems_ ,  Computer Physics Communications 183(8), 1760 (2012),  10.1016/j.cpc.2012.02.021.
  * [37] J. Johansson, P. Nation and F. Nori,  _Qutip 2: A python framework for the dynamics of open quantum systems_ ,  Computer Physics Communications 184(4), 1234 (2013),  10.1016/j.cpc.2012.11.019.
  * [38] J. V. V. Cassiano, A. de Lelis Araújo, P. E. F. Junior and G. J. Ferreira,  _DFT2kp: Effective kp models from ab-initio data_ ,  SciPost Phys. Codebases p. 25 (2024),  10.21468/SciPostPhysCodeb.25.