# Simulating scalar field theories on quantum computers with limited resources

Andy C. Y. Li, Alexandru Macridin, Stephen Mrenna, Panagiotis Spentzouris
Fermi National Accelerator Laboratory, Batavia, IL 60510

###### Abstract

We present a quantum algorithm for implementing $\phi^{4}$ lattice scalar
field theory on qubit computers. The field is represented in the discretized
field amplitude basis. The number of qubits and elementary gates required by
the implementation of the evolution operator is proportional to the lattice
size. The algorithm allows efficient $\phi^{4}$ state preparation for a large
range of input parameters in both the normal and broken symmetry phases. The
states are prepared using a combination of variational and adiabatic evolution
methods. First, the ground state of a local Hamiltonian, which includes the
$\phi^{4}$ self-interaction, is prepared using short variational circuits.
Next, this state is evolved by switching on the coupling between the lattice
sites adiabatically. The parameters defining the local Hamiltonian are
adjustable and constitute the input of our algorithm. We present a method to
optimize these parameters in order to reduce the adiabatic time required for
state preparation. For preparing broken symmetry states, the adiabatic
evolution problems caused by crossing the phase transition critical line and
by the degeneracy of the broken symmetry ground state can be addressed using
an auxiliary external field which gradually turns off during the adiabatic
process. We show that the time dependence of the external field during the
adiabatic evolution is important for addressing the broken symmetry ground
state degeneracy. The adiabatic time dependence on the inverse error tolerance
can be reduced from quadratic to linear by using a field strength that
decreases exponentially in time relative to one that decreases linearly.

††preprint: FERMILAB-PUB-22-757-QIS

##  I Introduction

Simulating highly entangled quantum systems is among the first applications of
quantum computers expected to show a practical advantage over classical
computers [1, 2]. The development and application of new quantum processors
[3, 4, 5] may allow for revolutionary calculations in quantum chemistry [6, 7,
8, 9, 10, 11], condensed-matter physics [12, 13, 14, 15, 16, 17], nuclear
physics [18, 19] and high-energy physics [20, 21, 22, 23, 24]. The simulation
of relativistic quantum field theory on quantum hardware [20, 21, 24, 25, 26,
27] has been an active research topic in recent years. In practice, the
mapping and preparation of continuous fields on near-future quantum hardware
of limited size and with limited control fidelity provides a number of
challenges. In this paper, we address the simulation of the $\phi^{4}$ scalar
field on digital quantum computers.

The $\phi^{4}$ scalar field model [28, 29] is a simplified model of the Higgs
field of the standard model of particle physics and has been intensively
studied over the years. Despite its apparent simplicity, it has rich physics.
For example, in $(1+1)$ and $(2+1)$ space-time dimensions, it exhibits a phase
transition to a broken-symmetry phase characterized by a finite vacuum
expectation value $\langle\phi\rangle$ [30, 31]. Perturbative methods based on
a diagrammatic expansion are only valid in the weak interaction regime. As a
result, the strong interacting regime has been studied numerically.
Calculations of the critical coupling and exponent for $(1+1)$ dimensions have
been performed using DMRG [32], tensor network methods [33, 34], Monte Carlo
methods [35, 36] and diagonaliztion methods [37, 38]. However, since the
Hilbert space of the $\phi^{4}$ model is exponentially large, the field
degrees of freedom must be truncated, making the extrapolation of the
numerical results to the continuous limit challenging and not always reliable.

Quantum simulations can overcome the size problems related to the Hilbert
space and, unlike most classical Monte Carlo methods, can calculate the real-
time correlations and non-equilibrium dynamics of the system. The bosonic
fields can be represented efficiently on qubits in a discretized field
amplitude basis [13, 12, 39]. However the preparation of field eigenstates on
qubits is not straightforward. For example, the method proposed in Refs [20,
26, 27] prepares an initial non-interacting multi-variate Gaussian state and
uses adiabatic continuation to reach the desired interacting state. However,
the construction of a multi-variate Gaussian wavefunction using the Kitaev-
Webb method [40] requires a very large number of qubits and is not feasible on
near-term quantum hardware. Moreover, the preparation of broken-symmetry
states using adiabatic continuation of non-interacting states is challenging
since the adiabatic path has to cross a critical region with a vanishing
excitation gap. Furthermore, the ground state of broken symmetry states is
degenerate, causing further complications.

Here, we present a quantum algorithm for lattice $\phi^{4}$ field evolution on
qubits and a method for initial state preparation suitable for near-term
quantum computers. A relatively small number of qubits per lattice site,
$n_{q}\approx 6\sim 8$, is sufficient to address even strong-coupling regimes.
The number of qubits and the number of gates scale proportionally to the
system size $N$. The most expensive part of the algorithm arises from the
implementation of the $\phi^{4}$ interaction, which requires
$\mathcal{O}(N\,n_{q}^{4})$ two-qubit gates, while the implementation of the
other terms in the Hamiltonian requires $\mathcal{O}(N\,n_{q}^{2})$ two-qubit
gates. The field state preparation combines variational and adiabatic
evolution approaches. The Hamiltonian is split into two parts, a local one
that sums contributions from each individual site and a non-local one
containing coupling between sites. The adiabatic process starts from the
ground state of the local Hamiltonian. Then the coupling between sites is turn
on adiabatically. Unlike previous approaches in the literature [20], our
method introduces self-interactions from the start. The ground state of the
local Hamiltonian is prepared accurately using short variational circuits.
Instead of preparing the full lattice states using variational ansatzes [41],
which would be difficult to scale up due to Barren plateaus [42], our
variational circuits prepare local states. The calculation of these circuits’
parameters is independent of the system size and can be done easily on
classical computers using various optimization methods. The input parameters
of the local Hamiltonian can be adjusted to minimize the time of the adiabatic
process. We find a direct correlation between the adiabatic time and the local
overlap of the initial wavefunction and the target wavefunction. We propose a
strategy to determine the optimal parameters of the local Hamiltonian by
maximizing this local overlap.

We also address the problems associated with the preparation of the broken
symmetry states, namely the crossing of the critical phase transition region
characterized by a vanishing excitation gap and the double degeneracy of the
broken symmetry state. Both of these problems can be mitigated by coupling the
scalar field to an external field. We propose a two step adiabatic process for
preparing broken symmetry states. The first adiabatic process starts from a
local state in the presence of a significant external field. Then,
adiabatically, the inter-site coupling term is turned on and the external
field is decreased. Due to the presence of the external field, the excitation
gap is robust during this process. The second adiabatic process starts from
the terminus of the first one. During this step the external field is
decreased to vanishing values. We find an reduction of the required adiabatic
time from $\mathcal{O}(\varepsilon^{-2})$ to
$\mathcal{O}(\varepsilon^{-1}\ln(\varepsilon^{-1}))$ with $\varepsilon$ being
the error bound when the external field decreases exponentially in time
compared to the case of linear decrease in time.

This paper is organized as follows. We review the $\phi^{4}$ model and its
lattice discretization in Section II. We then discuss the qubit encoding and
circuits to simulate the scalar field evolution on quantum computers in
Section III. In Section IV, we introduce our state preparation protocol
consisting of the variational local-state preparation (Section IV.1), and
adiabatic evolution for the normal phase (Section IV.2.1) and for the broken-
symmetry phase (Section IV.2.2) supported by numerical simulation of the
lattice $\phi^{4}$ model with up to four sites. Our summary and conclusions
are provided in Section V.

##  II The $\phi^{4}$ model

Figure 1:  The gap, $E_{10}=E_{1}-E_{0}$, and the energy difference between
the second and the first excited states, $E_{21}=E_{2}-E_{1}$, for $1$-site,
$2$-site and $4$-site $\phi^{4}$ lattices for negative values of $m_{0}^{2}$.
The system is nearly double degenerate when
$m_{0}^{2}/\lambda_{0}^{2/3}\ll-1$. The gap decreases exponentially with
increasing the magnitude of $|m_{0}^{2}|$. On the other hand, $E_{21}$
increases slightly with increasing the magnitude of $|m_{0}^{2}|$. With
increasing the number of sites, the gap $E_{10}$ decreases while $E_{21}$
increases.

The Hamiltonian density of the $\phi^{4}$ model can be written as
($\hbar=c=1$)

$\displaystyle\mathbb{H}$
$\displaystyle=\frac{1}{2}\pi^{2}+\frac{1}{2}m_{b}^{2}\phi^{2}+\frac{1}{2}\left(\nabla\phi\right)^{2}+\frac{\lambda_{b}}{4!}\phi^{4}+f_{b}\phi,$
(1)

where $m_{b}$ and $\lambda_{b}$ are the unrenormalized (bare) mass and
interaction strength, respectively. In order to investigate the broken
symmetry phase, it is convenient to consider a coupling term between the
scalar field and a external field $f_{b}$. The field operator $\phi$ and the
conjugate-field operator $\pi$ obey the commutation relation,

$\displaystyle\left[\pi(x),\phi(y)\right]=i\delta(x-y).$ (2)

For quantum simulation, we consider the lattice version of the $\phi^{4}$
model in $d+1$ spacetime dimensions given by

$\displaystyle H_{\mathrm{lat}}=a^{d}\sum_{j}\Bigg{[}$
$\displaystyle\frac{1}{2}\pi_{j}^{2}+\frac{1}{2}m_{b}^{2}\phi_{j}^{2}+\frac{1}{2a^{2}}\sum_{e=1}^{d}\left(\phi_{j+e}-\phi_{j}\right)^{2}$
$\displaystyle+\frac{\lambda_{b}}{4!}\phi_{j}^{4}+f_{b}\phi_{j}\Bigg{]}$ (3)

where $a$ is the lattice constant and $j$ labels lattice sites. The label
$j+e$ represents the next-nearest neighbor site of the site $j$ in the
direction $e$. Note that only the space dimension is discretized and on a
lattice, while the time is a continuous variable. This differs from most Monte
Carlo-based simulations, where both space and imaginary time dimensions are
discretized. The lattice field operators in Section II obey the commutation
relations

$\left[\phi_{i},\pi_{j}\right]=ia^{-d}\delta_{i,j}\text{\ \ and \ \
}\left[\phi_{i},\phi_{j}\right]=\left[\pi_{i},\pi_{j}\right]=0.$ (4)

The continuous limit is achieved by taking $a\rightarrow 0$. It is convenient
to introduce dimensionless field variables

$\displaystyle\mkern 1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern
1.5mu_{j}$ $\displaystyle=a^{\frac{d-1}{2}}\phi_{j}\text{\ \ and \ \ }\mkern
1.5mu\overline{\mkern-1.5mu\Pi\mkern-1.5mu}\mkern
1.5mu_{j}=a^{\frac{d+1}{2}}\pi_{j}$ (5)

which obey the canonical commutation relations

$\displaystyle\left[\mkern 1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern
1.5mu_{i},\mkern 1.5mu\overline{\mkern-1.5mu\Pi\mkern-1.5mu}\mkern
1.5mu_{j}\right]$ $\displaystyle=i\delta_{i,j}\text{\ \ and \ \ }\left[\mkern
1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern 1.5mu_{i},\mkern
1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern
1.5mu_{j}\right]=\left[\mkern
1.5mu\overline{\mkern-1.5mu\Pi\mkern-1.5mu}\mkern 1.5mu_{i},\mkern
1.5mu\overline{\mkern-1.5mu\Pi\mkern-1.5mu}\mkern 1.5mu_{j}\right]=0.$ (6)

Using these dimensionless operators, the Hamiltonian is given by

$\displaystyle\mkern 1.5mu\overline{\mkern-1.5muH\mkern-1.5mu}\mkern
1.5mu=\sum_{j}\Bigg{[}$ $\displaystyle\frac{1}{2}\mkern
1.5mu\overline{\mkern-1.5mu\Pi\mkern-1.5mu}\mkern
1.5mu_{j}^{2}+\frac{1}{2}m_{0}^{2}\mkern
1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern
1.5mu_{j}^{2}+\frac{1}{2}\sum_{e=1}^{d}\left(\mkern
1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern 1.5mu_{j+e}-\mkern
1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern 1.5mu_{j}\right)^{2}$
$\displaystyle+\frac{\lambda_{0}}{4!}\mkern
1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern 1.5mu_{j}^{4}+f_{0}\mkern
1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern 1.5mu_{j}\Bigg{]}.$ (7)

where $\mkern 1.5mu\overline{\mkern-1.5muH\mkern-1.5mu}\mkern 1.5mu\equiv
H_{\mathrm{lat}}a$, $m_{0}^{2}\equiv m_{b}^{2}a^{2}$,
$\lambda_{0}\equiv\lambda_{b}a^{3-d}$ and $f_{0}\equiv
f_{b}a^{\left(3+d\right)/2}$ are dimensionless. This Hamiltonian (with
$f_{0}=0$) was previously used in numerical simulations of the scalar field
model [33, 25]. It represents a set of coupled harmonic oscillators with an
anharmonic interaction.

The correlation length is a measurable parameter which determines how the
correlation between the field values at two separate points decays with the
distance between these points. To extrapolate the lattice results to a
meaningful continuous limit with a finite correlation length, it is necessary
to simulate large lattices for Hamiltonian parameters corresponding to large
lattice correlation lengths (measured in units of $a$). The physics of the
$\phi^{4}$ theory is extracted from the lattice results by taking
$a/\xi\rightarrow 0$ (continuous limit) and $L/\xi\rightarrow\infty$
(macroscopic limit), where $L$ is the lattice size and $\xi$ is the lattice
correlation length. Equivalently, this implies simulations with
$m_{p}a\rightarrow 0$ and $m_{p}a\gg 1/N$ where $N$ is the number of lattice
sites per dimension and $m_{p}\propto 1/\xi$ is the physical mass.

While the lattice physical parameters are needed for the extrapolation from
the lattice to the continuous theory, the bare parameters define the input of
the simulations. To be useful, a quantum algorithm should be able to prepare
efficiently quantum states for a large range of bare parameters, including
both negative and positive input parameter $m_{0}^{2}$. In principle, the
relation between the physical and the bare parameters can be established from
simulations, since the lattice physical parameters can be extracted from the
correlation functions. For the extrapolation to the continuous limit the bare
lattice parameters need to be chosen dependent on the lattice constant $a$.
This dependence is significant since the renormalization theory shows that, in
order to extrapolate to a continuous theory with finite physical observables,
the bare parameters diverge with $a\rightarrow 0$ in many cases. For example,
in $(1+1)$ and $(2+1)$ dimensions, the bare squared mass $m_{0}^{2}$ becomes
negative and proportional to $\ln(m_{p}a)$ and $-1/\left(m_{p}a\right)$,
respectively [43], for small $a$.

The $\phi^{4}$ model (with $f_{0}=0$) has a discrete $Z_{2}$ symmetry from the
transformation $\phi\rightarrow-\phi$. In $(3+1)$ dimensions, the theory is
”believed” to be trivial (i.e. the theory is actually non-interacting in the
continuum limit), although no rigorous proof exists [44, 45, 46, 47, 48]. For
higher dimensions, the triviality can be rigorously proven [28]. In $(1+1)$
and $(2+1)$ dimensions, the model exhibits a phase transition from a symmetric
state with $\langle\phi\rangle=0$ to a broken symmetry phase with finite
$\langle\phi\rangle$ [30, 31]. However, in finite size systems, like the ones
used for simulations, the ground state cannot have a broken-symmetry and there
is no phase transition. Nonetheless, the broken symmetry phase can be
investigated numerically by considering the coupling term $f_{0}\mkern
1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern 1.5mu$ in Section II and
extrapolating the results to the large lattice size ($L\rightarrow\infty$) and
the zero external field ($f_{0}\rightarrow 0$) limits.

It is interesting that, for negative values of $m_{0}^{2}$ and small
interaction strength (_i.e._ when $|m_{0}|^{3}/\lambda_{0}\gg 1$), the ground
state is nearly two-fold degenerate and exhibits properties characteristic of
the broken symmetry phase even for small lattices. A single site system
reduces to a double-well potential Hamiltonian for negative $m_{0}^{2}$. The
field distribution in the ground state is symmetric and double-peaked around
zero, with the two maxima located at the points $\pm\mkern
1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern 1.5mu_{\text{m}}$ which
minimize the potential. The gap is small, decreasing exponentially fast with
increasing magnitude of $|m_{0}|^{3}/\lambda_{0}$. A small external field
$f_{0}$ (of the order of the gap) coupled to the scalar field amplitude
produces a ground state with finite $\langle\mkern
1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern 1.5mu\rangle\approx\mkern
1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern 1.5mu_{\text{m}}$ (or
$\langle\mkern 1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern
1.5mu\rangle\approx-\mkern 1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern
1.5mu_{\text{m}}$, depending on $f_{0}$ sign). Numerical calculations of small
size systems show that the system remains nearly two-fold degenerate when the
number of sites is increased. The gap decreases with an increasing number of
sites, while the energy difference between the second and the first excited
states does not decrease, as can be seen in Fig. 1. This is a consequence of
the kinetic term in the Hamiltonian (the third term in Section II) which
favors similar field configurations at neighboring sites. These properties of
small size systems allows us to investigate quantum state preparation methods
for broken symmetry phase by using classical simulations of small lattices, as
discussed in Section IV.2.2.

##  III $\Phi^{4}$ field on qubits

This section describes the qubit encoding of the bosonic states and the
implementation of the evolution operator corresponding to the lattice
Hamiltonian $H_{\mathrm{lat}}$ in Section II. The representation of bosonic
fields on qubits was discussed in detail in [39]. We will review briefly the
general methodology in Section III.1, present the qubit encoding in Section
III.2, and address the evolution operator implementation in Section III.3.

###  III.1 Finite representation of bosonic fields

The lattice Hilbert space is a direct product of local Hilbert spaces (one at
each lattice site) such that
${\mathscr{H}}=\prod_{j=1}^{N}\otimes{\mathscr{H}}_{j}$, where $j$ labels the
lattice site and $N$ is the number of lattice sites. A local Hilbert space
${\mathscr{H}}_{j}$ is infinite dimensional. The number of bosons contributing
to the wavefunction is, in principle, unbounded. However, since we are
interested in the low-energy physics of the system, we postulate that, at
every lattice site, the number of bosons can be truncated with negligible
error by a cutoff number $N_{b}$.

The eigenvectors $\\{\left|\varphi\right\rangle_{j}\\}$ of the field operator,

$\mkern 1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern
1.5mu_{j}\left|\varphi\right\rangle_{j}=\varphi\left|\varphi\right\rangle_{j},$
(8)

form a convenient basis choice for representing the evolution operator since
the Hamiltonian interaction terms are diagonal in this basis. However, the
eigenvalues $\varphi\in\mathbb{R}$ are continuous and unbounded. Therefore,
discretization procedures are necessary to represent the truncated Hilbert
space in the field amplitude basis. We introduce the discretization procedure
below.

The low-energy subspace of the local Hilbert space ${\mathscr{H}}_{j}$ is
spanned by the states with a number of bosons below the cutoff $N_{b}$ and can
be represented with good accuracy by a finite Hilbert space ${\cal{H}}_{j}$ of
dimension $N_{\varphi}$, with $N_{\varphi}>N_{b}$, as described below. Let
$\\{\left|\varphi_{\alpha}\right\rangle_{j}\\}$ be a set of orthonormal
vectors belonging to ${\cal{H}}_{j}$ with $\alpha=\overline{0,N_{\varphi}-1}$.
We define the _discrete field operators_ ${\Phi}_{j}$ and ${\Pi}_{j}$ acting
on ${\cal{H}}_{j}$ as

$\displaystyle{\Phi}_{j}\left|\varphi_{\alpha}\right\rangle_{j}$
$\displaystyle=\varphi_{\alpha}\left|\varphi_{\alpha}\right\rangle_{j},$ (9)
$\displaystyle{\Pi}_{j}$
$\displaystyle=\mu{\cal{F}}_{j}{\Phi}_{j}{\cal{F}}_{j}^{-1}.$ (10)

where $\varphi_{\alpha}$ is the discrete eigenvalue,

$\displaystyle\varphi_{\alpha}$
$\displaystyle=\Delta_{\varphi}\left(\alpha-\frac{N_{\varphi}-1}{2}\right),~{}~{}\alpha\in\\{0,1,...,N_{\varphi}-1\\}$
(11) $\displaystyle\Delta_{\varphi}$
$\displaystyle=\sqrt{\frac{2\pi}{N_{\varphi}\mu}}.$ (12)

and ${\cal{F}}_{j}$ is the finite Fourier transform,

$\displaystyle{\cal{F}}_{j}=\frac{1}{\sqrt{N_{\varphi}}}\sum_{\alpha,\beta=0}^{N_{\varphi}-1}e^{i\frac{2\pi}{N_{\varphi}}\left(\alpha-\frac{N_{\varphi}-1}{2}\right)\left(\beta-\frac{N_{\varphi}-1}{2}\right)}\left|\varphi_{\alpha}\right\rangle_{j}\left\langle\varphi_{\beta}\right|_{j}.$
(13)

In Eq. 10, $\mu>0$ is the boson mass which is the parameter entering in the
definition of the lattice boson creation and annihilation operators,

$\displaystyle a_{j}^{\dagger}=\sqrt{\frac{\mu}{2}}\mkern
1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern
1.5mu_{j}-i\sqrt{\frac{1}{2\mu}}\mkern
1.5mu\overline{\mkern-1.5mu\Pi\mkern-1.5mu}\mkern
1.5mu_{j},~{}~{}a_{j}=\sqrt{\frac{\mu}{2}}\mkern
1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern
1.5mu_{j}+i\sqrt{\frac{1}{2\mu}}\mkern
1.5mu\overline{\mkern-1.5mu\Pi\mkern-1.5mu}\mkern 1.5mu_{j}.$ (14)

The definition of ${\Phi}_{j}$ given by Eqs. 9, 11 and 12 represents the
discretized and truncated version of Eq. 8. The set of states
$\\{\left|\kappa_{\beta}\right\rangle_{j}\\}_{\beta=\overline{0,N_{\varphi}-1}}$
obtained by applying the Fourier transform to the set
$\\{\left|\varphi_{\alpha}\right\rangle_{j}\\}$,

$\displaystyle\left|\kappa_{\beta}\right\rangle_{j}$
$\displaystyle\equiv{\cal{F}}_{j}\left|\varphi_{\beta}\right\rangle_{j}$ (15)
$\displaystyle=\frac{1}{\sqrt{N_{\varphi}}}\sum_{\alpha=0}^{N_{\varphi}-1}e^{i\frac{2\pi}{N_{\varphi}}\left(\alpha-\frac{N_{\varphi}-1}{2}\right)\left(\beta-\frac{N_{\varphi}-1}{2}\right)}\left|\varphi_{\alpha}\right\rangle_{j}$

are the eigenvectors of the discrete conjugate-field operator ${\Pi}_{j}$
defined by Eq. 10 such that

$\displaystyle{\Pi}_{j}\left|\kappa_{\beta}\right\rangle_{j}=\kappa_{\beta}\left|\kappa_{\beta}\right\rangle_{j}$
(16)

where

$\displaystyle\kappa_{\beta}$
$\displaystyle=\Delta_{\kappa}\left(\beta-\frac{N_{\varphi}-1}{2}\right),~{}~{}\beta\in\\{0,1,...,N_{\varphi}-1\\}$
(17) $\displaystyle\Delta_{\kappa}$
$\displaystyle=\sqrt{\frac{2\pi\mu}{N_{\varphi}}}.$ (18)

Equation 16 is the discretized version of the conjugate-field operator
eigenvalue equation, $\mkern 1.5mu\overline{\mkern-1.5mu\Pi\mkern-1.5mu}\mkern
1.5mu_{j}\left|\kappa\right\rangle_{j}=\kappa\left|\kappa\right\rangle_{j}$,
with continuous and unbounded $\kappa\in\mathbb{R}$.

Different representations corresponding to different values of $\mu$ can be
chosen to construct the finite representation. For a given problem and desired
accuracy, the cutoff $N_{b}$ depends on the boson mass $\mu$. In principle
$\mu$ should be optimized for the lowest possible cutoff $N_{b}$ to reduce the
computing resources. Moreover, as can be seen from Eqs. 12 and 18, the
discretization interval $\Delta_{\varphi}$ of the field amplitude variable and
the discretization interval $\Delta_{\kappa}$ of the conjugate-field variable
are also dependent on the boson mass parameter $\mu$. The parameter $\mu$ can
be tuned to adjust the accuracy of the discretization. Increasing $\mu$
decreases the field variable discretization interval and increases the
conjugate-field discretization interval. The discretized field and conjugate-
field variables are related by a finite Fourier transform, thus
$\Delta_{\varphi}\Delta_{\kappa}=2\pi/N_{\varphi}$. To decrease both
discretization intervals, $\Delta_{\varphi}$ and $\Delta_{\kappa}$, the number
of discretization points $N_{\varphi}$ should be increased. For quantum
simulations, tuning $\mu$ to increase the accuracy of the wavefunction’s
discretization is much easier than the process of optimizing $\mu$ to decrease
the boson number cutoff $N_{b}$, as discussed in [39].

On the subspace of ${\cal{H}}_{j}$ spanned by the first $N_{b}$ eigenstates of
the harmonic oscillator Hamiltonian
($H_{hj}=\frac{1}{2}{\Pi}_{j}^{2}+\frac{1}{2}\mu^{2}{\Phi}_{j}^{2}$), the
discrete field and conjugate-field operators obey, with ${\cal{O}}(\epsilon)$
accuracy, the canonical commutation relation,

$I_{Nb}\left[{\Phi}_{j},{\Pi}_{j}\right]I_{Nb}=iI_{Nb}+{\cal{O}}(\epsilon).$
(19)

Here, $I_{Nb}$ is the projector on the $N_{b}$ size low-energy subspace of the
harmonic oscillator. This is a consequence of the Nyquist-Shannon sampling
theorem applied to the fast decaying boson number wavefunctions, as discussed
in [39]. For a problem of interest, as long as $N_{b}$ is taken large enough
such that the contribution of states with more than $N_{b}$ bosons can be
neglected, the infinite Hilbert space ${\mathscr{H}}_{j}$ can be replaced by
the finite $N_{\varphi}$-size Hilbert space ${\cal{H}}_{j}$, and the lattice
field operators $\mkern 1.5mu\overline{\mkern-1.5mu\Phi\mkern-1.5mu}\mkern
1.5mu_{j}$ and $\mkern 1.5mu\overline{\mkern-1.5mu\Pi\mkern-1.5mu}\mkern
1.5mu_{j}$ [Eq. 5] can be replaced by the discrete operators ${\Phi}_{j}$ and
${\Pi}_{j}$ [Eqs. 9 and 10] with ${\cal{O}}(\epsilon)$ accuracy. For a fixed
$N_{b}$, the error ${\cal{O}}(\epsilon)$ decreases exponentially by increasing
$N_{\varphi}$. For practical purpose, we find numerically that a number of
discretization points $N_{\varphi}=2N_{b}$ yields an accuracy of order
$10^{-4}$.

The finite lattice representation is given by the finite Hilbert space
${\cal{H}}=\prod_{j=1}^{N}\otimes{\cal{H}}_{j}$ of dimension $N_{\varphi}^{N}$
and the set of local field and conjugate field-operators
$\\{{\Phi}_{j}\\}_{j=\overline{1,N}}$ and $\\{{\Pi}_{j}\\}_{j=\overline{1,N}}$
defined by Eq. 9 and Eq. 10, respectively. The discretized field amplitude
basis vectors are

$\displaystyle\left|\varphi_{\alpha}\right\rangle\equiv\left|\varphi_{\alpha
1}\right\rangle_{1}\left|\varphi_{\alpha
2}\right\rangle_{2}...\left|\varphi_{\alpha N}\right\rangle_{N}$ (20)

where

$\displaystyle\alpha=\\{\alpha_{1},\alpha_{2},...,\alpha_{N}\\}~{}~{}\text{with}~{}~{}\alpha_{j}\in\\{0,1,...,N_{\varphi}-1\\}.$
(21)

###  III.2 Qubit encoding of the finite representation

The discretized field amplitude basis
$\\{\left|\varphi_{\alpha}\right\rangle\\}$ [Eq. 20] can be encoded on qubits
using the binary representation of the label $\alpha$ [Eq. 21]. For each site,
a register of $n_{q}=\log_{2}(N_{\varphi})$ qubits is assigned. A local field
amplitude state $\left|\varphi_{\alpha_{j}}\right\rangle_{j}$ at site $j$ is
encoded as

$\displaystyle\left|\varphi_{\alpha_{j}}\right\rangle_{j}\equiv\left|\alpha_{0j}\right\rangle_{j}\left|\alpha_{1j}\right\rangle_{j}...\left|\alpha_{\left(n_{q}-1\right)j}\right\rangle_{j}$
(22)

where
$\left|\alpha_{qj}\right\rangle_{j}\in\\{\left|0\right\rangle,\left|1\right\rangle\\}$
is the $q$-th qubit-state from the register $j$ (i.e. allocated to represent
the field at the site $j$) such that

$\displaystyle\alpha_{j}=\sum_{q=0}^{n_{q}-1}\alpha_{qj}2^{n_{q}-1-q}.$ (23)

Note that the binary variables $\alpha_{qj}\in\\{0,1\\}$ defined by Eq. 23
yield the binary representation of the integer
$\alpha_{j}\equiv[\alpha_{0j}\alpha_{1j}...\alpha_{\left(n_{q}-1\right)j}]$. A
lattice state [Eq. 20] is encoded as a direct product of $N$ local states
encoded by Eq. 22. The lattice states require $N\log_{2}(N_{\varphi})$ qubits
for encoding.

The discrete field operator ${\Phi}_{j}$ acting on the $n_{q}$ qubits assigned
to encode the field at site $j$ can be written as

$\displaystyle{\Phi}_{j}$
$\displaystyle=-\Delta_{\varphi}\sum_{q=0}^{n_{q}-1}2^{n_{q}-1-q}\frac{\sigma^{z}_{qj}}{2}$
(24)

where $\sigma^{z}_{qj}=|0\vphantom{0}\rangle\langle
0\vphantom{0}|_{qj}-|1\vphantom{1}\rangle\langle 1\vphantom{1}|_{qj}$ is the
Pauli Z operator and $q$ is the qubit index. It can be directly checked that
${\Phi}_{j}$ defined here and the vector encoded as in Eq. 22 satisfy the
eigenvalue equation defined by Eqs. 9 and 11.

The definition of the conjugate-field operator ${\Pi}_{j}$ on the qubit space
requires first the qubit implementation of the Fourier transform
${\cal{F}}_{j}$ [see Eq. 10]. The implementation of Quantum Fourier transform
(QFT) on qubits is well known [49]. However, the Fourier transform
${\cal{F}}_{j}$ defined by Eq. 13 is centered, i.e. the summation index runs
from $-(N_{\varphi}-1)/2$ to $(N_{\varphi}-1)/2$, unlike the off-centered QFT
where the summation index runs from $0$ to $N_{\varphi}-1$. As shown in
Appendix A, the Fourier transform is related to the QFT by

$\displaystyle{\cal{F}}_{j}=$ $\displaystyle
e^{-i\frac{N_{\varphi}\delta^{2}}{2\pi}}\prod_{q=0}^{n_{q}-1}R^{z}_{qj}\left(2^{n_{q}-1-q}\delta\right)~{}\text{QFT}_{j}$
$\displaystyle\times\prod_{q=0}^{n_{q}-1}R^{z}_{qj}\left(2^{n_{q}-1-q}\delta\right)$
(25)

where $\delta=\pi\frac{N_{\varphi}-1}{N_{\varphi}}$ and $R^{z}_{qj}$ is a
single-qubit $z$ rotation acting on the qubit $q$ at site $j$ given by

$\displaystyle R^{z}_{qj}(\theta)\equiv
e^{-i\theta\frac{\sigma^{z}_{qj}}{2}}=e^{-i\frac{\theta}{2}}|0\vphantom{0}\rangle\langle
0\vphantom{0}|_{qj}+e^{i\frac{\theta}{2}}|1\vphantom{1}\rangle\langle
1\vphantom{1}|_{qj}.$ (26)

According to Eq. 10, the discrete conjugate-field operator is

$\displaystyle{\Pi}_{j}={\cal{F}}_{j}\left(-\Delta_{k}\sum_{q=0}^{n_{q}-1}2^{q}\frac{\sigma^{z}_{qj}}{2}\right){\cal{F}}_{j}^{-1}.$
(27)

Note that the factor before Pauli $\sigma^{z}_{qj}$ gate is $2^{q}$, unlike
the factor in Eq. 24 which is $2^{n_{q}-1-q}$. This is caused by the fact that
the qubit order is reversed after a QFT gate (unless additional swap
operations are performed to manually reverse the qubit order) [49].

###  III.3 Evolution operator

In order to implement the evolution operator we employ the Trotter-Suzuki [50,
51, 52] expansion. The evolution operator is written as a product of short-
time evolution operators corresponding to the different terms in the
Hamiltonian, called Trotter steps. Here, we present the qubit implementation
of the Trotter steps corresponding to the different terms present in the
$\phi^{4}$ Hamiltonian.

We start with the operator $e^{-i\theta\Phi_{j}}$, where $\theta$ is the time
interval of the Trotter step. This Trotter step implements the evolution of
the term $f_{0}\Phi_{j}$ in Section II. Employing Eq. 24 one has

$\displaystyle
e^{-i\theta{\Phi}_{j}}=\prod_{q=0}^{n_{q}-1}R^{z}_{qj}\left(-2^{n_{q}-1-q}\Delta_{\varphi}\theta\right).$
(28)

It reduces to $n_{q}$ single-qubit $z$ rotations.

The Trotter step $e^{-i\theta{\Phi}_{j}^{2}}$ can be written as

$\displaystyle
e^{-i\theta{\Phi}_{j}^{2}}=e^{-i\theta\Delta_{\varphi}^{2}\frac{N_{\varphi}^{2}-1}{12}}\prod_{p=0}^{n_{q}-1}\prod_{q=0}^{p-1}~{}\mathclap{Z}Z_{pj;qj}\left(\nu_{pq}\right)$
(29)

where

$\displaystyle~{}\mathclap{Z}Z_{pj;qj}(\nu)=$
$\displaystyle\,e^{-i\nu\sigma^{z}_{pj}\sigma^{z}_{qj}}$ (30)
$\displaystyle\nu_{pq}=$
$\displaystyle\,2^{2n_{q}-3-p-q}\Delta_{\varphi}^{2}\theta.$ (31)

The two-qubit gate $~{}\mathclap{Z}Z_{pj;qj}$ acts on the qubit $p$ at site
$j$ and on the qubit $q$ at site $j$. Note that one $~{}\mathclap{Z}Z_{pj;qj}$
can be decomposed into two CNOT gates and one $R^{z}$ gate [53]. Hence, the
Trotter step Eq. 29 consists of $n_{q}(n_{q}-1)$ CNOT gates.

The implementation of the Trotter step $e^{-i\theta{\Pi}_{j}^{2}}$ is given by

$\displaystyle e^{-i\theta{\Pi}_{j}^{2}}$
$\displaystyle={\cal{F}}_{j}e^{-i\mu\theta{\Phi}_{j}^{2}}{\cal{F}}_{j}^{-1}$
$\displaystyle=e^{-i\theta\Delta_{\kappa}^{2}\frac{N_{\varphi}^{2}-1}{12}}{\cal{F}}_{j}\left[\prod_{p=0}^{n_{q}-1}\prod_{q=0}^{p-1}~{}\mathclap{Z}Z_{pj;qj}\left(\nu^{\prime}_{pq}\right)\right]{\cal{F}}_{j}^{-1}$
(32)

where $\nu^{\prime}_{pq}=2^{p+q-1}\theta\Delta_{\kappa}^{2}$. The
$~{}\mathclap{Z}Z$ gate’s parameter $\nu^{\prime}_{pq}$ entering in Section
III.3 can be obtained from Eq. 30 by replacing $n_{q}-1-p\longrightarrow p$
and $n_{q}-1-q\longrightarrow q$, (consequence of reverse qubit order after
applying QFT) and $\Delta_{\varphi}\longrightarrow\Delta_{\kappa}$. Since QFT
requires $n_{q}(n_{q}-1)/2$ CNOT gates, this Trotter step consists of
$3n_{q}(n_{q}-1)$ CNOT gates.

The Trotter step $e^{-i\theta{\Phi}_{j}{\Phi}_{l}}$ corresponding to the
coupling term between the sites $j$ and $l$ is

$\displaystyle
e^{-i\theta{\Phi}_{j}{\Phi}_{l}}=\prod_{p=0}^{n_{q}-1}\prod_{q=0}^{n_{q}-1}~{}\mathclap{Z}Z_{pj;ql}\left(\nu^{\prime\prime}_{pq}\right)$
(33)

where $\nu^{\prime\prime}_{pq}=2^{2n_{q}-4-p-q}\theta\Delta_{\varphi}^{2}$.
This Trotter steps consists of $n_{q}^{2}$ $~{}\mathclap{Z}Z$ gates or
$2n_{q}^{2}$ CNOT gates. Since the interaction is not local, in this case the
$~{}\mathclap{Z}Z$ gates act on one qubit belonging to the qubit register
allocated for the field at site $j$ and on one qubit belonging to the qubit
register allocated for the field at site $l$.

The Trotter step corresponding to the $\phi^{4}$ interaction term is

$\displaystyle e^{-i\theta{\Phi}_{j}^{4}}=$
$\displaystyle\left[\prod_{p=0}^{n_{q}-1}\prod_{q=0}^{p-1}\prod_{r=0}^{q-1}\prod_{s=0}^{r-1}~{}\mathclap{Z}Z\mathclap{Z}Z_{pj;qj;rj;sj}(\rho_{qprs})\right]$
$\displaystyle\times\left[\prod_{p=0}^{n_{q}-1}\prod_{q=0}^{p-1}~{}\mathclap{Z}Z_{pj;qj}(\eta_{pq})\right]e^{i\xi}$
(34)

where

$\displaystyle~{}\mathclap{Z}Z\mathclap{Z}Z_{pj;qj;rj;sj}(\rho)=e^{-i\rho\sigma^{z}_{pj}\sigma^{z}_{qj}\sigma^{z}_{rj}\sigma^{z}_{sj}}$
(35)

and

$\displaystyle\rho_{pqrs}$
$\displaystyle=\frac{3N_{\varphi}^{4}}{32}\frac{1}{2^{p+q+r+s}}\theta\Delta_{\varphi}^{4}$
(36) $\displaystyle\eta_{pq}$
$\displaystyle=\frac{N_{\varphi}^{4}}{16}\frac{1}{2^{p+q}}\left(1-\frac{1}{N_{\varphi}^{2}}-\frac{1}{2^{2p+1}}-\frac{1}{2^{2q+1}}\right)\theta\Delta_{\varphi}^{4}$
(37) $\displaystyle\xi$
$\displaystyle=\left[\frac{(N_{\varphi}^{2}-1)^{2}}{48}-\frac{N_{\varphi}^{4}-1}{120}\right]\theta\Delta_{\varphi}^{4}$
(38)

This step requires $n_{q}(n_{q}-1)(n_{q}-3)(n_{q}-3)/24$ four-qubit
$~{}\mathclap{Z}Z\mathclap{Z}Z$ gates and $n_{q}(n_{q}-1)/2$ two-qubit
$~{}\mathclap{Z}Z$ gates.

| Operator
---
$e^{-i\Phi\theta}$ | $e^{-i\Phi^{2}\theta}$ | $e^{-i\Pi^{2}\theta}$ | $e^{-i\Phi_{j}\Phi_{k}\theta}$ | $e^{-i\Phi^{4}\theta}$
| Number of
---
CNOTs
0 | $n_{q}^{2}-n_{q}$ | $3n_{q}^{2}-3n_{q}$ | $2n_{q}^{2}$ |  | $\frac{1}{4}n_{q}^{4}-\frac{3}{2}n_{q}^{3}$
---
$+\frac{15}{4}n_{q}^{2}-\frac{5}{2}n_{q}$
Table 1: CNOT gates count of the Trotter steps required for the implementation
of the $\phi^{4}$ evolution operator. All to all qubit connectivity is
assumed.

Quantum simulation on near-term quantum devices is mainly limited by the two-
qubit gate fidelities. The implementation of the Trotter step corresponding to
the $\phi^{4}$ interaction term is computationally the most expensive one,
since it requires ${\cal{O}}(n_{q}^{4})$ of two-qubit gates. The number of
CNOT gates for all Trotter steps relevant for the implementation of the
$\phi^{4}$ evolution are summarized in Table 1. For comparison purpose, the
$~{}\mathclap{Z}Z$ and the $~{}\mathclap{Z}Z\mathclap{Z}Z$ gates are
decomposed in two-qubit $CX$ (CNOT) gates and single-qubit $R^{z}$ rotations
[53],

$\displaystyle\mathclap{Z}Z_{p;q}(\nu)=$ $\displaystyle
CX_{p;q}R^{z}_{q}(\nu)CX_{p;q},$ (39)
$\displaystyle~{}\mathclap{Z}Z\mathclap{Z}Z_{p;q;r;s}(\rho)=$ $\displaystyle
CX_{p;q}CX_{q;r}CX_{r;s}R^{z}_{s}(\rho)$ $\displaystyle\times
CX_{r;s}CX_{q;r}CX_{p;q}.$ (40)

Since the number of Trotter steps is proportional to the lattice size, the
computational cost of this algorithm scales linearly with $N$.

##  IV State preparation

Figure 2:  Ground state wavefunctions,
$\langle\varphi_{\alpha}|\psi^{loc}_{g}\rangle_{j}$, of the local Hamiltonian
$H_{loc,j}$ in Eq. 42 represented on $n_{q}=6$ qubits ($N_{\varphi}=64$)
versus the discretization index $\alpha=\overline{0,N_{\varphi}-1}$ [see Eq.
11]. The black circles illustrate the harmonic oscillator ground state
($\lambda_{I}=0$), which is a Gaussian. The red squares (green diamonds)
illustrate the anharmonic oscillator ground state with strong interaction,
$\lambda_{I}/m_{I}^{3}=100$, for a chosen discretization interval
$\Delta=\Delta_{0}\equiv\sqrt{2\pi/(N_{\varphi}m_{I})}$
($\Delta=\Delta_{0}/2$). The blue triangles illustrate the ground state for
the Hamiltonian with negative mass-squqred and small external field. The
magenta stars illustrate the ground state for the Hamiltonian with negative
mass-squared and significant external field. These states are obtained by
employing exact diagonalization. Figure 3:  (a) Fidelity of the quantum states
illustrated in Fig. 2 (same legend) prepared using an optimized quantum
circuit versus the number of circuit’s entanglement layers. (b) The same as in
(a) when the quantum states are represented on $n_{q}=7$ qubits
($N_{\varphi}=128$). (c) The same as in (a) when the quantum states are
represented on $n_{q}=8$ qubits ($N_{\varphi}=256$). The fidelity generally
increases with the number of entanglement layers of the variational circuit.
However the fidelity is not strictly monotonic since the circuit structure
with an odd number of entanglement layer is different than that with an even
number. For the cases with 6, 7 and 8 qubits, the variational method reaches a
fidelity $>0.9999$ with 6 entanglement layers.

This section addresses the preparation of the ground state on qubits for both
normal and broken-symmetry phases. Our method combines variational quantum
circuits and adiabatic evolution, and it is flexible enough to allow tuning of
different parameters to minimize circuit depth. To prepare broken symmetry
states, an interaction between the scalar field and an external field is
introduced that explicitly breaks the $Z_{2}$ symmetry. By properly choosing
the strength of the external field as a function of time during the adiabatic
process, the dual problems of degeneracy and broken symmetry in the ground
state are mitigated, as we discuss in Section IV.2.2.

The Hamiltonian $H$ employed for the quantum simulations of the $\phi^{4}$
model is given by (II) with the lattice field operators replaced by the
discrete field operators, as described in Section III.1. To prepare the ground
state, we divide $H$ in two parts,

$\displaystyle H=H_{loc}+H_{c},$ (41)

where

$\displaystyle H_{loc}$ $\displaystyle=\sum_{j=1}^{N}H_{loc,j}$ (42)
$\displaystyle=\sum_{j=1}^{N}\left(\frac{1}{2}\Pi_{j}^{2}+\frac{1}{2}m_{I}^{2}\Phi_{j}^{2}+\frac{\lambda_{I}}{4!}\Phi_{j}^{4}+f_{I}\Phi_{j}\right),$
$\displaystyle H_{c}$
$\displaystyle=\sum_{j=1}^{N}\Bigg{[}\frac{1}{2}\sum_{e=0}^{d}\left(\Phi_{j+e}-\Phi_{j}\right)^{2}+\frac{1}{2}\delta
m^{2}\Phi_{j}^{2}$ $\displaystyle\quad\quad\quad\
+\frac{\delta\lambda}{4!}\Phi_{j}^{4}+\delta f\Phi_{j}\Bigg{]},$ (43)

where $\delta m^{2}=m_{0}^{2}-m_{I}^{2}$,
$\delta\lambda=\lambda_{0}-\lambda_{I}$ and $\delta f=f_{0}-f_{I}$. The
Hamiltonian $H_{loc}$ is a sum of uncoupled local Hamiltonians $H_{loc,j}$
acting only at the lattice site $j$. The input parameters, $m_{I}^{2}$,
$\lambda_{I}$ and $f_{I}$ should be chosen to ensure that the adiabatic
evolution part of the state preparation is efficient, as we discuss in Section
IV.2. The first term in $H_{c}$ couples the fields at neighboring sites, while
the last three terms in $H_{c}$ are local.

Our state preparation protocol consists of two parts.

  1. 1.

The ground state of $H_{loc}$ is prepared using variational circuits, as we
describe in Section IV.1. It is a direct product of the ground state of
$H_{loc,j}$ at each lattice site $j$, $|\psi_{\textrm{g}}^{loc}\rangle_{j}$:

$\displaystyle|\psi_{\textrm{g}}^{loc}\rangle=\prod_{j=1}^{N}\otimes|\psi_{\textrm{g}}^{loc}\rangle_{j}.$
(44)

  2. 2.

The ground state of the full Hamiltonian is obtained by adiabatic evolution.
The Hamiltonian $H_{c}$ is turned on adiabatically. The system evolves under
the time dependent Hamiltonian,

$H(s)=H_{loc}+\alpha(s)H_{c},$ (45)

from $|\psi_{\textrm{g}}^{loc}\rangle$ to the ground state of the Hamiltonian
Section II. The time $t$ enters in Eq. 45 via the variable $s=t/T$, where $T$
is the total time of the adiabatic process and the function $\alpha(s)$ has
boundary conditions $\alpha(0)=0$ and $\alpha(1)=1$ in the time interval $T$.

###  IV.1 Variational preparation of local states

The ground state $|\psi_{\textrm{g}}^{loc}\rangle_{j}$ of the local
Hamiltonian $H_{loc,j}$ can be prepared accurately using short circuits on the
$n_{q}$ qubits assigned to represent the field at the site $j$. We propose a
hardware-efficient circuit ansatz to prepare the local ground state using one-
and two-qubit gates. The circuit parameters are determined using optimization
algorithms on classical computers, as discussed below.

First, the local wavefunction in the discrete field amplitude basis,
$\langle\varphi_{\alpha}|\psi_{\textrm{g}}^{loc}\rangle_{j}$ with
$\alpha=\overline{0,N_{\varphi}-1}$, is calculated on a classical computer.
This requires the diagonalization of a small size $N_{\varphi}\times
N_{\varphi}$ matrix corresponding to the Hamiltonian $H_{loc,j}$.

Second, parameterized quantum circuits are employed to produce $n_{q}$-qubit
quantum states. We consider circuit ansatzes made by successive one-qubit and
two-qubit layers. A one-qubit layer consists of one $R^{y}(\theta_{i})$
rotation followed by one $R^{z}(\theta_{j})$ rotation on every qubit. A two-
qubit layer, which is responsible for introducing entanglement, consists of CZ
gates acting on neighboring qubits. Qubit pairing in successive entanglement
layers differs from each other and alternates. The quantum state
$\left|\phi_{qc}\left(\bm{\theta}\right)\right\rangle$ depends on $M$ rotation
angles $\bm{\theta}=\left(\theta_{1},...,\theta_{M}\right)$ of the $R^{y}$ and
$R^{z}$ single-qubit gates in the circuit. Since, typically, $n_{q}$ is a
small number ($6\sim 8$), the state
$\left|\phi_{qc}\left(\bm{\theta}\right)\right\rangle$ can be computed on a
classical computer without memory limitation problems, using packages such as
Cirq [54] or Qiskit [55].

Third, the $M$ rotation angles, $\bm{\theta}$, that parametrize the circuit
are chosen so that the fidelity

$F(\bm{\theta})=\left|\langle\phi_{qc}(\bm{\theta})|\psi_{\textrm{g}}^{loc}\rangle_{j}\right|^{2}$
(46)

is as close to $1$ as possible. This can be accomplished, for example, by
using the Covariance Matrix Adaptation Evaluation Strategy (CMA-ES) [56] for
optimization on a classical computer. CMA-ES is an iterative, genetic
algorithm that generates a population of solutions at each iteration. The
covariance matrix, calculated from a population subset with the largest values
of $F(\bm{\theta})$, determines the population of solutions considered at the
next iteration. The algorithm terminates when the best $F(\bm{\theta})$ of the
population stops improving. The most difficult problem we encounter during the
optimization of $F(\bm{\theta})$ is trapping at points of local maxima. We
find that this problem can be avoided when CMA-ES runs with a large population
of solutions.

In Fig. 2, we show the ground states of the local Hamiltonian $H_{loc,j}$
represented on $n_{q}=6$ qubits for different Hamiltonian parameters
calculated using exact diagonalization. We are going to prepare these states
by the parameterized circuits to demonstrate the efficiency of the variational
preparation. For illustration, we have chosen parameters representing
different regimes, such as non-interacting, strong interacting with positive
squared mass, negative squared mass with small external field strength and
negative squared mass with a significant external field strength. Since, as
mentioned in Section III.1 and discussed at large in [39], the discretization
interval $\Delta_{\varphi}$ in an interacting model can be tuned to optimize
the performance of the algorithm, we present examples with
$\Delta_{\varphi}=\sqrt{2\pi/(N_{\varphi}m_{I})},~{}1/2\sqrt{2\pi/(N_{\varphi}m_{I})},~{}2\sqrt{2\pi/(N_{\varphi}m_{I})}$.
Note that the wavefunction representation on qubits depends significantly on
$\Delta_{\varphi}$. For example, the wavefunction plotted with red squares and
the one plotted with green diamonds both correspond to the same Hamiltonian
parameters $\lambda_{I}/m_{I}^{3}=100$ but the discretization interval of the
latter is a factor of $2$ smaller.

The state preparation fidelity of our parametrized circuits is shown in Fig. 3
as a function of entanglement layers, for $n_{q}=6$, $n_{q}=7$ and $n_{q}=8$
qubits. The target states are the one illustrated in Fig. 2. For all examples,
the fidelity of the local ground states is larger that $0.9999$ when at least
$6$ entanglement layers are used. This is sufficient to accurately simulate a
large lattice model. For example, the fidelity to prepare the local ground
state of a lattice with $N=100$ sites is estimated to be $0.9999^{100}\approx
0.99$, which is comparable to the typical two-qubit gate fidelity ($\sim
0.995$) on NISQ devices [3, 4, 5]. If a higher fidelity is needed, more
entanglement layers can be added to the variational circuit.

We do not encounter the difficulties seen in common variational quantum
approaches such as the variational quantum eigensolver (VQE) [6, 7] since the
calculation of the quantum circuit required for state preparation is done on
classical computers. We do not require a global minimum and any solution with
high enough fidelity (i.e. larger than the target accuracy) is acceptable. The
number of qubits $n_{q}$ needed to prepare the local wavefunction is not very
large. For example, $n_{q}=8$ is enough even for the strong interacting
regime, since it can accommodate $\approx 200$ bosons per lattice site with
great precision [39]. Hence, the barren plateau issue [42] is not a
significant concern in our optimization. The optimization problem does not
worsen when the system size is increased, since the optimized wavefunctions
are local. For the preparation of an $N$-site lattice wavefunction
$\left|\psi_{\textrm{g}}^{loc}\right\rangle$ [Eq. 44], a quantum circuit
running $N$ parallel $n_{q}$ quantum circuits should be used.

The variational circuit ansatz used here and constructed from
$\mathrm{R_{Y}}$, $\mathrm{R_{Z}}$ and $\mathrm{CZ}$ gates is just a
representative example. Different hardware-efficient circuit ansantzes can be
employed. For example, the results of this section will be similar if one
changes the $\mathrm{CZ}$ gate to the $\mathrm{CNOT}$ gate and the
$\mathrm{R_{Y}}$ to the $\mathrm{R_{X}}$.

###  IV.2 Adiabatic evolution

Our state preparation method relies on the adiabatic theorem [57], which
relates the ground state of the interacting $\phi^{4}$ Hamiltonian $H$ [Eq.
41] to the ground state of the local Hamiltonian $H_{loc}$ [Eq. 42] under the
action of the time dependent Hamiltonian $H(s=t/T)$ [Eq. 45] for a
sufficiently long time $T$.

There is a vast literature [58] addressing the necessary and sufficient
conditions the adiabatic time $T$ should fulfill. A necessary condition is
given by

$\displaystyle
T\gg\frac{1}{\epsilon}\max_{s\in[0,1]}\left|A_{m0}(s)\right|\text{\ \ for all
$m\neq 0$, with}$ (47) $\displaystyle A_{m0}(s)=\frac{\langle
E_{m}(s)|\dot{E}_{0}(s)\rangle}{E_{m}(s)-E_{0}(s)},$ (48)

where $\left|E_{m}(s)\right\rangle$ is the $m$-th instantaneous eigenstate of
$H(s)$ satisfying $H(s)|E_{m}(s)\rangle=E_{m}(s)|E_{m}(s)\rangle$. The system
starts evolving from $|E_{0}(0)\rangle\equiv|\psi_{\textrm{g}}^{loc}\rangle$.
The dot denotes the derivative with respect to the $s$ variable,
$|\dot{E}_{0}(s)\rangle\equiv\frac{d}{ds}|E_{0}(s)\rangle$, and
$\epsilon=\left|\left|U(1)\left|\psi_{\textrm{g}}^{loc}\right\rangle-\left|\psi_{\textrm{g}}^{trg}\right\rangle\right|\right|$
quantifies the difference between the state at the end of the the adiabatic
evolution $U(1)\left|\psi_{\textrm{g}}^{loc}\right\rangle$ and the target
state
$\left|\psi_{\textrm{g}}^{trg}\right\rangle\equiv\left|E_{0}(1)\right\rangle$.
Equation 47 is not a sufficient condition but provides a good estimate of $T$
for a large number of problems [59, 60]. The combination of this condition
with the relation

$\langle E_{m}(s)|\dot{E}_{0}(s)\rangle=\frac{\langle
E_{m}(s)\vphantom{E_{0}(s)}|\frac{dH(s)}{ds}|E_{0}(s)\vphantom{E_{m}(s)}\rangle}{E_{m}(s)-E_{0}(s)}~{}~{}\text{for}~{}m\neq
0$ (49)

implies that $T$ scales as the square of the minimum excitation gap. When the
excitation gap along the evolution path vanishes, as it does when the system
passes through a critical region, the adiabatic process fails.

The condition Eq. 47 is not always sufficient to ensure adiabatic evolution,
typical examples where it fails being Hamiltonians with oscillatory terms. A
further, necessary condition for the validity of the adiabatic approximation
[60] is given by

$\displaystyle
T\gg\frac{1}{\epsilon}\max_{s\in[0,1]}\left|\frac{d}{ds}A_{m0}(s)\right|\text{\
\ for all $m\neq 0$}.$ (50)

In our case, this second adiabatic condition is relevant for the preparation
of the broken symmetry state, as we will discuss in Section IV.2.2.

The adiabatic time is shown to be proportional to the changing rate of the
ground state wavefunction along the adiabatic path [61, 62], i.e.
$T\propto\int_{0}^{1}ds\left|\left|~{}\left|\dot{E}_{0}(s)\right\rangle~{}\right|\right|$.
Intuitively we expect that, the closer are the initial and the target
wavefunctions, the smaller is the overall changing rate of the ground state
along the path, and implicitly the required adiabatic time. In fact, for both
normal and broken symmetry phase preparations, we observe a direct correlation
between the adiabatic time and the local overlap of the initial and the final
wavefunctions. Namely, a larger overlap correlates with a shorter adiabatic
time, as we discuss in Sections IV.2.1 and IV.2.2.

The initial wavefunction is determined by the parameters $m_{I}$,
$\lambda_{I}$ and $f_{I}$. Next we will explore the influence of these
parameters on the adiabatic process. Since these parameters are adjustable in
our algorithm, we will make recommendations for their choices. For the
preparation of the broken symmetry states we will also address the
complications caused by the degeneracy of the ground state.

####  IV.2.1 Preparation of normal phase states

Figure 4: Adiabatic time $T$ required to prepare the $\phi^{4}$ normal state
with $0.97$ fidelity versus $m_{I}^{2}$ for different values of the parameters
$m_{0}^{2}$ and $\lambda_{0}$ for (a) two-site and (b) four-site lattices. $T$
sensitivity on $m_{I}^{2}$ increases with decreasing $m_{0}^{2}$. (c) and (d)
Local fidelity $F$, Eq. 52, measuring the overlap of the initial and the final
local density matrices, versus $m_{I}^{2}$ for two-site and, respectively,
four-site lattices. The values of $m_{I}^{2}$ yielding the smallest $T$ yield
the largest $F_{loc}$.  Figure 5: Adiabatic time $T$ as a function of
$\delta\lambda$ for $m_{I}^{2}=0.1$, different values of the parameters
$m_{0}^{2}$ and $\lambda_{0}$ for (a) two-site and (b) four-site lattices. The
adiabatic time $T$ dependence on $\delta\lambda=\lambda_{0}-\lambda_{I}$ is
weak.

Finding the optimal adiabatic path for adiabatic evolution is difficult
without the knowledge of the system’s excitation spectrum. However, our goal
in this section is less ambitious, and consists in investigating the effect of
the initial wavefunction $\left|\psi_{\textrm{g}}^{loc}\right\rangle$ (Eq. 44)
on the adiabatic process. For normal phase preparation we consider only
adiabatic paths with linear time dependence, i.e. we take $\alpha(s)=s$ in Eq.
45.

The normal phase of the $\phi^{4}$ model is characterized by a non-degenerate
ground state. The symmetry-breaking external field is unnecessary for ground
state preparation in this case, and we set $f_{I}=0$ and $\delta f=0$ in Eqs.
42 and IV. The dependence of $\left|\psi_{\textrm{g}}^{loc}\right\rangle$ on
$m_{I}$ and $\lambda_{I}$ can be understood by writing the local Hamiltonian
(Eq. 42) as

$\displaystyle\frac{H_{loc,j}}{\left|m_{I}\right|}$
$\displaystyle=\frac{1}{2}\left(\frac{\Pi_{j}}{\sqrt{\left|m_{I}\right|}}\right)^{2}+\text{sgn}(m_{I}^{2})\frac{1}{2}\left(\sqrt{\left|m_{I}\right|}{\Phi_{j}}\right)^{2}$
(51)
$\displaystyle+\frac{1}{4!}\frac{\lambda_{I}}{\left|m_{I}\right|^{3}}\left(\sqrt{\left|m_{I}\right|}\Phi_{j}\right)^{4}.$

Up to a field amplitude scaling factor $\sqrt{\left|m_{I}\right|}$ the
eigenfunctions of $H_{loc,j}$ are solely determined by the sign of $m_{I}^{2}$
and the ratio $\lambda_{I}/\left|m_{I}\right|^{3}$, while $\left|m_{I}\right|$
acts as a scaling factor for the energy.

The adiabatic evolution on small lattices is simulated on classical computers,
using the Trotterization method described in Section III.3. We observe that
$m_{I}$ has a strong influence on the adiabatic time and there is an optimal
value that minimizes the adiabatic time. Figure 4 –(a) and (b) shows the
adiabatic time $T$ needed to prepare the ground state with $0.97$ fidelity as
a function of $m_{I}^{2}$ for different values of the dimensionless bare
parameters $m_{0}^{2}$ and $\lambda_{0}$ for two- and four-site lattices. The
value of $T$ is especially sensitive to the mass parameter $m_{I}$ when
$m_{0}$ is small, a parameter regime relevant when taking the continuous limit
$a\rightarrow 0$ (remember that $m_{0}=m_{b}a$). We find that, in general, the
optimal value of $m_{I}^{2}$ is larger than the bare mass $m_{0}^{2}$.

The dependence of the adiabatic time on $m_{I}$ can be understood by
considering the strong influence of $m_{I}$ on the initial wavefunction, since
$\sqrt{\left|m_{I}\right|}$ acts as a scale factor for the field amplitude
variable. Two competing effects come into play to determine the optimal
$m_{I}$. On one hand, the initial gap increases with increasing $m_{I}^{2}$,
which favors the adiabatic process. On the other hand, increasing $m_{I}^{2}$
reduces the width of the initial wavefunction in the field amplitude basis.
When the initial wavefunction is too narrow, significant changes of the
wavefunction are required during evolution, which increases the adiabatic
time. In fact, for small lattices, we observe a direct correlation between the
optimal adiabatic time and the local fidelity defined as [49]

$\displaystyle
F_{loc}={}_{j}\langle\psi_{\textrm{g}}^{loc}\vphantom{\psi_{\textrm{g}}^{loc}}|\rho_{j}|\psi_{\textrm{g}}^{loc}\vphantom{\psi_{\textrm{g}}^{loc}}\rangle_{j},$
(52)

where $\rho_{j}=\prod_{k\neq
j}\mathrm{Tr}_{k}|\psi_{\textrm{g}}^{trg}\vphantom{\psi_{\textrm{g}}^{trg}}\rangle\langle\psi_{\textrm{g}}^{trg}\vphantom{\psi_{\textrm{g}}^{trg}}|$
is the reduced density matrix at site $j$ of the system final ground state
$\left|\psi_{\textrm{g}}^{trg}\right\rangle\equiv|E_{0}(s=1)\rangle$. In Eq.
52, $\mathrm{Tr}_{k}$ denotes the partial trace over a local basis at site
$k$. The smallest $T$ is obtained for the values of $m_{I}^{2}$ which yield
the largest values of $F_{loc}$, as can be seen by comparing the top panel and
bottom panel of Fig. 4.

Numerical simulations on small lattices reveal a weak dependence of the
adiabatic time on the coupling strength $\lambda_{I}$ used to prepare the
initial state. Fig. 5 shows the time $T$ required to prepare the two- and
four-site lattice ground state with fidelity $0.99$ as a function of
$\delta\lambda=\lambda_{0}-\lambda_{I}$. The parameters $m_{0}^{2}$ and
$\lambda_{0}$ chosen in these examples are small to accentuate the region
where $T$ is sensitive to $m_{I}^{2}$. The chosen value of $m_{I}^{2}$ is
close to the optimal one for $\delta\lambda=0$. Particularly for the four-site
lattice, the adiabatic time is nearly independent of $\lambda_{I}$.

####  IV.2.2 Preparation of broken-symmetry states

Two issues related to the vanishing of the excitation gap must be addressed
when preparing states in the broken-symmetry phase. One is the potential
crossing of the critical region during adiabatic evolution. The other is the
degeneracy of the broken symmetry ground state. These problems can be overcome
by introducing a time-dependent external field that couples to the scalar
field during the adiabatic evolution. The linear term in the Hamiltonian
breaks the symmetry $\phi\to-\phi$ and, hence, the degeneracy of the ground
state. The field strength also controls the energy gap between the ground and
first excited state. By dividing the adiabatic evolution into two steps, we
can focus on and discuss the gap and degeneracy problems independently. The
first adiabatic evolution starts from the ground state of the local
Hamiltonian [Eq. 42] and ends in the ground state of the full Hamiltonian with
a finite external field. The second adiabatic evolution starts from the
terminus of the first one and ends when the external field is brought to near
vanishing values (of the order of the desired error). Our notation is such
that the external field changes from $f_{I}$ to $f_{M}$ during the first stage
of the evolution and from $f_{M}$ to $f_{F}$ in the second stage.

##### First adiabatic path. Avoiding the critical region. ($f_{I}$
$\longrightarrow$ $f_{M}$)

Figure 6:  Diagnostics on the adiabatic evolution algorithm: (a) time $T$ to
reach a fidelity of 0.97 as a function of initial $m_{I}^{2}$ and external
field strength; (b) the local overlap in the same coordinate space; (c)
$m_{I}^{2}$ which minimizes the adiabatic time and maximizes the local overlap
for fixed $f_{I}$ versus $f_{I}$; (d) Adiabatic time $T$ and the local overlap
versus $f_{I}$ for $m_{I}^{2}(f_{I})$ minimizing $T$ and extracted from (c).
(Parameters: $m_{0}^{2}=-0.22,\lambda_{0}=\lambda_{I}=0.1,f_{M}=0.0011$)

Previous studies of quantum algorithms for $\phi^{4}$ field theory [20]
proposed state preparation via adiabatic evolution starting from the ground
state of a non-interacting Hamiltonian. Since the non-interacting ground state
belongs to the normal phase region of the phase diagram, preparing broken
symmetry states in this way implies crossing the critical region characterized
by a vanishing excitation gap. This is problematic, since the adiabatic
process requires a finite gap. Here we avoid crossing the critical region by
starting the adiabatic evolution from a broken symmetry state.

The first adiabatic path starts from a local state coupled to an external
field $f_{I}$ and ends in the ground state of the $\phi^{4}$ model coupled to
the external field $f_{M}$. The initial state is the ground state of Eq. 42
and is prepared variationally as described in Section IV.1. The first
adiabatic process here is described by Eq. 45 with linear time dependence,
$\alpha(s)=s$, and by Section IV with $\delta f=f_{M}-f_{I}$. At the end of
the first adiabatic path, the term containing the non-local coupling between
sites is fully switched on. The finite external field ensures a robust gap
during the evolution. Since the ground state of the broken symmetry phase of
the $\phi^{4}$ model in zero external field is doubly degenerate (or nearly
double degenerate for finite size lattices) and well separated from the rest
of the spectrum (as numerical simulations presented in Fig. 1 shows), $f_{M}$
can be chosen small enough such that the low-energy spectrum of the system at
the end of the first adiabatic path can be approximated by a coupled two-level
system (_i.e._ $f_{M}|\sum_{j}\langle 0|\Phi_{j}|1\rangle|\ll E_{2}-E_{1}$,
where $E_{1}$ and $E_{2}$ are the energies of the first and second excited
states, respectively). This choice of $f_{M}$, while providing a significant
gap during the first adiabatic path, will allow us to investigate analytically
the second adiabatic process where the external field is taken to vanishing
values.

Figure 7:  Initial and target wavefunctions. The probability distribution
along one of the field-amplitude coordinates (dimensionless) for different
initial wavefunctions and the target wavefunction being the ground state of
the full Hamiltonian. (a) $m_{I}^{2}$ dependence; (b) external field
dependence. (Parameters:
$m_{0}^{2}=-0.22,\lambda_{0}=\lambda_{I}=0.1,f_{M}=0.0011$)  Figure 8:
Adiabatic evolution time and local overlap as a function of $m_{I}^{2}$ and
$\lambda_{I}$ when $f_{I}$ is fixed to the value optimal for
$\lambda_{0}=\lambda_{I}$. (a) The time $T$ to reach a fidelity of 0.97 as a
function of input parameters $m_{I}^{2}$ and initial $\lambda_{I}$; (b) The
local overlap between the initial and the target wavefunction as a function of
input parameters $m_{I}^{2}$ and initial $\lambda_{I}$; (c) $m_{I}^{2}$
yielding the shortest $T$ (red triangles) and the largest local overlap (blue
triangles) at fixed $\lambda_{I}$ versus $\lambda_{I}/\lambda_{0}$. (d) $T$
and local overlap versus $\lambda_{I}/\lambda_{0}$ for
$m_{I}=m_{I}(\lambda_{I})$ yielding the shortest adiabatic time at fixed
$\lambda_{I}$ extracted from (c). The minimum adiabatic time and the maximum
local overlap are found around $\lambda_{I}=\lambda_{0}$. (Hamiltonian
parameters: $m_{0}^{2}=-0.8,\lambda_{0}=0.6,f_{M}=0.0014$)

Our goal is to develop an algorithm that minimizes the adiabatic time by
varying the algorithm input parameters, i.e., the initial magnitude field
$f_{I}$, the initial mass $m_{I}^{2}$ and the initial interaction strength
$\lambda_{I}$. Similar to the normal phase preparation, we find that the
adiabatic time remains correlated with the local overlap between the initial
state at $f_{I}$ and the intermediate state at $f_{M}$. In the following, we
investigate the adiabatic time dependence on input parameters through
numerical simulations of two, three and four site lattices.

Simulations on small size lattices show that the shortest adiabatic times are
correlated with the largest local overlap between the initial and final
wavefunctions, similar to the normal-phase case discussed in Section IV.2.1.
An example is shown in Fig. 6 for the Hamiltonian parameters
$m_{0}^{2}=-0.22$, $\lambda_{0}=0.1$, and $f_{M}=0.0011$. For this choice of
the parameters, the coupling between sites (which is $1/2$ in the
dimensionless units, see the first term in Section IV) is the dominant term in
the Hamiltonian. The parameter $\lambda_{I}$ is fixed and equal to
$\lambda_{0}$. In Fig. 6 (a) we show the adiabatic time required to reach the
fidelity $0.97$ while in Fig. 6 (b) we show the local overlap, both as a
function of $m_{I}^{2}$ and $f_{I}$. The correlation can be clearly observed
in Fig. 6(c), where we plot the value of $m_{I}^{2}$ that yields the smallest
$T$ (red upper triangles) and the largest local overlap (blue upside down
triangles) for a fixed $f_{I}$. In Fig. 6(d), we show the adiabatic time and
local overlap versus $f_{I}$, with $m_{I}^{2}=m_{I}^{2}(f_{I})$ extracted from
Fig. 6(c) to minimize $T$ at fixed $f_{I}$. We find that for the optimal
adiabatic time, the external field $f_{I}\gg f_{M}$ and the initial mass
parameter $m_{I}^{2}>0>\left|m_{0}^{2}\right|$. Similar investigations (with
fixed $\lambda_{I}$) for other parameters of the Hamiltonian, including
smaller $m_{0}^{2}$ and $\lambda_{0}$ lead to similar conclusions (not shown).

The relation between the local overlap (and implicitly of the adiabatic time)
and the parameters $m_{I}^{2}$ and $f_{I}$ can be understood by investigating
the effect of these parameters on the initial wavefunction. The initial
wavefunction’s main peak position and the peak’s width are dependent on both
$f_{I}$ and $m_{I}^{2}$. In Fig. 7, we show the initial wavefunction
distribution $|\langle\varphi_{\alpha}|\psi_{g}^{loc}\rangle_{j}|^{2}$ for
different input parameters $m_{I}^{2}$ and $f_{I}$ together with the local (at
site $j$) probability distribution of the target wavefunction, defined as
$p(\varphi)_{j}=\langle\varphi_{\alpha}\vphantom{\varphi_{\alpha}}|\rho_{j}|\varphi_{\alpha}\vphantom{\varphi_{\alpha}}\rangle$
[see Eq. 52 for the definition of $\rho_{j}$]. In the figure, the local field
amplitude $\left|\varphi_{\alpha}\right\rangle$ is discretized to
$N_{\varphi}=32$ points, corresponding to using five qubits for each site.
Figure 7(a) shows that the maximum overlap is obtained by choosing a value of
$f_{I}\gg f_{M}$, in agreement with the plot shown in Fig. 6(d). The initial
wavefunction with the largest local overlap (green crosses) with the target
wavefunction is centered at the same location and has a similar width as the
target wavefunction (black dots). By increasing (decreasing) $f_{I}$ while
keeping $m_{I}^{2}$ fixed, the wavefunction peak moves to the left (right). If
the initial mass $m_{I}^{2}$ is increased (decreased) such to keep the peak
aligned with the target wavefunction’s peak, the wavefunction distribution
becomes too narrow (wide) compared to the target one, as illustrated with red
squares (blue pluses). Consequently the overlap decreases. The case where the
external field $f_{I}$ is chosen small, comparable to $f_{M}$, is shown in
Fig. 7(b). In this case, the local overlap is not optimal because the initial
wavefunction exhibits a double-peak structure (green crosses and blue pluses).

An exhaustive, numerical search in a three dimensional space for the point
$(f_{I},m_{I},\lambda_{I})$ which minimizes the adiabatic time is infeasible.
Instead, we explore the adiabatic time dependence in the vicinity of the point
$\left({f_{I}}_{0},{m_{I}}_{0},\lambda_{0}\right)$, where ${f_{I}}_{0}$ and
${m_{I}}_{0}$ are the initial external field and initial mass parameter which
yield the shortest adiabatic time when $\lambda_{I}=\lambda_{0}$, as described
in the example shown in Fig. 6. We find that modifying $\lambda_{I}$ in the
vicinity this point does not reduce the adiabatic time. For example, in Fig.
8, we show results for Hamiltonian parameters $m_{0}^{2}=-0.8$,
$\lambda_{0}=0.6$, and $f_{M}=0.001$. We keep $f_{I}={f_{I}}_{0}=4.0008$ fixed
and vary $\lambda_{I}$ and $m_{I}^{2}$. The adiabatic time and the local
overlap as a function of $\lambda_{I}/\lambda_{0}$ and $m_{I}^{2}$ are shown
in Fig. 8(a) and (b), respectively. The minimum adiabatic time correlates with
the maximum local overlap as in the previous example, as can be seen from Fig.
8(c) where the $m_{I}^{2}$ which yields the shortest adiabatic time and
$m_{I}^{2}$ which yields the largest local overlap for fixed $\lambda_{I}$ are
shown. While the adiabatic time dependence on $\lambda_{I}$ is not negligible,
we find that the shortest adiabatic time occurs when
$\lambda_{I}\approx\lambda_{0}$, as can be seen in Fig. 8(d). Simulations for
other Hamiltonian parameters (not shown), lead to the same conclusion: varying
the input parameter $\lambda_{I}$ in the vicinity of $\lambda_{0}$ does not
decrease the adiabatic time.

In Figs. 6 and 8, it appears that there is a nearly linear relation between
$f_{I}$ and the corresponding best $m_{I}^{2}$ for both the adiabatic time and
local overlap. This could be explained by the fact that, in this case, the
local overlap is maximized when the potential energy minimum is the same for
the initial and final Hamiltonians. An approximation of that relation can be
obtained using the quasi-classical argument given in Appendix B and is shown
as black dashed lines.

##### Ground state degeneracy. ($f_{M}$ $\longrightarrow$ $f_{F}$)

The broken-symmetry phase of the $\phi^{4}$ model is characterized by a two-
fold degenerate ground state in the thermodynamic limit and a non-zero value
of the order parameter. In numerical simulations, the broken symmetry ground
state $\left|\psi_{a0}\right\rangle$ can be obtained from the ground state
$\left|\psi_{a}\left(L,f\right)\right\rangle$ of a system of finite size $L$
coupled to an external field $f$

$\displaystyle H(f)$ $\displaystyle=H+f\sum_{j=1}^{N}\Phi_{j}$ (53)

by taking the limits

$\left|\psi_{a0}\right\rangle=\lim_{f\to
0}\lim_{L\to\infty}\left|\psi_{a}\left(L,f\right)\right\rangle.$ (54)

Equation 53 is the same as Section II with $f\equiv f_{0}$ and with the field
operators replaced by the discrete field operators. The limiting order in Eq.
54 is important. When $L$ is finite, the system in the absence of the external
field is not truly degenerate. An arbitrary small external field can drive the
system to a broken symmetry state only after the limit $L\rightarrow\infty$.
When estimating the limit Eq. 54 numerically, $f$ should be decreased and $L$
increased subject to the condition that $|f\langle
0\vphantom{1}|\sum_{j}\Phi_{j}|1\vphantom{0}\rangle|\gg\Delta_{0}$, where the
ground state of the system in zero external field is $\left|0\right\rangle$,
the first excited state is $\left|1\right\rangle$, and the energy gap between
them is $\Delta_{0}=E_{1}-E_{0}$.

The degeneracy of the ground state in the broken-symmetry phase is a problem
for the adiabatic process, since the vanishing gap implies long adiabatic
times. However, the quadratic adiabatic time scaling of the vanishing gap
problem can be improved to linear by choosing an appropriate time dependence
of the external field during the adiabatic process, as discussed in this
section. We focus here on the adiabatic evolution in the vicinity of the final
state characterized by vanishing external field. In this region, we assume
that the nearly double degenerate ground state is well separated from the rest
of the spectrum. The external field during our adiabatic process is always
small such that the second term in Eq. 53 can be considered a small
perturbation.

We denote the two low-energy states in the presence of the external field $f$
by $\left|\psi_{a}\left(f\right)\right\rangle$ and
$\left|\psi_{b}\left(f\right)\right\rangle$. As described in Appendix C, a
perturbative analysis reveals that the difference between the ground state of
the system coupled to the external field $f$,
$\left|\psi_{a}(f)\right\rangle$, and the broken symmetry state,
$\left|\psi_{a0}\right\rangle$, is

$\left|\left|\left|\psi_{a}\left(f\right)\right\rangle-\left|\psi_{a0}\right\rangle\right|\right|=\frac{1}{2}fB+O(f^{2}),$
(55)

where $B$, explicitly derived in Eq. 111, is a quantity independent of $f$ and
dependent on the $\sum_{j}\Phi_{j}$ matrix elements coupling the low-energy
and the high-energy states. The gap dependence on $f$ is given by (see Eq.
112)

$\Delta_{ba}\equiv E_{b}-E_{a}=2fv+O(f^{2}),$ (56)

where

$v=\left|\langle\psi_{a0}\vphantom{\psi_{a0}}|\sum_{j}\Phi_{j}|\psi_{a0}\vphantom{\psi_{a0}}\rangle\right|$
(57)

is equal to the order parameter in the broken symmetry phase (see also Eq.
92). The wavefunction dependence on $f$ yields

$\displaystyle\langle\psi_{b}\left(f\right)\vphantom{\psi_{a}\left(f\right)}|\frac{d}{df}|\psi_{a}\left(f\right)\vphantom{\psi_{b}\left(f\right)}\rangle=\frac{B}{2}+O(f),$
(58)

as derived in Eq. 113.

The adiabatic process starts from $\left|\psi_{ini}\right\rangle$, the ground
state of Hamiltonian Eq. 53 with the external field $f(s=0)\equiv f_{M}$. At
the end of the adiabatic evolution the external field is $f(s=1)\equiv f_{F}$.
The error associated with the broken symmetry state preparation is

$\displaystyle\varepsilon=\left|\left|U(1)\left|\psi_{ini}\right\rangle-\left|\psi_{a0}\right\rangle\right|\right|,$
(59)

where $U(1)\left|\psi_{ini}\right\rangle$ is the adiabatically prepared state
and $\left|\psi_{a0}\right\rangle$ is the broken symmetry state (target
state). There are two contribution to $\varepsilon$. The first one is caused
by the finite value of the final external field, and using Eq. 55 is given by

$\displaystyle\varepsilon_{f}=\left|\left|\left|\psi_{a}(f_{F})\right\rangle-\left|\psi_{a0}\right\rangle\right|\right|=\frac{1}{2}f_{F}B,$
(60)

where $\left|\psi_{a}(f_{F})\right\rangle$ is the ground state of Hamiltonian
Eq. 53 with $f=f_{F}$. The second contribution is the error of the adiabatic
process defined as the difference between the adiabatically prepared state and
the ground state when the external field is $f_{F}$,

$\displaystyle\varepsilon_{ad}=\left|\left|U(1)\left|\psi_{ini}\right\rangle-\left|\psi_{a}(f_{F})\right\rangle\right|\right|,$
(61)

The triangle inequality implies that the total error (defined by Eq. 59) is
bounded by the sum of these two contributions

$\displaystyle\varepsilon\leq\varepsilon_{f}+\varepsilon_{ad}.$ (62)

The challenge to preparing the broken symmetry state through adiabatic
evolution can be understood by inspecting the adiabatic condition Eq. 47.
Employing Eqs. 56 and 58 the term $A_{10}$ reads

$\displaystyle A_{10}(s)=\frac{B}{4v}\frac{1}{f}\frac{df}{ds}.$ (63)

When $\frac{1}{f}\frac{df}{ds}$ is large, $\left|A_{10}\right|$ and,
implicitly, the adiabatic time become large.

For example, a large adiabatic time, with quadratic dependence on accuracy, is
required for an external field with linear time dependence

$\displaystyle f(s)=f_{M}-\delta_{f}s,$ (64)

where $\delta_{f}=f_{M}-f_{F}$. The magnitude of $A_{10}(s)$ at the end of the
adiabatic process becomes

$\displaystyle\left|A_{10}(s=1)\right|=\frac{B^{2}\delta_{f}}{8v}\frac{1}{\varepsilon_{f}}=\frac{\varepsilon_{s=0}\left(\varepsilon_{s=0}-\varepsilon_{f}\right)}{\Delta_{s=0}}\frac{1}{\varepsilon_{f}},$
(65)

where $\Delta_{s=0}=2f_{M}v$ is the initial gap of this adiabatic path and
$\varepsilon_{s=0}=Bf_{M}/2$ is the difference between the initial state and
the target state. If one naively employs Eq. 47 to estimate adiabatic time and
assumes $\varepsilon_{ad}\approx\varepsilon_{f}\approx\varepsilon/2$, the
adiabatic time reads
$T\approx\frac{B^{2}\delta_{f}}{8v}\frac{1}{\varepsilon_{ad}\varepsilon_{f}}\approx\frac{B^{2}\delta_{f}}{2v}\frac{1}{\varepsilon^{2}}$.
However, the second adiabatic condition Eq. 50, should also be considered for
assessing the validity of the adiabatic approximation. In this case it
predicts an even longer adiabatic time,
$T\approx\frac{1}{\varepsilon_{ad}}\left|\frac{dA_{10}}{ds}(s=1)\right|=\frac{B^{3}\delta_{f}^{2}}{16v\varepsilon_{ad}\varepsilon_{f}^{2}}\approx\frac{B^{3}\delta_{f}^{2}}{2v\varepsilon^{3}}$,
thus an adiabatic time scaling as ${\cal{O}}(\varepsilon^{-3})$. Fortunately,
in the perturbative region where the system can be reduced to a two level
effective model, one can do better than employing the adiabatic conditions
Eqs. 47 and 50 for estimating the adiabatic time. Analytical and numerical
calculations, presented in Section C.3.1, find that the best adiabatic time is
obtained when $f_{F}$ is taken small such that
$\varepsilon_{f}\ll\varepsilon_{ad}\approx\varepsilon$. In this case (see also
Eq. 138)

$\displaystyle T\approx\frac{\pi
B^{2}\delta_{f}}{16v}\frac{1}{\varepsilon^{2}}\approx\frac{\pi}{2}\frac{\varepsilon_{s=0}}{\Delta_{s=0}}\left(\varepsilon_{s=0}-\varepsilon_{f}\right)\frac{1}{\varepsilon^{2}}.$
(66)

Thus, for an external field with linear time dependence the required adiabatic
time scales as ${\cal{O}}(\varepsilon^{-2})$. This is an increase by a factor
of $\varepsilon^{-1}$ when compared to systems with non-vanishing excitation
gaps.

Since the excitation gap vanishes at the end of the adiabatic path, an
external field time dependence which slows towards the end of the adiabatic
evolution is expected to improve the adiabatic process. As Eq. 63 predicts,
the adiabatic time does not blow up when $\frac{1}{f}\frac{df}{ds}$ is kept in
bounds. For example, if the external field decreases exponentially, $T$ scales
as ${\cal{O}}(\varepsilon^{-1}\ln(\varepsilon^{-1}))$, which is better then
${\cal{O}}(\varepsilon^{-2})$ for a linearly decreasing field. Indeed, by
choosing

$\displaystyle f(s)=f_{M}\exp(-\gamma s),~{}\text{with
}\gamma=\ln(f_{M}/f_{F})=\ln(\varepsilon_{s=0}/\varepsilon_{f}),$ (67)

one has

$\displaystyle\left|A_{10}\right|=\gamma\frac{B}{4v}=\gamma\frac{\varepsilon_{s=0}}{\Delta_{s=0}},~{}~{}\text{
and }~{}~{}\frac{d}{ds}A_{10}=0.$ (68)

The adiabatic conditions Eqs. 47 and 50 predict an adiabatic time
$T\approx\varepsilon_{ad}^{-1}\ln(\varepsilon_{f}^{-1})$. Taking
$\varepsilon_{f}\approx\varepsilon_{ad}\approx\varepsilon/2$ this implies
$T\approx\varepsilon^{-1}\ln(\varepsilon^{-1})$. In fact, explicit analytical
and numerical calculations for a two-level systems, presented in Section
C.3.2, find that

$\displaystyle
T\approx\frac{B}{2v\varepsilon_{ad}}\ln\left(\frac{Bf_{F}}{2\varepsilon_{f}}\right)=2\frac{\varepsilon_{s=0}}{\Delta_{s=0}\varepsilon_{ad}}\ln\left(\frac{\varepsilon_{s=0}}{\varepsilon_{f}}\right).$
(69)

For a desired accuracy $\varepsilon$, one can show that (see Eqs. 150 and 151)

$\displaystyle\frac{B}{2v\varepsilon}\ln\left(\frac{Bf_{F}}{2\varepsilon}\right)\leq
T\leq\frac{B}{v\varepsilon}\ln\left(\frac{Bf_{F}}{\varepsilon}\right).$ (70)

Note that $T$ depends on the ratio $B/v=4\varepsilon_{s=0}/\Delta_{s=0}$. A
small ratio $B/v$ implies that the wavefunction’s dependence on the external
field is much weaker than the gap’s dependence on the external field. For the
small size lattices explored here, we find that $B/v$ is $\sim 10^{-2}$ (using
the dimensionless units defining the Hamiltonian Section II) close to the
critical region and is decreasing rapidly when moving further away from the
critical region. We conclude that, for the adiabatic preparation of the broken
symmetry states, the farther away from the critical region the states are the
smaller adiabatic time is needed.

####  IV.2.3 Adiabatic evolution input parameters selection

Finding the optimal input parameters $f_{I}$, $m_{I}^{2}$ and $\lambda_{I}$
for the preparation of large size lattice states is challenging. We propose a
strategy that avoids very long quantum circuits that are infeasible for
limited-coherence near-term quantum hardware. First, determine the optimal
input parameters for a small lattice, as discussed here. Then, use those
values as the starting point in a search for the optimal input parameters in
increasingly larger systems, guided by local overlap measuring. The local
overlap can be obtained by employing the SWAP test method [63]. This strategy
implies running multiple evolution circuits for systems smaller than the
target one. However, the evolution time for these runs is close to optimal,
implying relatively short circuits.

##  V Conclusion

In this paper, we present a circuit implementation of the evolution operator
of the $\phi^{4}$ lattice Hamiltonian on a qubit quantum computer and an
algorithm to prepare states in both the normal and broken-symmetry phases. The
implementation is efficient in its use of resources and leverages the
adjustable parameters of the problem to produce high fidelity states. It is
suitable for near-term quantum computers.

The evolution operator is implemented using the Trotterization method. The
scalar field is encoded on the qubits using the discretized field amplitude
representation [39] and requires a small number of qubits per site,
$n_{q}\approx 6\sim 8$. This number of qubits is adequate for exponential
precision in even the strongly interacting regimes. The required number of
qubits and gates per Trotter step scale with the lattice size $N$. The most
computationally expensive part of the evolution is the Trotter step associated
with the $\phi^{4}$ interaction, which requires a number of two-qubit gates
proportional to $Nn_{q}^{4}$.

Our state preparation combines a variational approach with adiabatic
evolution. The $\phi^{4}$ Hamiltonian is split into two parts: a local
Hamiltonian consisting of a sum of local terms (_i.e_ uncoupled terms
describing interaction at each site) and an inter-site coupling Hamiltonian
that is switched on adiabatically to restore the full Hamiltonian. The
adiabatic process starts from the ground state of the local Hamiltonian, which
is determined variationally.

The ground state of the local Hamiltonian is prepared with high fidelity using
short quantum circuits. These quantum circuits consists of a few ($\lessapprox
6$) two-qubit entangling layers ($CZ$ in our example) and parameterized
single-qubit gates. The circuit parameters are calculated on a classical
computer by optimizing the overlap of the circuit final state and the local
Hamiltonian ground state obtained from exact diagonalization. Since the local
Hamiltonian is a sum of uncoupled terms at each site, the circuit optimization
problem is reduced to the optimization of a circuit with a small number of
qubits ($n_{q}\approx 6\sim 8$) and is independent of the lattice size.

The local Hamiltonian contains the $\phi^{4}$ interaction and, for the
preparation of the broken symmetry phase, a coupling of the scalar field to an
external field. The parameters defining the local Hamiltonian, the initial
mass $m_{I}$, the initial interaction coupling $\lambda_{I}$ and the initial
external field $f_{I}$, constitute the input of our algorithm and can be
adjusted. The system’s Hamiltonian is restored by the adiabatic evolution. The
initial parameters can be optimized to decrease the adiabatic time necessary
to reach the full Hamiltonian ground state.

Our numerical investigation on small lattices finds a correlation between the
adiabatic time and the local overlap of the final wavefunction and the initial
wavefunction. To reduce the adiabatic time, the input parameters $m_{I}$,
$\lambda_{I}$ and $f_{I}$, should be chosen such that the final and the
initial states have a maximum local overlap. For state preparation in the
normal phase, we find a strong dependence of the adiabatic time on $m_{I}^{2}$
and a weak dependence on $\lambda_{I}$. In this case, the optimal $m_{I}^{2}$
is positive and larger than the $\phi^{4}$ Hamiltonian mass parameter
$\left|m_{0}\right|^{2}$. For state preparation in the broken symmetry phase,
we find that the optimal adiabatic time is achieved when the adiabatic process
starts from the ground state of a local Hamiltonian with significant external
field $f_{I}$ and a positive input mass parameter $m_{I}^{2}$.

The correlation between the adiabatic time and the local overlap allows us to
use the overlap as a tool to optimize the initial parameters. We propose an
iterative strategy for finding the optimal input parameters $m_{I}$,
$\lambda_{I}$ and $f_{I}$, starting with the optimal values for small lattices
and adjusting them in increasingly larger systems by maximizing the local
overlap.

There are two main challenges associated with the preparation of the broken
symmetry states that are addressed in this paper. The first is when the
adiabatic evolution crosses the critical phase transition region for an
initial state in the symmetric phase. The second is the vanishing gap of the
double degenerate broken symmetry phase. We avoid these challenges by coupling
the system to an external field during the adiabatic evolution. We propose an
adiabatic process consisting of two steps. In the first step, the adiabatic
evolution starts from a broken symmetry state prepared variationally by
coupling the local Hamiltonian to a finite external field. During this step,
the coupling term is switched on and the external field is decreased. During
the second step, the external field is decreased to vanishing values. The
error in the adiabatic process can be kept under control by choosing a linear
decrease of the external field in the first step and an exponential decrease
in the second step.

##  VI Acknowledgements

A. M. is partially supported by the DOE/HEP QuantISED program grant of the
theory consortium ”Intersections of QIS and Theoretical Particle Physics” at
Fermilab.” A.C.Y.L. and S.M. are partially supported by the DOE/HEP QuantISED
program grant ”HEP Machine Learning and Optimization Go Quantum”,
identification number 0000240323. This manuscript has been authored by Fermi
Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S.
Department of Energy, Office of Science, Office of High Energy Physics.

## References

  * Preskill [2018] J. Preskill, Quantum 2, 79 (2018).
  * Daley  _et al._ [2022] A. J. Daley, I. Bloch, C. Kokail, S. Flannigan, N. Pearson, M. Troyer, and P. Zoller, Nature 607, 667 (2022).
  * Reagor  _et al._ [2018] M. Reagor _et al._ , Sci. Adv. 4, eaao3603 (2018).
  * Arute  _et al._ [2019] F. Arute _et al._ , Nature 574, 505 (2019).
  * Jurcevic  _et al._ [2021] P. Jurcevic, A. Javadi-Abhari, L. S. Bishop, I. Lauer, D. F. Bogorin, M. Brink, L. Capelluto, O. Günlük, T. Itoko, N. Kanazawa, _et al._ , Quantum Sci. Technol. 6, 025020 (2021).
  * O’Malley  _et al._ [2016] P. J. J. O’Malley _et al._ , Phys. Rev. X 6, 031007 (2016).
  * Kandala  _et al._ [2017] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * Hempel  _et al._ [2018] C. Hempel, C. Maier, J. Romero, J. McClean, T. Monz, H. Shen, P. Jurcevic, B. P. Lanyon, P. Love, R. Babbush, A. Aspuru-Guzik, R. Blatt, and C. F. Roos, Phys. Rev. X 8, 031022 (2018).
  * Colless  _et al._ [2018] J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok, M. E. Kimchi-Schwartz, J. R. McClean, J. Carter, W. A. de Jong, and I. Siddiqi, Phys. Rev. X 8, 011021 (2018).
  * Kandala  _et al._ [2019] A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, Nature 567, 491 (2019).
  * Quantum  _et al._ [2020] G. A. Quantum _et al._ , Science 369, 1084 (2020).
  * Macridin  _et al._ [2018a] A. Macridin, P. Spentzouris, J. Amundson, and R. Harnik, Phys. Rev. Lett. 121, 110504 (2018a).
  * Macridin  _et al._ [2018b] A. Macridin, P. Spentzouris, J. Amundson, and R. Harnik, Phys. Rev. A 98, 042312 (2018b).
  * Gluza  _et al._ [2018] M. Gluza, M. Kliesch, J. Eisert, and L. Aolita, Phys. Rev. Lett. 120, 190501 (2018).
  * Ma  _et al._ [2020] H. Ma, M. Govoni, and G. Galli, npj Comput Mater 6, 1 (2020).
  * Uvarov  _et al._ [2020] A. Uvarov, J. D. Biamonte, and D. Yudin, Phys. Rev. B 102, 075104 (2020).
  * Ji  _et al._ [2020] W. Ji, L. Zhang, M. Wang, L. Zhang, Y. Guo, Z. Chai, X. Rong, F. Shi, X.-J. Liu, Y. Wang, and J. Du, Phys. Rev. Lett. 125, 020504 (2020).
  * Roggero  _et al._ [2020a] A. Roggero, A. C. Y. Li, J. Carlson, R. Gupta, and G. N. Perdue, Phys. Rev. D 101, 074038 (2020a).
  * Roggero  _et al._ [2020b] A. Roggero, C. Gu, A. Baroni, and T. Papenbrock, Phys. Rev. C 102, 064624 (2020b).
  * Jordan  _et al._ [2012] S. P. Jordan, K. S. M. Lee, and J. Preskill, Science 336, 1130 (2012).
  * Klco  _et al._ [2018] N. Klco, E. F. Dumitrescu, A. J. McCaskey, T. D. Morris, R. C. Pooser, M. Sanz, E. Solano, P. Lougovski, and M. J. Savage, Phys. Rev. A 98, 032331 (2018).
  * Pedernales  _et al._ [2018] J. S. Pedernales, M. Beau, S. M. Pittman, I. L. Egusquiza, L. Lamata, E. Solano, and A. del Campo, Phys. Rev. Lett. 120, 160403 (2018).
  * Hu  _et al._ [2019] J. Hu, L. Feng, Z. Zhang, and C. Chin, Nat. Phys. 15, 785 (2019).
  * Lamm  _et al._ [2019] H. Lamm, S. Lawrence, and Y. Yamauchi (NuQS Collaboration), Phys. Rev. D 100, 034518 (2019).
  * Farrelly and Streich [2020] T. Farrelly and J. Streich, Discretizing quantum field theories for quantum simulation (2020), arXiv:2002.02643 [quant-ph] .
  * Klco and Savage [2020] N. Klco and M. J. Savage, Phys. Rev. A 102, 012612 (2020).
  * Barata  _et al._ [2021] J. Barata, N. Mueller, A. Tarasov, and R. Venugopalan, Phys. Rev. A 103, 042410 (2021).
  * Aizenman [1981] M. Aizenman, Phys. Rev. Lett. 47, 1 (1981).
  * Fröhlich [1982] J. Fröhlich, Nucl. Phys. B 200, 281 (1982).
  * Chang [1976] S.-J. Chang, Phys. Rev. D 13, 2778 (1976).
  * Magruder [1976] S. F. Magruder, Phys. Rev. D 14, 1602 (1976).
  * Sugihara [2004] T. Sugihara, Journal of High Energy Physics 2004, 007 (2004).
  * Milsted  _et al._ [2013] A. Milsted, J. Haegeman, and T. J. Osborne, Phys. Rev. D 88, 085030 (2013).
  * Bañuls and Cichy [2020] M. C. Bañuls and K. Cichy, Rep. Prog. Phys 83, 024401 (2020).
  * Schaich and Loinaz [2009] D. Schaich and W. Loinaz, Phys. Rev. D 79, 056008 (2009).
  * Bosetti  _et al._ [2015] P. Bosetti, B. De Palma, and M. Guagnelli, Phys. Rev. D 92, 034509 (2015).
  * Lee  _et al._ [2001] D. Lee, N. Salwen, and D. Lee, Physics Letters B 503, 223 (2001).
  * Rychkov and Vitale [2015] S. Rychkov and L. G. Vitale, Phys. Rev. D 91, 085011 (2015).
  * Macridin  _et al._ [2022] A. Macridin, A. C. Y. Li, S. Mrenna, and P. Spentzouris, Phys. Rev. A 105, 052405 (2022).
  * Kitaev and Webb [2008] A. Kitaev and W. A. Webb, arXiv:0801.0342 (2008).
  * Liu  _et al._ [2021] J. Liu, Z. Li, H. Zheng, X. Yuan, and J. Sun, arXiv:2109.05547 10.48550/ARXIV.2109.05547 (2021).
  * McClean  _et al._ [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Nat. Commun. 9, 4812 (2018).
  * Jordan  _et al._ [2014] S. P. Jordan, K. S. M. Lee, and J. Preskill, Quant. Inf. Comput. 14, 1014 (2014), arXiv:1112.4833 [hep-th] .
  * Lüscher and Weisz [1987] M. Lüscher and P. Weisz, Nuclear Physics B 290, 25 (1987).
  * Freedman  _et al._ [1982] B. Freedman, P. Smolensky, and D. Weingarten, Physics Letters B 113, 481 (1982).
  * Fox and Halliday [1985] I. Fox and I. Halliday, Physics Letters B 159, 148 (1985).
  * Lang [1985] C. Lang, in  _Advances in lattice gauge theory_ (1985).
  * Drummond  _et al._ [1987] I. Drummond, S. Duane, and R. Horgan, Nuclear Physics B 280, 25 (1987).
  * Nielsen and Chuang [2002] M. A. Nielsen and I. Chuang, Quantum computation and quantum information (2002).
  * Trotter [1959] H. F. Trotter, Proc. Amer. Math. Soc. 10, 545 (1959).
  * Suzuki [1976] M. Suzuki, Communications in Mathematical Physics 51, 183 (1976).
  * Lloyd [1996] S. Lloyd, Science 273, 1073 (1996).
  * Welch  _et al._ [2014] J. Welch, D. Greenbaum, S. Mostame, and A. Aspuru-Guzik, New J. Phys. 16, 033040 (2014).
  * Developers [2021] C. Developers, Cirq (2021), See full list of authors on Github: https://github .com/quantumlib/Cirq/graphs/contributors.
  * ANIS  _et al._ [2021] M. S. ANIS _et al._ , Qiskit: An open-source framework for quantum computing (2021).
  * Hansen and Ostermeier [2001] N. Hansen and A. Ostermeier, Evol. Comput. 9, 159 (2001).
  * Born and Fock [1928] M. Born and V. Fock, Zeitschrift für Physik 51, 165 (1928).
  * Albash and Lidar [2018] T. Albash and D. A. Lidar, Rev. Mod. Phys. 90, 015002 (2018).
  * Amin [2009] M. H. S. Amin, Phys. Rev. Lett. 102, 220401 (2009).
  * Tong  _et al._ [2007] D. M. Tong, K. Singh, L. C. Kwek, and C. H. Oh, Phys. Rev. Lett. 98, 150402 (2007).
  * Boixo and Somma [2010] S. Boixo and R. D. Somma, Phys. Rev. A 81, 032308 (2010).
  * Cao  _et al._ [2013] H. Cao, Z. Guo, Z. Chen, and W. Wang, Science China Physics, Mechanics and Astronomy 56, 1401 (2013).
  * Buhrman  _et al._ [2001] H. Buhrman, R. Cleve, J. Watrous, and R. de Wolf, Phys. Rev. Lett. 87, 167902 (2001).
  * Shavitt and Redmon [1980] I. Shavitt and L. T. Redmon, The Journal of Chemical Physics 73, 5711 (1980).

##  Appendix A Fourier transform gate

In this appendix, we prove Section III.2, which expresses the discrete Fourier
transform gate ${\cal{F}}_{j}$ in terms of the QFT gate and single qubit
rotations. We begin from the definition of the discrete Fourier operator in
Eq. 13:

$\displaystyle{\cal{F}}_{j}=\frac{1}{\sqrt{N_{\varphi}}}\sum_{\alpha,\beta=0}^{N_{\varphi}-1}e^{i\frac{2\pi}{N_{\varphi}}\left(\alpha-\frac{N_{\varphi}-1}{2}\right)\left(\beta-\frac{N_{\varphi}-1}{2}\right)}\left|\varphi_{\alpha}\right\rangle_{j}\left\langle\varphi_{\beta}\right|_{j}.$
(71)

Expanding the phase factor and inserting the identity operator
$\openone=\sum_{\alpha=0}^{N_{\varphi}-1}\left|\varphi_{\alpha}\right\rangle_{j}\left\langle\varphi_{\alpha}\right|_{j}$
leads to

$\displaystyle{\cal{F}}_{j}=$
$\displaystyle\,e^{i\frac{N_{\varphi}\delta^{2}}{2\pi}}$
$\displaystyle\left(\sum_{\alpha=0}^{N_{\varphi}-1}e^{-i\delta\alpha}\left|\varphi_{\alpha}\right\rangle_{j}\left\langle\varphi_{\alpha}\right|_{j}\right)$
$\displaystyle\left(\frac{1}{\sqrt{N_{\varphi}}}\sum_{\mu,\nu=0}^{N_{\varphi}-1}e^{i\frac{2\pi}{N_{\varphi}}\mu\nu}\left|\varphi_{\mu}\right\rangle_{j}\left\langle\varphi_{\nu}\right|_{j}\right)$
$\displaystyle\left(\sum_{\beta=0}^{N_{\varphi}-1}e^{-i\delta\beta}\left|\varphi_{\beta}\right\rangle_{j}\left\langle\varphi_{\beta}\right|_{j}\right),$
(72)

where $\delta=\tfrac{(N_{\varphi}-1)\pi}{N_{\varphi}}$. The first line is a
phase factor, and it is relevant if we would like to implement a control
Fourier transform gate. The third line is a standard QFT gate [49]. The second
and the fourth lines can be implemented as single-qubit $z$ rotation gates
similar to Eq. 28 such that

$\displaystyle
e^{-i\delta\alpha}\left|\varphi_{\alpha}\right\rangle_{j}\left\langle\varphi_{\alpha}\right|_{j}$
$\displaystyle=$ $\displaystyle
e^{-i\delta\sum_{q=0}^{n_{q}-1}\alpha_{qj}2^{n_{q}-1-q}}\left|\varphi_{\alpha}\right\rangle_{j}\left\langle\varphi_{\alpha}\right|_{j}$
$\displaystyle=$ $\displaystyle
e^{-i\delta\sum_{q=0}^{n_{q}-1}\frac{2\alpha_{qj}-1}{2}2^{n_{q}-1-q}}e^{-i\delta\frac{N_{\varphi}-1}{2}}\left|\varphi_{\alpha}\right\rangle_{j}\left\langle\varphi_{\alpha}\right|_{j}$
$\displaystyle=$ $\displaystyle
e^{-i\frac{N_{\varphi}\delta^{2}}{2\pi}}\prod_{q=0}^{n_{q}-1}e^{-i\delta\frac{\sigma^{z}_{qj}}{2}2^{n_{q}-1-q}}$
$\displaystyle=$ $\displaystyle
e^{-i\delta^{2}\frac{N_{\varphi}\delta^{2}}{2\pi}}\prod_{q=0}^{n_{q}-1}R^{z}_{qj}\left(2^{n_{q}-1-q}\delta\right).$
(73)

In the second line of the above equation, we used the binary representation
$\alpha_{j}=\sum_{q=0}^{n_{q}-1}\alpha_{qj}2^{n_{q}-1-q}$ [Eq. 23]. Rewriting
Appendix A with QFT and $R^{z}$ gates, we arrive at Section III.2

$\displaystyle{\cal{F}}_{j}=$ $\displaystyle
e^{-i\frac{N_{\varphi}\delta^{2}}{2\pi}}\prod_{q=0}^{n_{q}-1}R^{z}_{qj}\left(2^{n_{q}-1-q}\delta\right)~{}\text{QFT}_{j}$
$\displaystyle\times\prod_{q=0}^{n_{q}-1}R^{z}_{qj}\left(2^{n_{q}-1-q}\delta\right).$
(74)

##  Appendix B Semi-classical derivation of the relation between the optimal
initial mass and the optimal initial external field

Here we present a semi-classical argument to explain the nearly linear
relation between the best $m_{I}^{2}$ and $f_{M}$ shown in Figs. 6 and 8. The
adiabatic preparation is optimized when the local overlap between the initial
ground state of $H(s=0)$ and the final ground state of $H(s=1)$ is close to
maximum. While we do not have an analytic form of the ground-state
wavefunction to determine the local overlap, we can determine the minimum of
the classical potential energy, which in our case approximately predicts the
center of the wavefunction. Hence, we find that, when the minima of the
initial potential energy and the final potential coincide, the adiabatic
protocol is close to optimal. The potential energy along the adiabatic path is
given by

$\displaystyle V(s)=\sum_{j=1}^{N}\bigg{[}$
$\displaystyle\frac{m_{I}^{2}+s\,\delta
m^{2}}{2}\Phi_{j}^{2}+\frac{s}{2}\sum_{e=0}^{d}\left(\Phi_{j+e}-\Phi_{j}\right)^{2}$
$\displaystyle+\frac{\lambda_{I}+s\,\delta\lambda}{4!}\Phi_{j}^{4}+\left(f_{I}+s\,\delta
f\right)\Phi_{j}\bigg{]}.$ (75)

Given $m_{0}^{2}<0$ and assuming translational symmetry, the minimum of the
potential energy is determined by

$\left.\frac{\partial{V(s)}}{\partial\Phi_{j}}\right|_{\Phi_{j}=\Phi_{0}(s)}=0,$
(76)

where $\Phi_{0}=\Phi_{1}=\cdots=\Phi_{N=1}=\Phi_{0}(s)$ is the location of the
minimum. By requiring the minimum to be at the same position in the beginning
(s=0) and at the end (s=1), i.e. $\Phi_{0}(0)=\Phi_{0}(1)$, we can find a
relation between $m_{I}^{2}$ and $f_{I}$ as shown in Figs. 6 and 8 as black
dashed lines. In particular the following linear relation,

$f_{I}=\mp\sqrt{\frac{3!\left|m_{0}^{2}\right|}{\lambda_{0}}}\left(m_{I}^{2}+\frac{\lambda_{I}}{\lambda_{0}}\left|m_{0}^{2}\right|\right)+\mathcal{O}(f_{M}),$
(77)

is a good approximation. We emphasize that while this approximated relation
aligns with the observations we made in the numerical simulation, this does
not mean that the system properties in the low-energy subspace can be
explained by a semi-classical argument. At best, the semi-classical argument
gives us an insight about where the wavefunction is confined by the semi-
classical potential in the field-amplitude basis. The distribution and most of
the properties of the wavefunction remain strongly influenced by quantum
effects and cannot be inferred from the semi-classical potential.

##  Appendix C Adiabatic preparation of a broken symmetry state in an
effective two-level system

In this appendix, we investigate the adiabatic condition for the preparation
of the broken-symmetry states near the end of the path in section IV.2.2.
Using an effective two-level system to approximate a nearly degenerate
subspace, we will derive the adiabatic condition and propose more efficient
adiabatic paths.

We consider a finite system of size $L$ in the parameter regime corresponding
to the broken symmetry phase. The two lowest energy states
$\left|0\right\rangle$ and $\left|1\right\rangle$ are nearly degenerate and
well separated from the rest of the spectrum,

$\displaystyle\Delta_{0}\equiv E_{1}-E_{0}\ll E_{2}-E_{1}\equiv\Delta_{1}.$
(78)

In the thermodynamic limit, the states $\left|0\right\rangle$ and
$\left|1\right\rangle$ are degenerate, i.e.

$\displaystyle\Delta_{0}(L)\xrightarrow[L\rightarrow\infty]{}0.$ (79)

The broken symmetry state can be reached by coupling the system to an external
field

$\displaystyle H_{f}$ $\displaystyle=H+f\Phi,~{}~{}\text{~{}~{}where~{}~{}}$
(80) $\displaystyle\Phi$ $\displaystyle=\sum_{j=1}^{N}\Phi_{j},$ (81)

and then taking the limits $L\rightarrow\infty$ followed by $f\rightarrow 0$.
Perturbation theory can be applied in this limit. We then apply a Schrieffer-
Wolff transformation [64] to obtain an effective two-level Hamiltonian.

###  C.1 Review of Schrieffer-Wolff transformation

In this subsection, we will provide a quick review of the Schrieffer-Wolff
transformation. In general, we can separate the full Hamiltonian into the
unperturbed part $H$ and the leading-order perturbation $\Phi$ such that

$H_{f}=H+f\Phi.$ (82)

The unperturbed part $H$ is two-fold degenerate due to a double-well
potential. We consider the nearly degenerate subspace to be spanned by the
states $|0\rangle$ and $|1\rangle$, each localized at one of the potential
wells with an energy $E_{0}(f=0)$.

We carry out the perturbative calculation in this subspace using a Schrieffer-
Wolff transformation, which constructs an effective Hamiltonian
$H_{\mathrm{eff}}=e^{-S}H_{f}e^{S}$ decoupling the two-fold degenerate
subspace from the higher-energy eigenstates up to an arbitrary order of $f$.
The generator $S$ can be constructed iteratively using a canonical van Vleck
formalism [64].

To derive our result, we separate $\Phi$ into a block diagonal piece $V_{D}$
and a piece $V_{X}$ that couples the degenerate subspace and the higher-energy
subspace. The generator $S$ has a series expansion

$S=\sum_{j=1}^{\infty}f^{j}S^{(j)}.$ (83)

Applying the Baker-Hausdorff lemma, we can expand the transformed Hamiltonian
into

$\displaystyle H_{\mathrm{eff}}=$
$\displaystyle\,e^{-S}\left(H+fV_{D}\right)e^{S}+e^{-S}fV_{X}e^{S}$
$\displaystyle=$ $\displaystyle\,H^{\prime}+[H+fV_{D},S]+\frac{1}{2}[[H,S],S]$
$\displaystyle+fV_{X}+f[V_{X},S]+\mathcal{O}(f^{3})$ $\displaystyle=$
$\displaystyle\,H+f\left(V_{D}+V_{X}+[H,S^{(1)}]\right)$
$\displaystyle+f^{2}\Big{(}\frac{1}{2}[[H,S^{(1)}],S^{(1)}]+[V_{X},S^{(1)}]+[H,S^{(2)}]$
$\displaystyle\ \ \ \ \ \ \ \ +[V_{d},S^{(1)}]\Big{)}+\mathcal{O}(f^{3}).$
(84)

To decouple the two subspaces, we pick

$\displaystyle[H,S^{(1)}]$ $\displaystyle=-V_{X},$ (85)
$\displaystyle[H,S^{(2)}]$ $\displaystyle=-[V_{d},S^{(1)}].$ (86)

We can use the resolvent operator technique to determine the generator such
that

$\displaystyle S^{(1)}\left|m\right\rangle$
$\displaystyle=-R_{m}[H,S^{(1)}]\left|m\right\rangle,$ (87)

where the resolvent operator is given by

$\displaystyle R_{m}$ $\displaystyle=\sum_{n\neq
m}\frac{|n\vphantom{n}\rangle\langle n\vphantom{n}|}{E_{m}-E_{n}}.$ (88)

This gives

$S^{(1)}=\sum_{k=0,1}\sum_{\gamma\geq
2}\left(\frac{\left<\gamma\left|V_{X}\right|k\right>}{E_{k}-E_{\gamma}}|\gamma\vphantom{k}\rangle\langle
k\vphantom{\gamma}|-\frac{\left<k\left|V_{X}\right|\gamma\right>}{E_{k}-E_{\gamma}}|k\vphantom{\gamma}\rangle\langle\gamma\vphantom{k}|\right).$

The effective Hamiltonian up to the second order is given by

$H_{\mathrm{eff}}=H+fV_{D}+f^{2}\frac{1}{2}[V_{X},S^{(1)}].$ (89)

In the nearly degenerate subspace, $H_{\mathrm{eff}}$ can be written as

$\displaystyle H_{\mathrm{eff}}=$
$\displaystyle\sum_{j=0,1}E_{j}|j\vphantom{j}\rangle\langle
j\vphantom{j}|+f\sum_{j,k=0,1}V_{jk}|j\vphantom{k}\rangle\langle
k\vphantom{j}|$
$\displaystyle+f^{2}\sum_{j,k=0,1}W_{jk}|j\vphantom{k}\rangle\langle
k\vphantom{j}|+\mathcal{O}(f^{3}),$ (90)

where $W_{jk}=\sum_{\gamma\geq 2}\frac{\langle
j|\Phi|\gamma\rangle\langle\gamma|\Phi|k\rangle}{2}\left(\frac{1}{E_{\gamma}-E_{0}}+\frac{1}{E_{\gamma}-E_{1}}\right)$
and $V_{jk}=\langle j|\Phi|k\rangle$. The first two terms in the equation are
simply the projection of $H+f\Phi$ in the nearly degenerate subspace. The last
term comes from the perturbative treatment and can be understood as the
virtual interaction between the two nearly degenerate states through the
higher-energy states.

###  C.2 Effective two-level model

Ignoring the higher-order terms in Section C.1, we get

$\displaystyle H_{\mathrm{eff}}$
$\displaystyle=\sum_{j=0,1}E_{j}\left|j\right\rangle\left\langle
j\right|+f\sum_{j,k=0,1}V_{jk}\left|j\right\rangle\left\langle k\right|$
$\displaystyle+f^{2}\sum_{j,k=0,1}W_{jk}\left|j\right\rangle\left\langle
k\right|.$ (91)

Note that $V_{00}=0$, and $V_{11}=0$ since $\left|0\right\rangle$ and
$\left|1\right\rangle$ have the full symmetry of the Hamiltonian $H$, while
$\Phi$ breaks the $Z_{2}$ symmetry. For convenience we will denote

$\displaystyle v=\left|V_{01}\right|=\left|V_{10}\right|.$ (92)

The implicit assumption made when applying the perturbation theory is that $f$
is small such that $fv\ll\Delta_{1}$. We also assume that the system is large
such that $\Delta_{0}\ll fv$. We ignore terms of ${\cal{O}}(\Delta_{0}/(fv))$
in the following.

The effective Hamiltonian (C.2) acts on the $2$-dimensional space spanned by
the states $\left|0\right\rangle$ and $\left|1\right\rangle$. The following
notation is convenient:

$\displaystyle s$
$\displaystyle=\frac{1}{2}\left[E_{0}+E_{1}+f^{2}\left(W_{00}+W_{11}\right)\right],$
(93) $\displaystyle\delta$
$\displaystyle=\frac{1}{2}\left[\Delta_{0}+f^{2}\left(W_{11}-W_{00}\right)\right],$
(94) $\displaystyle t$ $\displaystyle=fV_{01}+f^{2}W_{01},$ (95)
$\displaystyle D$ $\displaystyle=\frac{\delta}{\sqrt{\delta^{2}+|t|^{2}}},$
(96) $\displaystyle t$ $\displaystyle=\left|t\right|e^{-i2f},$ (97)
$\displaystyle\tan 2\theta$ $\displaystyle=-\frac{\left|t\right|}{\delta},$
(98) $\displaystyle\cos^{2}\theta$
$\displaystyle=\frac{1}{2}\left(1+D\right),~{}~{}\cos\theta=\sqrt{\frac{1+D}{2}},$
(99) $\displaystyle\sin^{2}\theta$
$\displaystyle=\frac{1}{2}\left(1-D\right),~{}~{}\sin\theta=-\sqrt{\frac{1-D}{2}}.$
(100)

The eigenstates of the Hamiltonian (C.2) can be written as

$\displaystyle\left|\psi_{a}\left(f\right)\right\rangle$
$\displaystyle=e^{-if}\cos\theta\left|0\right\rangle+e^{if}\sin\theta\left|1\right\rangle$
(101) $\displaystyle\left|\psi_{b}\left(f\right)\right\rangle$
$\displaystyle=-e^{-if}\sin\theta\left|0\right\rangle+e^{if}\cos\theta\left|1\right\rangle,$
(102)

while the corresponding energies are

$\displaystyle{\cal{E}}_{a}$ $\displaystyle=s-\sqrt{\delta^{2}+|t|^{2}}$ (103)
$\displaystyle{\cal{E}}_{b}$ $\displaystyle=s+\sqrt{\delta^{2}+|t|^{2}}.$
(104)

Defining the broken symmetry states as

$\displaystyle\left|\psi_{a0}\right\rangle$ $\displaystyle=\lim_{f\rightarrow
0}\left|\psi_{a}\left(f\right)\right\rangle=\frac{1}{\sqrt{2}}\left(e^{-if}\left|0\right\rangle-e^{if}\left|1\right\rangle\right),$
(105) $\displaystyle\left|\psi_{b0}\right\rangle$
$\displaystyle=\lim_{f\rightarrow
0}\left|\psi_{b}\left(f\right)\right\rangle=\frac{1}{\sqrt{2}}\left(e^{-if}\left|0\right\rangle+e^{if}\left|1\right\rangle\right),$
(106)

the eigenstates of the system can be written as

$\displaystyle\left|\psi_{a}(f)\right\rangle$
$\displaystyle=\frac{\left(\cos\theta-\sin\theta\right)}{\sqrt{2}}\left|\psi_{a0}\right\rangle+\frac{\left(\cos\theta+\sin\theta\right)}{\sqrt{2}}\left|\psi_{b0}\right\rangle$
(107) $\displaystyle\left|\psi_{b}(f)\right\rangle$
$\displaystyle=-\frac{\left(\cos\theta+\sin\theta\right)}{\sqrt{2}}\left|\psi_{a0}\right\rangle+\frac{\left(\cos\theta-\sin\theta\right)}{\sqrt{2}}\left|\psi_{b0}\right\rangle.$
(108)

The quantity $D$ in Eq. 96 can be written as

$\displaystyle D$
$\displaystyle=fB-f^{2}AB+{\cal{O}}(f^{3})+{\cal{O}}(\frac{\Delta_{0}}{fv}),$
(109)

where

$\displaystyle A$
$\displaystyle=\frac{\left(V_{01}W_{10}+V_{10}W_{01}\right)}{2v^{2}},$ (110)
$\displaystyle B$ $\displaystyle=\frac{\left(W_{11}-W_{00}\right)}{2v},$ (111)

are independent on the external field $f$. The following quantities, relevant
for the investigation of the adiabatic process, can be written up to
${\cal{O}}(f^{3})$ and ${\cal{O}}(\Delta_{0}/(fv))$ as

$\displaystyle\Delta_{ba}$
$\displaystyle={\cal{E}}_{b}-{\cal{E}}_{a}=2fv\left(1+fA\right)$ (112)
$\displaystyle\langle\psi_{b}\left(f\right)\vphantom{\psi_{a}\left(f\right)}|\frac{d}{df}|\psi_{a}\left(f\right)\vphantom{\psi_{b}\left(f\right)}\rangle$
$\displaystyle=-\sin\theta\frac{d}{df}\cos\theta+\cos\theta\frac{d}{df}\sin\theta$
(113) $\displaystyle=\frac{1}{2}B-fAB+\frac{1}{4}f^{2}B^{3}+{\cal{O}}(f^{3})$
$\displaystyle\langle\psi_{a}\left(f\right)\vphantom{\psi_{b}\left(f\right)}|\frac{d}{df}|\psi_{b}\left(f\right)\vphantom{\psi_{a}\left(f\right)}\rangle$
$\displaystyle=-\langle\psi_{b}\left(f\right)\vphantom{\psi_{a}\left(f\right)}|\frac{d}{df}|\psi_{a}\left(f\right)\vphantom{\psi_{b}\left(f\right)}\rangle$
(114)
$\displaystyle\left|\left|\left|\psi_{a}\left(f\right)\right\rangle-\left|\psi_{a0}\right\rangle\right|\right|$
$\displaystyle=\frac{1}{2}D=\frac{1}{2}fB-\frac{1}{2}f^{2}AB.$ (115)

###  C.3 Adiabatic evolution in a two-level system

Here, we investigate the two-level system’s evolution when the external field
$f(s)$, depending on the variable $s=t/T$, changes during a time interval $T$
from the initial value $f(s=0)=f_{i}$ to the final value $f(s=1)=f_{f}$. Note
that this evolution corresponds to the second stage of the adiabatic process
to prepare the broken symmetry states in Section IV.2.2 such that
$f_{i}=f_{M}$ and $f_{f}=f_{F}$. The initial wavefunction is the ground state
of the system when $f=f_{i}$, i.e.
$\left|\psi_{ini}\right\rangle\equiv\left|\psi(s=0)\right\rangle=\left|\psi_{a}(f_{i})\right\rangle$.
During the adiabatic evolution, the wavefunction is

$\displaystyle\left|\psi(s)\right\rangle=c_{a}(s)\left|\psi_{a}\left[f(s)\right]\right\rangle+c_{b}(s)\left|\psi_{b}\left[f(s)\right]\right\rangle.$
(116)

where $\left|\psi_{a}(f)\right\rangle$ and $\left|\psi_{b}(f)\right\rangle$
are the instantaneous eigenstates of the Hamiltonian (C.2) with external field
$f(s)$.

The state at the end of the adiabatic evolution is

$\displaystyle\left|\psi_{f}\right\rangle\equiv\left|\psi(s=1)\right\rangle=c_{af}\left|\psi_{a}(f_{f})\right\rangle+c_{bf}\left|\psi_{b}(f_{f})\right\rangle$
(117)

where we denote $c_{af}\equiv c_{a}(1)$ and $c_{bf}\equiv c_{b}(1)$. The
adiabatic error is given by the difference between the system’s state at the
end of the adiabatic evolution and the ground state when $f=f_{f}$,

$\displaystyle\epsilon_{ad}$
$\displaystyle=\left|\left|\left|\psi_{f}\right\rangle-\left|\psi_{a}(f_{f})\right\rangle\right|\right|=\sqrt{2\left(1-\left|c_{af}\right|\right)}$
(118)
$\displaystyle\approx\left|c_{bf}\right|+{\cal{O}}(\left|c_{bf}\right|^{4})$

The error caused by the finite final field, defined as the difference between
the eigenstate when $f=f_{f}$ and the broken symmetry state, is obtained by
applying Eq. 115:

$\displaystyle\epsilon_{f}\equiv\left|\left|\left|\psi_{a}(f_{f})\right\rangle-\left|\psi_{a0}\right\rangle\right|\right|=\frac{1}{2}Bf_{f}+{\cal{O}}(f^{2}).$
(119)

Employing Eqs. 107, 108 and 116, the adiabatically prepared state can be
written as

$\displaystyle\left|\psi_{f}\right\rangle$
$\displaystyle=\frac{\cos\theta_{f}\left(c_{af}-c_{bf}\right)-\sin\theta_{f}\left(c_{af}+c_{bf}\right)}{\sqrt{2}}\left|\psi_{a0}\right\rangle$
(120)
$\displaystyle+\frac{\cos\theta_{f}\left(c_{af}+c_{bf}\right)+\sin\theta_{f}\left(c_{af}-c_{bf}\right)}{\sqrt{2}}\left|\psi_{b0}\right\rangle,$

where $\cos\theta_{f}$ and $\sin\theta_{f}$ are given by Eqs. 99 and 100 when
$f=f_{f}$. The total error for preparing the broken symmetry state is

$\displaystyle\epsilon\equiv\left|\left|\left|\psi_{f}-\left|\psi_{a0}\right\rangle\right\rangle\right|\right|\approx\left|c_{af}\epsilon_{f}+c_{bf}\right|.$
(121)

The coefficients $c_{a}(s)$ and $c_{b}(s)$, which describe the adiabatic
evolution, obey the differential equations

$\displaystyle\frac{dc_{a}}{ds}$
$\displaystyle=-c_{b}(s)\frac{df}{ds}\langle\psi_{a}\left(f\right)\vphantom{\psi_{b}\left(f\right)}|\frac{d}{df}|\psi_{b}\left(f\right)\vphantom{\psi_{a}\left(f\right)}\rangle
e^{-iT\int_{0}^{s}\Delta_{ba}(u)du}$ (122) $\displaystyle\frac{dc_{b}}{ds}$
$\displaystyle=-c_{a}(s)\frac{df}{ds}\langle\psi_{b}\left(f\right)\vphantom{\psi_{a}\left(f\right)}|\frac{d}{df}|\psi_{a}\left(f\right)\vphantom{\psi_{b}\left(f\right)}\rangle
e^{iT\int_{0}^{s}\Delta_{ba}(u)du},$ (123)

with the initial conditions $c_{a}(0)=1$ and $c_{b}(0)=0$. These equations can
be solved numerically. Next we will present results for different choices of
the time dependence of $f$: (1) linear and (2) exponential.

####  C.3.1 Adiabatic evolution with linear dependence of the external field

Figure 9: Errors for preparing the broken symmetry state when the external
field has a linear time dependence in a two-level system with $B=0.2$, $v=10$
and $f_{i}=1$. a) Logarithmic scale. Errors $\epsilon$, $\epsilon_{f}$ and
$\epsilon_{ad}$ versus the adiabatic time $T$ when
$\epsilon_{f}=f_{f}B/2=10^{-3}$. The constants $a=B/4\sqrt{\pi\delta_{f}/v}$
and $b=\delta_{f}B^{2}/(8v\epsilon_{f})$, see Eq. (138) and, respectively, Eq.
(141). Inset: The error $\epsilon$ (red line) has an oscillatory component as
a function of $T$. The points of maximum are well approximated by
$\epsilon_{f}+\epsilon_{ad}$ (black line). b) The total error $\epsilon$
versus $T$ for different choices of the final external field
$f_{f}=2\epsilon_{f}/B$. The adiabatic time for a fixed $\epsilon$ decreases
with decreasing $\epsilon_{f}$. The optimal adiabatic time is obtained in the
parameter region where $\epsilon_{f}\ll\epsilon\approx\epsilon_{ad}$.

In this scenario, the field’s time dependence is

$\displaystyle f(s)=f_{i}+(f_{f}-f_{i})s=f_{i}-\delta_{f}s$ (124)

where $\delta_{f}=f_{i}-f_{f}$. Eqs. 122 and 123 reduce to

$\displaystyle\frac{dc_{a}}{ds}$
$\displaystyle=-\frac{B}{2}\delta_{f}e^{-iTv\left(2f_{i}s-\delta_{f}s^{2}\right)}c_{b}(s)$
(125) $\displaystyle\frac{dc_{b}}{ds}$
$\displaystyle=\frac{B}{2}\delta_{f}e^{iTv\left(2f_{i}s-\delta_{f}s^{2}\right)}c_{a}(s).$
(126)

At the end of this section we will present the numerical solution to these
equations. However, in order to have an analytical estimate of the error
dependence on the adiabatic time we proceed first by considering some
simplifying approximations. We will judge the accuracy of these approximations
by comparing the approximate results with the exact numerical solution.

First, taking into account that $\left|1-c_{a}(s)\right|\approx\epsilon_{ad}$,
Eq. 126 can be written as

$\displaystyle c_{b}(1)$
$\displaystyle=\frac{B}{2}\delta_{f}\int_{0}^{1}dse^{i2Tv\int_{0}^{s}\left(f_{i}-\delta_{f}u\right)du}+{\cal{O}}(\epsilon_{ad}^{2})$
(127)
$\displaystyle=\frac{B}{2}\delta_{f}\int_{0}^{1}dse^{-iTv\delta_{f}\left(-2\frac{f_{i}}{\delta_{f}}s+s^{2}\right)}+{\cal{O}}(\epsilon_{ad}^{2})$
$\displaystyle=\frac{B}{2}\delta_{f}e^{iTv\frac{f_{i}^{2}}{\delta_{f}}}\int_{0}^{1}dse^{-iTv\delta_{f}\left(s-\frac{f_{i}}{\delta_{f}}\right)^{2}}+{\cal{O}}(\epsilon_{ad}^{2})$
$\displaystyle=\frac{B}{2}\delta_{f}e^{iTv\frac{f_{i}^{2}}{\delta_{f}}}\frac{1}{\sqrt{Tv\delta_{f}}}\int_{\sqrt{Tv\delta_{f}}\frac{f_{f}}{\delta_{f}}}^{\sqrt{Tv\delta_{f}}\frac{f_{i}}{\delta_{f}}}dye^{-iy^{2}}+{\cal{O}}(\epsilon_{ad}^{2})$

The integral in Eq. (127) can be expressed in terms of the error function
$\mathrm{Erf}(x)$, as

$\displaystyle
c_{b}(1)=e^{i\frac{3}{4}\pi}\frac{B}{4}\sqrt{\frac{\pi\delta_{f}}{Tv}}e^{iTv\frac{f_{i}^{2}}{\delta_{f}}}\left[\mathrm{Erf}\left(e^{i\frac{\pi}{4}}\sqrt{\frac{Tv}{\delta_{f}}}f_{i}\right)-\mathrm{Erf}\left(e^{i\frac{\pi}{4}}\sqrt{\frac{Tv}{\delta_{f}}}f_{f}\right)\right]+{\cal{O}}(\epsilon^{2}).$
(128)

The error function expansion at small argument is

$\displaystyle\mathrm{Erf}(e^{i\frac{\pi}{4}}x)=\sqrt{\frac{2}{\pi}}\left(1+i\right)x-\frac{2}{\sqrt{\pi}}e^{i\frac{3\pi}{4}}\frac{x^{3}}{3}+{\cal{O}}(x^{5})$
(129)

while at large argument is

$\displaystyle\mathrm{Erf}(e^{i\frac{\pi}{4}}x)=1-e^{-ix^{2}}\frac{1-i}{\sqrt{2\pi}}\frac{1}{x}+{\cal{O}}(x^{-3}).$
(130)

We distinguish two cases.

##### Case $I$.

The adiabatic time satisfies

$\frac{\delta_{f}}{vf_{i}^{2}}\ll T\ll\frac{\delta_{f}}{vf_{f}^{2}}.$ (131)

In this case, the first error function in Eq. 128 has a large argument while
the second one has a small argument. The adiabatic error is

$\displaystyle\epsilon_{ad}$ $\displaystyle=\left|c_{b}(1)\right|$ (132)
$\displaystyle\approx
B\frac{\sqrt{\pi}}{4}\sqrt{\frac{\delta_{f}}{Tv}}\left(1-\sqrt{\frac{2}{\pi}}\sqrt{\frac{Tv}{\delta_{f}}}f_{f}\right)+{\cal{O}}(\frac{1}{T})$
$\displaystyle\approx
B\frac{\sqrt{\pi}}{4}\sqrt{\frac{\delta_{f}}{Tv}}-B\frac{f_{f}}{2\sqrt{2}}.$

This implies

$\displaystyle\left(\epsilon_{ad}+\frac{1}{\sqrt{2}}\epsilon_{f}\right)^{2}=\frac{\pi
B^{2}\delta_{f}}{16Tv}$ (133)

and an adiabatic time scaling as

$\displaystyle T\approx\frac{\pi
B^{2}\delta_{f}}{16v}\frac{1}{\left(\epsilon_{ad}+\frac{1}{\sqrt{2}}\epsilon_{f}\right)^{2}}.$
(134)

Equation 131 implies that this approximation is valid when

$\displaystyle\frac{1}{4\left(\epsilon_{ad}+\frac{1}{\sqrt{2}}\epsilon_{f}\right)^{2}}\ll\frac{1}{\pi\epsilon_{f}^{2}},$
(135)

or equivalently when

$\displaystyle\epsilon_{f}\ll\frac{2}{\sqrt{\pi}-\sqrt{2}}\epsilon_{ad}\approx
5.58\epsilon_{ad},$ (136)

i.e. when the final external field is chosen small enough that

$\displaystyle\epsilon\approx\epsilon_{ad}\gg\epsilon_{f}.$ (137)

In this case, the required adiabatic time scales inversely proportionally with
the squared accuracy,

$\displaystyle T\approx\frac{\pi B^{2}\delta_{f}}{16v}\frac{1}{\epsilon^{2}}.$
(138)

##### Case $II$.

The adiabatic time satisfies

$T\gg\frac{\delta_{f}}{vf_{f}^{2}}=\frac{\delta_{f}B^{2}}{4\epsilon_{f}^{2}}.$
(139)

In this case, both error functions appearing in Eq. (128) have a large
argument. The adiabatic error can be approximated by

$\displaystyle\epsilon_{ad}=\left|c_{b}(1)\right|\approx\frac{\delta_{f}B}{4Tvf_{f}}=\frac{\delta_{f}B^{2}}{8Tv\epsilon_{f}},$
(140)

which implies

$\displaystyle T\approx\frac{\delta_{f}B^{2}}{8v\epsilon_{f}\epsilon_{ad}}.$
(141)

Equation 139 implies that this approximation is valid when

$\displaystyle\epsilon_{f}\gg\epsilon_{ad}.$ (142)

Making the assumption (numerically verified) that both $\epsilon_{f}$ and
$\epsilon_{ad}$ are independent and the total error is the sum of these two
contributions, i.e. $\epsilon=\epsilon_{f}+\epsilon_{ad}$, one gets

$\displaystyle T>\frac{\delta_{f}B^{2}}{2v\epsilon^{2}}.$ (143)

By comparing the adiabatic time estimates for the two cases, Eqs. 138 and 141,
one concludes that, for a desired accuracy $\epsilon$, the choice of a large
final external field (when $\epsilon\approx\epsilon_{f}\gg\epsilon_{ad}$)
requires a larger $T$ than the choice of a small final external field (when
$\epsilon\approx\epsilon_{ad}\gg\epsilon_{f}$).

The errors calculated by solving Eqs. 125 and 126 numerically are shown in
Fig. 9. Since the right hand side of these differential equations contains
imaginary terms, the coefficients $c_{af}(T)$ and $c_{bf}(T)$ have an
oscillatory component. As a consequence, $\epsilon$ and $\epsilon_{ad}$
display an oscillatory behavior, as can be seen from the inset. For our
analysis of the numerical data, we consider the points where $\epsilon(T)$
reaches a local maxima. As shown in the inset,
$\epsilon\approx\epsilon_{f}+\epsilon_{ad}$ at these points, i.e. the
adiabatic error and the final field error can be considered as independent
contributions to the total error. In Fig. 9 (a), we show the errors
$\epsilon$, $\epsilon_{f}$ and $\epsilon_{ad}$ dependence on $T$ for a case
where the final external field yields $\epsilon_{f}=10^{-3}$. The numerical
results are in agreement with the analytical analysis discussed earlier. At
small $T$, where $\epsilon_{ad}\gg\epsilon_{f}$ the adiabatic time scales as
$T\propto\epsilon^{-2}$, as predicted by Eq. 138. At larger $T$ where
$\epsilon_{ad}\ll\epsilon_{f}$ the adiabatic time scales as
$T\propto\epsilon_{f}^{-1}\epsilon_{ad}^{-1}$, as predicted by Eq. 141.

In Fig. 9 (b) we show the total error $\epsilon$ versus $T$ for different
values of the final external field. As our analytical analysis predicts, for a
given $\epsilon$ (see for example the dashed line) the required adiabatic time
decreases with decreasing $\epsilon_{f}$.

To conclude, we find that for an adiabatic process with linear time dependence
of the external field, the adiabatic time scales as $\epsilon^{-2}$.

####  C.3.2 Adiabatic evolution with exponential dependence of the external
field

Figure 10: Adiabatic time to prepare the broken symmetry state in a two-level
system with $B=0.2$, $v=10$ and $f_{i}=1$. a) External field with exponential
time dependence. $T$ versus $\epsilon$ for different values of the final
external field $f=2\epsilon_{f}/B$. When $\epsilon\gg\epsilon_{f}$,
$T\approx\epsilon^{-1}\ln(\epsilon_{f}^{-1})$ in agreement with Eq. 149. The
approximation $\epsilon\approx\epsilon_{f}+\epsilon_{ad}$ does not hold well
when $\epsilon$ approaches $\epsilon_{f}$, but the bounds for the adiabatic
time given by Eq. 150 and Eq. 151 are valid (see dashed line where
$c_{1}=B/(2v)$ and $c_{2}=Bf_{i}/2$). b) Adiabatic time for the exponential
path (continuous lines, same legend as in (a)) and for the linear path (dotted
lines) versus accuracy. The adiabatic time for the exponential path scales as
$\epsilon^{-1}\ln\epsilon^{-1}$ which is better than the $\epsilon^{-2}$
scaling for the linear path.

In this scenario, the external field’s time dependence is

$\displaystyle f(s)=f_{i}e^{-\gamma s}~{}~{}\text{
with}~{}~{}\gamma=\ln\frac{f_{i}}{f_{f}}.$ (144)

Eqs. 122 and 123 reduce to

$\displaystyle\frac{dc_{a}}{ds}$ $\displaystyle=-\frac{1}{2}\gamma
fBe^{-i\frac{2Tv}{\gamma}\left(f_{i}-f\right)}c_{b}(s)$ (145)
$\displaystyle\frac{dc_{b}}{ds}$ $\displaystyle=\frac{1}{2}\gamma
fBe^{i\frac{2Tv}{\gamma}\left(f_{i}-f\right)}c_{a}(s)$ (146)

Employing the approximation $\left|1-c_{a}(s)\right|\approx\epsilon_{ad}$, one
has

$\displaystyle c_{b}(1)$
$\displaystyle\approx-\frac{1}{2}B\gamma\int_{0}^{1}dsf(s)e^{i2Tv\int_{0}^{s}f(u)du}+{\cal{O}}(\epsilon_{ad}^{2})$
(147)
$\displaystyle=\frac{iB\gamma}{4Tv}\left(e^{i2Tv\int_{0}^{1}f(s)ds}-1\right)+{\cal{O}}(\epsilon_{ad}^{2})$
$\displaystyle=\frac{iB\gamma}{4Tv}\left(e^{i\frac{2Tv\delta_{f}}{\gamma}}-1\right)+{\cal{O}}(\epsilon_{ad}^{2}).$

As for the linear adiabatic path case, $c_{b}(1)$ as a function of $T$ has an
oscillatory behavior. For the $T$ values which yield local maxima of
$\left|c_{b}(1)\right|$ (i.e. $T=(2k+1)\pi\gamma/(2v\delta_{f})$, with $k$
integer) the adiabatic error is

$\displaystyle\epsilon_{ad}=\left|c_{b}(1)\right|\approx\frac{B\gamma}{2Tv},$
(148)

implying

$\displaystyle
T\approx\frac{B}{2v}\frac{1}{\epsilon_{ad}}\ln\frac{f_{i}}{f_{f}}=\frac{B}{2v}\frac{1}{\epsilon_{ad}}\ln\frac{Bf_{i}}{2\epsilon_{f}}.$
(149)

Since $\epsilon_{f},\epsilon_{ad}\leq\epsilon$, a lower bound for $T$ is

$\displaystyle T>\frac{B}{2v\epsilon}\ln\frac{Bf_{i}}{2\epsilon}.$ (150)

For a desired accuracy $\epsilon$, it is possible to determine the optimal
$\epsilon_{f}$ and $\epsilon_{ad}$ which minimize $T$ numerically. To find an
asymptotic estimate for $T$, we assume $\epsilon_{f}$ and $\epsilon_{ad}$ are
independent contributions to the total error, i.e.
$\epsilon\approx\epsilon_{f}+\epsilon_{ad}$. Then the optimal $T$ is smaller
than or equal to the one obtained for
$\epsilon_{f}=\epsilon_{ad}\approx\epsilon/2$, i.e.

$\displaystyle T\leq\frac{B}{v\epsilon}\ln\frac{Bf_{i}}{\epsilon}.$ (151)

The bounds provided by Eqs. 150 and 151 show that the time required for
adiabatic evolution scales as

$\displaystyle T\propto\frac{1}{\epsilon}\ln\frac{1}{\epsilon},$ (152)

which is an improvement compared to using an external field with linear time
dependence, where $T\approx\epsilon^{-2}$ [see Eq. 138].

The results obtained by solving the differential equations (145) and (146)
numerically are shown in Fig. 10 (a). The adiabatic time satisfies Eq. 149.
The approximation $\epsilon\approx\epsilon_{f}+\epsilon_{ad}$ does not hold
well in the region where $\epsilon$ approaches $\epsilon_{f}$ (see the
continuous blue and the dot-dashed blue lines). However, the inequalities
(150) and (151) are true. In Fig. 10 (b) we compare the adiabatic time for the
linear adiabatic path (dotted lines) with the one corresponding to the
exponential adiabatic path (continuous lines). As Eqs. 138 and 152 predict,
the adiabatic time for the exponential case is much smaller than that for the
linear case.