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 j = 0 \sum M - 2 (c_{j}^{†} c_{j + 1} + c_{j + 1}^{†} c_{j}),$ 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 $j = 0 \sum M - 1 k = 0 \sum M - 1 (h_{0})_{jk} c_{j}^{†} c_{k},$ where $h_{0} \in C^{M \times M}$ is the, so called, hopping matrix.

If we define the matrix $h_{0}$ to be $h_{0} = 0 - t - t 0 - t - t ⋱ ⋱ ⋱ ⋱ - t - t 0,$ 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.

Note, for a fixed number of particles N the number of states Quantum Solver has to store is $(N M)$ . Hence, the number of states can grow rather quickly when increasing N while M is large. This problem gets worse when allowing the particle number to vary. Keep in mind, while Quantum Solver provides you with highly accurate computations, you will have to provide enough computational resources for it to shine.

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

#	$N = 1$	$N = 2$	$N = 3$
1	-1.97	-3.83	-5.53
2	-1.86	-3.67	-5.31
3	-1.70	-3.57	-5.14
4	-1.48	-3.44	-5.04
5	-1.21	-3.34	-5.04

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} ⟩ = ⟨ ψ ∣ n_{j} ∣ ψ ⟩ = ⟨ ψ ∣ c_{j}^{†} c_{j} ∣ ψ ⟩ = (c_{j} ∣ ψ ⟩)^{†} c_{j} ∣ ψ ⟩ .$ Hence, we just need to compute $c_{j} ∣ ψ ⟩$ 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")

The site occupation. — The site occupations for `M=16`, `N=1`, `E=0`.

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} = \frac{2}{M + 1} sin (\frac{( j + 1 ) ( k + 1 ) π}{M + 1}) 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} = j = 0 \sum M - 1 Q_{jk}^{*} c_{j}$ we can write the Hamiltonian $H$ as $H = k \sum ε_{k} d_{k}^{†} d_{k}$ where $ε_{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} ⟩ = ⟨ ψ ∣ d_{k}^{†} d_{k} ∣ ψ ⟩ .$

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.

The energy occupation. — The energy occupations for `M=16`, `N=2`, `E=0`.

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()