I am trying to simulate the Hamiltonian evolution of the 1+1D $\lambda\phi^4$ scalar field theory by digitising it and encoding on a quantum computer. The process of digitising is taken from this paper: https://arxiv.org/abs/1808.10378.
The evolution operator is: $e^{-iHt} = exp(-it(\frac{1}{2}m^2\phi^2+\frac{1}{2}\Pi^2 + \frac{\lambda}{4!}\phi^4+\frac{1}{2}(\nabla\phi(x))^2))$
The term that brings trouble is: $(\nabla\phi(x))^2 = (\phi(x+i)-\phi(x))^2= \phi(x+i)^2 + \phi(x)^2 -2\underline{\phi(x+1)\phi(x)}$
In short terms, space is discredited into sites, each site is assigned a set number of qubits (3 in my case) to represent the value of the field at that site. $\phi(x)$ is the field operator that acts on a site $x$ ($\phi(x+i)$ is the field operator acting on the next adjacent site of the latice), which in the computational basis is a diagonal matrix that returns the value of the field of a targeted site: (in the 3 qubit case)
$\phi(x) = \frac{\phi_{max}}{7}\begin{pmatrix}7&0&0&0&0&0&0&0\\0&5&0&0&0&0&0&0\\0&0&3&0&0&0&0&0\\0&0&0&1&0&0&0&0\\0&0&0&0&-1&0&0&0\\0&0&0&0&0&-3&0&0\\0&0&0&0&0&0&-5&0\\0&0&0&0&0&0&0&-7\end{pmatrix}$
Since the operator targets 3 qubits, its matrix representation is a $n_s \cdot n_s = 8\cdot8$ matrix. Implementing $e^{-it\phi(x)}$ and $e^{-it\phi(x)^2}$ is straightforward, as it is easy to square and to take an exponent of a diagonal matrix. However, how would one implement a $e^{it\phi(x+i)\phi(x)}$ operator? And what would its matrix size be?
Currently, I assumed that this is not a standard matrix multiplication like in the $\phi(x)^2$ case and is instead a tensor product, which would yield a 64x64 matrix that one can apply to the 6 qubits representing two sites. However, this assumption doesn't have any basis except the fact that one would expect this size.
I am working with a quantum simulator (CIRQ), so for the time being am implementing custom gates directly as big matrices without decomposing them as pauli matrices.
