Gaussian Hamiltonians

Section author: Josh Izaac <josh@xanadu.ai>

Note

This section assumes a familiarity with the operators and conventions used in quantum optics and continuous-variable (CV) quantum computation. For an introduction to these concepts, see the Strawberry Fields documentation.

Quadrature operators

A Gaussian Hamiltonian is a Hamiltonian that can be written as a quadratic polynomial of the quadrature operators \(\q_i\) (canonical position, also written as \(\x_i\)) and \(\p_j\) (canonical momentum) [2]:

\[\hat{H} = \frac{1}{2}\r^T A\r + \r^T \mathbf{d}\]

where:

  • \(\r=(\q_1,\dots,\q_{N},\p_1,\dots,\p_N)\) is a vector of the canonical position and momentum operators acting on qumode \(i\), in \(\q\p\)-ordering,
  • \(A\in\mathbb{R}^{2N\times 2N}\) is a symmetric matrix containing the quadratic coefficients of terms of the form \(\hat{x}_i\hat{y}_j\),
  • \(\mathbf{d}\in\mathbb{R}^{2N}\) is a vector containing the linear coefficients corresponding to terms of the form \(\hat{x}_i\).

For example, the following are all quadratic Hamiltonians:

  • \(\hat{H} = \q_0^2 -q_0\)
  • \(\hat{H} = \q_0 \q_1 - \p_0\p_1\)

while \(\hat{H}=\q_0^2\p_1\), with a third-order term in the quadrature operators, is not.

Ladder operators

In terms of the bosonic annihilation and creation operators

\[\a = \sqrt{\frac{1}{2 \hbar}} (\q +i\p), ~~~~ \ad = \sqrt{\frac{1}{2 \hbar}} (\q -i\p),\]

this corresponds to a Hamiltonian of the form [2]

\[\hat{H} = i\left(\av^\dagger \mathbf{w} - \mathbf{w}^\dagger\av +\av^\dagger F \av - \av F^\dagger \av^\dagger +\av^\dagger G {\av^\dagger}^T -{\av}^T G^\dagger \av\right)\]

where \(\av = (\a_1, \a_2,\dots,\a_N)\), \(\mathbf{w}\in\mathbb{C}^N\) and \(F,G\in\mathbb{C}^{N\times N}\). This corresponds to a linear unitary Bogoliubov transform of the annihilation and creation operators:

\[e^{-i\hat{H}t/\hbar}\a e^{i\hat{H}t/\hbar} = A\a + B\ad + \mathbf{w}\]

where \(AB^T=BA^T\) and \(AA^\dagger = BB^\dagger+\I\).

Time propagation

Once the matrix of quadratic coefficients \(A\) and vector of linear coefficients \(\mathbf{d}\) have been found for the Gaussian Hamiltonian, it can be shown that the Hamiltonian can always be written in the following form, up to a local phase factor [3]:

\[\hat{H} = \frac{1}{2}(\r-\mathbf{d}')^T A(\r-\mathbf{d}')\]

where \(\mathbf{d}'=-A^{-1}\mathbf{d}\). Let’s consider the case of zero linear coefficients and non-zero linear coefficients separately.

Zero linear coefficients

In the Heisenberg picture, with \(\hat{H}=\frac{1}{2}\r^T A\r\), the time-evolution of the quadrature operators must satisfy the Heisenberg equations of motion:

\[\frac{d}{dt}\r_j = \frac{i}{\hbar}[\hat{H},\r_j] ~~\Leftrightarrow ~~ \frac{d}{dt}\r = \Omega A \r ,\]

where

\[\begin{split}\Omega = \begin{bmatrix} 0 & \I_N \\-\I_N & 0 \\\end{bmatrix}\end{split}\]

is the symplectic matrix, coming from the definition of the canonical commutation relation \([\r_i,\r_j]=i\hbar \Omega_{ij}\).

Solving this differential equation, we find that the symplectic Gaussian transformation describing the time-evolution of the Hamiltonian \(\hat{H}\) acting on the quadrature operators, for time \(t\), is given by:

\[S = \exp{\left(\Omega A t\right)}\]

Non-zero linear coefficients

If we have non-zero linear coefficients, we need to take this into account by performing the required displacement operation. Consider the following Gaussian Hamiltonian:

\[\hat{H}(\mathbf{d}') = \frac{1}{2}(\r-\mathbf{d}')^TA(\r-\mathbf{d}')\]

where \(\mathbf{d}'=-A^{-1}\mathbf{d}\), as before. This corresponds to the action of a Weyl operator or displacement of the form \(\hat{D}(\mathbf{s})\) with \(\mathbf{s}=-\mathbf{d}'/\sqrt{2\hbar}\):

\[\hat{H}(\mathbf{d}') = \frac{1}{2}\hat{D}(\mathbf{s})\r^T A\r \hat{D}(\mathbf{s})^\dagger = \hat{D}(\mathbf{s})\hat{H}(0)\hat{D}(\mathbf{s})^\dagger.\]

Calculating the time-evolution operator,

\[\hat{U}(t) = e^{-i\hat{H}(d) t/\hbar} = e^{-i\hat{D}(\mathbf{s})\hat{H}(0)\hat{D}(\mathbf{s})^\dagger t} = \hat{D}(\mathbf{s})e^{-i\hat{H}(0) t}\hat{D}(\mathbf{s})^\dagger.\]

In order to write this as a symplectic matrix transformation, we need to move all displacement operators to the left. To do this, we can post-multiply by \(\I=e^{i\hat{H}(0)t}e^{-i\hat{H}(0)t}\):

\[\hat{U}(t) = \hat{D}(\mathbf{s})\left[e^{-i\hat{H}(0) t}\hat{D}(\mathbf{s})^\dagger e^{i\hat{H}(0)t}\right]e^{-i\hat{H}(0)t}\]

Finally, we can rewrite this as a symplectic transformation, by making the substitution \(e^{-i\hat{H}(0)t}\rightarrow e^{\Omega A t}\) and by noting that the bracketed term is simply a displacement by \(-\mathbf{s}\), evolved under \(\hat{H}(0)\) for time \(t\):

\[S = \hat{D}(\mathbf{s} -\left(e^{\Omega A t}\right)^T \mathbf{s}) e^{\Omega A \hbar t}\]

Definition

For a quadratic Hamiltonian of the form \(\hat{H} = \frac{1}{2}\r^T A\r + \r^T \mathbf{d}\), the symplectic transformation \(S\in\mathbb{R}^{2N\times 2N}\) characterizing the time-evolution unitary operator \(\hat{U}(t) = e^{-i\hat{H}t/\hbar}\) is given by

\[S = \hat{D}(\mathbf{s} -{e^{-\Omega A t}} \mathbf{s}) e^{\Omega A t}\]

where \(\Omega\) is the symplectic matrix, \(\hat{D}\) the displacement operation, and \(\mathbf{s} = -A^{-1}\mathbf{d}/\sqrt{2\hbar}\).

Tip

Implemented in SFOpenBoson as a quantum operation by sfopenboson.ops.GaussianPropagation

Warning

In the case where the quadratic coefficient matrix \(A\) is singular, for example the Hamltonian \(\hat{H}=\frac{1}{2}p_0^2+q_0\), in order to calculate \(A^{-1}\) to determine the resulting displacement, the matrix \(A\) is perturbed by \(\epsilon\ll 1\):

\[\mathbf{s} = -\frac{(A+\epsilon)^{-1}\mathbf{d}}{\sqrt{2\hbar}}, ~~~S = \hat{D}\left(\left(\I -{e^{-\Omega (A+\epsilon) t}}\right) \frac{\mathbf{s}}{\epsilon}\right) e^{\Omega A t}\]