# How to implement the cross term (multi-qubit) in the square of the finite difference operator?

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.

After some looking into this, the two operator product is a tensor product instead of normal matrix multiplication. For a 2 qubit system, operator $$\phi(1)$$ is represented as a matrix acting on the full system as $$\phi(1) \otimes I(2)$$.