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. [21]

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 [22][23].

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 [24][25]. Efficient implementations of Hamiltonian simulation also exist in the CV formulation [26], with specific application to Bose-Hubbard Hamiltonians (describing a system of interacting bosonic particles on a lattice of orthogonal position states [27]). As such, this method is ideally suited to photonic quantum computation.

CV 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. [26] for more details.

Blackbird code

The Hamiltonian simulation circuit displayed above for the 2-node lattice, can be implemented using the Blackbird quantum circuit language:

# prepare the initial state
Fock(2) | q[0]

# Two node tight-binding
# Hamiltonian simulation

for i in range(k):
    BSgate(theta, 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). After constructing the circuit and running the engine,

results = eng.run(ham_simulation)

the site occupation probabilities can be calculated via

>>> state = results.state
>>> state.fock_prob([2,0])
>>> state.fock_prob([1,1])
>>> state.fock_prob([0,2])

As Hamiltonian simulation is particle preserving, these probabilities should add up to one; indeed, summing them results in a value \(\sim1.0000000000000053\).

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])
>>> np.abs(np.dot(expm(-1j*t*H), init_state))**2
[ 0.52249102,  0.23516277,  0.24234621]

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


A fully functional Strawberry Fields simulation containing the above Blackbird code is included at examples/hamiltonian_simulation.py.