Quantum computation of phase transition in interacting scalar quantum field theory

Abstract

It has been demonstrated that the critical point of the phase transition in scalar quantum field theory with a quartic interaction in one space dimension can be approximated via a Gaussian Effective Potential (GEP). We discuss how this critical point can be estimated using quantum hardware. Performing quantum computations with various lattice sizes, we obtain evidence of a transition from a symmetric to a symmetry-broken phase using both discrete- and continuous-variable quantum computation. The ten-site case is implemented on IBM quantum hardware using the Variational Quantum Eigensolver algorithm to minimize the GEP and identify lattice level crossings. These are extrapolated via simulations to find the continuum critical point.

Appendices

Appendix A: Critical coupling

In this Appendix, we provide details of the calculation of minima of the GEP. In Sect. 3, we asserted that the effective potential, as a function of the mass parameter $\Omega$ (Eq. (28)), has a second minimum at some $\Omega _1>m$. We can use Eq. (29) to eliminate the bare mass $m_0$ in favor of the renormalized mass $m = \Omega _0$ in (28) to obtain

$$\begin{aligned} V_{\text{ G }} (\Omega )= & {} \frac{1}{2} \left[ \Omega _0^2 - \frac{\lambda }{2} I_0 (\Omega _0) \right] \left[ \frac{2(\Omega ^2 -\Omega _0^2 )}{\lambda } +I_0 (\Omega _0)- I_0 \right] \nonumber \\{} & {} + \frac{\lambda }{4!} \left[ \frac{2(\Omega ^2 -\Omega _0^2 )}{\lambda } + I_0 (\Omega _0)- I_0 \right] ^2 + I_1 - \frac{\lambda }{8} I_0^2 \end{aligned}$$

(A1)

The value of the potential at the renormalized mass is

$$\begin{aligned} V_{\text{ G }} (\Omega _0) = I_1 (\Omega _0) - \frac{\lambda }{8} I_0^2 (\Omega _0) \end{aligned}$$

(A2)

The derivative of (A1) is

$$\begin{aligned} \frac{{\text{ d }}V_{\text{ G }}}{{\text{ d }}\Omega ^2} = \frac{1}{3\lambda } \left[ 1 -\frac{\lambda }{2} \frac{dI_0}{d\Omega ^2} \right] \left[ \Omega ^2 +2\Omega _0^2 + \lambda I_0- \lambda I_0 (\Omega _0) \right] \end{aligned}$$

(A3)

It vanishes for $\Omega = \Omega _1$ such that

$$\begin{aligned} \lambda = \frac{\Omega _1^2 + 2\Omega _0^2}{I_0 (\Omega _0) -I_0 (\Omega _1)} \end{aligned}$$

(A4)

For the phase transition, we need to know the sign of $\Delta V_{\text{ G }}$ (Eq. (31)). We obtain

$$\begin{aligned}{} & {} \Delta V_{\text{ G }} =\nonumber \\{} & {} \quad \frac{[-\Omega _1^4-4\Omega _1^2\Omega _0^2 +2\Omega _0^4] I_0 (\Omega _1) +[-\Omega _1^4+2\Omega _1^2\Omega _0^2 +2\Omega _0^4] I_0 (\Omega _0) + 4(\Omega _1^2 + 2\Omega _0^2) ( I_1(\Omega _1) -I_1 (\Omega _0)) }{4(\Omega _1^2 + 2\Omega _0^2)}\nonumber \\ \end{aligned}$$

(A5)

where we used (A4). Given $m= \Omega _0$, the solution $\Omega = \Omega _1$ to $\Delta V_{\text{ G }}=0$ with $\Omega _1 > \Omega _0$ is the critical value $\Omega _c$. These values are plotted for different lattice sizes on panel (a) in Fig. 3. Substituting in (A4), we obtain the critical coupling $\lambda _c$ (Eq. (32)). Their values for different lattice sizes are shown on panel (b) in Fig. 3.

Appendix B: VQE

In this Appendix, we provide details of the terms in the Hamiltonian that are relevant to our VQE algorithm and the CV quantum circuits.

To take advantage of the form of the Ansatz state, we express the Hamiltonian in terms of the quadrature operators q and p and then compute expectation values of the various terms in the Hamiltonian using a quantum computer. Notice that the Gaussian Ansatz in Fig. 4 is even under the reflection operators for nonzero modes, $e^{\pi i \left[ N (k)+ N (L-k)\right] }$ ($1\le k < \frac{L}{2}$) and $e^{\pi i N (\frac{L}{2})}$. Therefore, the only Hamiltonian terms contributing are those with an even number of quadrature operators acting on the $k=\frac{L}{2}$ mode or an even total number of quadrature operators acting on a $(k,L-k)$ mode pair, since each quadrature operator is odd under reflection. Since the Hamiltonian only contains an even number of quadrature operators per term, this implies that only an even number can act on the zero mode as well. We can also eliminate terms with an odd number of quadratures p(k) acting on a single mode. This follows from the fact that such an operator is purely imaginary while the Ansatz unitaries are real. These observations help us reduce the number of terms to be computed.

It is convenient to express the fields (Eq. (21)) as

$$\begin{aligned} \phi (x)= & {} \frac{1}{\sqrt{L}}\left[ \frac{q(0)}{\sqrt{\omega (0)}}+\frac{q (\frac{L}{2})}{\sqrt{\omega (\frac{L}{2})}}(-)^x+\sum _{1\le k< \frac{L}{2}} \sqrt{\frac{2}{\omega (k)}}\left( q_+(k)\cos \frac{2\pi kx}{L}\right. \right. \nonumber \\{} & {} \left. \left. -p_-(k)\sin \frac{2\pi kx}{L}\right) \right] \nonumber \\ \pi (x)= & {} \frac{1}{\sqrt{L}}\left[ \sqrt{\omega (0)} p(0)+\sqrt{\omega \left( \frac{L}{2}\right) } p\left( \frac{L}{2}\right) (-)^x+\sum _{1 \le k < \frac{L}{2}} \sqrt{2\omega (k)} \right. \nonumber \\{} & {} \left. \left( p_+(k)\cos \frac{2\pi kx}{L}+q_-(k)\sin \frac{2\pi kx}{L}\right) \right] \end{aligned}$$

(B1)

in terms of the linear combinations of quadratures,

$$\begin{aligned} q_\pm (k)=\frac{1}{\sqrt{2}}\left( q(k)\pm q(L-k)\right) \ ,\ \ p_\pm (k)=\frac{1}{\sqrt{2}}\left( p(k)\pm p(L-k)\right) \ ,\ \ 1 \le k < \frac{L}{2} \end{aligned}$$

(B2)

which can be implemented with beam splitters.

Fig. 15

Alternative quantum circuit to the one in Fig. 5

The Hamiltonian can be written as

$$\begin{aligned} H= & {} \frac{\omega (0)}{2} \left[ q^2(0)+p^2(0)\right] + \frac{\omega (\frac{L}{2})}{2} \left[ q^2(\frac{L}{2}) + p^2( \frac{L}{2}) \right] \nonumber \\{} & {} + \sum _{1\le k< \frac{L}{2}} \frac{\omega (k)}{2} \left[ q_+^2(k)+q_-^2(k)+p_+^2(k)+p_-^2(k)\right] ] \nonumber \\{} & {} + \frac{\Omega ^2 - m^2}{2}\left[ \frac{q^2(0)}{\omega (0)}+\frac{q^2(L/2)}{\omega (\frac{L}{2})} + \sum _{1 \le k< \frac{L}{2}} \frac{1}{\omega (k)} \left( q_+^2(k)+p_-^2(k)\right) \right] \nonumber \\{} & {} + \frac{\lambda }{24L}\left[ \frac{q^4(0)}{\omega ^2 (0)}+\frac{q^4(\frac{L}{2})}{\omega ^2 (\frac{L}{2})}+6\frac{q^2(0)q^2( \frac{L}{2})}{\omega (0) \omega (\frac{L}{2})} + 6\left( \frac{q^2(0)}{\omega (0)} + \frac{q^2(\frac{L}{2})}{\omega (\frac{L}{2})} \right) \right. \nonumber \\{} & {} \sum _{1\le k< \frac{L}{2}} \frac{1}{\omega (k)} \left( q_+^2(k)+p_-^2(k)\right) \nonumber \\{} & {} + \frac{1}{2}\sum _{1\le k< \frac{L}{2}}\frac{1}{\omega ^2(k)} \left[ 3\left( q_+^2(k)+p_-^2(k)\right) ^2\left( 1-\delta _{k,\frac{L}{4}} \right) +4\left( q_+^4(k)+p_-^4(k)\right) \delta _{k,\frac{L}{4}} \right] \nonumber \\{} & {} + 3\sum _{1\le k< k' < \frac{L}{2}} \frac{1}{\omega (k) \omega (k')} \left[ 2\left( q_+^2(k)+p_-^2(k)\right) \left( q_+^2(k')+p_-^2(k')\right) \right. \nonumber \\{} & {} \left. \left. +\left( q_+^2(k)-p_-^2(k)\right) \left( q_+^2(k')-p_-^2(k')\right) \delta _{k+k',\frac{L}{2}} \right] \right] \end{aligned}$$

(B3)

In the VQE circuit shown in Fig. 5, we used a two-mode squeezer to generate all of the squeezed states, including the $k=0, \frac{L}{2}$ modes. We also applied 50/50 beamsplitters to all pairs of qumodes in order to implement the transformation (B2). For simulations, the circuit in Fig. 5 results in significant sampling error from computing derivatives with respect to $\Gamma$, even when using the parameter shift rule. This sampling error can be mitigated on actual CV hardware by increasing the magnitude of the shifts on the $\Gamma$ parameter (up to the limit the platform can handle). However, in a finite Hilbert space necessitated by limitations of simulations, this greatly increases truncation error. We therefore implemented the equivalent circuit shown in Fig. 15, instead.

We calculated expectation values of quadratures that contribute to the expectation value of the Hamiltonian using photon number measurements. The expectation value $\langle q^2(k)\rangle$ was computed in (39). More generally, we have

$$\begin{aligned} \langle q^{2n}(k)\rangle = \frac{2^n}{(2n)!}\frac{d^{2n}}{d\Gamma ^{2n}}\left\langle {\Phi (\Gamma )}\right| N_\text {anc}^{n}\left| {\Phi (\Gamma )}\right\rangle \end{aligned}$$

(B4)

where $\left| {\Phi (\Gamma )}\right\rangle$ is given in Eq. (38). This follows from the fact that the CX gate $e^{i\Gamma p_\text {anc}\otimes q(k)}$ transforms $q_\text {anc} \rightarrow q_\text {anc}+\Gamma q(k)$, $p_\text {anc} \rightarrow p_\text {anc}$.

We re-expressed these derivatives as linear combinations of expectation values using CX gates with different parameters $\Gamma$; see, Eq. (40) for $\langle q^2(k)\rangle$. For the quartic term, $\langle q^4(k)\rangle$, we used a five-term parameter shift:

$$\begin{aligned} \langle q^4(k)\rangle= & {} \frac{1}{6s^4}\Big [ \left\langle {\Phi (2s)}\right| N_\text {anc} \left| {\Phi (2s)}\right\rangle + \left\langle {\Phi (-2s)}\right| N_\text {anc} \left| {\Phi (-2s)}\right\rangle - 4\left\langle {\Phi (s)}\right| N_\text {anc} \left| {\Phi (s)}\right\rangle \nonumber \\- & {} 4\left\langle {\Phi (-s)}\right| N_\text {anc} \left| {\Phi (-s)}\right\rangle +6 \left\langle {\Phi (0)}\right| N_\text {anc} \left| {\Phi (0)}\right\rangle \Big ] \end{aligned}$$

(B5)

For modes with $1\le k < \frac{L}{2}$, we considered pairs $(k,L-k)$ as one another’s ancilla by making the replacements $N_\text {anc}\rightarrow N(k)+N(L-k)$, and $e^{-i\Gamma p_\text {anc}\otimes q(k)}\rightarrow e^{-i\Gamma p(k)\otimes q(L-k)}$.

Appendix C: Derivatives for physical parameters

Here, we discuss further the implementation of derivatives of expectation values with respect to physical parameters which can be computed exactly in quantum circuits. As we discussed, the introduction of parameter shift rules allowed us to compute derivatives with respect to the squeezing parameters r(k) (Eq. (34)). This was enabled by the Gaussian form of our Ansatz. It is of great interest to extend results of exact computation of derivatives to cases where parameter shift rules are challenging to implement, such as in the presence of non-Gaussian gates.

To this end, one can leverage the tools of bosonic qiskit [30] which allow coupling between qubits and qumodes. Such couplings exist in hybrid quantum systems such as circuit QED [45]. In particular, one can take advantage of the Selective Number-dependent Arbitrary Phase (SNAP) gate,

$$\begin{aligned} \text {SNAP}_n\left( \theta \right) =e^{-i\theta Z\otimes \left| {n}\right\rangle \left\langle {n}\right| } \end{aligned}$$

(C1)

which involves a phase that is conditioned on the photon number state $\left| {n}\right\rangle$ of the physical qumode.

Fig. 16

Quantum circuit employing two ancilla qubits for the calculation of derivatives with respect to squeezing parameters for the modes with $k=0,\frac{L}{2}$. It features single-mode squeezers, an unspecified qumode unitary U, Hadamard gates on the qubits, and SNAP gates entangling qumodes and qubits. At the end, qubits are measured in the computational basis whereas the (unspecified) Hermitian operator A is measured on the qumode

To compute derivatives with respect to a squeezing parameter r(k), notice that for the expectation value of a Hermitian operator A in a state constructed using a squeezer $S(r) = e^{\frac{i}{2} (a^{\dagger 2} - a^2)}$ and an unspecified, possibly multi-modal, unitary U, taking a derivative with respect to r, we obtain

$$\begin{aligned} \frac{{\text{ d }}}{{\text{ d }}r} \left\langle {0}\right| S^\dagger (r)U^\dagger AUS(r) \left| {0}\right\rangle = {\sqrt{2}} \Re { \left\langle {2}\right| S^\dagger (r)U^\dagger AUS(r)\left| {0}\right\rangle } \end{aligned}$$

(C2)

involving number states $\left| {0}\right\rangle$ and $\left| {2}\right\rangle$. For a single qumode, this is implemented with the quantum circuit in Fig. 16. This circuit allows us to compute derivatives using an ancilla. We first apply a squeezer of arbitrary squeezing t that generates all the even number states, and then, we use SNAP gates to change the states of the qubits. The qubit states allow us to label the number states for the physical qumode. By projecting the ancilla qubit labeled 1 onto the state $\left| {1}\right\rangle$, we pick out the states of the qumode that we are interested in: $\left| {0}\right\rangle ,\left| {2}\right\rangle$. Finally, we measure $Z_0\otimes A$. To obtain (C2), we need divide by the appropriate factor which contains the matrix elements $\left\langle {0}\right| S(t)\left| {0}\right\rangle$ and $\left\langle {2}\right| S(t)\left| {0}\right\rangle$. This method is an extension of the Hadamard test [46] to hybrid systems. We used the circuit of Fig. 16 for the $k = 0, \frac{L}{2}$ modes.

Fig. 17

Quantum circuit employing two ancilla qubits for the calculation of derivatives with respect to squeezing parameters for pairs of modes $(k,L-k)$ with $1\le k <\frac{L}{2}$. It features single-mode squeezers, an unspecified qumode unitary U, Hadamard gates on the qubits, and controlled beam splitter, rotation, and SNAP gates entangling qumodes and qubits. At the end, qubits are measured in the computational basis whereas the photon number is measured on the qumodes

For modes with $1\le k < \frac{L}{2}$, we treat them as pairs $(k,L-k)$. The quantum circuit for derivatives with respect to squeezing parameters in this case is slightly modified. We need include additional circuit elements that incorporate some of bosonic qiskit’s controlled Gaussian operations, namely a controlled beam splitter and rotation gate. The resulting circuit is shown in Fig. 17.

Turning to derivatives with respect to $\phi _C$, they can be computed by removing the displacement D from the circuit, taking $H\rightarrow D^\dagger H D$, and expressing it as a polynomial in $\phi _C$. In general, derivatives with respect to $\phi _C$ are straightforward. The only terms affected by the displacement are those that contain the zero-mode quadrature. Thus, we need only compute $\langle D^\dagger q^2(0) D\rangle$ and $\langle D^\dagger q^4(0) D \rangle$ in a state that lies in the even parity subspace (squeezed vacuum). Dropping terms that are odd in the quadrature q(0), we easily obtain

$$\begin{aligned} \frac{{\text{ d }}}{{\text{ d }}c}\langle D^\dagger (c)q^2(0)D(c) \rangle =2c \ . \end{aligned}$$

(C3)

For the quartic term, we obtain

$$\begin{aligned} \frac{d}{dc} \langle D^\dagger (c)q^4(0)D(c)\rangle = 12\langle D^\dagger (c)q^2(0)D(c) \rangle -8c^3\ , \end{aligned}$$

(C4)

expressing the derivative in terms of the expectation value of $q^2(0)$ which has already been computed (Eq. (40)).