Introduction

HQS Quantum Solver enables you to run precise simulations of quantum systems by accounting for all interactions between the particles in a system.

HQS Quantum Solver gives you easy access to high-performance routines and solvers needed to work with quantum systems.

Providing a consistent and extensible interface, it offers routines for the construction of Hamiltonians and operators in either full configurational space or subspaces characterized by conserved quantities such as the spin or the particle number.

Interested in high-quality time evolution? We offer an optimized Krylov implementation of the time-evolution propagator. By using Krylov methods, we target exactly the space of interest and preserve the unitarity of the wavefunction at a fraction of the cost of conventional methods. See the theory section for more details.

In addition to providing a framework, multiple backends have been implemented for the efficient evaluation of properties depending on the user's needs.

We provide a consistent protocol allowing for user-driven extensibility; users can write their own backends, tools, and algorithms and interface with ours.

Installation

HQS Quantum Solver

Visit cloud.quantumsimulations.de/software, search for "HQS Quantum Solver" in the package list and select "Download". Follow the instructions on that page to install HQS Quantum Solver.

To use HQS Quantum Solver, you furthermore need to install the Intel Math Kernel Library (MKL), which is discussed in the next section.

MKL

If you want, you can manually provide a version of the MKL by making sure that the file libmkl_rt.so is found by the dynamic linker. That means, that a system-wide installation should be found automatically.

If you use HQStage to install HQS Quantum Solver, you can install the MKL into an HQStage managed environment by running the following command.

hqstage envs install-mkl

When using your own virtual environment you can install the MKL via pip. However, you need to create a symlink for the file libmkl_rt.so. The commands below perform the necessary steps.

pip install mkl
PLATLIB_DIR="$(python -c "import sysconfig; print(sysconfig.get_path('platlib'))")"
ln -fs libmkl_rt.so.2 "${PLATLIB_DIR}/../../libmkl_rt.so"

Getting Started

Since quantum mechanical systems are described by operators acting on state vectors living in Hilbert spaces, this library is mainly concerned with providing such operators and functions that manipulate state vectors. In this chapter, you will learn the essentials of using the HQS Quantum Solver software by considering an example computation.

In this and the following examples, we shall consider fermion systems, but other types of systems are available as well.

Let us consider the tight-binding Hamiltonian defined on a chain lattice, i.e., the Hamiltonian $H$ is given by $$H=-t\sum _{j=0}^{M-2}\left(c_j^{†}{c}_{j+1}+c_{j+1}^{†}{c}_j\right),$$ where $M$ is the number of lattice sites, $t$ the strength of the hopping term and $c_j$ the annihilation operator at site $j$.

Imports and Parameters

First, we need to import the functions, classes and modules that are required for the script. If you want to follow along, copy the following imports into your script file. We will explain the classes and functions as we go along.

# Title    : HQS Quantum Solver Getting Started Guide
# Filename : getting_started.py

# Python, NumPy, SciPy, and Matplotlib imports
from math import pi, sqrt
import numpy as np
from scipy.linalg import toeplitz
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt

# HQS Quantum Solver imports
from hqs_quantum_solver.spinless_fermions import VectorSpace, Operator, hopping, annihilation
from hqs_quantum_solver.util import unit_vector

Next, we define the parameters that we will refer to in the script.

M = 16  # site count
N = 1   # particle count
E = 0   # energy level

With that out of the way, we can start setting up the tight-binding Hamiltonian.

Setting up the Hamiltonian

We want to look at a system given by the tight-binding Hamiltonian. Thus, the first step is to construct said Hamiltonian. Quantum Solver provides a generic interface for constructing quantum mechanical operators, as well as several convenience functions.

The first step in constructing an operator is to define the vector space the operator will be acting on. In quantum mechanics we consider vector spaces where each dimension corresponds to a state of the system. Hence, to define a Quantum Solver vector space, we need to specify the states of the system. For this particular example, we want the states that have exactly N particles distributed across M sites. The corresponding vector space is constructed by creating an instance of the VectorSpace in the following way.

v = VectorSpace(sites=M, particle_numbers=N)

Note, that Quantum Solver is also able to create vector spaces for a varying number of particles (see spinless_fermions.VectorSpace, HQS Quantum Solver Reference). Fixing the number of particles, however, reduces the amount of memory and computation time Quantum Solver needs.

The second step, is to define the mathematical terms that define the operator. In Quantum Solver, operators are built by combining multiple mathematical terms. For the current example, we use the function called hopping, which is part of the spinless_fermions module, and which implements the term $$\sum _{j=0}^{M-1}\sum _{k=0}^{M-1}(h_0{)}_{jk}{c}_j^{†}{c}_k,$$ where $h_0\in {\mathbb{C}}^{M\times M}$ is the, so called, hopping matrix.

If we define the matrix $h_0$ to be $$h_0=\left[\begin{array}{ccccc}0& -t& & & \\ -t& 0& -t& & \\ & -t& \ddots & \ddots & \\ & & \ddots & \ddots & -t\\ & & & -t& 0\end{array}\right],$$ then the term in the previous formula becomes the tight-binding Hamiltonian defined at the beginning of this chapter. This matrix can be easily constructed using the toeplitz function from SciPy.

t = 1
h0 = toeplitz(-t * unit_vector(v.sites, 1))

With this matrix at hand, we can construct the Hamiltonian, by passing the term and the vector space defining the operator to the constructor of the Operator class.

H = Operator(hopping(h0), domain=v)

The Quantum Solver operators primarily provide matrix-vector and matrix-matrix routines. More precisely, most operators, like the one we have constructed here, implement the methods defined in the OperatorProtocol class. As a first example for using the operator, we want to look at its spectrum. For that, we can use the eigsh function from SciPy. Quantum Solver operators implement the SciPy LinearOperator interface, and are therefore directly compatible with many SciPy functions.

k = min(12, H.shape[0] - 1)  # Number of eigenvalues to compute
evals, evecs = eigsh(H, k=k, which="SA")

print(evals)

plt.figure("eigenvalues", clear=True)
plt.title(f"Eigenvalues (M = {M}, N = {N})")
plt.plot(evals, ".")

Now, we have enough code to run the script for the first time. Doing so should produce the plot shown below, which shows the ten smallest eigenvalues of the Hamiltonian, and hence, the lowest energy levels of the system.

Before continuing, modify the values of $M$ and $N$, run the script, and look at the result a few times.

Below, we have listed the 5 smallest eigenvalues for M = 16 and different values of N.

We can make our first observation: Adding the first two values of the $N=1$ column gives you the first value in the $N=2$ column, and adding the first three values of the $N=1$ column will give you the first value in the $N=3$ column. The reason for this is that in a two- and three-particle system of non-interacting paricles, the lowest energy to obtain is the one of the particles occupying the lowest two and three energy levels, respectively. Since we are dealing with a fermionic system, the particles are not allowed to occupy the same state, which explains this behavior.

We shall later come back to this observation.

Note how few lines of code were necessary, to produce physically meaningful results.

Site Occupations

To take this example a little bit further, let us look at the computation of site occupation expectation values. More precisely, we want to compute $⟨n_j⟩$ for $j=0,\dots ,M-1$, where $n_j=c_j^{†}{c}_j$. Note that $$⟨n_j⟩=⟨\psi |n_j|\psi ⟩=⟨\psi |c_j^{†}{c}_j|\psi ⟩=(c_j|\psi ⟩{)}^{†}{c}_j|\psi ⟩.$$ Hence, we just need to compute $c_j|\psi ⟩$ and then compute the square norm of this vector.

First, we have to construct the annihilation operator $c_j$. The annihilation operator reduces the number of particles in the system by one. Since, Quantum Solver assumes by default that the domain and the codomain are identical, we need to specify the codomain explicitly. To construct the codomain vector space, we can use the copy method on the existing vector space to obtain a modified copy.

v_minus = v.copy(particle_number_change=-1)

The term that describes the annihilation operator is given by the annihilation function, which takes the keyword argument site to specify the site at which to apply the annihilation operator. Passing this term, the domain, and the codomain to the constructor of the Operator class constructs the desired operator.

Once we have constructed the operator, we can use the dot method, to apply it to a state vector. Below, we define a function site_occupation which computes the expectation value for a given state vector and site. We then plot the expectation value for every site and for the state vector corresponding to the energy level E.

def site_occupation(psi, j):
    annihilation_operator = Operator(annihilation(site=j), domain=v, codomain=v_minus)
    c_psi = annihilation_operator.dot(psi)
    return c_psi.conj() @ c_psi


occupations = np.array([site_occupation(evecs[:, E], j) for j in range(M)])
plt.figure("site_occupation", clear=True)
plt.title(f"Site Occupation (M = {M}, N = {N}, E = {E})")
plt.plot(occupations, "x-")
plt.axis((0, M - 1, 0, 1))
plt.xlabel("j")
plt.ylabel("Expectation Value")

Energy occupations

While looking a site occupations is certainly interesting, for most applications it is more useful to look at energy occupations instead.

Consider the matrix $Q$ given by $$Q_{jk}=\sqrt{\frac{2}{M+1}}\sin \left(\frac{(j+1)(k+1)\pi }{M+1}\right)\text{for}j,k=0,\dots ,M-1.$$ The columns of this matrix are the eigenvectors of the matrix $h_0$. It is a well-known fact that if we define $$d_k=\sum _{j=0}^{M-1}{Q}_{jk}^{\ast }{c}_j$$ we can write the Hamiltonian $H$ as $$H=\sum _k{\epsilon }_k{d}_k^{†}{d}_k$$ where ${\epsilon }_k$ are the energy levels of the system. Consequently, we can measure the expectation value of the occupation of the $k$-th energy level, by computing $$⟨d_k^{†}{d}_k⟩=⟨\psi |d_k^{†}{d}_k|\psi ⟩.$$

We can carry out this computation in a similar way to the computation of the site occupation. All we have to do is to set the annihilation_operator variable to a linear combination of annihilation operators, which is given by the definition of $d_k$. Luckily, this is such a common operation that the Quantum Solver includes a convenience method for exactly this case. Instead of specifying the site argument of the annihilation function we can use the coef argument to specify a linear combination of annihilation operators.

J, K = np.meshgrid(np.arange(M), np.arange(M), indexing='ij')
Q = sqrt(2) / sqrt(M+1) * np.sin((J+1) * (K+1) * pi / (M+1))

def energy_occupation(psi, k):
    annihilation_operator = Operator(
        annihilation(coef=Q[:,k].conj()), domain=v, codomain=v_minus
    )
    c_psi = annihilation_operator.dot(psi)
    return (c_psi.conj() @ c_psi).real


ks = np.arange(M)
expectation = np.array([energy_occupation(evecs[:, E], k) for k in ks])
plt.figure("energy", clear=True)
plt.title(f"Energy Occupation (M = {M}, N = {N}, E = {E})")
plt.plot(ks, expectation, "x-")
plt.xlabel("p")
plt.ylabel("Expectation Value")
plt.show()

Below, you can find the plots that this script produces for different parameter configurations.

In these figures, you can see how the fermions arrange in the different one-particle energy states to achieve the energy levels we have observed in the table in the Setting up the Hamiltonian section.

Note that if you do not know the eigenvectors of the matrix $h_0$, you can easily obtain them by using the eigh function of NumPy.

Summary

To conclude this introduction, let us sum up the main takeaways of this chapter.

• Quantum Solver works by providing classes for working with quantum mechanical operators.
• Operators are constructed by specifying a term and a vector space.
• For operators that change the particle number, the codomain must be given as well.
• Operators provide operations on vectors and matrices (as described in the OperatorProtocol class), like the dot method, which applies the operator to a state vector.

Complete Code

# Title    : HQS Quantum Solver Getting Started Guide
# Filename : getting_started.py

# Python, NumPy, SciPy, and Matplotlib imports
from math import pi, sqrt
import numpy as np
from scipy.linalg import toeplitz
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt

# HQS Quantum Solver imports
from hqs_quantum_solver.spinless_fermions import VectorSpace, Operator, hopping, annihilation
from hqs_quantum_solver.util import unit_vector

# ===== Parameters =====

M = 16  # site count
N = 1   # particle count
E = 0   # energy level

# ===== Creation of the Hamiltonian =====

v = VectorSpace(sites=M, particle_numbers=N)

t = 1
h0 = toeplitz(-t * unit_vector(v.sites, 1))

H = Operator(hopping(h0), domain=v)

# ===== Computing the Eigenvalues and Eigenvectors =====

k = min(12, H.shape[0] - 1)  # Number of eigenvalues to compute
evals, evecs = eigsh(H, k=k, which="SA")

print(evals)

plt.figure("eigenvalues", clear=True)
plt.title(f"Eigenvalues (M = {M}, N = {N})")
plt.plot(evals, ".")

# ===== Plotting the Site Occupations =====

v_minus = v.copy(particle_number_change=-1)

def site_occupation(psi, j):
    annihilation_operator = Operator(annihilation(site=j), domain=v, codomain=v_minus)
    c_psi = annihilation_operator.dot(psi)
    return c_psi.conj() @ c_psi


occupations = np.array([site_occupation(evecs[:, E], j) for j in range(M)])
plt.figure("site_occupation", clear=True)
plt.title(f"Site Occupation (M = {M}, N = {N}, E = {E})")
plt.plot(occupations, "x-")
plt.axis((0, M - 1, 0, 1))
plt.xlabel("j")
plt.ylabel("Expectation Value")

# ===== Plotting the Energy Occupations =====

J, K = np.meshgrid(np.arange(M), np.arange(M), indexing='ij')
Q = sqrt(2) / sqrt(M+1) * np.sin((J+1) * (K+1) * pi / (M+1))

def energy_occupation(psi, k):
    annihilation_operator = Operator(
        annihilation(coef=Q[:,k].conj()), domain=v, codomain=v_minus
    )
    c_psi = annihilation_operator.dot(psi)
    return (c_psi.conj() @ c_psi).real


ks = np.arange(M)
expectation = np.array([energy_occupation(evecs[:, E], k) for k in ks])
plt.figure("energy", clear=True)
plt.title(f"Energy Occupation (M = {M}, N = {N}, E = {E})")
plt.plot(ks, expectation, "x-")
plt.xlabel("p")
plt.ylabel("Expectation Value")
plt.show()

Building Operators

Quantum Solver allows constructing quantum mechanical operators in a flexible and convenient way. As we have seen in the "Getting Started" section, quantum mechanical operators are constructed by providing a mathematical term and a vector space. The constructor of the Operator class, however, does not only accept single terms, but arbitrary linear combinations of terms, which makes this interface flexible.

To demonstrate this feature, let us consider the following example. Assume we want to construct the Hamiltonian $H$, $$H=-t\sum _{j=0}^{M-2}(c_j^{†}{c}_{j+1}+c_{j+1}^{†}{c}_j)+U\sum _{j=0}^{M-2}(n_j-\frac{1}{2})(n_{j+1}-\frac{1}{2}),$$ where $n_j=c_j^{†}{c}_j$. We start by creating a VectorSpace object, like in the "Getting Started" section, i.e.:

v = VectorSpace(sites=M, particle_numbers=N)

We can now consider the individual parts of the formula for the Hamiltonian. The term ${\sum }_{j=0}^{M-2}(c_j^{†}{c}_{j+1}+c_{j+1}^{†}{c}_j)$ can be obtained by passing the matrix $$\left[\begin{array}{ccccc}0& 1& & & \\ 1& 0& 1& & \\ & 1& \ddots & \ddots & \\ & & \ddots & \ddots & 1\\ & & & 1& 0\end{array}\right]$$ to the hopping function, and this matrix can be constructed by:

neighbors = toeplitz(unit_vector(v.sites, 1))

The term ${\sum }_{j=0}^{M-2}(n_j-\frac{1}{2})(n_{j+1}-\frac{1}{2})$ can be obtained by passing the matrix $$\left[\begin{array}{ccccc}0& 1& & & \\ & 0& 1& & \\ & & \ddots & \ddots & \\ & & & 0& 1\\ & & & & 0\end{array}\right]$$ to the interaction function, and this matrix can be constructed by:

right_neighbors = toeplitz(np.zeros(v.sites), unit_vector(v.sites, 1))

We can now combine the individual terms, to construct the desired Hamiltonian.

H = Operator(-t * hopping(neighbors) + U * interaction(right_neighbors), domain=v)

Using the same code as discussed in the "Getting Started" section, we can compute the expectation value of the occupation for different values of $U$.

Complete Code

# Title    : HQS Quantum Solver Hubbard Interaction
# Filename : hubbard_interaction.py

# NumPy, SciPy, and Matplotlib imports
import numpy as np
from scipy.linalg import toeplitz
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt

# HQS Quantum Solver imports
from hqs_quantum_solver.spinless_fermions import (
    VectorSpace,
    Operator,
    hopping,
    interaction,
    annihilation,
)
from hqs_quantum_solver.util import unit_vector

# ===== Parameters =====

M = 16  # site count
N = 2  # particle count
t = 1  # hopping strength
U = 10  # interaction strength

# ===== Creation of the Hamiltonian =====

v = VectorSpace(sites=M, particle_numbers=N)

neighbors = toeplitz(unit_vector(v.sites, 1))
right_neighbors = toeplitz(np.zeros(v.sites), unit_vector(v.sites, 1))

H = Operator(-t * hopping(neighbors) + U * interaction(right_neighbors), domain=v)

# ===== Computing the Eigenvalues and Eigenvectors =====

k = min(12, H.shape[0] - 1)  # Number of eigenvalues to compute
evals, evecs = eigsh(H, k=k, which="SA")
groundstate = evecs[:, 0]

# ===== Plotting the Site Occupations =====

v_minus = v.copy(particle_number_change=-1)


def site_occupation(psi, j):
    annihilation_operator = Operator(annihilation(site=j), domain=v, codomain=v_minus)
    c_psi = annihilation_operator.dot(psi)
    return c_psi.conj() @ c_psi


occupations = np.array([site_occupation(groundstate, j) for j in range(M)])
plt.figure("site_occupation", clear=True)
plt.title(f"Site Occupation (M = {M}, N = {N}, t={t}, U={U})")
plt.plot(occupations, "x-")
plt.axis((0, M - 1, 0, 1))
plt.xlabel("j")
plt.ylabel("Expectation Value")
plt.show()

Lattice Builder

When building Hamiltonians with Quantum Solver, it is often required to create adjacency matrices that represent the lattice structure of the system under consideration, e.g., for the hopping term. Building those terms in 1d is relatively straight forward. For higher dimensions and more complicated lattices, however, this becomes more difficult. The HQS Lattice Builder and Lattice Validator are libraries that simplify this process and can directly be used with Quantum Solver.

As an example, we consider the computation of the groundstate of the tight-binding Hamiltonian, but this time on a rectangular lattice. The Hamiltonian is given by $$H=-t\sum _{(j,k)\in \mathcal{N}}{c}_j^{†}{c}_k,$$ where $\mathcal{N}$ is the set of all tuples $(j,k)$ such that $j$ is a neighbor of $k$.

We start the example by specifying the required imports.

# NumPy, SciPy, and Matplotlib imports
import numpy as np
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt

# HQS Quantum Solver imports
from hqs_quantum_solver.spinless_fermions import (
    VectorSpace, Operator, lattice_terms, annihilation
)

# HQS Lattice Builder/Validator imports
import lattice_builder
import lattice_validator

Then, we define the parameters of the simulation.

Mx = 5  # sites in x-direction
My = 4  # sitex in y-direction
N = 1  # particle count
t = 1  # hopping strengh

The next step is to create a Lattice Builder/Validator configuration that specifies the Hamiltonian on an Mx $\times$ My rectangular lattice. In this example we have just one type of atom (lattice site) that is repeated in the direction of the vectors in "lattice_vectors". The hopping term is defined in the bonds variable, which adds an entry of value -t into the hopping matrix for neighboring lattice sites.

For more details on the Lattice Builder configuration, take a look at the Hamiltonians in "Physics", HQS Lattice Documentation documentation, and at "Schema definitions", Lattice Validator Documentation.

system = {
    "site_type": "spinless_fermions",
    "cluster_size": [Mx, My, 1],
    "system_boundary_conditions": ["hard-wall", "hard-wall", "hard-wall"]
}
atoms = [{"id": 0}]
bonds = [
    {"id_from": 0, "id_to": 0, "translation": [1, 0, 0], "t": -t},
    {"id_from": 0, "id_to": 0, "translation": [0, 1, 0], "t": -t},
]
unitcell = {
    "atoms": atoms,
    "bonds": bonds,
    "lattice_vectors": [[1, 0, 0], [0, 1, 0]]
}
conf = {
    "system": system,
    "unitcell": unitcell
}

The configuration needs to be validated and normalized, which is done by using the Lattice Validator. Then, a Lattice Builder instance can be created. This instance is then passed to the lattice_terms function, to convert it into an operator definition.

Note, if you do not need access to the Lattice Builder instance, you can also directly pass the conf dict to lattice_terms, which will take care of the validation and builder construction.

is_valid, conf = lattice_validator.Validator().validateConfiguration(conf)
if not is_valid:
    raise ValueError(f"Configuration not valid: {conf}")
builder = lattice_builder.Builder(conf)

v = VectorSpace(sites=Mx * My, particle_numbers=N)
H = Operator(lattice_terms(builder), domain=v)

After having obtained the Hamiltonian operator, the computation of the groundstate and the definition of the site_occupation function are the same as in the "Getting Started" section.

k = min(12, H.shape[0] - 1)  # Number of eigenvalues to compute
evals, evecs = eigsh(H, k=k, which="SA")
groundstate = evecs[:, 0]

v_minus = v.copy(particle_number_change=-1)

def site_occupation(psi, j):
    annihilation_operator = Operator(annihilation(site=j), domain=v, codomain=v_minus)
    c_psi = annihilation_operator.dot(psi)
    return c_psi.conj() @ c_psi

To plot the result we use the get_atom_positions method from the Lattice Builder to show the groundstate occupation on a 2d grid.

occupations = np.zeros((My, Mx))
for j, position in enumerate(builder.get_atom_positions().tolist()):
    x, y = round(position[0]), round(position[1])
    # Set value for row=y and column=x.
    occupations[y, x] = site_occupation(groundstate, j)

plt.figure("site_occupation", clear=True)
plt.title(f"Site Occupation (lattice: {Mx}x{My}, particles: {N})")
plt.imshow(occupations, origin="lower", aspect="equal")
plt.xlabel("x")
plt.ylabel("y")
plt.xticks(np.arange(Mx))
plt.yticks(np.arange(My))
plt.colorbar()
plt.show()

Finally, we can run the code, to obtain a plot of the occupation, as shown below.

Complete Code

# Title    : HQS Quantum Solver using the Lattice Builder
# Filename : getting_started.py

# NumPy, SciPy, and Matplotlib imports
import numpy as np
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt

# HQS Quantum Solver imports
from hqs_quantum_solver.spinless_fermions import (
    VectorSpace, Operator, lattice_terms, annihilation
)

# HQS Lattice Builder/Validator imports
import lattice_builder
import lattice_validator

# ===== Parameters =====

Mx = 5  # sites in x-direction
My = 4  # sitex in y-direction
N = 1  # particle count
t = 1  # hopping strengh

# ===== Definition of the Lattice =====

system = {
    "site_type": "spinless_fermions",
    "cluster_size": [Mx, My, 1],
    "system_boundary_conditions": ["hard-wall", "hard-wall", "hard-wall"]
}
atoms = [{"id": 0}]
bonds = [
    {"id_from": 0, "id_to": 0, "translation": [1, 0, 0], "t": -t},
    {"id_from": 0, "id_to": 0, "translation": [0, 1, 0], "t": -t},
]
unitcell = {
    "atoms": atoms,
    "bonds": bonds,
    "lattice_vectors": [[1, 0, 0], [0, 1, 0]]
}
conf = {
    "system": system,
    "unitcell": unitcell
}

# ===== Creation of the Hamiltonian =====

is_valid, conf = lattice_validator.Validator().validateConfiguration(conf)
if not is_valid:
    raise ValueError(f"Configuration not valid: {conf}")
builder = lattice_builder.Builder(conf)

v = VectorSpace(sites=Mx * My, particle_numbers=N)
H = Operator(lattice_terms(builder), domain=v)

# ===== Computing the Eigenvalues and Eigenvectors =====

k = min(12, H.shape[0] - 1)  # Number of eigenvalues to compute
evals, evecs = eigsh(H, k=k, which="SA")
groundstate = evecs[:, 0]

# ===== Plotting the Site Occupations =====

v_minus = v.copy(particle_number_change=-1)

def site_occupation(psi, j):
    annihilation_operator = Operator(annihilation(site=j), domain=v, codomain=v_minus)
    c_psi = annihilation_operator.dot(psi)
    return c_psi.conj() @ c_psi

occupations = np.zeros((My, Mx))
for j, position in enumerate(builder.get_atom_positions().tolist()):
    x, y = round(position[0]), round(position[1])
    # Set value for row=y and column=x.
    occupations[y, x] = site_occupation(groundstate, j)

plt.figure("site_occupation", clear=True)
plt.title(f"Site Occupation (lattice: {Mx}x{My}, particles: {N})")
plt.imshow(occupations, origin="lower", aspect="equal")
plt.xlabel("x")
plt.ylabel("y")
plt.xticks(np.arange(Mx))
plt.yticks(np.arange(My))
plt.colorbar()
plt.show()

System Types

Up until this point, we have only considered quantum mechanical systems composed of "spinless fermions". HQS Quantum Solver, however, provides further types of quantum systems. Each of these systems is implemented in its own module, and all these modules are structured in a similar way, meaning that they all provide a VectorSpace, an Operator, and a set of functions providing operator terms.

Spins

The spins module implements systems consisting of $M$ spin-$\frac{1}{2}$ particles. The VectorSpace is spanned by vectors of the form $$\underset{j=0}{\overset{M-1}{⨂}}|{\sigma }_j⟩\text{where}{\sigma }_j\in \{\uparrow ,\downarrow \}.$$

Spinless Fermions

The spinless_fermions module implements systems of fermions having a site index but no (dedicated) spin index. The VectorSpace is spanned by vectors of the form $$\prod _{j=0}^{M-1}(c_j^{†}{)}^{n_j}|0⟩\text{where}{n}_j\in \{0,1\},$$ and $c_j^{†}$ is the creation operator for site $j$.

Spinful Fermions

The spinful_fermions module implements system of fermions having a site index and a spin index. The VectorSpace is spanned by vectors of the form $$\prod _{\sigma \in \{\uparrow ,\downarrow \}}\prod _{j=0}^{M-1}(c_{j\sigma }^{†}{)}^{n_{j\sigma }}|0⟩\text{where}{n}_{j\sigma }\in \{0,1\},$$ and $c_{j\sigma }^{†}$ is the creation operator for site $j$ and spin polarization $\sigma$.

Bosons

The bosons module implements systems of bosons. The VectorSpace is spanned by vectors of the form $$\prod _{j=0}^{M-1}\sqrt{n_j!}(b_j^{†}{)}^{n_j}|0⟩\text{where}{n}_j\in \{0,1,\dots ,N_j\},$$ and $b_j^{†}$ is the creation operator for site $j$. The variable $N_j$ is the maximal occupation of site $j$.

Examples

Before looking at the examples, it is advised that you read the Getting Started Guide.

The examples described below can be accessed via HQStage.

1. spinful_fermions_groundstate.ipynb

   In this example, we show how to compute the ground state of the Hamiltonian of a many-particle system.

2. spin_waves.ipynb

   This example introduces spin systems and time evolution. We consider the time evolution of a Heisenberg spin model in an excited state.

3. nonequilibrium_electron_transport.ipynb

   This example simulates the time evolution of a system of electrons starting in a nonequilibrium state.

4. spinless_fermions_greens_function.ipynb

   Example for the Calculation of the Spectral function using HQS Quantum Solver and Lattice Functions.

   This method combines operators, Chebyshev expansion, cluster perturbation theory and clever Green's function mathematics to create the reciprocal-space resolved spectral function (a. k. a. the interacting band structure) in a honeycomb lattice causing an opening of the Dirac point.

5. spinful_fermions_greens_function.ipynb

   Prototype example for the Calculation of the Spectral function of the Hubbard model using hqs_quantum_solver. In this approach we represent the wavefunction as a dyadic product of up and down states.

   Similar to the previous example, this example computes the interacting band structure of a square lattice under the influence of interaction. We retrieve the expected behaviour of opening a gap at half-filling.

6. bose_hubbard_groundstate.ipynb

   In this example we compute the groundstates of Bose-Hubbard systems, which is a simple bosonic system.

7. lidar.ipynb

   In this scientific example we study properties of light emitted by multiple atoms, a fundamental model in quantum optics. We use the HQS Quantum Solver to time evolve the system according to the Lindblad equation of motion and calculate various expectation values of system operators describing properties of the emitted light.

Reference

You can find the API reference here.

Theory

Time Evolution using Krylov Methods

In this section, we will discuss some of the mathematics behind time-evolution using Krylov methods.

First, a refresher on a Krylov subspace. Given an initial state $v$, and matrix $A$, the $n^{\text{th}}$ order Krylov subspace is given by: $$\begin{array}{c}{\mathcal{K}}_n\equiv \text{span}\left\{v,Av,A^{2}v,...,A^{n-1}v\right\}\end{array}$$

where span refers to the space spanned by the vectors. In practice, this list of vectors is converted into an orthonormal basis by (modified-)Gram-Schmidt methods, Givens rotations or Householder reflections.

For the evaluation of the time-propagator, we seek approximations of the form:

$$\begin{array}{c}{e}^{A}v\approx p_{n-1}(A)v\end{array}$$

that is, that we can expand the effect of our evolution operator as a polynomial of order $n-1$ (in $A$). To do so, first we define our orthonormal basis using the Arnoldi method:

compute $v_0:=v/||v|{|}_2$
for $k=1,2,\dots ,n-1$:
    compute $v_k:=A{v}_{k-1}$
    for $j$ in $0,\dots ,k-1$:
        $h_{j,k-1}:=v_j^{†}{v}_k$
        $v_k:=v_k-h_{j,k-1}$
    $h_{k,k-1}:=||v_k|{|}_2$
    $v_k:=v_k/h_{k,k-1}$

then our basis $V_n$ is comprised of the vectors $v_k$ and the subscript $n$ refers to the fact that this is the basis formed from the $n^{\text{th}}$ order Krylov subspace.

The Arnoldi method has the useful property, that upon applying the basis vectors to the original matrix, the resulting smaller matrix $H_n\in {\mathbb{C}}^{n\times n}$, $$\begin{array}{c}{V}_n^{†}A{V}_n=H_n,\end{array}$$ is tridiagonal when $A$ is Hermitian, and is upper-Hessenberg otherwise. Both are more easily exponentiated, but the tridiagonal case especially so.

Now that we have our basis, we seek to find an approximate solution $e^{A}v\approx V_n{y}_0\in {\mathcal{K}}_n$. More precisely, we are looking for a solution of the minimization problem $$\begin{array}{c}{y}_0=\underset{y}{\mathrm{argmin}}{\Vert e^{A}v-V_{n}y\Vert }_2.\end{array}$$

Since the columns of $V_n$ are orthonormal, we can split the identity into the projection onto the range of $V$ and its orthogonal complement via $I=V_n{V}_n^{†}+(I-V_n{V}_n^{†})$. Thus,

$$\begin{array}{c}{\Vert e^{A}v-V_{n}y\Vert }_2={\Vert (V_n{V}_n^{†}+(I-V_n{V}_n^{†}))(e^{A}v-V_{n}y)\Vert }_2\end{array}$$

Expanding, using the pythagorean theorem, and simplifying gives

$$\begin{array}{c}{\Vert e^{A}v-V_{n}y\Vert }_2={\Vert V_n{V}_n^{†}{e}^{A}v-V_{n}y{\Vert }_2+\Vert e^{A}v-V_n{V}_n^{†}{e}^{A}v\Vert }_2.\end{array}$$

Note, that the second norm on the right does not include the variable $y$, while choosing

$$\begin{array}{c}y=V_n^{†}{e}^{A}v,\end{array}$$

makes the first norm on the right equal to zero. Since norms are non-negative, this expression is, therefore, the desired minimizing argument $y_0$.

This equation, however, still involves the original matrix exponential. By noting in the orthogonalization process of the Krylov subspace, that $$\begin{array}{c}v\equiv V_n{e}_0\end{array}$$ ($e_0$ being the unit vector at index $0$), we rewrite our solution as $$y_0=V_n^{†}{e}^A{V}_n{e}_0.$$

Hence

$$\begin{array}{c}{y}_0=V_n^{†}{e}^A{V}_n{e}_0=V_n^{†}\sum _{k=0}^{\infty }\frac{A^k}{k!}{V}_n{e}_0=\sum _{k=0}^{\infty }\frac{V_n^{†}{A}^k{V}_n}{k!}{e}_0.\end{array}$$

We would like to move the matrice $V_n$ into the exponentiation, because then we would be able to compute $y_0$ by

$$\begin{array}{c}\sum _{k=0}^{\infty }\frac{(V_n^{†}A{V}_n{)}^k}{k!}{e}_0=\sum _{k=0}^{\infty }\frac{H_n^k}{k!}{e}_0=e^{H_n}{e}_0,\end{array}$$

which only involves the small matrix $H_n$. However,

$$\begin{array}{c}(V_n^{†}A{V}_n{)}^k=V_n^{†}A{V}_n{V}_n^{†}A{V}_n\cdots V_n^{†}A{V}_n\ne V_n^{†}(A{)}^k{V}_n,\end{array}$$ since $V_n{V}_n^{†}\ne I$. Let us assume for a moment that the range of $V_n$ is invariant under the operation of $A$, i.e., $Av\in \mathcal{R}(V_n)$ for $v\in \mathcal{R}(V_n)$. Then, $V_n{V}_n^{†}A{V}_n=A{V}_n$, since $V_n{V}_n^{†}$ is the orthogonal projector onto $\mathcal{R}(V_n)$. Thus, under this assumption the desired equation does hold, and if we assume that the range of $V_n$ is almost invariant under the action of $A$ then the equation is still a good approximation.

$$\begin{array}{c}{e}^{A}v\approx V_n{y}_0\approx V_n{e}^H{e}_0.\end{array}$$

As a final point, if we want to compute the time evolution of a quantum system defined by a Hamiltonian $A$, we want to compute $e^{-\mathrm{i}A\tau }v$ for some value of $\tau \in \mathbb{R}$. We note that as the decomposition is invariant with respect to (non-zero) scalar multiples, we may modify $A\to -\mathrm{i}A\tau$ to give $$\begin{array}{c}{e}^{-\mathrm{i}A\tau }v\approx V_n{e}^{-\mathrm{i}{H}_n\tau }{e}_0.\end{array}$$

For more details on related methods see the following literature:

E. Gallopoulos, Y. Saad "On the parallel solution of parabolic equations." Proceedings of the 3rd international conference on Supercomputing, 1989 17-28. doi: 10.1145/318789.318793

E. Gallopoulos, Y. Saad "Efficient solution of parabolic equations by Krylov approximation methods." SIAM journal on scientific and statistical computing 13(5) 1992: 1236-1264. doi: 10.1137/0913071

Y. Saad, "Analysis of some Krylov subspace approximations to the matrix exponential operator." SIAM Journal on Numerical Analysis 29(1) 1992: 209-228. doi: 10.1137/0729014

Implementation

Protocols

HQS Quantum Solver makes extensive use of Protocols, which means that it is easy to write code interfacing with Quantum Solver. Whenever you encounter a protocol class in the type hints of a function or method, the protocol class defines which methods and properties the corresponding argument object is expected to have and the corresponding return type is going to provide, respectively.

For example, the evolve_using_krylov function expects a hamiltonian argument which adheres to the SimpleOperatorProtocol. To be allowed to pass an object as hamiltonian, the object just needs to implement all the methods and properties listed in SimpleOperatorProtocol, which is in contrast to abstract base classes, where the object would need to be explicitly derived from the abstract base class.

StateMaps

HQS Quantum Solver is able to exploit various conserved properties by implementing, so-called, StateMaps. A StateMap provides a list of the Fock states that are used to build an operator; it defines the basis set of our Hilbert space. By specifically generating Hamiltonians and generic operators, with custom StateMaps, you can efficiently exploit these properties. This is typically done automatically when specifying things like conserved particle number, by mod_N = 0, or conserved Sz quantum number, by mod_Sz = 0.

Dyadic representation

We also provide a special kind of operator that exploits the so-called dyadic representation. Specifically, this methodology can be used for systems where the state exists in a product of spaces, i.e., ${\psi }_{ab}={\psi }_a\otimes {\psi }_b$ and the operators acting on it are either diagonal in the extended space or may be described by Kronecker sums/products of operators in the individual spaces, implemented by DyadicOperatorSum and DyadicOperatorProduct, respectively.

We demonstrate the usefulness of this concept by evaluating the Hubbard Model. Starting with a spinful fermion model, we assume that the wavefunction can be expressed as a Kronecker product in the spin-up and spin-down degrees of freedom $\Psi ={\psi }_{\uparrow }\otimes {\psi }_{\downarrow }$.

To do this we use a well-known matrix vectorization trick in reverse, let $\text{vec}(X)$ denote the so-called row-major vectorization operator acting on the matrix $X$, then the Kronecker product has the property that $$(A\otimes B)\text{vec}(X)=\text{vec}(AX{B}^T).$$

Similarly, if the operator is a Kronecker sum, $$A\oplus B\equiv (A\otimes I_B)+(I_A\otimes B).$$ Using the first property, we get that $$(A\otimes I_B)\text{vec}(X)=\text{vec}(AX{I}_B)=\text{vec}(AX)$$ and $$(I_A\otimes B)\text{vec}(X)=\text{vec}(I_{A}X{B}^T)=\text{vec}(X{B}^T)$$ and therefore, combining with the second property, $$(A\otimes I_B+I_A\otimes B)\text{vec}(X)=\text{vec}(AX+X{B}^T).$$ This expression is the representation that is preferable to directly representing the operator on the product space, because it requires less storage and fewer operations to execute when applying the operator.