Hamiltonian simulation

“The problem of simulating the dynamics of quantum systems was the original motivation for quantum computers and remains one of their major potential applications.” - Berry et al. [1]

The simulation of atoms, molecules and other biochemical systems is another application uniquely suited to quantum computation. For example, the ground state energy of large systems, the dynamical behaviour of an ensemble of molecules, or complex molecular behaviour such as protein folding, are often computationally hard or downright impossible to determine via classical computation or experimentation [2][3].

In the discrete-variable qubit model, efficient methods of Hamiltonian simulation have been discussed at-length, providing several implementations depending on properties of the Hamiltonian, and resulting in a linear simulation time [4][5]. Efficient implementations of Hamiltonian simulation also exist in the CV formulation [6], with specific application to Bose-Hubbard Hamiltonians (describing a system of interacting bosonic particles on a lattice of orthogonal position states [7]). As such, this method is ideally suited to photonic quantum computation.

Implementation

For a quick example, consider a lattice composed of two adjacent nodes:

This graph is represented by the \(2\times 2\) adjacency matrix \(A=\begin{bmatrix}0&1\\1&0\end{bmatrix}\). Here, each node in the graph represents a qumode, so we can model the dynamics of Bosons on this structure via a 2-qumode CV circuit.

The Bose-Hubbard Hamiltonian with on-site interactions is given by

\[H = J\sum_{i}\sum_j A_{ij} \ad_i\a_j + \frac{1}{2}U\sum_i \hat{n}_i(\hat{n}_i-1)= J(\ad_1 \a_2 + \ad_2\a_1) + \frac{1}{2}U ( \hat{n}_1^2 - \hat{n}_1 + \hat{n}_2^2 - \hat{n}_2)\]

where \(J\) represents the transfer integral or hopping term of the boson between nodes, and \(U\) is the on-site interaction potential. Here, \(\ad_1 \a_2\) represents a boson transitioning from node 1 to node 2, while \(\ad_2\a_1\) represents a boson transitioning from node 2 to node 1, and \(\hat{n}_i=\ad_i\a_i\) is the number operator applied to mode \(i\). Applying the Lie-product formula, we find that

\[e^{-iHt} = \left[\exp\left({-i\frac{ J t}{k}(\ad_1 \a_2 + \ad_2\a_1)}\right)\exp\left(-i\frac{Ut}{2k}\hat{n}_1^2\right)\exp\left(-i\frac{Ut}{2k}\hat{n}_2^2\right) \exp\left(i\frac{Ut}{2k}\hat{n}_1\right)\exp\left(i\frac{Ut}{2k}\hat{n}_2\right)\right]^k+\mathcal{O}\left(t^2/k\right),\]

where \(\mathcal{O}\left(t^2/k\right)\) is the order of the error term, derived from the Lie product formula. Comparing this to the form of various CV gates, we can write this as the product of symmetric beamsplitters (BSgate), Kerr gates (Kgate), and rotation gates (Rgate):

\[e^{iHt} = \left[BS\left(\theta,\phi\right)\left(K(r)R(-r)\otimes K(r)R(-r)\right)\right]^k+\mathcal{O}\left(t^2/k\right).\]

where \(\theta=-Jt/k\), \(\phi=\pi/2\), and \(r=-Ut/2k\).

For the case \(k=2\), this can be drawn as the circuit diagram

For more complex CV decompositions, including those with interactions, see Kalajdzievski et al. [6] for more details.

Code

The Hamiltonian simulation circuit displayed above for the 2-node lattice, can be implemented using Strawberry Fields:

import numpy as np

# set the random seed
np.random.seed(42)

import strawberryfields as sf
from strawberryfields.ops import *

ham_simulation = sf.Program(2)

# set the Hamiltonian parameters
J = 1           # hopping transition
U = 1.5         # on-site interaction
k = 20          # Lie product decomposition terms
t = 1.086       # timestep
theta = -J*t/k
r = -U*t/(2*k)

with ham_simulation.context as q:
    # prepare the initial state
    Fock(2) | q[0]

    # Two node tight-binding
    # Hamiltonian simulation

    for i in range(k):
        BSgate(theta, np.pi/2) | (q[0], q[1])
        Kgate(r)  | q[0]
        Rgate(-r) | q[0]
        Kgate(r)  | q[1]
        Rgate(-r) | q[1]

where, for this example, we have set J=1, U=1.5, k=20, t=1.086, theta = -J*t/k, and r = -U*t/(2*k).

Since are initial state is non-Gaussian, we will initialize the "fock" backend, with a cutoff of 5:

eng = sf.Engine(backend="fock", backend_options={"cutoff_dim": 5})

After constructing the circuit and running the engine,

results = eng.run(ham_simulation)

the site occupation probabilities can be calculated via the fock_prob state method, which returns the probability of measuring some output Fock state:

state = results.state
print("P(|0, 2>) = ", state.fock_prob([0, 2]))
print("P(|1, 1>) = ", state.fock_prob([1, 1]))
print("P(|2, 0>) = ", state.fock_prob([2, 0]))

Out:

P(|0, 2>) =  0.5224012457200213
P(|1, 1>) =  0.23565287685672467
P(|2, 0>) =  0.24194587742325993

As Hamiltonian simulation is particle preserving, these probabilities should add up to one:

result = [state.fock_prob([0, 2]), state.fock_prob([1, 1]), state.fock_prob([2, 0])]
print(np.sum(result))

Out:

1.000000000000006

We can compare this result to the analytic matrix exponential \(e^{-iHt}\), where the matrix elements of \(H\) can be computed in the Fock basis. Considering the diagonal interaction terms,

\[\begin{split}& \braketT{0,2}{H}{0,2} = \frac{1}{2}U\braketT{0}{(\hat{n}^2-\hat{n})}{0} + \frac{1}{2}U\braketT{2}{(\hat{n}^2-\hat{n})}{2} = \frac{1}{2}U(2^2-2) = U\\[7pt] & \braketT{1,1}{H}{1,1} = 0\\[7pt] & \braketT{2,0}{H}{2,0} = U\end{split}\]

as well as the off-diagonal hopping terms,

\[\begin{split}& \braketT{1,1}{H}{0,2} = J\braketT{1,1}{\left(\ad_1\a_2 + \a_1\ad_2\right)}{0,2} = J(\sqrt{1}\sqrt{2} + \sqrt{0}\sqrt{3}) = J\sqrt{2}\\[7pt] & \braketT{1,1}{H}{2,0} = J\sqrt{2}\end{split}\]

and taking into account the Hermiticity of the system, we arrive at

\[\begin{split}H = \begin{bmatrix}U&J\sqrt{2}&0\\ J\sqrt{2} & 0 & J\sqrt{2}\\ 0 & J\sqrt{2} & U\end{bmatrix}\end{split}\]

which acts on the Fock basis \(\{\ket{0,2},\ket{1,1},\ket{2,0}\}\). Using the SciPy matrix exponential function scipy.linalg.expm:

from scipy.linalg import expm

H = J*np.sqrt(2)*np.array([[0,1,0],[1,0,1],[0,1,0]]) + U*np.diag([1,0,1])
init_state = np.array([0,0,1])

theoretical_result = np.abs(np.dot(expm(-1j*t*H), init_state))**2
print(theoretical_result)

Out:

[0.52249102 0.23516277 0.24234621]

which agrees, within the expected error margin, with our Strawberry Fields Hamiltonian simulation:

print(np.all(np.abs(theoretical_result - result) < 1e-2))

Out:

True

References

1

Dominic W. Berry, Andrew M. Childs, and Robin Kothari. Hamiltonian simulation with nearly optimal dependence on all parameters. In 2015 IEEE 56th Annual Symposium on Foundations of Computer Science. IEEE, Oct 2015. doi:10.1109/focs.2015.54.

2

Alán Aspuru-Guzik, Anthony D. Dutoi, Peter J. Love, and Martin Head-Gordon. Simulated quantum computation of molecular energies. Science, 309(5741):1704–1707, 2005. doi:10.1126/science.1113479.

3

James D. Whitfield, Jacob Biamonte, and Alán Aspuru-Guzik. Simulation of electronic structure Hamiltonians using quantum computers. Molecular Physics, 109(5):735–750, 2011. doi:10.1080/00268976.2011.552441.

4

Andrew M. Childs and Nathan Wiebe. Hamiltonian simulation using linear combinations of unitary operations. Quantum Information and Computation, 12(11-12):901–924, 2012. arXiv:1202.5822.

5

Dominic W. Berry, Graeme Ahokas, Richard Cleve, and Barry C. Sanders. Efficient quantum algorithms for simulating sparse Hamiltonians. Communications in Mathematical Physics, 270(2):359–371, 2006. doi:10.1007/s00220-006-0150-x.

6(1,2)

Timjan Kalajdzievski, Christian Weedbrook, and Patrick Rebentrost. Continuous-variable gate decomposition for the Bose-Hubbard model. 2018. arXiv:1801.06565.

7

Tomasz Sowiński, Omjyoti Dutta, Philipp Hauke, Luca Tagliacozzo, and Maciej Lewenstein. Dipolar molecules in optical lattices. Physical Review Letters, 108:115301, 2012. doi:10.1103/PhysRevLett.108.115301.

demos/run_hamiltonian_simulation

 

Download Python script

 

Download Notebook