# 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. 

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 .

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 . Efficient implementations of Hamiltonian simulation also exist in the CV formulation , with specific application to Bose-Hubbard Hamiltonians (describing a system of interacting bosonic particles on a lattice of orthogonal position states ). 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.  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

# Two node tight-binding
# Hamiltonian simulation

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


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