# Numerical Methods for Quantum Spin Dynamics

Danny Goodacre
Supervisor: Dr. Pranav Singh

(April 2022)

## Abstract

This report is concerned with the efficiency of numerical methods for
simulating quantum spin systems. We aim to implement an improved method for
simulation with a time-dependent Hamiltonian and behaviour which displays
chirped pulses at a high frequency.

Working in the density matrix formulation of quantum systems, we study
evolution under the Liouville-von Neumann equation in the Liouvillian picture,
presenting analysis of and benchmarking current numerical methods. The
accuracy of existing techniques is assessed in the presence of chirped pulses.

We discuss the Magnus expansion and detail how a truncation of it is used to
solve differential equations. The results of this work are implemented in the
Python package MagPy to provide a better error-to-cost ratio than current
approaches allow for time-dependent Hamiltonians. Within this method we also
consider how to approximate the computation of the matrix exponential.

### Motivation

The spin of a particle is an entirely quantum mechanical quantity that has no
classical analogue. Accurately predicting the dynamics of spins systems is of
great importance in applications such as quantum computing and nuclear
magnetic resonance imaging (NMR/MRI) [5, 16]. Our work is focused on systems
displaying highly oscillatory chirped pulses (HOCP), which are prevalent in
nuclear spin dynamics.

Simulating larger systems (2-5 spins) is common in spectroscopy and quickly
becomes prohibitively expensive using current techniques, motivating the
development of a method which is more computationally efficient.

### Outline

The report is organised as follows: the first section provides background in
the mathematics and theoretical concepts, presenting relevant equations and
proving theorems which are used in our application throughout.

Within the first section we introduce the formalism that we use to represent
quantum spin systems and discuss the equation underpinning their evolution
over time. We also detail the structure of the Hamiltonian which we
consider–for both single and multiple spins, interacting and non-
interacting—and introduce notation for concisely writing such Hamiltonians.

Next, we discuss the Liouvillian picture of quantum mechanics, giving us a
form to which we can readily apply various numerical methods. We prove some
useful results of the Liouvillian superoperator which will be later used in
our implementation.

Finally, we discuss how to measure the expectation values and spin components
of spin systems and visualise these in 3D space where possible.

In the second section, we discuss the Python package QuTiP and review current
numerical methods for solving initial value problems for ordinary differential
equations and both polynomial and iterative techniques for approximating the
matrix exponential.

For each method in this section we discuss their shortcomings regarding highly
oscillatory systems, such as what properties they do or do not conserve, and
how the error in the approximation is bounded depending on the qualities of
the system.

The third section begins with defining the Magnus expansion for time-dependent
first-order homogeneous linear differential equations. Truncating the
expansion, we reformulate the remaining terms to be applied more effectively
to the Hamiltonians with which we are concerned. An example of our
implementation in Python is provided with instructions for its use.

The fourth section is primarily concerned with comparing the effectiveness of
our Magnus-based approach against current techniques. We analyse the rates of
convergence for both our approach and others, and also consider the error-to-
cost ratio for different approximations for evaluating the integals involved
in the Magnus expansion.

###### Contents

  1. 1 Background
    1. 1.1 Density Matrices
    2. 1.2 The Hamiltonian
      1. 1.2.1 Single Spin Hamiltonian
      2. 1.2.2 Multi-Spin Hamiltonian
      3. 1.2.3 Highly Oscillatory Chirped Pulse Systems
      4. 1.2.4 Notation
    3. 1.3 The Liouville-von Neumann Equation
      1. 1.3.1 The Liouvillian Picture
      2. 1.3.2 Properties of the Liouvillian
    4. 1.4 Spin Components and Expectation Values
      1. 1.4.1 The Bloch Sphere
  2. 2 Current Methods
    1. 2.1 QuTiP
    2. 2.2 Standard ODE Methods
      1. 2.2.1 Euler’s Method
      2. 2.2.2 The Trapezoidal Rule
      3. 2.2.3 RK4
    3. 2.3 Approximating the Matrix Exponential
      1. 2.3.1 Taylor Series
      2. 2.3.2 Padé Approximants
      3. 2.3.3 Krylov Subspace Methods
  3. 3 The Magnus Expansion
    1. 3.1 Definition
    2. 3.2 Reformulation
      1. 3.2.1 Derivation of the Second Term
    3. 3.3 Implementation in Python
  4. 4 Numerical Comparison
    1. 4.1 One-term Magnus
    2. 4.2 Two-term Magnus
    3. 4.3 Larger Systems
  5. 5 Discussion
    1. 5.1 Summary
    2. 5.2 Further Work
  6. Appendix
    1. The Kronecker Product
    2. An Error Bound of Scaling and Squaring
    3. The Lanczos Algorithm

##  1 Background

In this section we discuss the mathematical representation and properties of
quantum spin systems; we define the Hamiltonian of a system, which—alongside a
given initial condition—determines the evolution of the system over time; and
we discuss the resulting differential equation describing the dynamics of a
quantum spin system and reformulate to the Liouvillian picture to apply
existing numerical methods. Furthermore, throughout this section we prove
properties of the Liouvilillan superoperator and introduce some new notation.

###  1.1 Density Matrices

We represent quantum states by density matrices, denoted $\rho(t)$. For a $n$
spin system, it is a complex-valued $2^{n}\times 2^{n}$ matrix. Density
matrices are Hermitian, positive semi-definite, and have unit trace

$\rho^{\dagger}=\rho,\;\;\;\rho\geq 0,\;\;\;\operatorname{Tr}\\{\rho\\}=1.$

The initial condition of the density matrix is denoted as $\rho_{0}$. Later in
this section we discuss how to use the density matrix to determine properties
of the system.

For a system of $n$ spins that are not interacting, the density matrix can be
written as

$\rho(t)=\rho_{1}(t)\otimes\cdots\otimes\rho_{n}(t),$ (1)

where $\rho_{k}(t)$ is the density matrix for each respective spin and
$\otimes$ is the Kronecker product. When there is interaction between the
spins, the density matrix cannot be written as such.

###  1.2 The Hamiltonian

The Hamiltonian for an $n$ spin system is an $2^{n}\times 2^{n}$ complex-
valued Hermitian (self-adjoint) matrix. For such a Hamiltonian,

$iH\in\mathfrak{su}(2^{n}),$

where $\mathfrak{su}(2^{n})$ is the Lie algebra of $\mathrm{SU}(2^{n})$ and
consists of $n\times n$ skew-Hermitian matrices with trace zero, with the
operator

$[A,B]=AB-BA,\;\;\;\;A,\,B\in\mathfrak{su}(2^{n}),$

called the commutator [15]. The space $\mathfrak{su}(2^{n})$ is closed under
commutation of any two of its elements, meaning the commutator as a matrix is
also skew-Hermitian

$[A,B]^{\dagger}=-[A,B].$

####  1.2.1 Single Spin Hamiltonian

The space $\mathfrak{su}(2)$, has dimension 3 and a basis formed of the skew-
Hermitian matrices

$i\sigma_{x},\;\;i\sigma_{y},\;\;i\sigma_{z},$

where

$\sigma_{x}=\begin{pmatrix}0&1\\\
1&0\end{pmatrix},\;\;\sigma_{y}=\begin{pmatrix}0&-i\\\
i&0\end{pmatrix},\;\;\sigma_{z}=\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix},$

called the Pauli spin matrices [15]. These matrices satisfy the following
commutation relations,

$[\sigma_{x},\sigma_{y}]=2i\sigma_{z},\;\;[\sigma_{y},\sigma_{z}]=2i\sigma_{x},\;\;[\sigma_{z},\sigma_{x}]=2i\sigma_{y}.$

Given two real functions and a constant,

$f\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}\to\mathbb{R},\;\;\;g\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}\to\mathbb{R},\;\;\;\Omega\in\mathbb{R},$

we define

$H\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}\to\mathbb{C}^{2\times 2}$
$H(t)=f(t)\sigma_{x}+g(t)\sigma_{y}+\Omega\sigma_{z},$

which is the form of a single spin Hamiltonian used throughout this report. As
the Pauli matrices are Hermitian and $\mathfrak{su}(2)$ is a linear subspace,
this form of Hamiltonian is still Hermitian. Note that for a time-independent
Hamiltonian we simply take $f$ and $g$ to be constant functions.

####  1.2.2 Multi-Spin Hamiltonian

When the system contains more than one spin, it can be either interacting or
non-interacting. Suppose the system contains $n$ spins. For each spin, we
define a Hamiltonian

$H_{j}(t)=f_{j}(t)\sigma_{x}+g_{j}(t)\sigma_{y}+\Omega_{j}\sigma_{z},$

for $j\in\\{1,\ldots,n\\}$. Then, for the non-interacting case, the
Hamiltonian is defined as

$\displaystyle H_{0}(t)=\,\,$ $\displaystyle H_{1}(t)\otimes
I\otimes\cdots\otimes I$ $\displaystyle+I\otimes H_{2}(t)\otimes\cdots\otimes
I$ $\displaystyle+\cdots$ $\displaystyle+I\otimes\cdots\otimes I\otimes
H_{n}(t),$

where $I$ is the $2\times 2$ identity matrix. If the spins interact, there is
an additional term in the Hamiltonian,

$H_{J}\in\mathbb{C}^{2^{n}\times 2^{n}},$

which is time-independent and Hermitian. Thus the Hamiltonian for an
interacting multi-spin system is

$H(t)=H_{0}(t)+H_{J}.$

Therefore, in order to specify an $n$ spin Hamiltonian, we require

$\left(f_{j}(t),\,g_{j}(t),\,\Omega_{j}\right),\;\;\;j\in\\{1,\ldots,n\\},$

and the constant matrix $H_{J}$ if the system is interacting.

####  1.2.3 Highly Oscillatory Chirped Pulse Systems

A system with which we are primarily concerned is one where each spin displays
highly oscillatory chirped pulses (HOCP) in its evolution. We define such a
system for a single spin as follows: let

$\displaystyle e(t)$
$\displaystyle=\beta\exp\bigg{(}\frac{-(t-10)^{8}}{10^{7}}\bigg{)},$
$\displaystyle w(t)$
$\displaystyle=\exp\bigg{(}i\,\gamma\,(t-10)^{2}\bigg{)}.$

Then,

$\displaystyle f(t)$
$\displaystyle=e(t)\,\operatorname{Re}\big{\\{}w(t)\big{\\}},$ $\displaystyle
g(t)$ $\displaystyle=e(t)\,\operatorname{Im}\big{\\{}w(t)\big{\\}},$

where $\operatorname{Re}$ and $\operatorname{Im}$ denote the real and
imaginary parts, respectively. Throughout this report we define such a system
by specifying the values of $\beta$, $\gamma$, and $\Omega$ for each spin’s
Hamiltonian, an initial condition $\rho_{0}$, and a constant matrix $H_{J}$ if
the spins are interacting.

Figure 1: The function $f$ with $\beta=10$ and $\gamma=2$.

####  1.2.4 Notation

Since the notation for a Hamiltonian is laborious to write, we introduce the
following operators $\mathbf{I}$ and $\mathbf{S}$.

######  Definition 1.1.

$\mathbf{I}_{j}\\{A\\}\mathrel{\mathop{\mathchar
58\relax}}=I\otimes\cdots\otimes I\otimes A\otimes I\otimes\cdots\otimes I,$

where $A$ is in the $j$th position. The number of terms in the product is
inferred from the number of spins in the system.

This operator is linear due to the linearity of the Kronecker product. Using
this, we write the Hamiltonian as

$H(t)=\sum_{j=0}^{n}\mathbf{I}_{j}\left\\{H_{j}(t)\right\\}+H_{J}$

When $A$ is a Pauli matrix, we simplify further,

$\mathbf{I}_{j}^{a}\mathrel{\mathop{\mathchar
58\relax}}=\mathbf{I}_{j}\\{\sigma_{a}\\},\;\;\;a\in\\{x,y,z\\}.$

Noting that
$H_{j}(t)=f_{j}(t)\,\sigma_{x}+g_{j}(t)\,\sigma_{y}+\Omega_{j}\,\sigma_{z},$
we apply the linearity of the Kronecker product to get

$\mathbf{I}_{j}\left\\{H_{j}(t)\right\\}=f_{j}(t)\,\mathbf{I}_{j}^{x}+g_{j}(t)\,\mathbf{I}_{j}^{y}+\Omega_{j}\,\mathbf{I}_{j}^{z}.$

Therefore the full interacting $n$ spin Hamiltonian can be written as

$H(t)=\sum_{j=1}^{n}\left(f_{j}(t)\,\mathbf{I}_{j}^{x}+g_{j}(t)\,\mathbf{I}_{j}^{y}+\Omega_{j}\,\mathbf{I}_{j}^{z}\right)+H_{J}.$
(2)

This leads us to the following,

######  Lemma 1.2.

$\int_{a}^{b}H(t)\,dt=\sum_{j=1}^{n}\left(\int_{a}^{b}f_{j}(t)\,dt\,\mathbf{I}_{j}^{x}+\int_{a}^{b}g_{j}(t)\,dt\,\mathbf{I}_{j}^{y}+\Omega_{j}(b-a)\,\mathbf{I}_{j}^{z}\right)+H_{J}(b-a)$

###### Proof.

Follows directly from definition of $H(t)$ and linearity of integration. ∎

######  Lemma 1.3 (Generalised mixed-product).

Given matrices $A_{j},B_{j}\in\mathbb{C}^{n\times n},\;\\\
j\in\\{1,\ldots,k\\},$

$(A_{1}\otimes A_{2}\otimes\cdots\otimes A_{k})(B_{1}\otimes
B_{2}\otimes\cdots\otimes B_{k})=A_{1}B_{1}\otimes
A_{2}B_{2}\otimes\cdots\otimes A_{k}B_{k}.$

###### Proof.

Induct on $k$—c.f. Broxson [3]. ∎

######  Definition 1.4.

Given matrices $A,B\in\mathbb{C}^{n\times n}$, let

$\mathbf{S}_{i,\,j}\left\\{A,B\right\\}\mathrel{\mathop{\mathchar
58\relax}}=I\otimes\cdots\otimes I\otimes A\otimes I\otimes\cdots\otimes
I\otimes B\otimes I\otimes\cdots\otimes I,$

where $A$ is in the $i$th position and $B$ is in the $j$th position. When
$i=j$ we write

$\displaystyle\mathbf{S}_{j,\,j}\left\\{A,B\right\\}$
$\displaystyle=I\otimes\cdots\otimes I\otimes AB\otimes I\otimes\cdots\otimes
I$ $\displaystyle=\mathbf{I}_{j}\left\\{AB\right\\}.$

######  Lemma 1.5.

Given matrices $A,B\in\mathbb{C}^{n\times n}$,

  1. 1.

$\mathbf{S}_{i,\,j}\left\\{A,B\right\\}=\mathbf{I}_{i}\left\\{A\right\\}\,\mathbf{I}_{j}\left\\{B\right\\},$

  2. 2.

When $i\neq j$,
$\mathbf{S}_{i,\,j}\left\\{A,B\right\\}=\mathbf{S}_{j,\,i}\left\\{B,A\right\\},$

  3. 3.

$\mathbf{S}_{j,\,j}\\{A,B\\}=\mathbf{S}_{j,\,j}\\{B,A\\}\Longleftrightarrow
AB=BA.$

###### Proof.

  1. 1.

Apply the generalised Kronecker product to
$\mathbf{I}_{i}\left\\{A\right\\}\,\mathbf{I}_{j}\left\\{B\right\\},$

  2. 2.

This follows directly from the definition.

  3. 3.

This also follows directly from the definition.

∎

To conclude, the following conventions will be maintained throughout:

  * •

$H_{j}(t),\;j\in\mathbb{N}$ refers to a function of the form

$f_{j}(t)\,\sigma_{x}+g_{j}(t)\,\sigma_{y}+\Omega_{j}\,\sigma_{z},$

  * •

$H_{0}(t)$ is

$\sum_{j=1}^{n}\mathbf{I}_{j}\left\\{H_{j}(t)\right\\},$

  * •

and $H(t)$ is

$H_{0}(t)+H_{J}.$

###  1.3 The Liouville-von Neumann Equation

Given a Hamiltonian, $H(t)$, the corresponding dynamics of the spin system’s
density matrix, $\rho(t)$, are determined by the following ordinary
differential equation,

$\frac{\partial\rho(t)}{\partial t}=-i\left[H(t),\,\rho(t)\right],$

called the Liouville-von Neumann equation [14]. Along with an initial
condition, $\rho_{0}$, we get a matrix-valued initial value problem. Using our
reduced notation, the initial density matrix for an $n$ particle non-
interacting system is written as

$\rho_{0}=\sum_{j=1}^{n}\mathbf{I}_{j}\left\\{\rho_{j}\right\\}.$

####  1.3.1 The Liouvillian Picture

In order to more readily apply existing numerical methods to the Liouville-von
Neumann equation, we change from a Hamiltonian framework to a Liouvillian
framework. This is achieved through use of the Liouvillian superoperator,

$L\mathrel{\mathop{\mathchar 58\relax}}\mathbb{C}^{n\times
n}\to\mathbb{C}^{n^{2}\times n^{2}}$ $L\\{X\\}=I_{n}\otimes X-X^{T}\otimes
I_{n},$

and column-major vectorisation

$\text{vec}\mathrel{\mathop{\mathchar 58\relax}}\mathbb{C}^{n\times
n}\to\mathbb{C}^{n\times 1},$

denoted $\text{vec}(X)$. From [12] we have

$\text{vec}(XY)=(I\otimes X)\text{vec}(Y)=(Y^{T}\otimes I)\text{vec}(X),$

and so, denoting $\text{vec}(\rho(t))$ by $\mathbf{r}(t)$, the Liouvillian-von
Neumann equation is transformed into

$\frac{\partial\mathbf{r}(t)}{\partial t}=-iL\\{H(t)\\}\,\mathbf{r}(t),$ (3)

which is in the form of a standard vector-valued initial value problem:

$\mathbf{r}^{\prime}(t)=A(t)\,\mathbf{r}(t),\;\;\mathbf{r}(t_{0})=\mathbf{r}_{0},$

where $\mathbf{r}_{0}=\text{vec}(\rho_{0}).$ It is in this form to which we
apply our numerical methods in later sections.

####  1.3.2 Properties of the Liouvillian

The main result here is Theorem 1.9 and is vital in the later implementation
of the Magnus expansion. Theorem 1.6 is also useful as it shows that the
Liouvillian preserves a useful property of the Hamiltonian.

######  Theorem 1.6 (Convservation of Hermitivity).

For a Hermitian matrix $H$, $L\left\\{H\right\\}$ is Hermitian.

###### Proof.

For a Hermitian matrix $H$, consider

$L\left\\{H\right\\}=I\otimes H-H^{T}\otimes I.$

Taking the conjugate transpose of both sides gives

$L\left\\{H\right\\}^{\dagger}=\left(I\otimes H\right)^{\dagger}-(H^{T}\otimes
I)^{\dagger}.$

Since conjugate transposition is distributive over the Kronecker product and
$H$ is Hermitian, we then get

$L\left\\{H\right\\}^{\dagger}=I\otimes H-H^{T}\otimes I.$

Thus $L\left\\{H\right\\}$ is Hermitian. ∎

######  Lemma 1.7.

For a matrix $X\in\mathbb{C}^{n\times n},$

$\int L\left\\{X\right\\}=L\left\\{\int X\right\\}.$

###### Proof.

This follows from the linearity of integration and of the Kronecker product. ∎

######  Lemma 1.8.

For $A(t)=-iL\left\\{H(t)\right\\},$

$\big{[}A(s),\,A(r)\big{]}=L\left\\{\big{[}H(r),\,H(s)\big{]}\right\\}.$

###### Proof.

We first note that by the mixed-product property of the Kronecker product,

$\displaystyle L\left\\{H(s)\right\\}\,L\left\\{H(r)\right\\}$
$\displaystyle=\left(I\otimes H(s)-H(s)^{T}\otimes I\right)\left(I\otimes
H(r)-H(r)^{T}\otimes I\right)$ $\displaystyle=I\otimes
H(s)H(r)-H(s)^{T}\otimes H(r)$ $\displaystyle-H(r)^{T}\otimes
H(s)+H(s)^{T}H(r)^{T}\otimes I.$

This implies that

$\displaystyle\big{[}-iL\left\\{H(s)\right\\},\,-iL\left\\{H(r)\right\\}\big{]}$
$\displaystyle=L\left\\{H(r)\right\\}L\left\\{H(s)\right\\}-L\left\\{H(s)\right\\}L\left\\{H(r)\right\\}$
$\displaystyle=I\otimes H(r)H(s)-I\otimes H(s)H(r)$
$\displaystyle+H(r)^{T}H(s)^{T}\otimes I-H(s)^{T}H(r)^{T}\otimes I.$

By linearity of the Kronecker product and properties of transposition,

$\displaystyle\big{[}-iL\left\\{H(s)\right\\},\,-iL\left\\{H(r)\right\\}\big{]}$
$\displaystyle=I\otimes\big{(}H(r)H(s)-H(s)H(r)\big{)}$
$\displaystyle-\big{(}H(r)H(s)-H(s)H(r)\big{)}^{T}\otimes I$
$\displaystyle=L\left\\{\big{[}H(r),\,H(s)\big{]}\right\\}.$

∎

######  Theorem 1.9.

$\int_{0}^{t}\int_{0}^{s}\big{[}-iL\left\\{H(s)\right\\},\,-iL\left\\{H(r)\right\\}\big{]}\,dr\,ds=L\left\\{\int_{0}^{t}\int_{0}^{s}\big{[}H(r),\,H(s)\big{]}\,dr\,ds\right\\}.$

###### Proof.

Apply Lemma 1.8 and then apply Lemma 1.7 twice. ∎

###  1.4 Spin Components and Expectation Values

An important aspect of a spin system is being able to discern properties and
analyse how these evolve over time. From the density matrix we can determine
the expectation value for a given operator. Operators are represented by
$n\times n$ Hermitian matrices for system of $n$ spins. We denote the
expectation value for an operator $A$ by $\langle A\rangle.$

######  Definition 1.10 (Frobenius Inner Product).

For $A,B\in\mathbb{C}^{n\times n}$,

$\langle A,B\rangle_{F}\mathrel{\mathop{\mathchar
58\relax}}=\operatorname{Tr}\\{A^{\dagger}B\\}.$

######  Theorem 1.11.

For any two Hermitian matrices $A,B\in\mathbb{C}^{n\times n},$ their Frobenius
inner product is real.

###### Proof.

For a Hermitian matrix $A\in\mathbb{C}^{n\times n},$ we have $A=A^{\dagger}$
and can write $A^{T}=\overline{A}.$ Applying these and standard properties of
conjugation and transposition, we get

$\displaystyle\operatorname{Tr}\\{A^{\dagger}B\\}$
$\displaystyle=\operatorname{Tr}\\{AB\\}$
$\displaystyle=\overline{\operatorname{Tr}\\{AB\\}}$
$\displaystyle=\overline{\operatorname{Tr}\\{A^{\dagger}B\\}},$

meaning

$\operatorname{Tr}\\{A^{\dagger}B\\}\in\mathbb{R}.$

∎

All single spin density matrices can be expressed as a linear combination of
the three Pauli matrices. When the operator is a Pauli matrix, the resulting
expectation value is the component of spin in the respective direction. For
example, consider the following HOCP system for time $t\in[0,5]$,

$\begin{cases}\beta_{1}=10,\;\;\gamma_{1}=2,\;\;\Omega_{1}=1,\\\
\rho_{0}=\sigma_{x}.\end{cases}$ (4)

Figure 2: Normalised component of (4) in $\sigma_{x}$ direction over time.

We calculate the $x$-component of spin—equivalently the expectation of the
operator $\sigma_{x}$—using the Frobenius inner product,

$\langle\sigma_{x}\rangle(t)=\langle{\rho(t),\sigma_{x}\rangle}_{F}.$

The normalised $x$-component is then

$d_{x}(t)=\frac{1}{2}\langle{\rho(t),\sigma_{x}\rangle}_{F}.$

And $d_{y}$ and $d_{z}$ are defined similarly. These components are such that

$d_{x}(t)^{2}+d_{y}(t)^{2}+d_{z}(t)^{2}=1.$

In Figure 1 at time $t=0$, we see that the $x$-component is $1$, meaning that
the spin is wholly in the $x$-direction. Thus $d_{y}$ and $d_{z}$ must be
zero.

For a $n$-spin system that is not interacting, with a density matrix of the
form (1), we can measure the $\sigma_{x}$ component of the $k$th spin using
the operator $\mathbf{I}_{k}^{x}$. The normalisation constant is $2^{n}$ for
$n$ particles. For example, the normalised $\sigma_{x}$ component of the first
of two non-interacting particles is

$\frac{1}{4}\langle{\rho(t),\sigma_{x}\otimes I\rangle}_{F}.$

####  1.4.1 The Bloch Sphere

It is often difficult to visualise the spin of a quantum system, especially
when there is more than one spin in the system. However, for a single spin we
can use the fact that an analogy can be drawn between the normalised
components and 3D space.

In particular, we can treat the three Pauli matrices as the canonical basis of
$\mathbb{R}^{3}$ and the three components as a vector projected from the
origin onto the surface of the unit sphere. This is called the Bloch sphere.

Figure 3: Evolution of (4) over a small interval visualised on the Bloch
sphere.

##  2 Current Methods

Here we consider different approaches to solving ordinary differential
equations. We first look at some important standard methods applied to the
relevant ODE and then move onto more advanced approaches using approximations
of the matrix exponential, discussing their rates of convergence, error
bounds, and computational efficiency.

For the initial value problem

$\mathbf{x}^{\prime}(t)=A(t)\mathbf{x}(t),\;\;\;\mathbf{x}(t_{0})=\mathbf{x}_{0},\;\;\;t\in[t_{0},\,t_{f}],$
(5)

we discretise the time interval with a time-step $h=0.5^{k}$, where
$k\in\mathbb{N}$ is specified, and approximate the solution at each such point
in time, beginning with $\mathbf{x}_{0}$ at $t_{0}$. From now on we give the
step-size by specifying the value of $k$. The general solution this approach
gives is

$\mathbf{x}(t_{n+1})=P(t_{n})\,\mathbf{x}(t_{n}),$

for some matrix-valued function $P(t)$. This inductively gives an explicit
form for $\mathbf{x}(t_{n}),$

$\mathbf{x}(t_{n})=P(t_{n})^{n}\,\mathbf{x}_{0}.$

Note how equation (5) is analogous to the vectorised Liouville-von Neumann
equation (3).

###  2.1 QuTiP

A current implementation for solving the Liouville-von Neumann equation is the
open-source Python package QuTiP, which provides methods and classes for
representing and manipulating quantum objects [9]. The function with which we
are concerned is the mesolve function, which evolves a density matrix using a
given Hamiltonian in the Liouville-von Neumann equation, evaulating the
density matrix at the given times.

As opposed to the other methods we consider, where the given times are used to
approximate the density matrix, QuTiP sub-samples the density matrix at the
given times from a much more accurate solution. This enables mesolve to
maintain an accuracy of $10^{-4}$ regardless of how large a time-step is
taken.

Figure 4: QuTiP sub-sampling for three solutions to (4).

###  2.2 Standard ODE Methods

####  2.2.1 Euler’s Method

We begin with the most basic implicit methods for solving ordinary
differential equations. There are two forms of Euler’s method: explicit and
implicit.

  * •

Explicit Euler:

$\mathbf{x}(t_{n+1})=\big{(}I+hA(t_{n})\big{)}\,\mathbf{x}(t_{n}),$

  * •

Implicit Euler:

$\mathbf{x}(t_{n+1})=\big{(}I-hA(t_{n+1})\big{)}^{-1}\mathbf{x}(t_{n}).$

Implicit methods typically require more computations, but result in more
stable solutions for larger time-steps. For Implicit Euler we see that a
matrix inversion is needed, which is incredibly costly and warrants further
numerical methods.

Consider the following system:

$\begin{cases}H(t)=\sigma_{x}+\sigma_{y}+\sigma_{z},\\\
\rho_{0}=\sigma_{x}.\end{cases}$ (6)

Plotting the solution of (6) for these methods and QuTiP’s mesolve for $k=5$,
we see how both Euler methods violate the conservation of energy: the
normalised $\sigma_{x}$ component oscillates between $1$ and $-1$, but these
oscillations increase and decrease over time for the two methods,
respectively.

Figure 5: Normalised component of system $(2)$ in $\sigma_{x}$ direction over
time for $k=5$.

####  2.2.2 The Trapezoidal Rule

Continuing with the same ODE introduced previously, the trapezoidal rule is as
follows:

$\mathbf{x}(t_{n+1})=\bigg{(}\big{(}I-\frac{h}{2}A(t_{n})\big{)}\big{(}I+\frac{h}{2}A(t_{n})\big{)}\bigg{)}^{-1}\mathbf{x}(t_{n}).$

Being a second-order implicit method, the trapezoidal rule is more stable and
converges to our reference solution faster than both of the Euler methods.

Although we won’t employ this method for larger systems, it serves to
highlight how midpoint methods converge faster than their initial point
counterparts. To change the previous methods from initial point to midpoint,
we replace $t_{n}$ with $t_{n}+h/2$.

Figure 6: log-log plot of maximum error against time-step $h$ for the
trapezoidal rule.

From Figure 6 we see that for system $(1)$, initial point trapezoidal rule has
a rate of convergence $1.255\approx 1$, whereas midpoint has a rate of
convergence $1.965\approx 2$.

####  2.2.3 RK4

The most widely known Runge-Kutta method—commonly referred to as RK4—is an
explicit fourth-order method that uses the previous value of $\mathbf{x}$ and
the weighted average of four increments to determine the following value of
$\mathbf{x}$. The system is evolved as follows:

$\mathbf{x}(t_{n+1})=\mathbf{x}(t_{n})+\frac{1}{6}h(k_{1}+2k_{2}+2k_{3}+k_{4}),$

where

$\displaystyle k_{1}$ $\displaystyle=A(t_{n})\,\mathbf{x}(t_{n}),$
$\displaystyle k_{2}$
$\displaystyle=A\left(t_{n}+\frac{h}{2}\right)\left(\mathbf{x}(t_{n})+h\frac{k_{1}}{2}\right),$
$\displaystyle k_{3}$
$\displaystyle=A\left(t_{n}+\frac{h}{2}\right)\left(\mathbf{x}(t_{n})+h\frac{k_{2}}{2}\right),$
$\displaystyle k_{4}$
$\displaystyle=A\left(t_{n}+h\right)\left(\mathbf{x}(t_{n})+hk_{3}\right).$

When the right-hand side of the ODE is independent of $\mathbf{x}$, RK4
reduces to Simpson’s rule [17]. Interestingly, when $A(t)$ is independent of
$t$, the function $P(t)$ for RK4 becomes

$I+hA+\frac{h^{2}}{2!}A^{2}+\frac{h^{3}}{3!}A^{3}+\frac{h^{4}}{4!}A^{4},$

which is the first four terms of the Taylor expansion of the matrix
exponential.

Figure 7: The true solution and the RK4 approximation of the system (4) for
$k=4$.

###  2.3 Approximating the Matrix Exponential

An important function in the solution of differential equations is the matrix
exponential. Here we discuss some more sophisticated methods that can be
implemented in numerically solving ODEs.

####  2.3.1 Taylor Series

This method directly calculates the matrix exponential using the definition

$e^{A}=I+A+A^{2}/2!+\dotsc,$

which is truncated to the first $K$ terms, giving

$T_{K}(A)=\sum_{j=0}^{K}A^{j}/j!$

as the approximation for $e^{A}$.

An arbitrarily high degree of accuracy can be achieved given a large enough
$K$ is used, but even when approximating the scalar exponential this method is
substandard. One may also find example matrices which cause ”catastrophic
cancellation” in floating point arithmetic when using ”short arithmetic”. This
is not, however, a direct consequence of the approximation, but a result of
using insufficiently high relative accuracy in the calculation [19].

Liou [11] describes an approach to calculating $e^{tA}$ via power series,
writing

$e^{tA}=M+R,$

where $M$ is the $K$-term approximation

$\sum_{j=0}^{K}\frac{t^{j}A^{j}}{j!}$

and $R$ the remainder

$\sum_{j=K+1}^{\infty}\frac{t^{j}A^{j}}{j!}.$

In order the achieve an accuracy of $d$ significant digits, we require

$|r_{ij}|\leq 10^{-d}\,|m_{ij}|,$ (7)

for the elements of $R$ and $M$. The upper bound derived by Liou for
$|r_{ij}|$, given $K$ terms in the approximation, was sharpened by Everling
[4] to

$|r_{ij}|\leq\mathinner{\\!\left\lVert\frac{t^{K}A^{K}}{K!}\right\rVert}\,\frac{t\mathinner{\\!\left\lVert
A\right\rVert}}{K+1}\;\frac{1}{1-\epsilon},\;\;\;\epsilon=\frac{t\mathinner{\\!\left\lVert
A\right\rVert}}{K+2},$

which allows for checking whether the desired accuracy defined by (7) has been
achieved.

This approach is much worse than some of the methods addressed later, but is
worth covering briefly for completeness-sake. One issue is that it does not
take advantage of any of the properties of $A$, such as its eigenvalues or
Hermitivity, as other methods do.

####  2.3.2 Padé Approximants

When approximating a function, the Taylor series expansion is limiting since
polynomials are not a good class of functions to use if the approximated
function has singularities. Rational functions are the simplest functions with
singularities, making them an obvious first choice for a better approach [18].

The Padé approximation is a rational function with numerator of degree $\leq
p$ and denominator of degree $\leq q$. We may adapt this approach to
calculating the matrix exponential as follows [19]:

$e^{A}\approx R_{pq}(A)=[D_{pq}(A)]^{-1}N_{pq}(A),$

where

$N_{pq}(A)=\sum_{j=0}^{q}\frac{(p+q-j)!p!}{(p+q)!j!(p-j)!}A^{j}$

and

$D_{pq}(A)=\sum_{j=0}^{q}\frac{(p+q-j)!q!}{(p+q)!j!(q-j)!]}(-A)^{j}.$

Rational approximations of the matrix exponential are preferred over
polynomial approximations due to their better stability properties. When
$p\leq q$ the Padé approximation gives rise to unconditionally stable methods
[20].

When the norm of A is not too large, $p=q$ is a preferred choice for the order
of the numerator and denominator. However, the computational costs and
rounding errors for Padé approximants increase as $\mathinner{\\!\left\lVert
A\right\rVert}$ increases, where $\mathinner{\\!\left\lVert A\right\rVert}$
denotes the matrix operator norm. We may mitigate this by scaling and
squaring, which takes advantage of the following property:

$e^{A}=(e^{A/{2^{j}}})^{2^{j}}.$

This allows a $j$ to be chosen such that the Padé approximation is not too
costly. This is the algorithm used by SciPy’s expm function [1].

Given an error tolerance, we may determine an appropriate $q$ and $j$ to
approximate $e^{A/{2^{j}}}$ using Padé approximants, and multiply to get
$e^{A}$. For these values of $q$ and $j$, we find that Padé approximants are
more efficient than a Taylor approximant. Specifically, for a small
$\mathinner{\\!\left\lVert A\right\rVert}$, the Taylor expansion requires
about twice as much work for the same level of accuracy [19]. The appendices
contains a table of values for which a diagonal Padé approximant is optimal,
given a norm $\mathinner{\\!\left\lVert A\right\rVert}$ and error tolerance
$\epsilon$.

####  2.3.3 Krylov Subspace Methods

It is not always necessary to calculate $e^{A}$ in its entirety. Sometimes we
only need the product of the matrix exponential with a vector. Even if $A$ is
a sparse matrix $e^{A}$ will typically still be dense, so when dealing with
large matrix systems it may be computationally infeasible to calculate $e^{A}$
completely. This is easily seen with quantum systems, where a ten spin system
results in a Liouvillian with dimension $2^{20}$.

Given a matrix $A\in\mathbb{C}^{n\times n}$ and a compatible vector
$\mathbf{b}$, we can approximate $e^{A}\mathbf{b}$ using a Krylov subspace
method. Being iterative methods, they produce a sequence of improving
approximations, and achieve the actual solution (up to arithmetic accuracy)
with $n$ iterations.

We consider an order-$m$ Krylov subspace,

$\mathcal{K}_{m}(A,\mathbf{v})=\text{span}\\{\mathbf{v},A\mathbf{v},\dotsc,A^{m-1}\mathbf{v}\\},$

where
$\mathbf{v}=\mathbf{b}/\mathinner{\\!\left\lVert\mathbf{b}\right\rVert}$, and
wish to find an orthonormal basis. The approximation to $e^{A}\mathbf{b}$ is
an element of this space. For Hermitian matrices, we employ the Lanczos
algorithm with $m\leq n$ iterations [10]. The algorithm outputs an $n\times m$
matrix $V_{m}$ with orthonormal columns and an $m\times m$ real, symmetric,
tridiagonal matrix $T_{m}$ such that $T_{m}=V_{m}^{\dagger}\,A\,V_{m}.$ For
$m=n$, $V_{m}$ is a unitary matrix and so $A\approx
V_{m}\,T_{m}\,V_{m}^{\dagger}.$ Then,

$e^{A}\,\mathbf{b}\,\approx\,\mathinner{\\!\left\lVert\mathbf{b}\right\rVert}V_{m}\,e^{T_{m}}\,e_{1},$

where $e_{1}$ is the first column of the $m$-dimensional identity matrix.
Since $T_{m}$ is tridiagonal, Padé approximants can be used to calculate
$e^{T_{m}}$ efficiently [6].

Krylov methods are preferred over polynomial or rational function
approximations for large and/or sparse matrices, since it replaces calculating
$e^{A}$ with calculating the much easier $e^{T_{m}}$. This is achieved through
approximating the eigenvalues of $A$ with the eigenvalues of $H_{m}$.

For a negative semi-definite Hermitian matrix $A$, with eigenvalues in the
interval $\left[-4\rho,0\right]$, a scalar $t$, and a unit vector $v$,
Hochbruck and Lubich [7] derived a bound for the error,

$\varepsilon_{m}\mathrel{\mathop{\mathchar
58\relax}}=\mathinner{\\!\left\lVert
e^{tA}v-V_{m}e^{tT_{m}}e_{1}\right\rVert}.$

The error is bounded as follows:

$\varepsilon_{m}\leq\begin{cases}10\,e^{-m^{2}/(5\rho
t)},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\sqrt{4\rho t\vphantom{|}}\leq m\leq
2\rho t,\\\ 10(\rho t)^{-1}e^{-\rho t}\Big{(}\dfrac{e\rho
t}{m}\Big{)}^{m},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;m\geq 2\rho
t.\end{cases}$ (8)

Figure 8: The error bound (8) and the experimental error of our example.

Following the approach of [7], we consider a diagonal matrix with eigenvalues
in the interval $[-40,0]$ and a 1001-dimensional unit vector. This error bound
and the error in the approximation from our implementation—see appendix—is
shown in Figure 8.

##  3 The Magnus Expansion

In this section we aim to implement a method which converges faster than
current approaches by applying the Magnus expansion to the differential
equation (5), providing an exponential representation of the solution to such
an equation [13]. Solutions derived from the Magnus expansion give
approximations which preserve certain qualitative properties of the exact
solution, making it a desirable method to consider [2].

We first discuss the theory behind the approach, detailing the components
relevant to our implementation, and then optimise the solution to best solve
ODEs driven by Hamiltonians in the form of (2). In addition to this, we
describe how our code was implemented in the Python package MagPy and give an
overview of how to use it.

###  3.1 Definition

Given an $n\times n$ time-dependent matrix $A(t)$, the Magnus expansion
provides a solution to the initial value problem

$\mathbf{x}^{\prime}(t)=A(t)\,\mathbf{x}(t),\;\;\mathbf{x}(t_{0})=\mathbf{x}_{0}.$

for a vector $\mathbf{x}(t)$. The solution takes the form

$\mathbf{x}(t)=\exp{\left(\Omega(t,t_{0})\right)}\,\mathbf{x_{0}},$

where $\Omega(t,t_{0})$ is the infinite series,

$\Omega(t,t_{0})=\sum_{k=1}^{\infty}\Omega_{k}(t),$

where $\Omega_{k}(t)$ are matrix-valued functions. We will apply this
approximation iteratively, using the same approach as described at the
beginning of section 4. For $t_{0}=0$ we write $\Omega(t,0)$ as $\Omega(t)$.
Arbitrarily high degrees of accuracy can be achieved given a sufficient number
of terms in the summation, but we will truncate to the first two terms,

$\displaystyle\Omega_{1}(t)$ $\displaystyle=\int_{0}^{t}A(s)\,ds,$ (9)
$\displaystyle\Omega_{2}(t)$
$\displaystyle=\int_{0}^{t}\int_{0}^{s}\left[A(s),A(r)\right]dr\,ds.$ (10)

When $A(s)$ and $A(r)$ commute $\Omega(t)$ reduces to $\Omega_{1}(t)$, and for
a time independent $A$ the commutator terms cancel and this reduces to the
well-known solution

$\mathbf{x}(t)=\exp\\{\left(t-t_{0})A\right)\\}\mathbf{x_{0}}.$

###  3.2 Reformulation

Since the Liouvilian takes the form of a matrix-valued function, it is
difficult to implement the integral terms in the Magnus expansion directly.
The core idea here is to write the integrals in terms of scalar integrals of
$f_{j}(t)$, $g_{j}(t)$, and $\Omega_{j}$, reducing the complexity of the code
and the number of computations needed. We consider the two terms separately,
detailing the necessary properties used in implementing the Magnus method in
Python. For the first term the following theorem is sufficient to write the
integral in the required form.

######  Theorem 3.1 (1st Term of Magnus Expansion).

For $A(t)=-iL\left\\{H(t)\right\\},$

$\int_{0}^{t}A(s)\,ds=-iL\left\\{\sum_{j=1}^{n}K_{j}(s)+H_{J}\,t\right\\},$

where

$K_{j}(s)=\int_{0}^{t}f_{j}(s)\,ds\,\mathbf{I}_{j}^{x}+\int_{0}^{t}g_{j}(s)\,ds\,\mathbf{I}_{j}^{y}+\Omega_{j}\,t\,\mathbf{I}_{j}^{z}.$

###### Proof.

Lemma 1.7 gives us

$\int_{0}^{t}A(s)\,ds=-iL\left\\{\int_{0}^{t}H(s)\,ds\right\\},$

and applying Lemma 1.2 to $\int_{0}^{t}H(s)\,ds$ gives the result. ∎

####  3.2.1 Derivation of the Second Term

For the second term of the expansion we are concerned with calculating a
double integral of the commutator,

$\big{[}A(s),A(r)\big{]},$

where

$A(t)=-iL\left\\{H(t)\right\\}.$

We build up to this through the following results. Theorem 1.9 reduces this to
first calculating the commutator of

$H(t)=\sum_{j=1}^{n}\mathbf{I}_{j}\left\\{H_{j}(t)\right\\}+H_{J},$

integrating, and then applying the Liouvillian superoperator.

######  Lemma 3.2.

For
$H_{j}(t)=f_{j}(t)\,\sigma_{x}+g_{j}(t)\,\sigma_{y}+\Omega_{j}\,\sigma_{z},\;j\in\mathbb{N},$

$\mathbf{I}_{j}\left\\{\big{[}H_{j}(r),\,H_{j}(s)\big{]}\right\\}=2\,i\,\Omega_{j}\big{(}g_{j}(r)-g_{j}(s)\big{)}\mathbf{I}_{j}^{x}\\\
+2\,i\,\Omega_{j}\big{(}f_{j}(s)-f_{j}(r)\big{)}\mathbf{I}_{j}^{y}\\\
+2\,i\,\big{(}f_{j}(r)g_{j}(s)-g_{j}(r)f_{j}(s)\big{)}\mathbf{I}_{j}^{z}.$

###### Proof.

Expanding $\big{[}H_{j}(r),\,H_{j}(s)\big{]}$, applying the commutator
properties of the Pauli matrices, and simplifying gives

$\big{[}H_{j}(r),\,H_{j}(s)\big{]}=2\,i\,\Omega_{j}\big{(}g_{j}(r)-g_{j}(s)\big{)}\sigma_{x}+2\,i\,\Omega_{j}\big{(}f_{j}(s)-f_{j}(r)\big{)}\sigma_{y}\\\
+2\,i\,\big{(}f_{j}(r)g_{j}(s)-g_{j}(r)f_{j}(s)\big{)}\sigma_{z}.$

Then the linearity of $\mathbf{I}$ yields the result. ∎

######  Corollary 3.2.1.

$\displaystyle\int_{0}^{t}\int_{0}^{s}\mathbf{I}_{j}\left\\{\big{[}H_{j}(r),\,H_{j}(s)\big{]}\right\\}dr\,ds$
$\displaystyle=2\,i$
$\displaystyle\,\Omega_{j}\int_{0}^{t}\int_{0}^{s}g_{j}(r)-g_{j}(s)\,dr\,ds\,\mathbf{I}_{j}^{x}$
$\displaystyle+$ $\displaystyle
2\,i\,\Omega_{j}\int_{0}^{t}\int_{0}^{s}f_{j}(s)-f_{j}(r)\,dr\,ds\,\mathbf{I}_{j}^{y}$
$\displaystyle+$ $\displaystyle
2\,i\int_{0}^{t}\int_{0}^{s}f_{j}(r)g_{j}(s)-g_{j}(r)f_{j}(s)\,dr\,ds\,\mathbf{I}_{j}^{z}.$

###### Proof.

This follows directly from Lemma 3.2 and from linearity of integration. ∎

Now we look at the commutator of
$H_{0}(t)=\sum_{j=1}^{n}\mathbf{I}_{j}\left\\{H_{j}(t)\right\\}$ and prove a
useful result

######  Lemma 3.3.

Given $H_{0}(t)=\sum_{j=1}^{n}\mathbf{I}_{j}\left\\{H_{j}(t)\right\\},$

$\big{[}H_{0}(r),\,H_{0}(s)\big{]}=\sum_{j=1}^{n}\mathbf{I}_{j}\left\\{\left[H_{j}(r),H_{j}(s)\right]\right\\}.$

###### Proof.

By linearity of summation and Lemma 1.5.1,

$\displaystyle\big{[}H_{0}(r),\,H_{0}(s)\big{]}$
$\displaystyle=\sum_{i,\,j=1}^{n}\left(\mathbf{I}_{i}\big{\\{}H_{i}(r)\big{\\}}\,\mathbf{I}_{j}\big{\\{}H_{j}(s)\big{\\}}-\mathbf{I}_{j}\big{\\{}H_{j}(s)\big{\\}}\,\mathbf{I}_{i}\big{\\{}H_{i}(r)\big{\\}}\right),$
$\displaystyle=\sum_{i,\,j=1}^{n}\left(\mathbf{S}_{i,\,j}\big{\\{}H_{i}(r),H_{j}(s)\big{\\}}-\mathbf{S}_{j,\,i}\big{\\{}H_{j}(s),H_{i}(r)\big{\\}}\right).$

Applying Lemma 1.5.2,

$\mathbf{S}_{i,\,j}\big{\\{}H_{i}(r),H_{j}(s)\big{\\}}=\mathbf{S}_{j,\,i}\big{\\{}H_{j}(s),H_{i}(r)\big{\\}},$

for $i\neq j$. And so, by Definition 1.4,

$\displaystyle\big{[}H_{0}(r),\,H_{0}(s)\big{]}$
$\displaystyle=\sum_{j=1}^{n}\left(\mathbf{S}_{j,\,j}\big{\\{}H_{j}(r),H_{j}(s)\big{\\}}-\mathbf{S}_{j,\,j}\big{\\{}H_{j}(s),H_{j}(r)\big{\\}}\right)$
$\displaystyle=\sum_{j=1}^{n}\left(\mathbf{I}_{j}\big{\\{}H_{j}(r)\,H_{j}(s)\big{\\}}-\mathbf{I}_{j}\big{\\{}H_{j}(s)\,H_{j}(r)\big{\\}}\right)$
$\displaystyle=\sum_{j=1}^{n}\mathbf{I}_{j}\left\\{\big{[}H_{j}(r),\,H_{j}(s)\big{]}\right\\}.$

∎

######  Corollary 3.3.1.

$\int_{0}^{t}\int_{0}^{s}\big{[}H_{0}(r),\,H_{0}(s)\big{]}\,dr\,ds\\\
=\sum_{j=1}^{n}\bigg{(}2\,i\,\Omega_{j}\int_{0}^{t}\int_{0}^{s}g_{j}(r)-g_{j}(s)\,dr\,ds\,\mathbf{I}_{j}^{x}\\\
+2\,i\,\Omega_{j}\int_{0}^{t}\int_{0}^{s}f_{j}(s)-f_{j}(r)\,dr\,ds\,\mathbf{I}_{j}^{y}\\\
+2\,i\int_{0}^{t}\int_{0}^{s}f_{j}(r)g_{j}(s)-g_{j}(r)f_{j}(s)\,dr\,ds\,\mathbf{I}_{j}^{z}\bigg{)}.$

###### Proof.

Integrate Lemma 3.3 and apply Corollary 3.2.1. ∎

Next, we consider the commutator of $H(t)=H_{0}(t)+H_{J}.$ By expanding and
simplifying we can write

$\big{[}H(r),H(s)\big{]}=\big{[}H_{0}(r),H_{0}(s)\big{]}+\big{[}H_{J},H_{0}(s)\big{]}+\big{[}H_{0}(r),H_{J}\big{]}.$

We have already dealt with the first term on the right hand side, and now deal
with the other two terms. The following lemma and corollary tell us how.

######  Lemma 3.4.

Given $H_{0}(t)=\sum_{j=1}^{n}\mathbf{I}_{j}\left\\{H_{j}(t)\right\\}$ and
$H_{J}\in\mathbb{C}^{n^{2}\times n^{2}},$

$\big{[}H_{J},H_{0}(s)\big{]}+\big{[}H_{0}(r),H_{J}\big{]}=\sum_{j=1}^{n}\bigg{(}\left(f_{j}(r)-f_{j}(s)\right)\big{[}\mathbf{I}_{j}^{x},H_{J}\big{]}\\\
+\left(g_{j}(r)-g_{j}(s)\right)\big{[}\mathbf{I}_{j}^{y},H_{J}\big{]}\bigg{)}.$

###### Proof.

Note that

$\displaystyle\big{[}H_{J},H_{0}(s)\big{]}$
$\displaystyle=\sum_{j=1}^{n}\bigg{(}H_{J}\big{(}f_{j}(s)\mathbf{I}_{j}^{x}+g_{j}(s)\mathbf{I}_{j}^{y}+\Omega_{j}\mathbf{I}_{j}^{z}\big{)}$
$\displaystyle-\big{(}f_{j}(s)\mathbf{I}_{j}^{x}+g_{j}(s)\mathbf{I}_{j}^{y}+\Omega_{j}\mathbf{I}_{j}^{z}\big{)}H_{J}\bigg{)},$
$\displaystyle=\sum_{j=1}^{n}f_{j}(s)\big{[}H_{J},\mathbf{I}_{j}^{x}\big{]}+g_{j}(s)\big{[}H_{J},\mathbf{I}_{j}^{y}\big{]}+\Omega_{j}\big{[}H_{J},\mathbf{I}_{j}^{z}\big{]},$

and similarly that

$\displaystyle\big{[}H_{0}(r),H_{J}\big{]}$
$\displaystyle=\sum_{j=1}^{n}f_{j}(r)\big{[}\mathbf{I}_{j}^{x},H_{J}\big{]}+g_{j}(r)\big{[}\mathbf{I}_{j}^{y},H_{J}\big{]}+\Omega_{j}\big{[}\mathbf{I}_{j}^{z},H_{J}\big{]}.$

Summing these two equations gives the result. ∎

######  Corollary 3.4.1.

$\int_{0}^{t}\int_{0}^{s}\big{[}H_{J},H_{0}(s)\big{]}+\big{[}H_{0}(r),H_{J}\big{]}\,dr\,ds\\\
=\sum_{j=1}^{n}\bigg{(}\int_{0}^{t}\int_{0}^{s}f_{j}(r)-f_{j}(s)\,dr\,ds\big{[}\mathbf{I}_{j}^{x},H_{J}\big{]}\\\
+\int_{0}^{t}\int_{0}^{s}g_{j}(r)-g_{j}(s)\,dr\,ds\big{[}\mathbf{I}_{j}^{y},H_{J}\big{]}\bigg{)}.$

###### Proof.

This follows directly from Lemma 3.4 and linearity of integration. ∎

We now have all we need to write the main result for the second term,

$\Omega_{2}(t)=\frac{1}{2}\int_{0}^{t}\int_{0}^{s}\left[A(s),A(r)\right]dr\,ds.$

######  Theorem 3.5 (2nd Term of Magnus Expansion).

Given

$H(t)=\sum_{j=1}^{n}f_{j}(t)\mathbf{I}_{j}^{x}+g_{j}(t)\mathbf{I}_{j}^{y}+\Omega_{j}\mathbf{I}_{j}^{z}+H_{J},$

we can rewrite

$\Omega_{2}(t)=\frac{1}{2}\int_{0}^{t}\int_{0}^{s}\left[A(s),A(r)\right]dr\,ds$

as

$\Omega_{2}(t)=\frac{1}{2}L\Bigg{\\{}\sum_{j=1}^{n}\bigg{(}\int_{0}^{t}\int_{0}^{s}g_{j}(r)-g_{j}(s)\,dr\,ds\,\big{(}2\,i\,\Omega_{j}\,\mathbf{I}_{j}^{x}+\big{[}\mathbf{I}_{j}^{y},H_{J}\big{]}\big{)}\\\
-\int_{0}^{t}\int_{0}^{s}f_{j}(r)-f_{j}(s)\,dr\,ds\,\big{(}2\,i\,\Omega_{j}\mathbf{I}_{j}^{y}+\big{[}H_{J},\mathbf{I}_{j}^{x}\big{]}\big{)}\\\
+2\,i\int_{0}^{t}\int_{0}^{s}f_{j}(r)g_{j}(s)-g_{j}(r)f_{j}(s)\,dr\,ds\,\mathbf{I}_{j}^{z}\bigg{)}\Bigg{\\}}.$

###### Proof.

We begin with

$\Omega_{2}(t)=\frac{1}{2}\int_{0}^{t}\int_{0}^{s}\left[A(s),A(r)\right]dr\,ds.$

Theorem 1.9 gives

$\frac{1}{2}L\big{\\{}\int_{0}^{t}\int_{0}^{s}\left[H(r),H(s)\right]dr\,ds\big{\\}}.$

Writing

$\left[H(r),H(s)\right]=\left[H_{0}(r),H_{0}(s)\right]+\left[H_{J},H_{0}(s)\right]+\left[H_{0}(r),H_{J}\right],$

we apply Corollary 3.3.1 and Corollary 3.4.1 and simplify to get the result. ∎

###  3.3 Implementation in Python

The results of Theorems 3.1 and 3.5 are applied in the MagPy function
lvnsolve. Given a Hamiltonian of the form (2), an initial density matrix
$\rho_{0}$, and a list of times, this function evolves a density matrix under
the Liouville-von Neumann equation and returns the density matrix evaluated
across the specified times, starting at $\rho_{0}$. Currently, the integrals
in the Magnus expansion are calculated using SciPy’s integrate package.

We introduce the function linspace to construct the list of times. This is a
simple wrapper of the NumPy linspace function, using a given step-size instead
of a specified number of points. MagPy also has built-in definitions of the
Pauli matrices, which are $2\times 2$ NumPy ndarrays: $\sigma_{x}$ is denoted
sigmax and the rest are defined analogously. The following code snippet shows
how to apply lvnsolve to a two spin time-dependent Hamiltonian.

⬇

import numpy as np

import magpy as mp

def f1(t): return t

def g1(t): return np.sqrt(t)

omega1 = 1

def f2(t): return -t**2

def g2(t): return np.sin(t)

omega2 = 2

H_coeffs = [[f1, g1, omega1], [f2, g2, omega2]]

HJ = np.kron(mp.sigmax(), mp.sigmay())

rho0 = np.kron(mp.sigmax(), mp.sigmax())

tlist = mp.linspace(0, 20, 0.5**10)

density_matrices = mp.lvnsolve(H_coeffs, rho0, tlist, HJ)

Figure 9: Example of the lvnsolve function.

##  4 Numerical Comparison

In this section we compare the rate of convergence for different variations on
Magnus-based techniques for simulating spin systems displaying HOCP. As a
reference solution we use one-term Magnus with the integrals approximated
using the midpoint method:

$\int_{a}^{b}f(t)\,dt=(b-a)\,f\left(\frac{a+b}{2}\right).$

Our interval of simulation is $[0,20]$ unless otherwise specified, with a
time-step defined by $k=20$. We present our analysis using log-log plots for
the maximum error of our approximation from the chosen reference against the
time-step.

###  4.1 One-term Magnus

Recall for a single spin density matrix driven by the Liouville-von Neumann
equation with a Hamiltonian of the form (2), one-term Magnus approximates the
density matrix as follows:

$\mathbf{x}(t_{n+1})=\exp{\left\\{\Omega(t_{n+1},t_{n})\right\\}}\,\mathbf{x}(t_{n}),$

where

$\Omega(t_{n+1},t_{n})=\int_{t_{n}}^{t_{n+1}}f_{1}(s)\,ds\,\sigma_{x}+\int_{t_{n}}^{t_{n+1}}g_{1}(s)\,ds\,\sigma_{y}+\Omega_{1}(t_{n+1}-t_{n})\,\sigma_{z}.$

Figure 10: Functions $f_{1}(t)$ and $g_{1}(t)$ for (11).

We apply our analysis of one-term Magnus to the HOCP system defined as
follows:

$\begin{cases}\beta_{1}=10,\;\;\gamma_{1}=2,\;\;\Omega_{1}=1,\\\
\rho_{0}=\sigma_{x}.\end{cases}$ (11)

The integral approximations that we consider are initial point, midpoint, and
SciPy’s integrate.quad, which uses a Clenshaw-Curtis method with Chebyshev
moments [21]. SciPy’s method will be treated as the exact integral. We compare
the error in these approximations for (11) with $\beta_{1}$ and $\Omega_{1}$
scaled by 1, 2, $\frac{1}{10}$, and $\frac{1}{100}$. Physically this is
equivalent to changing the Hamiltonian’s field strength.

We can clearly see in Figure 11 the advantage of using a midpoint
approximation over initial point. The midpoint approximation maintains second-
order convergence, whereas initial point displays approximately first-order
convergence for all simulations.

Figure 11: Error against time-step for (11) with the Hamiltonian’s field
strength scaled four different ways.

Interestingly, midpoint approximations for the integral display second order
convergence, which is equivalent to the rate of convergence when using the
exact integral. However, midpoint begins to converge at a slightly larger
time-step than when using the exact integral. For scales $\frac{1}{10}$ and
$\frac{1}{100}$ we see that the exact integral method initially converges
faster, but is overtaken by midpoint at a time-step of size $10^{-1}$.

It must be emphasised how surprising these results are. We see that it is not
necessary to use more accurate approximation techniques for the integrals when
the field strength is stronger. Plus, the midpoint approximations actually
perform better than the exact integral. The complete opposite behaviour is
expected.

We also see that the smaller the scale of the system, the smaller the error is
in the approximation, with a maximum error of $10^{-6}$ being achieved when
the system is scaled by $\frac{1}{100}.$

###  4.2 Two-term Magnus

For two-term of Magnus, there are multiple ways to approximate the integrals
in the calculation. The double integrals will all be approximated using
SciPy’s integrate.quad function being treated as the exact integral, and the
single integrals will be approximated using the following:

  1. 1.

SciPy’s integrate.quad function,

  2. 2.

Midpoint,

  3. 3.

Gauss-Legendre quadrature of order 3.

Gauss-Legendre quadrature approximates integrals using quadrature weights and
roots of Legendre polynomials [8]. Order 3 gives the following approximation:

$\int_{a}^{b}f(x)\,dx\approx\frac{b-a}{2}\sum_{i=1}^{3}w_{i}\,f\left(\frac{b-a}{2}\xi_{i}+\frac{a+b}{2}\right),$

where

$\displaystyle\xi_{1}$
$\displaystyle=0,\;\;\;\;\;\;\;\;\;\,w_{1}=\frac{8}{9},$
$\displaystyle\xi_{2}$
$\displaystyle=\sqrt{\frac{3}{5}},\;\;\;\;\;\,w_{2}=\frac{5}{9},$
$\displaystyle\xi_{3}$
$\displaystyle=-\sqrt{\frac{3}{5}},\;\;\;w_{3}=\frac{5}{9}.$

We compare the convergence of these methods against one-term Magnus using
midpoint integral approximations, using midpoint as a reference solution.
System (11)

In figure 11 we see that both two-term Magnus with SciPy and with order 3
Gauss-Legendre Quadrature display fourth-order convergence. When the system is
scaled by $\frac{1}{10}$ and $\frac{1}{100}$, we achieve the accuracy of one-
term Magnus with midpoint with a larger time-step of $10^{-2}$. The log-log
plots plateau at a maximum error of $10^{-9}$.

Up to a time-step of $10^{-1}$ the exact integral method converges faster than
Gauss-Legendre quadrature, but then they display an equal rate of convergence
from that point onward. This behaviour is not seen for the system scaled by 1
or 2.

Using the midpoint method to approximate the first term in the two-term Magnus
expansion is not any better than one-term Magnus using a midpoint
approximation.

Figure 12: Error against time-step for (11) with the Hamiltonian’s field
strength scaled four different ways.

###  4.3 Larger Systems

Consider the following system of two spins:

$i$ | $\beta_{i}$ | $\gamma_{i}$ | $\Omega_{i}$
---|---|---|---
1 | 10 | 2 | 5
2 | -40 | 25 | -12

The interacting component and initial condition are

$H_{J}=\sigma_{x}\otimes\sigma_{y},\;\;\;\;\;\rho_{0}=\sigma_{x}\otimes
I+I\otimes\sigma_{y}.$

For the non-interacting case, we can measure the z-component of the first spin
without any interference from the second spin.

Figure 13: The normalised z-component of the first spin of system of the two-
spin system.

It is not possible to isolate a single spin in an interacting multi-spin
system, but we can use a compound operator to visualise the evolution of the
system. Specifically, we measure the proportion of the total density matrix in
the direction $\sigma_{x}\otimes\sigma_{x}.$ Figure 14 highlights the highly-
oscillatory nature of the system, showing how an insufficiently small time-
step will cause these oscillations to be missed in the simulation.

Figure 14: A snippet of the evolution of the compound spin direction of the
two particle system.

##  5 Discussion

###  5.1 Summary

In this project we have implemented a more efficient method for simulation of
quantum spin systems displaying highly oscillatory chirped pulses using a
truncated Magnus expansion. The approach has been optimised to handle a
specific Hamiltonian common in the work of NMR and MRI. We opted for the
Liouvillian picture in which to evolve our spin systems, and thus proved
properties of the Liouvillian superoperator to write the matrix-valued
integrals of the Magnus expansion in terms of scalar-valued integrals.

Using the midpoint integral approximation as a reference, we compared the
effectiveness of current numerical methods and approximations of the matrix
exponential, looking at their rates of convergence and error bounds. This was
then compared against our implementation of the Magnus expansion, resulting in
us achieve fourth-order convergence—twice as fast as current implementations.
We looked at the same system scaled four different ways and found that for
systems with smaller oscillations our implementation achieves higher accuracy.
We also discussed the Python package QuTiP and its approach to simulating the
same spin systems.

This has been collated—along with other relevant functions pertaining o spin
system simulation—in the Python package MagPy. The packages main function,
lvnsolve, can evolve a density matrix given an initial condition and a
Hamiltonian over a specified list of points in time.

###  5.2 Further Work

A next step in the expansion of the work of this report would be to adapt the
function lvnsolve to handle any matrix-valued function as a Hamiltonian. This
would enable more systems than those relevant to NMR or MRI to take advantage
of the increased accuracy and efficiency of MagPy. In addition to this, we
would like to test scaling and splitting methods in conjunction with the
Magnus expansion to further increase the effectiveness of our implementation.

It is also worth investigating further the convergence of different integral
approximations. Specifically that more accurate techniques for integration do
not necessarily result in faster convergence to our reference solution. We do
not expect the midpoint approximation to perform better than the SciPy
package’s more sophisticated approximation techniques.

Testing of systems with more spins, both interacting and not, would aid in
further verifying the effectiveness of our implementation, and potentially
find more systems to which our method can be adapted. Alongside this, testing
of the speed and number of computations required by lvnsolve would determine
the complexity of the function.

## Appendix

### The Kronecker Product

The Kronecker product of $A\in\mathbb{C}^{m\times n}$ and
$B\in\mathbb{C}^{p\times q}$ is

$A\otimes B=\begin{bmatrix}a_{11}B&\cdots&a_{1n}B\\\ \vdots&\ddots&\vdots\\\
a_{m1}B&\cdots&a_{mn}B\end{bmatrix}\in\mathbb{C}^{pm\times qn}.$

The mixed-product property of the Kronecker product is that for matrices
$A,B,C,D$ such that the products $AC$ and $BD$ can be formed, we have

$\left(A\otimes B\right)\left(C\otimes D\right)=AC\otimes BD.$

### An Error Bound of Scaling and Squaring

If $||A||\leq 2^{j-1},$ then

$[R_{qq}(A)]^{2j}=e^{A+E},$

where $E$ is a matrix such that

$\frac{\mathinner{\\!\left\lVert E\right\rVert}}{\mathinner{\\!\left\lVert
A\right\rVert}}\leq 8\left[\frac{\mathinner{\\!\left\lVert
A\right\rVert}}{2^{j}}\right]^{2q}\left(\frac{(q!)^{2}}{(2q)!(2q+1)!}\right).$

And so, we may choose a $(q,j)$ pair such that

$\frac{\mathinner{\\!\left\lVert E\right\rVert}}{\mathinner{\\!\left\lVert
A\right\rVert}}\leq\epsilon,$

for a chosen error tolerance $\epsilon.$

$\mathinner{\\!\left\lVert A\right\rVert}$ $\epsilon$ | $10^{-3}$ | $10^{-6}$ | $10^{-9}$ | $10^{-12}$ | $10^{-15}$
---|---|---|---|---|---
$10^{-2}$ | $(1,0)$ | $(1,0)$ | $(2,0)$ | $(3,0)$ | $(3,0)$
$10^{-1}$ | $(1,0)$ | $(2,0)$ | $(3,0)$ | $(4,0)$ | $(4,0)$
$10^{0}$ | $(2,1)$ | $(3,1)$ | $(4,1)$ | $(5,1)$ | $(6,1)$
$10^{1}$ | $(2,5)$ | $(3,5)$ | $(4,5)$ | $(5,5)$ | $(6,5)$
$10^{2}$ | $(2,8)$ | $(3,8)$ | $(4,8)$ | $(5,8)$ | $(6,8)$
$10^{3}$ | $(2,11)$ | $(3,11)$ | $(4,11)$ | $(5,11)$ | $(6,11)$
Table 1: Optimal values for $(q,j)$ with diagonal Padé.

Since $[R_{qq}(A/{2^{j}})]^{2^{j}}$ requires about $(q+j+\frac{1}{3})n^{3}$
flops to evaluate, it makes sense to choose a $(q,j)$ pair such that $q+j$ is
minimum [19]. Table (1) shows the optimal scaling and squaring parameters for
a fixed $\mathinner{\\!\left\lVert A\right\rVert}$ and error tolerance
$\epsilon$.

### The Lanczos Algorithm

Given a Hermitian matrix $A$ of size $n\times n$, a column vector $b$ of size
$n$, and a number of iterations $m$,

  1. 1.

$v_{1}=b/\mathinner{\\!\left\lVert b\right\rVert}_{2}$.

  2. 2.

$w_{1}=Av_{1}-\langle Av_{1},v_{1}\rangle v_{1}$.

  3. 3.

For $j=2,\ldots,m$, do

    1. a)

$\beta_{j}=\mathinner{\\!\left\lVert w_{j-1}\right\rVert}$.

    2. b)

$v_{j}=w_{j-1}/\beta_{j}$.

    3. c)

$\alpha_{j}=\langle Av_{j},v_{j}\rangle$.

    4. d)

$w_{j}=Av_{j}-\alpha_{j}v_{j}-\beta_{j}v_{j-1}$.

  4. 4.

V is the matrix with columns $v_{1},\dots,v_{m}$, and

$T=\begin{pmatrix}\alpha_{1}&\beta_{2}&&&&0\\\
\beta_{2}&\alpha_{2}&\beta_{3}&&&\\\ &\beta_{3}&\alpha_{3}&\ddots&&\\\
&&\ddots&\ddots&\beta_{m-1}&\\\ &&&\beta_{m-1}&\alpha_{m-1}&\beta_{m}\\\
0&&&&\beta_{m}&\alpha_{m}\end{pmatrix}.$

⬇

import numpy as np

import scipy as sp

def lanczos(A, b, m=None):

n = A.shape[0]

if m is None:

m = n

V = np.zeros((n, m), dtype=’complex’)

W = np.zeros((n, m), dtype=’complex’)

alpha = np.zeros((m, 1), dtype=’complex’)

beta = np.zeros((m, 1), dtype=’complex’)

V[:, 0] = b / np.linalg.norm(b)

W[:, 0] = A @ V[:, 0]

alpha[0] = np.vdot(W[:, 0], V[:, 0])

W[:, 0] = W[:, 0] - alpha[0]*V[:, 0]

for j in range(1, m):

beta[j] = np.linalg.norm(W[:, j-1])

V[:, j] = W[:, j-1] / beta[j]

W[:, j] = A @ V[:, j]

alpha[j] = np.vdot(W[:, j], V[:, j])

W[:, j] = W[:, j] - alpha[j]*V[:, j] - beta[j]*V[:, j-1]

return V, np.diagflat(alpha) + np.diagflat(beta[1:], 1) +
np.diagflat(beta[1:], -1)

def krylov_expm(A, b, m=None):

V, T = lanczos(A, b, m)

return np.linalg.norm(b) * V @ sp.linalg.expm(T) @ np.eye(1, T.shape[0])[0]

Figure 15: Lanczos algorithm and Krylov subspace approximation.

## Bibliography

  * [1] A. H. Al-Mohy and N. J. Higham.  A New Scaling and Squaring Algorithm for the Matrix Exponential.  The University of Manchester, 2009.
  * [2] S. Blanes, F. Casas, J. A. Oteo, and J. Ros.  The Magnus expansion and some of its applications.  Physics Reports, 2008.
  * [3] B. J. Broxson.  The Kronecker Product.  University of North Florida, 2006.
  * [4] W. Everling and M. L. Liou.  On the evaluation of $e^{AT}$ by power series.  IEEE, 1967.
  * [5] M. Foroozandeh.  Spin dynamics during chirped pulses: applications to homonuclear decoupling and broadband excitation.  Journal of Magnetic Resonance, 2020.
  * [6] E. Gallopoulous and Y. Saad.  On the paralell solution of parabolic equations.  RIACS, 1989.
  * [7] M. Hochbruck and C. Lubich.  On Krylov subspace approximations to the matrix exponential operator.  SIAM, 1997.
  * [8] C. G. J. Jacobi.  Ueber Gauß neue Methode, die Werthe der Integrale näherungsweise zu finden.  Journal für Reine und Angewandte Mathematik, 1826.
  * [9] J. R. Johansson, P. D. Nation, and F. Nori.  QuTiP 2: A Python framework for the dynamics of open quantum systems.  Computer Physics Communications, 2013.
  * [10] C. Lanczos.  An Iteration Method for the Solution of the Eigenvalue Problem of Linear Differential and Integral Operators.  Journal of Research of the National Bureau of Standards, 1950.
  * [11] M. L. Liou.  A Novel Method of Evaluating Transient Response.  IEEE, 1966.
  * [12] H. D. Macedo and J. N. Oliveira.  Typing Linear Algebra: A Biproduct-oriented Approach.  Science of Computer Programming, 2013.
  * [13] W. Magnus.  On the Exponential Solution of Differential Equations for a Linear Operator.  Communications on Pure and Applied Mathematics, 1954.
  * [14] G. Mazzi.  Numerical Treatment of the Liouville-von Neumann Equation for Quantum Spin Dynamics.  University of Edinburgh, 2010.
  * [15] W. Pfeifer.  The Lie Algebra su(N).  Birkhäuser, Basel, 2003.
  * [16] J. Randall, A. M. Lawrence, S. C. Webster, S. Weidt, N. V. Vitanov, and W. K. Hensinger.  Generation of high-fidelity quantum control methods for multilevel systems.  Physical Review A, 2018.
  * [17] E. Süli and D. F. Mayers.  An Introduction to Numerical Analysis.  Cambridge University Press, 2003.
  * [18] W. Van Assche.  Padé and Hermite-Padé approximation and orthogonality.  Surveys in Approximation Theory, 2006.
  * [19] C. Van Loan and C. Moler.  Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later.  SIAM Review, 2003.
  * [20] R. S. Varga.  On higher order stable implicit methods for solving parabolic partial differential equations.  Journal of Mathematical Physics, 1961.
  * [21] P. Virtanen et al.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 2020.