r""" .. role:: html(raw) :format: html .. _ham_sim: 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: .. raw:: html
.. image:: /tutorials/images/graph.svg :align: center :width: 40% :target: javascript:void(0); .. raw:: html
This graph is represented by the :math:2\times 2 adjacency matrix :math: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 .. math:: 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 :math:J represents the transfer integral or hopping term of the boson between nodes, and :math:U is the on-site interaction potential. Here, :math:\ad_1 \a_2 represents a boson transitioning from node 1 to node 2, while :math:\ad_2\a_1 represents a boson transitioning from node 2 to node 1, and :math:\hat{n}_i=\ad_i\a_i is the number operator applied to mode :math:i. Applying the Lie-product formula, we find that .. math:: 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 :math:\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 :ref:symmetric beamsplitters  (:class:~strawberryfields.ops.BSgate), :ref:Kerr gates  (:class:~strawberryfields.ops.Kgate), and :ref:rotation gates  (:class:~strawberryfields.ops.Rgate): .. math:: 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 :math:\theta=-Jt/k, :math:\phi=\pi/2, and :math:r=-Ut/2k. For the case :math:k=2, this can be drawn as the circuit diagram | .. image:: /tutorials/images/hamsim.svg :align: center :width: 70% :target: javascript:void(0); | 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])) ###################################################################### # 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)) ###################################################################### # We can compare this result to the analytic matrix exponential :math:e^{-iHt}, where the matrix # elements of :math:H can be computed in the Fock basis. Considering the diagonal interaction # terms, # # .. math:: # & \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 # # as well as the off-diagonal hopping terms, # # .. math:: # & \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} # # and taking into account the Hermiticity of the system, we arrive at # # .. math:: # H = \begin{bmatrix}U&J\sqrt{2}&0\\ J\sqrt{2} & 0 & J\sqrt{2}\\ 0 & J\sqrt{2} & U\end{bmatrix} # # which acts on the Fock basis :math:\{\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) ###################################################################### # which agrees, within the expected error margin, with our Strawberry Fields Hamiltonian simulation: print(np.all(np.abs(theoretical_result - result) < 1e-2)) ###################################################################### # References # ---------- # # ..  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. # # ..  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. # # ..  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. # # ..  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. # # ..  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. # # ..  Timjan Kalajdzievski, Christian Weedbrook, and Patrick Rebentrost. Continuous-variable gate # decomposition for the Bose-Hubbard model. 2018. arXiv:1801.06565. # # ..  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.