# A Quantum N-Queens Solver

Valentin Torggler  Institute for Theoretical Physics, University of Innsbruck,
A-6020 Innsbruck, Austria  Philipp Aumann  Institute for Theoretical Physics,
University of Innsbruck, A-6020 Innsbruck, Austria  Helmut Ritsch  Institute
for Theoretical Physics, University of Innsbruck, A-6020 Innsbruck, Austria
Wolfgang Lechner<EMAIL_ADDRESS>Institute for Theoretical
Physics, University of Innsbruck, A-6020 Innsbruck, Austria  Institute for
Quantum Optics and Quantum Information of the Austrian Academy of Sciences,
A-6020 Innsbruck, Austria

(August 27, 2024)

###### Abstract

The $N$-queens problem is to find the position of $N$ queens on an $N$ by $N$
chess board such that no queens attack each other. The excluded diagonals
$N$-queens problem is a variation where queens cannot be placed on some
predefined fields along diagonals. This variation is proven NP-complete and
the parameter regime to generate hard instances that are intractable with
current classical algorithms is known. We propose a special purpose quantum
simulator that implements the excluded diagonals $N$-queens completion problem
using atoms in an optical lattice and cavity-mediated long-range interactions.
Our implementation has no overhead from the embedding allowing to directly
probe for a possible quantum advantage in near term devices for optimization
problems.

##  1 Introduction

Quantum technology with its current rapid advances in number, quality and
controllability of quantum bits (qubits) is approaching a new era with
computational quantum advantage for numerical tasks in reach [1, 2, 3, 4, 5,
6, 7, 8, 9, 10, 11]. While building a universal gate-based quantum computer
with error-correction is a long-term goal, the requirements on control and
fidelity to perform algorithms with such a universal device that outperform
their classical counterparts are still elusive. Building special purpose
quantum computers with near-term technology and proving computational
advantage compared to classical algorithms is thus a goal of the physics
community world wide [12]. Quantum simulation with the aim to solve
Hamiltonian systems may serve as a building block of such a special purpose
quantum computer [13, 14, 15]. In particular, adiabatic quantum computing [16,
17, 18] has been proposed to solve computationally hard problems by finding
the ground state of Ising spin glasses [19]. Despite considerable theoretical
[18] and experimental [20] efforts, quantum speedup in adiabatic quantum
computing has not been demonstrated in an experiment yet [21]. Thus,
demonstrating quantum advantage by solving optimization problems using quantum
simulation tools is a crucial step towards the development of general
programmable quantum optimizers [22, 23].

Figure 1: Sketch of the setup. \- (a) The $N$-queens problem is to place $N$
non-attacking queens on an $N$ by $N$ board. A variation thereof is the
N-queens completion problem where some of the queens are already placed
(yellow). In addition, some excluded diagonals are introduced (dashed blue
lines) on which no queen can be placed. (b) Each queen is represented by an
atom which is trapped in an anisotropic optical potential (blue) allowing for
tunneling in $x$-direction only. Collective scattering of pump laser light
(green arrows) into an optical resonator induces atom-atom interactions,
preventing atoms from aligning along the $y$-axis and along the diagonals.
After initial preparation in superposition states delocalized in $x$-direction
(red tubes), increasing the interactions transfers the system into its solid
phase, which is the solution of the queens problem (black balls).

Here we present a scheme that aims at solving the $N$-queens problem, and
variations of it, using atoms with cavity-mediated long-range interactions
[24, 25, 26, 27, 28]. We note that the $N$-queens problem is not just of
mathematical interest but also has some applications in computer science [29].
In this work, variations of the problem are used as a testbed [30] to study a
possible quantum advantage in solving classical combinatorial problems in near
term quantum experiments.

Our proposed setup consists of $N$ ultracold atoms in an optical lattice that
represent the queens on the chess board [31]. The non-attacking conditions are
enforced by a combination of restricted hopping [32] and interactions between
the atoms stemming from collective scattering of pump laser light into a
multi-mode cavity [33, 34, 35, 36, 37, 38, 39, 40] (see Fig. 1). For the
excluded diagonals variation of the $N$-queens problem, additional repulsive
optical potentials are introduced. The solution of the problem (or the ground
state of the many-body quantum system) is attained via a superfluid-to-solid
transition. From the measurement of photons that leave the cavity [41] it can
be determined if a state is a solution to the $N$-queens problem. The position
of the atoms can in addition be read out with single site resolved
measurement. The final solution is a classical configuration and thus easy to
verify. We show that a full quantum description of the dynamics is required to
find this solution.

Following Ref. [1], we identify a combination of several unique features of
the proposed model that makes it a viable candidate to test quantum advantage
in near term devices. (a) The completion and excluded diagonals problem is
proven to be NP-complete and hard instances for the excluded diagonals variant
are known from computer science literature [30], (b) the problem maps
naturally to the available toolbox of atoms in cavities and thus can be
implemented without intermediate embedding and no qubit overhead, (c) the
verification is computationally simple and (d) the number of qubits required
to solve problems which are hard for classical computers ($N>21$ for the
solvers used in Ref. [30]) is available in the lab.

Methods such as minor embedding [23, 42], LHZ [22, 43, 44] or nested embedding
[45] always cause a qubit overhead. Here the intermediate step of embedding
the optimization problem in an Ising model is removed by implementing the
infinite-range interactions with cavity-mediated forces tailored to the
problem’s geometry in combination with constrained tunneling [32]. Hence there
is no qubit overhead and the mode resources scale linearly with $N$. The
required number of qubits is reduced from several hundreds to below 50, which
is available in current experiments. By implementing our scheme with less than
50 atoms the problem is already hard to tackle with current classical
algorithms [30].

Light-mediated coupled tunneling gives rise to non-local quantum fluctuations
across the whole lattice in the intermediate stage of the transition [46, 47].
Their non-uniform signs stemming from the relative phases of the cavity fields
from site to site indicate that the system’s Hamiltonian is non-stoquastic and
can thus not be efficiently simulated with path integral Monte Carlo methods
on a classical computer [48]. The question of a quantum speed-up is thus open
unless a local transformation to a stoquastic Hamiltonian is found [49, 50].

In our model implementation the non-local qubit interactions are mediated via
the field modes of an optical resonator, which will attain non-classical atom-
field superposition states during the parameter sweep. This appears to be an
essential asset of the system as we find that the ground-state is reached only
with a very low probability, when the full quantum dynamics of the fields is
replaced by a classical mean-field approximation.

As a final feature let us point out here that the verification of a solution
is computationally trivial as the final state is classical and no quantum
tomography is needed. In principle the convergence to a solution can be simply
deduced from the cavity outputs at the end of the sweep. With this, the
proposed setup may serve as a platform to demonstrate combinatorial quantum
advantage in near-term experiments.

This work is organized as follows: In Sec. 2 we introduce a quantum model
based on coupled quantum harmonic oscillators simulating the $N$ queens
problem. A proposed physical implementation using ultracold atoms in optical
lattices and light-mediated atom-atom interaction is described in Sec. 3. In
Sec. 4 we present a numerical comparison between model and implementation
including photon loss. Finally we discuss in Sec. 5 how light leaking out of
the cavity can be used for read-out and we conclude in Sec. 6.

##  2 Quantum simulation of the $N$-queens problem

Following the idea of adiabatic quantum computation [16, 17, 18], we construct
a classical problem Hamiltonian $H_{\mathrm{pr}}$ such that its ground state
corresponds to the solutions of the $N$-queens problem. In order to find this
ground state, the system is evolved with the time-dependent Hamiltonian

$H(t)=H_{\mathrm{kin}}+\frac{t}{\tau}H_{\mathrm{pr}}$ (1)

from $t=0$ to $t=\tau$. Initially at $t=0$, the system is prepared in the
ground state of $H(0)=H_{\mathrm{kin}}$. During the time evolution, the second
term is slowly switched on. If this parameter sweep is slow enough, the system
stays in the instantaneous ground state and finally assumes the ground state
of $H(\tau)$ at $t=\tau$. If the lowest energy gap of $H_{\mathrm{pr}}$ is
much larger than the one of $H_{\mathrm{kin}}$, this state is close to the
ground state of $H_{\mathrm{pr}}$ and thus the solution of the optimization
problem.

In the following we construct the problem Hamiltonian $H_{\mathrm{pr}}$ and
the driver Hamiltonian $H_{\mathrm{kin}}$. The system is modeled as a 2D Bose-
Hubbard model with annihilation (creation) operators $b_{ij}$
($b_{ij}^{\dagger}$) on the sites $(i,j)$. A position of a queen is
represented by the position of an atom in an optical lattice with the total
number of atoms being fixed to $N$. The non-attacking condition between
queens, which amounts to interactions between two sites $(i,j)$ and $(k,l)$,
is implemented with four constraints: There can not be two queens on the same
line along (i) the $x$-direction $j=l$, (ii) the $y$-direction $i=k$, the
diagonals (iii) $i+j=k+l$ and (iv) $i-j=k-l$.

Condition (i) is implemented by using an initial state with one atom in each
horizontal line at $y_{j}$ and restricting the atomic movement to the
$x$-direction [see Fig. 1(b)]. Thereby we use the a priori knowledge that a
solution has one queen in a row, which reduces the accessible configuration
space size from ${N^{2}}\choose{N}$ to $N^{N}$ configurations. In this vein
the restricted tunneling Hamiltonian [32] is given by

$H_{\mathrm{kin}}=-J\sum_{i,j=1}^{N}\hat{B}_{ij},$ (2)

where $J$ is the tunneling amplitude and
$\hat{B}_{ij}=b_{i,j}^{\dagger}b_{i+1,j}+b_{i+1,j}^{\dagger}b_{i,j}$ with
$\hat{B}_{Nj}=0$ are the tunneling operators.

Constraints (ii), (iii) and (iv) are enforced by infinite range interactions
between the atoms with

$H_{\mathrm{Q}}=U_{\mathrm{Q}}\sum_{ijkl=1}^{N}A_{ijkl}\
\hat{n}_{ij}\hat{n}_{kl},$ (3)

where $\hat{n}_{i,j}=b_{i,j}^{\dagger}b_{i,j}$ and $U_{\mathrm{Q}}>0$. The
interaction matrix is

$A_{ijkl}=\begin{cases}3&\text{if }(i,j)=(k,l)\\\ 1&\text{if }i=k\lor
i+j=k+l\lor i-j=k-l\\\ 0&\text{otherwise},\end{cases}$ (4)

where in the first case all three constraints are broken.

In order to implement variations of the $N$-queens problem, we need to exclude
diagonals (for the excluded diagonals problem) and pin certain queens (for the
completion problem). These additional conditions are implemented by local
energy offsets of the desired lattice sites

$H_{\mathrm{pot}}=U_{\mathrm{D}}\sum_{i,j=1}^{N}D_{ij}\
\hat{n}_{ij}-U_{\mathrm{T}}\sum_{i,j=1}^{N}T_{ij}\ \hat{n}_{ij}.$ (5)

For $U_{\mathrm{D}}>0$ the first term renders occupations of sites on chosen
diagonals energetically unfavorable. Each diagonal (in $+$ and $-$ direction)
has an index summarized in the sets $\mathcal{D}_{+}$ and $\mathcal{D}_{-}$,
respectively, and the coefficients are

$D_{ij}=\begin{cases}2&\text{if }i+j-1\in\mathcal{D}_{+}\land
i-j+N\in\mathcal{D}_{-}\\\ 1&\text{if }i+j-1\in\mathcal{D}_{+}\lor
i-j+N\in\mathcal{D}_{-}\\\ 0&\text{otherwise}.\end{cases}$ (6)

For $U_{\mathrm{T}}>0$, the second term favors occupations of certain sites.
The sites where queens should be pinned to are pooled in the set $\mathcal{T}$
and therefore the coefficients are given by

$T_{ij}=\begin{cases}1&\text{if }(i,j)\in\mathcal{T}\\\
0&\text{otherwise}.\end{cases}$ (7)

The problem Hamiltonian of the $N$-queens problem with excluded diagonals is
then

$H_{\mathrm{pr}}=H_{\mathrm{Q}}+H_{\mathrm{pot}}.$ (8)

Note that due to the initial condition atoms never meet and sites are occupied
by zero or one atom only. Hence the system can be effectively described by
spin operators [31, 40], also without large contact interactions.

Let us illustrate the parameter sweep in Eq. (1) for a specific example
instance with $N=5$ queens (see Fig. 1). The excluded diagonals chosen here
restrict the ground state manifold to two solutions, and by biasing site
$(3,5)$ one of these solutions is singled out. The time evolution of the site
occupations $\langle\hat{n}_{ij}\rangle$ from numerically solving the time-
dependent Schrödinger equation is shown in Fig. 2. Initially, the atoms are
spread out in $x$-direction since the ground state of $H(0)=H_{\mathrm{kin}}$
is a superposition of excitations along each tube. After evolving for a
sufficiently large time $J\tau/\hbar=49$, the system is in the ground state of
$H_{\mathrm{pr}}$ and thus assumed the solution of the optimization problem.

Figure 2: Time evolution of site occupations. \- Each subplot shows a snapshot
of the site occupations $\langle\hat{n}_{ij}\rangle$ for the parameter sweep
in Eq. (1) and a sweep time $J\tau/\hbar=49$. This is the instance shown in
Fig. 1a, where the excluded diagonals are indexed by
$\mathcal{D}_{+}=\\{2,3,6,9\\}$ and $\mathcal{D}_{-}=\\{1,2,8,9\\}$ and one
queen is pinned at site $(3,5)$, i.e. $\mathcal{T}=\\{(3,5)\\}$ (see text).
The final values of the sweep are $U_{\mathrm{Q}}=J$, $U_{\mathrm{D}}=5J$,
$U_{\mathrm{T}}=2J$ while $J$ is kept constant. Since the sweep time is large
enough, the state of the system adiabatically converges to the unique solution
of the problem, which can be easily verified. Figure 3: Energy spectrum. \-
Eigenvalue spectrum of $H(t)$ for (a) the model Hamiltonian [Eq. (1)] and (b)
the Hamiltonian created by the field-atom interactions [Eq. (10)] for the
instance described in Fig. 2. The orange and green lines show ground state and
first excited state, respectively. The blue lines depict higher energy
eigenstates. The low energy sector, and especially the minimal gap, are
qualitatively similar.

The energy spectrum of the given instance is shown in Fig. 3(a). The minimal
gap between ground state (orange) and first excited state (green) determines
the minimum sweep time $\tau$ to remain in the ground state according to the
Landau-Zener formula. At the end of the sweep, the ground state closely
resembles the solution to the excluded diagonals problem shown in Fig. 1.

The Hilbert space for the atomic state corresponding to the configuration
space mentioned above grows exponentially as $N^{N}$ and thus, as usual for
quantum systems, the computational costs get large for rather small systems.
Simulations with significantly larger systems are hence not easily tractable.

##  3 Implementation

Here we propose a specific and at least conceptionally simple and
straightforward physical implementation of Eq. (1), where the $N$ queens are
directly represented by $N$ ultracold atoms in an $N$x$N$ two-dimensional
optical lattice. Assuming tight binding conditions the atoms are confined to
the lowest band of the lattice. They can coherently tunnel between sites [51]
and interact via collective light scattering within an optical resonator. To
implement the required queens interactions via light scattering we make use of
a set of optical field modes in a multi-mode standing-wave resonator [see Fig.
1(b)]. It was shown before that this configuration in principle allows to
implement arbitrary site to site interactions [40] if a sufficient number of
modes is used.

To reduce the necessary Hilbert space without loss of generality, it is
sufficient to enable tunneling only along one dimension ($x$), e.g. the rows
of the lattice and increase the lattice depth in the column ($y$) direction.
This confines the atoms (queens) to move only along tubes in the $x$-direction
forming a parallel array of $N$ 1D optical lattices as routinely used to study
1D physics with cold atoms [52].

Luckily, for the special case of the infinite range queens interactions, one
can find a strikingly simple and intuitive example configuration requiring
only few field modes for each of the three remaining queens interaction
directions. For this we consider the optical lattice to be placed in the
center symmetry plane of the optical resonator, where one has a common anti-
node of all symmetric eigenmodes. The trapped atoms are then illuminated by
running plane wave laser beams from three different directions within the
lattice plane with frequencies matched to different longitudinal cavity modes.

Depending on the size of the problem we need to add more laser frequencies in
each direction to avoid a periodic recurrence of the interaction within the
lattice. Note that all frequencies can be easily derived and simultaneously
stabilized from a single frequency comb with a spacing matched to the cavity
length. Since all light frequencies are well separated compared to the cavity
line width, each is scattered into a distinct cavity mode and as the cavity
modes are not directly coupled no relative phase stability is needed [40]. The
resulting position-dependent collective scattering into the cavity then
introduces the desired infinite-range interactions between the atoms. Choosing
certain frequencies and varying their relative pump strengths allows for
tailoring these interactions to simulate the three non-attacking conditions in
the queens problem discussed in the previous section. Additional light sheets
and local optical tweezers can be used to make certain diagonals energetically
unfavorable (for the excluded diagonal problem) or pin certain queens.
Alternatively pinning can be achieved by selective cavity mode injection.

In the first part (Sec. 3.1) we will now introduce a Bose-Hubbard-type
Hamiltonian for the trapped atoms interacting via collective scattering for
general illumination fields in more detail. For this we first consider the
full system of coupled atoms and cavity modes and later adiabatically
eliminate the cavity fields to obtain effective light-mediated atom-atom
interactions.

In the second part (Sec. 3.2) we discuss how to implement the non-attacking
conditions of the $N$-queens problem with such light-mediated atom-atom
interactions. Specifically, we consider the limit of a deep optical lattice
leading to simple analytical expressions. These formulas allow us to find a
specific pump configuration of wave numbers and pump strengths leading to the
ideal queens interaction Hamiltonian in Eq. (3). As a deep lattice depth slows
down atomic tunneling it requires long annealing times. Luckily, it turns out
that the pump configuration derived for a deep lattice still gives the same
ground state for moderate lattice depth with faster tunneling. We show this by
numerical simulations in Sec. 4. While the details of the interaction are
altered it still sufficiently well approximates the $N$-queens interaction.

###  3.1 Tight-binding model for atoms interacting via light

Driving far from any atomic resonance the internal degrees of freedom of the
atoms can be eliminated. In this so-called dispersive limit the resulting
effective Hamiltonian couples the atomic motion to the light fields [53].

Single-particle Hamiltonian. For a single particle of mass $m_{A}$, the motion
of the atoms in the $x$-$y$-plane is described by [25, 54]

$\displaystyle H_{1}=$
$\displaystyle\frac{\hat{p}_{x}^{2}+\hat{p}_{y}^{2}}{2m_{A}}+V_{\mathrm{L}}^{x}\cos^{2}(k_{\mathrm{L}}\hat{x})+V_{\mathrm{L}}^{y}\cos^{2}(k_{\mathrm{L}}\hat{y})$
(9)
$\displaystyle+V_{\mathrm{bias}}|\mathcal{F}(\hat{x},\hat{y})|^{2}-\hbar\sum_{m=1}^{M_{\mathrm{tot}}}\tilde{\Delta}_{c,m}a_{m}^{\dagger}a_{m}$
$\displaystyle+\hbar\sum_{m=1}^{M_{\mathrm{tot}}}\eta_{m}\left(h_{m}^{*}(\hat{x},\hat{y})a_{m}+a_{m}^{\dagger}h_{m}(\hat{x},\hat{y})\right).$

The first line contains the kinetic term with the momentum operators
$\hat{p}_{x}$ and $\hat{p}_{y}$. Classical electric fields create optical
potentials with depths $V_{\mathrm{L}}^{x}$, $V_{\mathrm{L}}^{y}$ and
$V_{\mathrm{bias}}$. The first two create the optical lattice with wave number
$k_{\mathrm{L}}$ and lattice spacing $a=\pi/k_{\mathrm{L}}$, while
$V_{\mathrm{bias}}$ is much smaller and only responsible for a bias field on
certain sites, for instance for excluding diagonals. Thereby
$\mathcal{F}(x,y)$ is an electric field distribution whose maximum is
normalized to one.

The last two terms describe the free evolution of the cavity fields and atom-
state-dependent scattering of the pump fields into the cavity, the atom-light
interaction. The quantized electric cavity fields are described by $a_{m}$
($a_{m}^{\dagger}$), the annihilation (creation) operators of a photon in the
$m$-th mode. These fields are coupled to the classical pump fields with mode
functions $h_{m}(x,y)$ via the effective scattering amplitudes
$\eta_{m}=g_{m}\Omega_{m}/\Delta_{a,m}$, with the pump laser Rabi frequencies
$\Omega_{m}$, the atom-cavity couplings $g_{m}$ and the detunings between pump
lasers and atomic resonance frequency $\Delta_{a,m}$. The effective cavity
detunings $\tilde{\Delta}_{c,m}=\Delta_{c,m}-NU_{0,m}$ are given by the
detunings between pump laser and cavity mode frequencies $\Delta_{c,m}$ and
the dispersive shifts of the cavity resonance due to the presence of the atoms
in the cavity $NU_{0,m}$ [24].

Note that, for example by placing the optical lattice ($x$-$y$-plane) in a
common anti-node of the standing wave cavity modes and exciting only TEM00
modes, the atom-cavity coupling is uniform in space in our model. Thus the
only spatial dependence in the cavity term is due to the pump fields.

Generalized Bose-Hubbard Hamiltonian. The atom-atom interactions are taken
into account by introducing bosonic field operators
$\hat{\Psi}(\boldsymbol{x})$ with
$\tilde{H}=\int\mathrm{d}^{2}x\hat{\Psi}^{\dagger}(\boldsymbol{x})H_{1}\Psi(\boldsymbol{x})$
[55, 25]. Note that we do not include contact interactions since atoms never
meet due to the initial condition we will use. We assume that the optical
lattice with depths $V_{\mathrm{L}}^{x}$ and $V_{\mathrm{L}}^{y}$ is so deep,
that the atoms are tightly bound at the potential minima and only the lowest
vibrational state (Bloch band) is occupied. Moreover, the optical potential
created by bias and cavity fields is comparably small, such that the form of
the Bloch wave functions only depends on the optical lattice [54]. In this
limit we can expand the bosonic field operators in a localized Wannier basis
$\hat{\Psi}(\boldsymbol{x})=\sum_{i,j}w_{2D}(\boldsymbol{x}-\boldsymbol{x}_{ij})b_{ij}$
with the lowest-band Wannier functions $w_{2D}(\boldsymbol{x})$ coming from
Bloch wave functions of the lattice [56]. We split the resulting Hamiltonian
in three terms

$\tilde{H}=H_{\mathrm{kin}}+H_{\mathrm{cav}}+H_{\mathrm{pot}},$ (10)

which will be explained in the following.

As in the standard Bose-Hubbard model, one obtains a tunneling term
$H_{\mathrm{kin}}$ as in Eq. (2). Tunneling in $y$-direction is frozen out by
ensuring $V_{\mathrm{L}}^{y}\gg V_{\mathrm{L}}^{x}$. The other terms
$H_{\mathrm{cav}}$ and $H_{\mathrm{pot}}$ originate from the weak cavity-pump
interference fields and the bias fields introduced above and should resemble
$H_{\mathrm{pr}}$ [Eq. (8)]. In order to realize the sweep Eq. (1), the
relative strength of these terms and the kinetic term has to be tuned, e.g. by
ramping up the pump laser and bias field intensity (make $H_{\mathrm{cav}}$
and $H_{\mathrm{pot}}$ larger) or the lattice depth (make $H_{\mathrm{kin}}$
smaller).

The cavity-related terms in Eq. (9) give rise to

$\displaystyle H_{\mathrm{cav}}=$
$\displaystyle-\hbar\sum_{m}\tilde{\Delta}_{c,m}a_{m}^{\dagger}a_{m}$ (11)
$\displaystyle+\hbar\sum_{m}N\eta_{m}\left(\hat{\Theta}_{m}^{\dagger}a_{m}+a_{m}^{\dagger}\hat{\Theta}_{m}\right)$

with the order operator of cavity mode $m$

$\hat{\Theta}_{m}=\frac{1}{N}\sum_{i,j=1}^{N}\left(v_{m}^{ij}\hat{n}_{ij}+u_{m}^{ij}\hat{B}_{ij}\right).$
(12)

The structure of the fields enters in the on-site and nearest neighbor atom-
mode overlaps

$\displaystyle v_{m}^{ij}$ $\displaystyle=\int
dx\,w^{2}(x-x_{i})h_{m}(x,y_{j})$ (13) $\displaystyle u_{m}^{ij}$
$\displaystyle=\int dx\,w(x-x_{i})h_{m}(x,y_{j})w(x-x_{i+1})$ (14)

where $y_{j}=ja$ with $j=0,...,N-1$ are the tube positions and $w(x)$ the one-
dimensional Wannier functions in $x$-direction. This is because for
$V_{\mathrm{L}}^{y}\gg V_{\mathrm{L}}^{x}$ we can approximate the
$y$-dependence of the Wannier functions by a Dirac delta:
$w_{2D}(\boldsymbol{x})=w(x)\delta(y)$.

The last term $H_{\mathrm{pot}}$ describes all extra fields responsible for
local energy off-sets at certain sites that stem from the weak classical
fields with the distribution $\mathcal{F}(x,y)$ ($V_{\mathrm{bias}}\ll
V_{\mathrm{L}}^{x}$). The off-sets ought to be calculated from the overlaps of
fields and Wannier functions, analogously to $v_{m}^{ij}$. Thus the fields
have to be chosen such that the resulting Hamiltonian resembles Eq. (5). We do
not detail the derivation further here and use $H_{\mathrm{pot}}$ for
numerical simulations.

Atom-atom interaction Hamiltonian. The main focus of this work is to show how
to create the tailored all-to-all particle interactions via collective
scattering. We derive this interaction by eliminating the cavity fields
introduced in the previous section [57, 58, 54]. This can be done because the
cavity fields decay through the mirrors with the rates $\kappa_{m}$ and thus
end up in a particular steady-state for each atomic configuration. Assuming
that the atomic motion is much slower than the cavity field dynamics, i.e.
$J/\hbar\ll|\tilde{\Delta}_{c,m}+i\kappa_{m}|$, this steady-state is a good
approximation at all times. The stationary cavity field amplitudes are given
by

$a_{m}^{\mathrm{st}}\equiv\frac{\eta_{m}}{\tilde{\Delta}_{c,m}+i\kappa_{m}}N\hat{\Theta}_{m}$
(15)

and thus replaced by atomic operators (see Appendix B).

In the coherent regime $|\tilde{\Delta}_{c,m}|\gg\kappa_{m}$, the atom-light
interaction is then described by an effective interaction Hamiltonian for the
atoms [28]

$H_{\mathrm{cav}}^{\mathrm{eff}}=\hbar\sum_{m}\frac{\tilde{\Delta}_{c,m}\eta_{m}^{2}}{\tilde{\Delta}_{c,m}^{2}+\kappa_{m}^{2}}N^{2}\hat{\Theta}_{m}^{\dagger}\hat{\Theta}_{m}.$
(16)

The collective, state-dependent scattering induces interactions between each
pair of sites $(i,j)$ and $(k,l)$: Density-density interactions due to the
terms containing $\hat{n}_{ij}\hat{n}_{kl}$ and a modified tunneling amplitude
due an occupation or a tunneling event somewhere else in the lattice due to
$\hat{n}_{ij}\hat{B}_{kl}$, $\hat{B}_{ij}\hat{n}_{kl}$ and
$\hat{B}_{ij}\hat{B}_{kl}$. While the density-density interactions constitute
the problem Hamiltonian [see Eq. (3)], the latter cavity-induced tunneling
terms lead to non-local fluctuations which might help to speed up the
annealing process [46].

Since the Wannier functions are localized at the lattice sites, the on-site
overlaps $v_{m}^{ij}$ tend to be much larger than the nearest-neighbor
overlaps $u_{m}^{ij}$. Thus density-density interactions are expected to be
the dominant contribution to $H_{\mathrm{cav}}^{\mathrm{eff}}$. As intuitively
expected, the atoms localize stronger for deeper lattices, where analytical
expressions for the overlaps can be obtained within a harmonic approximation
of the potential wells leading to Gaussian Wannier functions with a width
$\propto(V_{\mathrm{L}}^{x})^{-1/4}$. Apart from a correction factor due to
this width, the on-site overlaps are given by the pump fields at the lattice
sites. The nearest-neighbor overlaps correspond to the pump fields in between
the lattice sites, but are exponentially suppressed (see Appendix C).
Consequently, in the deep lattice limit (large $V_{\mathrm{L}}^{x}$) when the
width tends to zero we get $v_{m}^{ij}=h_{m}(x_{i},y_{j})$ and $u_{m}^{ij}=0$
[28], and an interaction Hamiltonian [from Eq. (16)]

$\displaystyle H_{\mathrm{cav}}^{\mathrm{dl}}=\hbar\sum_{m}$
$\displaystyle\frac{\tilde{\Delta}_{c,m}\eta_{m}^{2}}{\tilde{\Delta}_{c,m}^{2}+\kappa_{m}^{2}}$
(17) $\displaystyle\times$
$\displaystyle\sum_{ijkl}h_{m}^{*}(x_{i},y_{j})h_{m}(x_{k},y_{l})\hat{n}_{ij}\hat{n}_{kl},$

which only depends on density operators, and hence does not include cavity
induced-tunneling.

###  3.2 $N$-queens interaction

In this section we aim to find pump fields $h_{m}(x,y)$ such that the
interaction Hamiltonian in the deep lattice limit Eq. (17) corresponds to the
desired queens Hamiltonian $H_{\mathrm{Q}}$ [Eq. (3)] containing the non-
attacking conditions. Using these pump fields we later show numerically in
Sec. 4, that the atom-atom interaction for realistic lattice depths [Eq.
(16)], although slightly altered, still well resembles the queens interaction.

We consider three sets of $M$ parallel running wave laser beams with different
propagation directions, each of which could be created by a frequency comb.
The three directions are perpendicular to the lines along which queens should
not align, that is along the $x$-direction and along the diagonals. We denote
the corresponding wave vectors with
$\boldsymbol{k}_{m}^{x}=(k_{m}^{0},0)^{T}$,
$\boldsymbol{k}_{m}^{+}=(k_{m}^{0},k_{m}^{0})^{T}$ and
$\boldsymbol{k}_{m}^{-}=(k_{m}^{0},-k_{m}^{0})^{T}$, respectively, with the
wave numbers $k_{m}^{0}$. Therefore the pump fields are given by

$h_{m}(x,y)=e^{i\boldsymbol{k}_{m}\boldsymbol{x}},$ (18)

where $\boldsymbol{x}=(x,y)^{T}$ is the position vector and
$\boldsymbol{k}_{m}$ is a wave vector in any of the three directions.

With running wave pump fields, Eq. (17) can be written as

$H_{\mathrm{cav}}^{\mathrm{dl}}=U_{\mathrm{Q}}\sum_{ijkl}\tilde{A}_{ijkl}\hat{n}_{ij}\hat{n}_{kl}.$
(19)

This formally corresponds to $H_{\mathrm{Q}}$ [Eq. (3)], where the quantities
now have a physical meaning: The interaction matrix is given by

$\tilde{A}_{ijkl}=\sum_{m}f_{m}\cos(\boldsymbol{k}_{m}(\boldsymbol{x}_{ij}-\boldsymbol{x}_{kl}))$
(20)

with lattice site connection vectors $\boldsymbol{x}_{ij}-\boldsymbol{x}_{kl}$
and

$f_{m}U_{\mathrm{Q}}=\hbar\frac{\tilde{\Delta}_{c,m}\eta_{m}^{2}}{\tilde{\Delta}_{c,m}^{2}+\kappa_{m}^{2}}$
(21)

with $\sum_{m=0}^{M-1}f_{m}=1$. The dimensionless parameters $f_{m}$ capture
the relative strengths of the modes, determining the shape of the interaction.
They have to be chosen such that $\tilde{A}$ approximates $A$ [Eq. (4)]. The
overall strength of the interaction term is captured by the energy
$U_{\mathrm{Q}}$, which can be easily tuned by the cavity detunings or the
pump intensities to implement the parameter sweep in Eq. (1). For the
following discussion we define an interaction function

$\tilde{\mathcal{A}}(\boldsymbol{r})=\sum_{m}f_{m}\cos(\boldsymbol{k}_{m}\boldsymbol{r})$
(22)

which returns the interaction matrix when evaluated at lattice site connection
vectors
$\tilde{A}_{ijkl}=\tilde{\mathcal{A}}(\boldsymbol{x}_{ij}-\boldsymbol{x}_{kl})$.

We note that one set of parallel $\boldsymbol{k}_{m}^{\mu}$
($\mu\in\\{x,+,-\\}$) creates an interaction $\tilde{\mathcal{A}}$ which is
constant and infinite range (only limited by the laser beam waist) in the
direction perpendicular to the propagation direction
$\boldsymbol{r}\perp\boldsymbol{k}_{m}^{\mu}$. Along the propagation direction
$\boldsymbol{r}\parallel\boldsymbol{k}_{m}^{\mu}$ instead, the interaction is
shaped according to the sum of cosines, and can be modified by the choice of
wave numbers $k_{m}^{\mu}=|\boldsymbol{k}_{m}^{\mu}|$ and their relative
strengths $f_{m}$.

In the following we will use the example wave numbers

$k_{m}^{0}=k_{\mathrm{L}}\left(1+\frac{2m+1}{2M}\right)$ (23)

with $m=0,...,M-1$ and uniform $f_{m}=1/M$. Taking into account
$\boldsymbol{k}_{m}^{x}$ only, the interaction along the $x$-direction
$\boldsymbol{r}\parallel\boldsymbol{k}_{m}^{x}$ at lattice site distances
$r_{j}=j\pi/k_{\mathrm{L}}=ja$ has the values

$\tilde{\mathcal{A}}(r_{j})=\begin{cases}(-1)^{l}\;&\text{for}\;j=2Ml,\;l\in\mathbb{Z}\\\
0\;&\text{otherwise},\end{cases}$ (24)

as shown in Appendix D. If we guarantee that $-2M<j<2M$ this results in an
interaction which is zero everywhere apart from $j=0$, i.e. at zero distance.
So for repulsive interactions ($U_{\mathrm{Q}}>0$ and thus
$\tilde{\Delta}_{c,m}>0$), the wave vectors $\boldsymbol{k}_{m}^{x}$ create
the non-attacking interaction along the $y$-direction ($\tilde{A}_{ijkl}=1$ if
$i=k$ and $0$ otherwise) as long as $N\leq 2M$. This is illustrated in Fig.
4(a) for $M=N=5$. Analogously, $\boldsymbol{k}_{m}^{\pm}$ cause the non-
attacking interactions along the diagonals. In a square lattice the diagonals
have the distance $r_{j}/\sqrt{2}$, which is compensated by
$k_{m}^{\pm}=|\boldsymbol{k}_{m}^{\pm}|=\sqrt{2}k_{m}^{0}$. Since there are
$2N-1$ diagonals, one has to make sure that $2N-1\leq 2M$. Upon combining all
wave vectors from three directions we finally obtain the full queens
interaction, as shown in Fig. 4(b), which is realized with
$M_{\mathrm{tot}}=3M=3N$ frequencies in our example.

Figure 4: Energy penalty created by one atom. \- These density plots show the
energy penalty $\tilde{\mathcal{A}}(\boldsymbol{x}-\boldsymbol{x}_{a})$ for an
atom at position $\boldsymbol{x}=(x,y)^{T}$ created by an atom at position
$\boldsymbol{x}_{a}=(a,a)^{T}$ (red dot) for $N=M=5$ [see Eq. (22)]. The dots
indicate lattice site positions and the solid (dashed) contour lines indicate
where $\tilde{\mathcal{A}}(\boldsymbol{x}-\boldsymbol{x}_{a})=0$
($\tilde{\mathcal{A}}(\boldsymbol{x}-\boldsymbol{x}_{a})=1$) is fulfilled. (a)
Pumping along the $x$-axis with the wave vectors
$\boldsymbol{k}_{m}^{x}=(k_{m}^{0},0)^{T}$ with
$k_{m}^{0}/k_{\mathrm{L}}=1.1,1.3,1.5,1.7,1.9$ [according to Eq. (23)] creates
interactions along $y$. (b) Additionally including diagonal pump lasers
$\boldsymbol{k}_{m}^{+}=(k_{m}^{0},k_{m}^{0})^{T}$ and
$\boldsymbol{k}_{m}^{-}=(k_{m}^{0},-k_{m}^{0})^{T}$ implements the full queens
interaction along diagonals and vertical lines [Eq. (3)].

Note that there are several combinations of wave numbers and mode strengths
which, at least approximately, create the desired line-shaped interactions
perpendicular to the light propagation direction. For this it is insightful to
reformulate the interaction as a Fourier transform. To deal with continuous
functions, we define an envelope $f(k)$ with $f(k_{m})=f_{m}$, which is
sampled at the wave numbers $n\Delta k$ with $n\in\mathbb{Z}$ containing all
$k_{m}^{\mu}$. Considering one illumination direction for simplicity, the
interaction [Eq. (22)] along $\boldsymbol{r}\parallel\boldsymbol{k}_{m}^{\mu}$
with $r=|\boldsymbol{r}|$ can be written as

$\displaystyle\tilde{\mathcal{A}}(r)$
$\displaystyle=\mathrm{Re}\left[\sqrt{2\pi}\mathcal{F}\left\\{f(k)\sum_{n=-\infty}^{\infty}\delta(k-n\Delta
k)\right\\}(r)\right]$ (25)
$\displaystyle=\sum_{l=-\infty}^{\infty}\mathrm{Re}\left[\sqrt{2\pi}\mathcal{F}\\{f\\}\left(r-l\frac{2\pi}{\Delta
k}\right)\right],$

where
$\mathcal{F}\\{f\\}(r)=\int_{-\infty}^{\infty}\mathrm{d}kf(k)e^{ikr}/\sqrt{2\pi}$
is the Fourier transform of $f(k)$ and $\delta(x)$ is a Dirac delta at $x=0$.
See Appendix D for a detailed derivation.

The last line allows for a simple interpretation: The interaction consists of
peaks repeating with a spatial period $R=2\pi/\Delta k$. Each of these peaks
has the shape of the real part of the Fourier transform of the envelope
function $f(k)$ with a width corresponding to the inverse of the mode
bandwidth $\sigma\sim 2\pi/\Delta k_{\mathrm{BW}}$.

For the (approximate) non-attacking condition ($\tilde{\mathcal{A}}(ja)\approx
1$ for $j=0$ and $|\tilde{\mathcal{A}}(ja)|\ll 1$ otherwise) there are two
conditions. Firstly, at most one peak should be within the region of the
atoms. Thus the period has to be larger than the (diagonal) size of the
optical lattice $R\geq Na$ ($R\geq\sqrt{2}Na$). Secondly, the width of one
peak has to be smaller than the lattice spacing $\sigma\lesssim a$. Combining
these conditions to $R\gtrsim N\sigma$, we see that the minimum number of
modes per direction $M$ scales linearly with $N$

$M\approx\Delta k_{\mathrm{BW}}/\Delta k\gtrsim N.$ (26)

Therefore, with only $\sim N$ modes this quite generically allows for creating
an interaction along lines perpendicular to the light propagation. Note
however, that the second condition also implies that the spatial frequency
spread has to be at least on the order of the lattice wave number $\Delta
k_{\mathrm{BW}}\gtrsim k_{\mathrm{L}}$.

##  4 Numerical justification of assumptions

We compare the ideal model Hamiltonian described in Sec. 2 [Eq. (1)] to the
physically motivated tight-binding Hamiltonian for finite lattice depths
introduced in Sec. 3 [Eq. (10)]. As for the ideal model in Fig. 2, we consider
the time evolution during a slow linear sweep of $U_{\mathrm{Q}}$,
$U_{\mathrm{T}}$ and $U_{\mathrm{D}}$ by numerically integrating the time-
dependent Schrödinger equation. Physically, this sweep can be realized by
ramping up the pump and the bias field intensities. Moreover, we show that
evolving the system using a classical approximation for the cavity mode fields
does not result in a solution to the $N$-queens problem, and address the
effect of dephasing by photon loss with open system simulations.

In the following we use a realistic lattice depth of
$V_{\mathrm{L}}^{x}=10E_{\mathrm{R}}$ with the recoil energy
$E_{\mathrm{R}}=\hbar^{2}k_{\mathrm{L}}^{2}/(2m_{\mathrm{A}})$. For example,
for rubidium 87Rb and $\lambda_{\mathrm{L}}=785.3\,\mathrm{nm}$ it is
$E_{\mathrm{R}}/\hbar=23.4\,\mathrm{kHz}$ [26]. The chosen lattice depth leads
to a tunneling amplitude $J\approx 0.02E_{\mathrm{R}}$, which can be obtained
from the band structure of the lattice. We consider our cavity model in Eq.
(10) for $N=5$. The pump modes are as in Sec. 3.2 and Fig. 4. While in the
limit of a deep lattice this would result in the ideal model interactions,
here they depend on the overlaps between Wannier functions and pump modes [Eq.
(13)] and are thus altered. In the following the overlaps are calculated with
Wannier functions which where numerically obtained from the band structure of
the lattice. It turns out, that the deviation from the ideal overlaps does not
qualitatively change the interaction for the realistic parameters used.

###  4.1 Coherent dynamics

The energy spectrum of the Hamiltonian for finite lattice depths in Eq. (16)
is shown in Fig. 3(b) and is qualitatively of the same form as for the ideal
model in Fig. 3(a). In comparison the eigenvalue gaps tend to be smaller at
the end of the sweep. This is because the on-site atom-mode overlaps decrease
for shallower lattices and less localized atoms due to a smoothing of the mode
functions by the finite width Wannier functions (see Appendix C). Moreover, we
consider the time evolution during the nearly adiabatic sweep for
$J\tau/\hbar=49$ for the same parameters by integrating the time-dependent
Schrödinger equation. Snapshots of the site occupations
$\langle\hat{n}_{ij}\rangle$ for several times are shown in Fig. 5(a), where
we observe a similar behavior as for the model in Fig. 2. This suggests that
the system is robust against the errors introduced by the moderate lattice
depth. Note that the physical Hamiltonian used here is non-stoquastic in the
occupation basis due to cavity-induced tunneling.

Figure 5: Comparison of the time evolution with quantum and classical fields.
\- Lattice site occupations for a time evolution during a nearly adiabatic
sweep using the cavity Hamiltonian including cavity-assisted tunneling for the
same parameters as in Fig. 2. Subplot (a) shows the dynamics using the full
quantum interaction Hamiltonian [Eq. (16)]. It closely resembles the results
from the model Hamiltonian in Fig. 2. Subplot (b) shows the time evolution
with the classical approximation of the cavity fields [Eq. (27)]. The state
does not converge to the solution, also not for much larger sweep times. The
modes for both cases where chosen as in Fig. 4(b).

###  4.2 Classical cavity fields

In the following we show that quantum correlations are crucial for the
efficiency of the sweep. In particular, if we substitute the field operators
$a_{m}^{\mathrm{st}}$ by its expectation values representing classical cavity
fields the solution is not found. We consider the semi-classical Hamiltonian

$\displaystyle H_{\mathrm{cav}}^{\mathrm{class}}=$
$\displaystyle\hbar\sum_{m}\frac{N^{2}\tilde{\Delta}_{c,m}\eta_{m}^{2}}{\tilde{\Delta}_{c,m}^{2}+\kappa_{m}^{2}}$
(27)
$\displaystyle\times\left(\hat{\Theta}_{m}^{\dagger}\langle\hat{\Theta}_{m}\rangle+\langle\hat{\Theta}_{m}^{\dagger}\rangle\hat{\Theta}_{m}-|\langle\hat{\Theta}_{m}\rangle|^{2}\right),$

where the expectation values have to be calculated self-consistently with the
current atom state vector. This substitution amounts to considering only first
order fluctuations around the mean of $\hat{\Theta}_{m}$ in Eq. (16).

Consequently, the dynamics are described by a differential equation which is
non-linear in the state vector $|\psi\rangle$. We numerically solve this
equation by self-consistently updating the expectation value in each time
step. It turns out that even for very long sweep times, using classical fields
does not lead to a solution of the queens problem. The time evolution for
$J\tau/\hbar=49$ is depicted in Fig. 5(b). The discrepancy shows the necessity
of entangled light-matter states in our procedure.

###  4.3 Dephasing due to cavity field loss

Finally, motivated by experimental considerations, we consider the open system
including photon decay through the cavity mirrors. This system can be
described by a master equation for the atoms [54]

$\displaystyle\dot{\rho}=$
$\displaystyle-\frac{i}{\hbar}[H_{\mathrm{kin}}+H_{\mathrm{cav}}^{\mathrm{eff}},\rho]$
(28)
$\displaystyle+\sum_{m}\frac{N^{2}\eta_{m}^{2}\kappa_{m}}{\Delta_{c,m}^{2}+\kappa_{m}^{2}}\left(2\hat{\Theta}_{m}\rho\hat{\Theta}_{m}^{\dagger}-\\{\hat{\Theta}_{m}^{\dagger}\hat{\Theta}_{m},\rho\\}\right),$

where curly brackets denote the anti-commutator (see Appendix B). The model is
suitable for analyzing the dephasing close to the coherent regime before any
steady state is reached.

More insight can be gained by rewriting it in the basis of scattering
eigenstates $|\nu\rangle$ with
$\hat{\Theta}_{m}|\nu\rangle=\theta_{m}^{\nu}|\nu\rangle$, which scatter a
field
$\alpha_{m}^{\nu}=\frac{N\eta_{m}}{\tilde{\Delta}_{c,m}+i\kappa_{m}}\theta_{m}^{\nu}$.
These states converge to the occupation states in the deep lattice limit. The
time evolution for the matrix elements
$\rho_{\mu\nu}=\langle\mu|\rho|\nu\rangle$ reads

$\displaystyle\dot{\rho}_{\mu\nu}=$
$\displaystyle(-\Gamma_{\mu\nu}-i\Omega_{\mu\nu})\rho_{\mu\nu}$ (29)
$\displaystyle-i\frac{J}{\hbar}\sum_{k}\left(\langle\mu|\hat{B}|k\rangle\rho_{k\nu}-\langle
k|\hat{B}|\nu\rangle\rho_{\mu k}\right)$

with the rates

$\displaystyle\Omega_{\mu\nu}=$
$\displaystyle\sum_{m}\tilde{\Delta}_{c,m}(|\alpha_{m}^{\mu}|^{2}-|\alpha_{m}^{\nu}|^{2})$
(30) $\displaystyle\Gamma_{\mu\nu}=$
$\displaystyle\sum_{m}\kappa_{m}|\alpha_{m}^{\mu}-\alpha_{m}^{\nu}|^{2}.$ (31)

While the energy gaps $\hbar\Omega_{\mu\nu}$ describe the coherent dynamics we
considered up to now, the dephasing rates $\Gamma_{\mu\nu}$ stem from photon
loss. Hence a superposition of two scattering eigenstates $|\mu\rangle$ and
$|\nu\rangle$ looses its coherence depending on the difference of the
scattered fields, or how distinguishable the states are by field measurement.

We now return to the example we had before, using uniform mode strengths
$f_{m}=1/M$, decay rates $\kappa$ and detunings $\tilde{\Delta}_{c}$. Figure 6
compares the coherent Schrödinger time evolution from Sec. 4.1, the mean-field
approximation with classical cavity fields from Sec. 4.2 and the open system
dynamics with dephasing. For the latter the master equation is approximated by
using Monte-Carlo wave function simulations. As a measure of similarity
between states we use the fidelity. For two mixed states it is defined as
$F(\rho,\sigma)=\mathrm{Tr}(\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}})$ and reduces
to the overlap $|\langle\psi|\phi\rangle|$ for pure states.

Figure 6: Comparison to open system dynamics. \- (a) The figure shows the
fidelity $F_{\mathrm{sol}}$ between the solution and the instantaneous
eigenstate (black dotted), and between solution and dynamical state from the
time evolution using the Schrödinger equation (blue), classical fields (green)
and the open system (red). The open system fidelity depends on
$\tilde{\Delta}_{c}/\kappa$, where we present results for the values
$5,50,100,200,500,1000$ (from light red to dark red). (b) The fidelity between
solution and the final state after the sweep as a function of
$\tilde{\Delta}_{c}/\kappa$ for the curves of (a), where the colors are the
same as in (a). The parameters used here are as in Fig. 5.

The open system dynamics is depicted for different detunings
$\tilde{\Delta}_{c}/\kappa$ while keeping $U_{\mathrm{Q}}$ fixed. This can be
achieved by adjusting the pump strength $\eta$ correspondingly. In this case
the dephasing rates

$\frac{\hbar\Gamma_{\mu\nu}}{U_{\mathrm{Q}}}=\frac{\kappa}{\tilde{\Delta}_{c}}N^{2}\sum_{m}|\theta_{m}^{\mu}-\theta_{m}^{\nu}|^{2}$
(32)

go to zero for $\tilde{\Delta}_{c}/\kappa\gg 1$. Thus, as expected, the open
system converges to the coherent Schrödinger dynamics in this limit (see Fig.
6(b)).

Note that the coherence between states creating similar fields is preserved
much longer than for other states, which is expected to be important at the
late stage of the sweep. For states with fixed similarity (e.g. one atom
moved), $|\theta_{m}^{\mu}-\theta_{m}^{\nu}|^{2}$ is on the order of $N^{-2}$,
and thus the dephasing rates do not scale with $N$ for such states.

##  5 Read-out

After the parameter sweep we need to determine if the obtained state is a
solution or not. This can in principle be done by reading out the final atomic
state with single site resolution using a quantum gas microscope [59, 60].
However, as we consider an open system with the cavity output fields readily
available, we will show that by proper measurements on the output light we can
directly answer this question without further additions. Note that after the
sweep at the stage of the read-out, quantum coherences do not have to be
preserved since the solution is a classical state. This gives the freedom to
increase the lattice depth to some high value in the deep lattice regime to
suppress further tunneling, and to increase the pump power or decrease the
detunings in order to get a stronger signal at the detector.

###  5.1 Intensity measurement

For uniform cavity detunings, a state corresponding to the solution of the
$N$-queens problem scatters less photons than all other states. Thus the
measurement of the total intensity in principle allows one to distinguish a
solution from other states. To illustrate this we consider the total rate of
photons impinging on a detector scattered by an atomic state $|\psi\rangle$

$\displaystyle P(|\psi\rangle)=$
$\displaystyle\sum_{m}2\kappa_{m}\langle(a_{m}^{\mathrm{st}})^{\dagger}a_{m}^{\mathrm{st}}\rangle$
(33) $\displaystyle=$
$\displaystyle\sum_{m}\frac{2\kappa_{m}}{\tilde{\Delta}_{c,m}}\frac{\tilde{\Delta}_{c,m}\eta_{m}^{2}N^{2}}{\tilde{\Delta}_{c,m}^{2}+\kappa_{m}^{2}}\langle\hat{\Theta}_{m}^{\dagger}\hat{\Theta}_{m}\rangle$
$\displaystyle\approx$ $\displaystyle
U_{\mathrm{Q}}\frac{\zeta}{\hbar}\sum_{ijkl}\tilde{A}_{ijkl}\langle\hat{n}_{ij}\hat{n}_{kl}\rangle=\frac{\zeta}{\hbar}\langle
H_{\mathrm{cav}}^{\mathrm{dl}}\rangle.$

In the last line we assumed that $\zeta=2\kappa_{m}/\tilde{\Delta}_{c,m}$ does
not depend on $m$ and a deep lattice.

Since $P$ is proportional to the energy expectation value, the ground state,
i.e. the solution of the queens problem, causes a minimal photon flux at the
detector $P_{0}=3NU_{\mathrm{Q}}\zeta$. It stems from the on-site terms
$(i,j)=(k,l)$, where the factor $3$ comes from the three pump directions. In
contrast, each pair of queens violating the non-attacking condition in
$\tilde{A}$ leads to an increase of the photon flux by $\Delta
P=2U_{\mathrm{Q}}\zeta$. The two atoms create an energy penalty for one
another, explaining the factor $2$. The relative difference of the photon flux
due to a state with $L$ attacking pairs and a solution is given by

$\frac{L\Delta P}{P_{0}}=\frac{2L}{3N}.$ (34)

As this scales with $1/N$ it is difficult to distinguish solutions from other
states via measurement of the intensity for large $N$. Note that for non-
uniform $\kappa_{m}/\tilde{\Delta}_{c,m}$, photons from different modes have
to be distinguished.

Figure 7: Signature of an atomic state in the cavity output. \- In left and
right column the signature of two different atomic (pure) states are compared,
whose occupations are shown on the top (a). On the left there is a solution to
the $N$-queens problem, while on the right one atom was moved. (b) The polar
plot shows the cavity field expectation value $\langle
a_{m}^{\mathrm{st}}\rangle$ in the complex plane. Measuring the fields by
homodyne detection yields a certain quadrature of the field depending on the
phase angle $\phi$. An example is illustrated by the black line. (c) The
measurable quadratures are shown for the example angle. The sequences can be
Fourier transformed to obtain the occupations of columns and of the diagonals.
The parameters used here are $U_{\mathrm{Q}}=5J$ and
$\tilde{\Delta}_{c}/\kappa=10$. We used $M=2N-1=9$ modes per direction with
wave numbers as in Eq. (23). Different colors encode the different directions
of the pump modes: $\boldsymbol{k}_{m}^{x}$ (blue), $\boldsymbol{k}_{m}^{+}$
(green) and $\boldsymbol{k}_{m}^{-}$ (red). The lighter the color the larger
the wave number $|\boldsymbol{k}_{m}^{x,\pm}|$.

###  5.2 Field measurement

Measurement of the output field quadratures, for example by homodyne
detection, gives insight about the absolute position of the atoms projected
onto the pump laser propagation direction. For the three directions used in
our setup, this yields the occupations of each column
$\langle\hat{N}_{i}^{x}\rangle$ and each diagonal
$\langle\hat{N}_{i}^{+}\rangle$ and $\langle\hat{N}_{i}^{-}\rangle$. Since a
solution of the queens problem has maximally one atom on each diagonal and
exactly one atom on each column, it must fulfill

$\langle\hat{N}_{i}^{x}\rangle=1\land\langle\hat{N}_{i}^{+}\rangle\leq
1\land\langle\hat{N}_{i}^{-}\rangle\leq 1.$ (35)

The output field quadratures can thus be used to determine if a classical
final state is a solution or not, which is the answer to the blocked diagonals
decision problem we aim to solve. The signatures of two example states in the
cavity fields are depicted in Fig. 7.

Let us illustrate the measurement by considering only light scattered from the
$x$-direction with incident wave vectors $\boldsymbol{k}_{m}^{x}$. Since these
plane waves are constant in $y$-direction, the atom-field overlaps do not
depend on $j$. Neglecting cavity-induced tunneling, the field quadratures for
a phase difference $\phi$ are

$\displaystyle\mathrm{Re}(\langle a_{m}^{\mathrm{st}}\rangle e^{-i\phi})$
$\displaystyle=\sum_{i}\mathrm{Re}\left(\frac{\eta_{m}e^{-i\phi}}{\tilde{\Delta}_{c,m}+i\kappa_{m}}v_{m}^{i1}\right)\langle\hat{N}_{i}^{x}\rangle,$
(36)

revealing that cavity fields are determined by the total occupations of the
columns $\hat{N}_{i}^{x}=\sum_{j}\hat{n}_{ij}$.

For at least $N$ modes ($M\geq N$) this system of equations can be inverted
yielding the column occupations $\langle\hat{N}_{i}^{x}\rangle$. By measuring
the cavity output field quadratures scattered from the diagonal pump light we
obtain the occupations of each diagonal $\langle\hat{N}_{i}^{+}\rangle$ and
$\langle\hat{N}_{i}^{-}\rangle$. Inverting the system of equations for
diagonals demands at least as many pump modes as diagonals, that is $2N-1$.
Inversion can also be done efficiently and intuitively by using a discrete
Fourier transform and its inverse. To reveal the Fourier relation, one has to
express the $v_{m}^{ij}$’s in Eq. (36) within the harmonic (or in the deep
lattice) approximation (see Appendix C). The so obtained approximate inversion
formula also works well for realistic lattice depths.

Strictly speaking the condition in Eq. (35) is sufficient only for classical
configurations, like occupation number basis states $|\phi_{\nu}\rangle$. Some
superpositions $|\psi\rangle=\sum_{\nu}c_{\nu}|\phi_{\nu}\rangle$ which are no
solutions might also fulfill the above criterion, because summands in the
field expectation values $\langle
a_{m}^{\mathrm{st}}\rangle=\sum_{\nu}|c_{\nu}|^{2}\langle\phi_{\nu}|a_{m}^{\mathrm{st}}|\phi_{\nu}\rangle$
can cancel each other. For instance, for $U_{\mathrm{T}}=0$ the solution from
our example in Fig. 2 $|\psi_{\mathrm{sol}}\rangle=|1,4,2,5,3\rangle$ scatters
the same fields $\langle a_{m}^{\mathrm{st}}\rangle$ as the superposition
$|\psi_{\mathrm{nosol}}\rangle=(|\psi^{1}_{\mathrm{nosol}}\rangle+|\psi^{2}_{\mathrm{nosol}}\rangle)/\sqrt{2}$
with $|\psi^{1}_{\mathrm{nosol}}\rangle=|1,3,2,5,4\rangle$ and
$|\psi^{2}_{\mathrm{nosol}}\rangle=|1,4,5,2,3\rangle$, both of which are no
solution. In this notation the state $|i_{1},i_{2},...,i_{N}\rangle$ has one
atom on each site $(i_{j},j)$. However, these macroscopic superpositions are
highly unstable. Even theoretically the measurement back-action [61, 41]
projects superpositions of states scattering different fields (such as
$|\psi_{\mathrm{nosol}}\rangle$) to one of its constituents. The inclusion of
measurement back-action due to continuous measurement might thus lead to
intriguing phenomena beyond those presented here and is subject to future
work.

We emphasize again that the measurements described above answer the question
if we found a solution or not, which is the answer to the combinatorial
decision problem. The exact configuration of the final state can be measured
with single site resolution as demonstrated in several experiments [59, 60].

##  6 Conclusions

We present a special purpose quantum simulator with the aim to solve
variations of the $N$-queens problem based on atoms in a cavity. This
combinatorial problem may serve as a benchmark to study a possible quantum
advantage in intermediate size near term quantum experiments. From the
algorithmic point of view, the problem is interesting for quantum advantage as
it is proven NP-hard and instances can be found that are not solvable with
current state-of-the-art algorithms. From the implementation point of view,
the proposed quantum simulator implements the queens problem without overhead
and thus a few tens of atoms are sufficient to enter the classically
intractable regime. The proposed setup of atoms in a cavity fits the queens
problem naturally as the required infinite range interactions arise there
inherently. We find that by treating the light field classically the
simulation does not find the solutions suggesting that quantum effects like
atom-field entanglement cannot be neglected. Moreover, we investigate the
influence of photon loss on the coherence time.

The queens problem is formulated as a decision problem, asking whether there
is a valid configuration of queens or not given the excluded diagonals and
fixed queens. Remarkably, to answer the decision problem, a read-out of the
atom positions is not required as the necessary information is encoded in the
light that leaves the cavity. To determine the position of the queens requires
single site resolved read-out, which is also available in several current
experimental setups [59].

In this work we concentrated on the coherent regime. The driven-dissipative
nature of the system provides additional features which can be exploited for
obtaining the ground state. For certain regimes, cavity cooling [24, 62] can
help to further reduce sweep times and implement error correction. Moreover,
the back action of the field measurement onto the atomic state can be used for
preparing states [41].

Note that an implementation of the $N$-queens problem for a gate-based quantum
computer was proposed in Ref. [63] aiming to find a solution of the
unconstrained $N$-queens problem. Our work in contrast employs an adiabatic
protocol and intends to answer the question if a solution exists given
constraints of blocked diagonals or already placed queens, which was shown to
be NP-complete and numerically hard [30].

Acknowledgments. We thank I. Gent, C. Jefferson and P. Nightingale for
fruitful discussions. Simulations were performed using the open source
QuantumOptics.jl framework in Julia [64] and we thank D. Plankensteiner for
related discussions. V. T. and H. R. are supported by Austrian Science Fund
Project No. I1697-N27. W. L. acknowledges funding by the Austrian Science Fund
(FWF) through a START grant under Project No. Y1067-N27 and the SFB BeyondC
Project No. F7108-N38, the Hauser-Raspe foundation, and the European Union’s
Horizon 2020 research and innovation program under grant agreement No. 817482
PasQuanS.

## References

  * Harrow and Montanaro [2017] Aram W Harrow and Ashley Montanaro.  Quantum computational supremacy.  _Nature_ , 549(7671):203–209, 2017.  doi: 10.1038/nature23458.
  * Cirac and Zoller [2012] J Ignacio Cirac and Peter Zoller.  Goals and opportunities in quantum simulation.  _Nat. Phys._ , 8(4):264–266, apr 2012.  ISSN 1745-2473.  doi: 10.1038/nphys2275.
  * Blatt and Roos [2012] Rainer Blatt and Christian F Roos.  Quantum simulations with trapped ions.  _Nature Physics_ , 8(4):277–284, 2012.  doi: 10.1038/nphys2252.
  * Bloch et al. [2012] Immanuel Bloch, Jean Dalibard, and Sylvain Nascimbene.  Quantum simulations with ultracold quantum gases.  _Nat. Phys._ , 8(4):267–276, apr 2012.  ISSN 1745-2473.  doi: 10.1038/nphys2259.
  * Aspuru-Guzik and Walther [2012] Alán Aspuru-Guzik and Philip Walther.  Photonic quantum simulators.  _Nature Physics_ , 8(4):285–291, 2012.  doi: 10.1038/nphys2253.
  * De Léséleuc et al. [2018] Sylvain De Léséleuc, Sebastian Weber, Vincent Lienhard, Daniel Barredo, Hans Peter Büchler, Thierry Lahaye, and Antoine Browaeys.  Accurate mapping of multilevel rydberg atoms on interacting spin-1/2 particles for the quantum simulation of ising models.  _Physical review letters_ , 120(11):113602, 2018\.  doi: 10.1103/PhysRevLett.120.113602.
  * Houck et al. [2012] Andrew A Houck, Hakan E Türeci, and Jens Koch.  On-chip quantum simulation with superconducting circuits.  _Nature Physics_ , 8(4):292–299, 2012.  doi: 10.1038/nphys2251.
  * Georgescu et al. [2014] I. M. Georgescu, S. Ashhab, and Franco Nori.  Quantum simulation.  _Rev. Mod. Phys._ , 86:153–185, Mar 2014.  doi: 10.1103/RevModPhys.86.153.
  * Saffman [2016] Mark Saffman.  Quantum computing with atomic qubits and rydberg interactions: progress and challenges.  _Journal of Physics B: Atomic, Molecular and Optical Physics_ , 49(20):202001, 2016.  doi: 10.1088/0953-4075/49/20/202001.
  * Preskill [2018] John Preskill.  Quantum computing in the nisq era and beyond.  _Quantum_ , 2:79, 2018.  doi: 10.22331/q-2018-08-06-79.
  * Boixo et al. [2018] Sergio Boixo, Sergei V Isakov, Vadim N Smelyanskiy, Ryan Babbush, Nan Ding, Zhang Jiang, Michael J Bremner, John M Martinis, and Hartmut Neven.  Characterizing quantum supremacy in near-term devices.  _Nature Physics_ , 14(6):595, 2018.  doi: 10.1038/s41567-018-0124-x.
  * Preskill [2012] John Preskill.  Quantum computing and the entanglement frontier.  _arXiv:1203.5813 [quant-ph]_ , 2012.
  * Bernien et al. [2017] Hannes Bernien, Sylvain Schwartz, Alexander Keesling, Harry Levine, Ahmed Omran, Hannes Pichler, Soonwon Choi, Alexander S Zibrov, Manuel Endres, Markus Greiner, et al.  Probing many-body dynamics on a 51-atom quantum simulator.  _Nature_ , 551(7682):579, 2017.  doi: 10.1038/nature24622.
  * Kim et al. [2010] Kihwan Kim, M-S Chang, Simcha Korenblit, Rajibul Islam, Emily E Edwards, James K Freericks, G-D Lin, L-M Duan, and Christopher Monroe.  Quantum simulation of frustrated ising spins with trapped ions.  _Nature_ , 465(7298):590–593, 2010.  doi: 10.1038/nature09071.
  * Schauß et al. [2012] Peter Schauß, Marc Cheneau, Manuel Endres, Takeshi Fukuhara, Sebastian Hild, Ahmed Omran, Thomas Pohl, Christian Gross, Stefan Kuhr, and Immanuel Bloch.  Observation of spatially ordered structures in a two-dimensional rydberg gas.  _Nature_ , 491(7422):87–91, 2012.  doi: 10.1038/nature11596.
  * Kadowaki and Nishimori [1998] Tadashi Kadowaki and Hidetoshi Nishimori.  Quantum annealing in the transverse ising model.  _Phys. Rev. E_ , 58:5355–5363, Nov 1998.  doi: 10.1103/PhysRevE.58.5355.
  * Farhi et al. [2000] Edward Farhi, Jeffrey Goldstone, Sam Gutmann, and Michael Sipser.  Quantum computation by adiabatic evolution.  _arXiv:quant-ph/0001106_ , 2000.
  * Albash and Lidar [2018] Tameem Albash and Daniel A Lidar.  Adiabatic quantum computation.  _Reviews of Modern Physics_ , 90(1):015002, 2018\.  doi: 10.1103/RevModPhys.90.015002.
  * Lucas [2014] Andrew Lucas.  Ising formulations of many np problems.  _Frontiers in Physics_ , 2:5, 2014.  doi: 10.3389/fphy.2014.00005.
  * Boixo et al. [2014] Sergio Boixo, Troels F Rønnow, Sergei V Isakov, Zhihui Wang, David Wecker, Daniel A Lidar, John M Martinis, and Matthias Troyer.  Evidence for quantum annealing with more than one hundred qubits.  _Nature Physics_ , 10(3):218, 2014.  doi: 10.1038/nphys2900.
  * Hauke et al. [2019] Philipp Hauke, Helmut G Katzgraber, Wolfgang Lechner, Hidetoshi Nishimori, and William D Oliver.  Perspectives of quantum annealing: Methods and implementations.  _arXiv:1903.06559 [quant-ph]_ , 2019.
  * Lechner et al. [2015] Wolfgang Lechner, Philipp Hauke, and Peter Zoller.  A quantum annealing architecture with all-to-all connectivity from local interactions.  _Science Advances_ , 1(9), 2015.  doi: 10.1126/sciadv.1500838.
  * Choi [2008] Vicky Choi.  Minor-embedding in adiabatic quantum computation: I. the parameter setting problem.  _Quantum Information Processing_ , 7(5):193–209, 2008.  doi: 10.1007/s11128-008-0082-9.
  * Ritsch et al. [2013] Helmut Ritsch, Peter Domokos, Ferdinand Brennecke, and Tilman Esslinger.  Cold atoms in cavity-generated dynamical optical potentials.  _Reviews of Modern Physics_ , 85(2):553, 2013\.  doi: 10.1103/RevModPhys.85.553.
  * Maschler and Ritsch [2005] Christoph Maschler and Helmut Ritsch.  Cold atom dynamics in a quantum optical lattice potential.  _Physical review letters_ , 95(26):260401, 2005\.  doi: 10.1103/PhysRevLett.95.260401.
  * Landig et al. [2016] Renate Landig, Lorenz Hruby, Nishant Dogra, Manuele Landini, Rafael Mottl, Tobias Donner, and Tilman Esslinger.  Quantum phases from competing short-and long-range interactions in an optical lattice.  _Nature_ , 532(7600):476, 2016.  doi: 10.1038/nature17409.
  * Mekhov and Ritsch [2012] Igor B Mekhov and Helmut Ritsch.  Quantum optics with ultracold quantum gases: towards the full quantum regime of the light–matter interaction.  _Journal of Physics B: Atomic, Molecular and Optical Physics_ , 45(10):102001, 2012.  doi: 10.1088/0953-4075/45/10/102001.
  * Caballero-Benitez et al. [2016] Santiago F Caballero-Benitez, Gabriel Mazzucchi, and Igor B Mekhov.  Quantum simulators based on the global collective light-matter interaction.  _Physical Review A_ , 93(6):063632, 2016.  doi: 10.1103/PhysRevA.93.063632.
  * Bell and Stevens [2009] Jordan Bell and Brett Stevens.  A survey of known results and research areas for n-queens.  _Discrete Mathematics_ , 309(1):1–31, 2009.  doi: 10.1016/j.disc.2007.12.043.
  * Gent et al. [2017] Ian P Gent, Christopher Jefferson, and Peter Nightingale.  Complexity of n-queens completion.  _Journal of Artificial Intelligence Research_ , 59:815–848, 2017.  doi: 10.1613/jair.5512.
  * Hen [2017] Itay Hen.  Realizable quantum adiabatic search.  _EPL (Europhysics Letters)_ , 118(3):30003, 2017\.  doi: 10.1209/0295-5075/118/30003.
  * Hen and Spedalieri [2016] Itay Hen and Federico M Spedalieri.  Quantum annealing for constrained optimization.  _Physical Review Applied_ , 5(3):034007, 2016\.  doi: 10.1103/PhysRevApplied.5.034007.
  * Domokos and Ritsch [2002] Peter Domokos and Helmut Ritsch.  Collective cooling and self-organization of atoms in a cavity.  _Physical review letters_ , 89(25):253003, 2002\.  doi: 10.1103/PhysRevLett.89.253003.
  * Black et al. [2003] Adam T Black, Hilton W Chan, and Vladan Vuletić.  Observation of collective friction forces due to spatial self-organization of atoms: from rayleigh to bragg scattering.  _Physical review letters_ , 91(20):203001, 2003\.  doi: 10.1103/PhysRevLett.91.203001.
  * Vaidya et al. [2018] Varun D Vaidya, Yudan Guo, Ronen M Kroeze, Kyle E Ballantine, Alicia J Kollár, Jonathan Keeling, and Benjamin L Lev.  Tunable-range, photon-mediated atomic interactions in multimode cavity qed.  _Physical Review X_ , 8(1):011002, 2018.  doi: 10.1103/PhysRevX.8.011002.
  * Gupta et al. [2007] Subhadeep Gupta, Kevin L Moore, Kater W Murch, and Dan M Stamper-Kurn.  Cavity nonlinear optics at low photon numbers from collective atomic motion.  _Physical review letters_ , 99(21):213601, 2007\.  doi: 10.1103/PhysRevLett.99.213601.
  * Gopalakrishnan et al. [2011] Sarang Gopalakrishnan, Benjamin L Lev, and Paul M Goldbart.  Frustration and glassiness in spin models with cavity-mediated interactions.  _Physical review letters_ , 107(27):277201, 2011\.  doi: 10.1103/PhysRevLett.107.277201.
  * Gopalakrishnan et al. [2012] Sarang Gopalakrishnan, Benjamin L Lev, and Paul M Goldbart.  Exploring models of associative memory via cavity quantum electrodynamics.  _Philosophical Magazine_ , 92(1-3):353–361, 2012\.  doi: 10.1080/14786435.2011.637980.
  * Krämer and Ritsch [2014] Sebastian Krämer and Helmut Ritsch.  Self-ordering dynamics of ultracold atoms in multicolored cavity fields.  _Physical Review A_ , 90(3):033833, 2014.  doi: 10.1103/PhysRevA.90.033833.
  * Torggler et al. [2017] Valentin Torggler, Sebastian Krämer, and Helmut Ritsch.  Quantum annealing with ultracold atoms in a multimode optical resonator.  _Physical Review A_ , 95(3):032310, 2017.  doi: 10.1103/PhysRevA.95.032310.
  * Mekhov and Ritsch [2009] Igor B Mekhov and Helmut Ritsch.  Quantum nondemolition measurements and state preparation in quantum gases by light detection.  _Physical review letters_ , 102(2):020403, 2009\.  doi: 10.1103/PhysRevLett.102.020403.
  * Choi [2011] Vicky Choi.  Minor-embedding in adiabatic quantum computation: Ii. minor-universal graph design.  _Quantum Information Processing_ , 10(3):343–353, 2011.  doi: 10.1007/s11128-010-0200-3.
  * Rocchetto et al. [2016] Andrea Rocchetto, Simon C Benjamin, and Ying Li.  Stabilizers as a design tool for new forms of the lechner-hauke-zoller annealer.  _Science advances_ , 2(10):e1601246, 2016.  doi: 10.1126/sciadv.1601246.
  * Glaetzle et al. [2017] Alexander W Glaetzle, Rick MW van Bijnen, Peter Zoller, and Wolfgang Lechner.  A coherent quantum annealer with rydberg atoms.  _Nature Communications_ , 8:15813, 2017.  doi: 10.1038/ncomms15813.
  * Vinci and Lidar [2018] Walter Vinci and Daniel A Lidar.  Scalable effective-temperature reduction for quantum annealers via nested quantum annealing correction.  _Physical Review A_ , 97(2):022308, 2018.  doi: 10.1103/PhysRevA.97.022308.
  * Hormozi et al. [2017] Layla Hormozi, Ethan W Brown, Giuseppe Carleo, and Matthias Troyer.  Nonstoquastic hamiltonians and quantum annealing of an ising spin glass.  _Physical Review B_ , 95(18):184416, 2017.  doi: 10.1103/PhysRevB.95.184416.
  * Albash [2019] Tameem Albash.  Role of nonstoquastic catalysts in quantum adiabatic optimization.  _Physical Review A_ , 99(4):042334, 2019.  doi: 10.1103/PhysRevA.99.042334.
  * Isakov et al. [2016] Sergei V Isakov, Guglielmo Mazzola, Vadim N Smelyanskiy, Zhang Jiang, Sergio Boixo, Hartmut Neven, and Matthias Troyer.  Understanding quantum tunneling through quantum monte carlo simulations.  _Physical review letters_ , 117(18):180402, 2016\.  doi: 10.1103/PhysRevLett.117.180402.
  * Klassen and Terhal [2019] Joel Klassen and Barbara M. Terhal.  Two-local qubit Hamiltonians: when are they stoquastic?  _Quantum_ , 3:139, May 2019.  ISSN 2521-327X.  doi: 10.22331/q-2019-05-06-139.
  * Marvian et al. [2019] Milad Marvian, Daniel A Lidar, and Itay Hen.  On the computational complexity of curing non-stoquastic hamiltonians.  _Nature communications_ , 10(1):1571, 2019.  doi: 10.1038/s41467-019-09501-6.
  * Bloch [2005] Immanuel Bloch.  Ultracold quantum gases in optical lattices.  _Nature Physics_ , 1(1):23–30, 2005.  doi: doi.org/10.1038/nphys138.
  * Meinert et al. [2014] Florian Meinert, Manfred J Mark, Emil Kirilov, Katharina Lauber, Philipp Weinmann, Michael Gröbner, Andrew J Daley, and Hanns-Christoph Nägerl.  Observation of many-body dynamics in long-range tunneling after a quantum quench.  _Science_ , 344(6189):1259–1262, 2014.  doi: 10.1126/science.1248402.
  * Domokos and Ritsch [2003] Peter Domokos and Helmut Ritsch.  Mechanical effects of light in optical resonators.  _JOSA B_ , 20(5):1098–1130, 2003.  doi: 10.1364/JOSAB.20.001098.
  * Maschler et al. [2008] Christoph Maschler, Igor B Mekhov, and Helmut Ritsch.  Ultracold atoms in optical lattices generated by quantized light fields.  _The European Physical Journal D_ , 46(3):545–560, 2008.  doi: 10.1140/epjd/e2008-00016-4.
  * Jaksch et al. [1998] Dieter Jaksch, Christoph Bruder, Juan Ignacio Cirac, Crispin W Gardiner, and Peter Zoller.  Cold bosonic atoms in optical lattices.  _Physical Review Letters_ , 81(15):3108, 1998\.  doi: 10.1103/PhysRevLett.81.3108.
  * Kohn [1959] Walter Kohn.  Analytic properties of bloch waves and wannier functions.  _Physical Review_ , 115(4):809, 1959.  doi: 10.1103/PhysRev.115.809.
  * Nagy et al. [2009] D Nagy, P Domokos, A Vukics, and H Ritsch.  Nonlinear quantum dynamics of two bec modes dispersively coupled by an optical cavity.  _The European Physical Journal D_ , 55(3):659, 2009.  doi: 10.1140/epjd/e2009-00265-7.
  * Habibian et al. [2013] Hessam Habibian, André Winter, Simone Paganelli, Heiko Rieger, and Giovanna Morigi.  Bose-glass phases of ultracold atoms due to cavity backaction.  _Physical review letters_ , 110(7):075304, 2013\.  doi: 10.1103/PhysRevLett.110.075304.
  * Bakr et al. [2009] Waseem S Bakr, Jonathon I Gillen, Amy Peng, Simon Fölling, and Markus Greiner.  A quantum gas microscope for detecting single atoms in a hubbard-regime optical lattice.  _Nature_ , 462(7269):74, 2009.  doi: 10.1038/nature08482.
  * Sherson et al. [2010] Jacob F Sherson, Christof Weitenberg, Manuel Endres, Marc Cheneau, Immanuel Bloch, and Stefan Kuhr.  Single-atom-resolved fluorescence imaging of an atomic mott insulator.  _Nature_ , 467(7311):68, 2010.  doi: 10.1038/nature09378.
  * Carmichael [1993] H. Carmichael.  _An Open Systems Approach to Quantum Optics: Lectures Presented at the Université Libre de Bruxelles, October 28 to November 4, 1991_.  Number Bd. 18 in An Open Systems Approach to Quantum Optics: Lectures Presented at the Université Libre de Bruxelles, October 28 to November 4, 1991\. Springer Berlin Heidelberg, 1993.  ISBN 9783540566342.
  * Wolke et al. [2012] Matthias Wolke, Julian Klinner, Hans Keßler, and Andreas Hemmerich.  Cavity cooling below the recoil limit.  _Science_ , 337(6090):75–78, 2012.  doi: 10.1126/science.1219166.
  * Jha et al. [2018] Rounak Jha, Debaiudh Das, Avinash Dash, Sandhya Jayaraman, Bikash K Behera, and Prasanta K Panigrahi.  A novel quantum n-queens solver algorithm and its simulation and application to satellite communication using ibm quantum experience.  _arXiv:1806.10221 [quant-ph]_ , 2018.
  * Krämer et al. [2018] Sebastian Krämer, David Plankensteiner, Laurin Ostermann, and Helmut Ritsch.  Quantumoptics.jl: A julia framework for simulating open quantum systems.  _Computer Physics Communications_ , pages –, 2018.  ISSN 0010-4655.  doi: 10.1016/j.cpc.2018.02.004.

## Appendix

##  Appendix A Instance parameters

Table 1 provides an overview of the chosen parameters for the exemplary linear
parameter sweep in the main text (Figs. 2, 3 and 5).

Parameter | Symbol | Value
---|---|---
System size | $N$ | $5$
Final queens interaction energy | $U_{\mathrm{Q}}$ | $J$
Final excluded diagonals penalty | $U_{\mathrm{D}}$ | $5J$
Final trapping energy | $U_{\mathrm{T}}$ | $2J$
Sweep time | $\tau$ | $49\hbar/J$
Excluded sum-diagonals |  | $\\{2,3,6,9\\}$
Excluded difference-diagonals |  | $\\{1,2,8,9\\}$
Trapping sites |  | $\\{(3,5)\\}$
Number of modes per direction | $M$ | 5
Table 1: Parameters of the exemplary instance used in the figures in the main
text.

We now describe how we choose the parameters used in our example. For this we
calculate the minimal gap and the overlap with the final solution for several
parameters to find a region with large minimal gap and large overlap. Note
that this is only done to find good parameters for our small example, where we
already know the solution. For large systems such a calculation would beyond
classical numerical capabilities, which is why the problem poses a potential
application for a quantum simulator.

The minimal gap in the spectrum (e.g. the one shown in Fig. 3) depends on the
final queens interaction energy $U_{\mathrm{Q}}$, the final trapping energy
$U_{\mathrm{T}}$, the tunneling amplitude $J$ and the final excluded diagonals
penalty $U_{\mathrm{D}}$. To find proper values for these parameters we
determine the minimal gap in a wide parameter range. In order to get the
minimal gap, some of the Hamiltonian’s lowest eigenenergies are calculated for
discrete time steps during the sweep. Subsequently, the minimum of the
difference between the groundstate and the first exited state at all time
steps is taken to be the minimal gap. The values of the minimal gap have to be
scrutinized carefully since its accuracy depends on the resolution of the
discrete time steps. Therefore a more detailed analysis of the minimal gap
might require a more careful analysis, especially for high interaction
strengths.

To analyze how well the quantum system reproduces the solution of the
$N$-queens problem we study the overlap

$F=|\braket{\phi}{\psi}|$ (37)

between the state $\ket{\phi}$ that corresponds to the solution of the chosen
instance of the queens problem introduced in Fig. 1 and the state at the end
of an adiabatic sweep $\ket{\psi}$ (i.e. the ground state of our spectrum on
the right side). This is necessary because we do not switch off the kinetic
Hamiltonian in our example, and thus the "perfect" solution is only obtained
in the limit of large energy penalties $U_{\mathrm{Q}}$, $U_{\mathrm{D}}$ and
$U_{\mathrm{T}}$.

Figure 8(a) suggests that in order to increase the minimal gap the ratio
$U_{\mathrm{Q}}/U_{\mathrm{D}}$ has to be chosen as small as possible. We vary
the ratio by fixing $U_{\mathrm{D}}$ and varying $U_{\mathrm{Q}}$. Therewith,
Fig. 8(a) indicates that $U_{\mathrm{Q}}$ should be as small as possible.
However, as it can be seen in Fig. 8(b), a small $U_{\mathrm{Q}}$ also
decreases the overlap with the solution and the physical system does not
resemble the desired solution of the queens problem anymore. We therefore have
to make a compromise between a reasonably large overlap and an optimized
minimal gap.

If we set $U_{\mathrm{D}}$ to $5J$ and $U_{\mathrm{T}}$ to $2J$ we find that
for $U_{\mathrm{Q}}=1J$ the overlap is $F\approx 0.93$ and the minimal gap is
around $0.44J$. These values were used for Figs. 2 and 3.

Figure 8: Discussion of the impact of parameter values. \- (a) The minimal gap
of the energy spectrum for the sweep as a function of the final parameters.
$U_{\mathrm{D}}$ is fixed to $5J$ and $U_{\mathrm{Q}}$ is varied. For the
discussion on fidelity [subplot (b)] we choose $U_{\mathrm{T}}=2J$ which is
indicated by the black line. The values of $U_{\mathrm{T}}$ are non-zero to
avoid degenerate ground states. We choose an instance with parameter values
denoted by the red dot. (b) The overlap defined in Eq. (37). The black line
indicates the fixed value of $U_{\mathrm{D}}$ for subplot (a) and the red dot
indicates the chosen set of parameters as before.

##  Appendix B Derivation of the effective Hamiltonian

The derivation essentially follows App. C in Ref. [40], generalized to two
dimensions and additionally including cavity-induced tunneling (see also e.g.
Ref. [54]). Including the decay of the cavity fields through the mirrors with
the rates $\kappa_{m}$, the full open-system dynamics in the tight binding
limit is given by the Lindblad equation

$\dot{\rho}=-\frac{i}{\hbar}[\tilde{H},\rho]+\sum_{m}\kappa_{m}(2a_{m}\rho
a_{m}^{\dagger}-a_{m}^{\dagger}a_{m}\rho-\rho a_{m}^{\dagger}a_{m}),$ (38)

where $\tilde{H}$ is specified in Eq. (10). The important term for the
following discussion discussion is the atom-light interaction
$H_{\mathrm{cav}}$ given in Eq. (11).

From this master equation we obtain the Heisenberg-Langevin equations of the
cavity fields

$\dot{a}_{m}=(i\tilde{\Delta}_{c,m}-\kappa_{m})a_{m}-iN\eta_{m}\hat{\Theta}_{m}+\sqrt{2\kappa_{m}}\hat{\xi}_{m}$
(39)

with the quantum noise operator $\hat{\xi}_{m}$ obeying
$[\hat{\xi}_{m}(t),\hat{\xi}_{m}^{\dagger}(t^{\prime})]=\delta(t-t^{\prime})$.

Assuming that the cavity mode fields evolve on a much faster time scale than
the atomic motion ($J/(\hbar|\tilde{\Delta}_{c,m}+i\kappa_{m}|)\ll 1$), they
can be approximated by their steady state on a course grained time scale [40,
57, 58, 54]. From the Heisenberg-Langevin equation, to zeroth order in
$J/(\hbar|\tilde{\Delta}_{c,m}+i\kappa_{m}|)$, we get

$a_{m}^{\mathrm{st}}\equiv\frac{\eta_{m}}{\tilde{\Delta}_{c,m}+i\kappa_{m}}N\hat{\Theta}_{m}+\sqrt{\frac{2\kappa_{m}}{\tilde{\Delta}_{c,m}^{2}+\kappa_{m}^{2}}}\hat{\xi}_{m}$
(40)

with
$\hat{\Theta}_{m}=\frac{1}{N}\sum_{i,j=1}^{N}\left(v_{m}^{ij}\hat{n}_{ij}+u_{m}^{ij}\hat{B}_{ij}\right)$.
That is, at steady state the effect of the field can be expressed by atomic
operators only.

We now substitute the cavity field operators by their corresponding steady-
state approximations in the Heisenberg equation of the atomic annihilation
operators

$\displaystyle\dot{b}_{ij}=\frac{1}{i\hbar}[b_{ij},H_{\mathrm{cav}}]+...=$
$\displaystyle-i\sum_{m}\frac{(N\eta_{m})^{2}}{\tilde{\Delta}_{c,m}^{2}+\kappa_{m}^{2}}\left[\tilde{\Delta}_{c,m}\left([b_{ij},\hat{\Theta}_{m}^{\dagger}]\hat{\Theta}_{m}+\hat{\Theta}_{m}^{\dagger}[b_{ij},\hat{\Theta}_{m}]\right)\right.$
(41) $\displaystyle\hskip
91.04872pt\left.-i\kappa_{m}\left([b_{ij},\hat{\Theta}_{m}^{\dagger}]\hat{\Theta}_{m}-\hat{\Theta}_{m}^{\dagger}[b_{ij},\hat{\Theta}_{m}]\right)\right]$
$\displaystyle-i\sum_{m}\frac{N\eta_{m}\sqrt{2\kappa_{m}}}{\sqrt{\tilde{\Delta}_{c,m}^{2}+\kappa_{m}^{2}}}\left([b_{ij},\hat{\Theta}_{m}^{\dagger}]\hat{\xi}_{m}+\hat{\xi}_{m}^{\dagger}[b_{ij},\hat{\Theta}_{m}]\right)$

where we only report terms including the cavity. At this point, ordering of
atomic and field operators becomes important, since
$a_{m}^{\mathrm{st}}\propto\hat{\Theta}_{m}$ as opposed to $a_{m}$ does not
necessarily commute with atomic operators. Here we choose normal ordering, as
already done in Eq. (11). The expression contains coherent terms proportional
to $\tilde{\Delta}_{c,m}$ and incoherent terms proportional to $\kappa_{m}$.

For $|\tilde{\Delta}_{c,m}|\gg\kappa_{m}$ we can neglect the incoherent part
and the Heisenberg equation can be obtained from

$\dot{b}_{ij}=\frac{1}{i\hbar}[b_{ij},H_{\mathrm{cav}}^{\mathrm{eff}}]+...$
(42)

Thus the dynamics in the coherent regime is described by the effective
Hamiltonian $H_{\mathrm{cav}}^{\mathrm{eff}}$ given in Eq. (16). Otherwise the
Heisenberg equation is equivalent to the master equation (28).

The results can also be obtained by naively substituting $a_{m}$ with
$a_{m}^{\mathrm{st}}$ directly in the Hamiltonian Eq. (11) or the master
equation (38) with the given ordering.

Note that the effective Hamiltonian can also be written in the form

$H_{\mathrm{cav}}^{\mathrm{eff}}=\hbar\sum_{m}\tilde{\Delta}_{c,m}(a_{m}^{\mathrm{st}})^{\dagger}a_{m}^{\mathrm{st}},$
(43)

which allows for a simple interpretation: For $\tilde{\Delta}_{c,m}>0$ the
lowest energy states tend to minimize the intensity of the cavity fields
$\langle(a_{m}^{\mathrm{st}})^{\dagger}a_{m}^{\mathrm{st}}\rangle$.

##  Appendix C Harmonic approximation of potential wells

In this section we investigate the limit of a deep lattice in more detail. In
Section 3.1 we presented results in the "infinitely" deep lattice limit, where
the Wannier functions become delta functions. To gain more insight to deep but
finite lattice depths, we use a harmonic approximation for the potential
wells. The ground state wave function is then an approximation to the lowest-
band Wannier function

$w_{\mathrm{har}}(x)=\pi^{-\frac{1}{4}}a_{0}^{-\frac{1}{2}}e^{-\frac{x^{2}}{2a_{0}^{2}}}$
(44)

with the size $a_{0}=(E_{\mathrm{R}}/V_{\mathrm{L}})^{1/4}/k_{\mathrm{L}}$
[55].

With this the atom-mode overlap integrals [Eq. (13)] can be calculated
analytically using running wave mode functions [Eq. (18)]. For the on-site
term we obtain

$v_{m}^{ij}=h_{m}(x_{i},y_{j})e^{-\left(\frac{k_{m}^{x}}{2k_{\mathrm{L}}}\right)^{2}\sqrt{\frac{E_{\mathrm{R}}}{V_{\mathrm{L}}}}}.$
(45)

It consists of the mode function at the lattice site and an exponential which
reduces the overlap due to Gaussian smoothing of the mode function. As
intuitively expected, the smoothing has a stronger effect for large mode wave
numbers $k_{m}^{x}$. For $V_{\mathrm{L}}/E_{\mathrm{R}}\gg 1$, we obtain
$v_{m}^{i,j}=h_{m}(x_{i},y_{j})$, as in the main text.

For the off-site overlaps we obtain

$\displaystyle u_{m}^{ij}=$ $\displaystyle
h_{m}((x_{i}+x_{i+1})/2,y_{j})e^{-\left(\frac{k_{m}^{x}}{2k_{\mathrm{L}}}\right)^{2}\sqrt{\frac{E_{\mathrm{R}}}{V_{\mathrm{L}}}}}e^{-\frac{\pi^{2}}{4}\sqrt{\frac{V_{\mathrm{L}}}{E_{\mathrm{R}}}}}.$
(46)

The overlap consists of three terms: First, it is the mode function evaluated
in between the lattice sites. Second, there is again the Gaussian smoothing
term as for the on-site overlap. Lastly, there is an exponential independent
of the modes, which comes from the overlap of the two Gaussians. It goes to
zero for $V_{\mathrm{L}}/E_{\mathrm{R}}\gg 1$, leading to $u_{m}^{ij}=0$.

The order operator is then

$\hat{\Theta}_{m}^{\mathrm{har}}=\frac{1}{N}e^{-\left(\frac{k_{m}^{x}}{2k_{\mathrm{L}}}\right)^{2}\sqrt{\frac{E_{\mathrm{R}}}{V_{\mathrm{L}}}}}\sum_{i,j=1}^{N}\left(h_{m}(x_{i},y_{j})\hat{n}_{ij}+h_{m}((x_{i}+x_{i+1})/2,y_{j})\hat{B}_{ij}e^{-\frac{\pi^{2}}{4}\sqrt{\frac{V_{\mathrm{L}}}{E_{\mathrm{R}}}}}\right)$
(47)

leading to an interaction Hamiltonian [Eq. (16)] given by

$\displaystyle H_{\mathrm{cav}}^{\mathrm{har}}=$ $\displaystyle
U_{Q}\sum_{m}f_{m}e^{-2\left(\frac{k_{m}^{x}}{2k_{\mathrm{L}}}\right)^{2}\sqrt{\frac{E_{\mathrm{R}}}{V_{\mathrm{L}}}}}\sum_{ijkl}\left(h_{m}^{*}(x_{i},y_{j})\hat{n}_{ij}+h_{m}^{*}((x_{i}+x_{i+1})/2,y_{j})\hat{B}_{ij}e^{-\frac{\pi^{2}}{4}\sqrt{\frac{V_{\mathrm{L}}}{E_{\mathrm{R}}}}}\right)$
(48) $\displaystyle\hskip
118.07875pt\times\left(h_{m}(x_{k},y_{l})\hat{n}_{kl}+h_{m}((x_{k}+x_{k+1})/2,y_{l})\hat{B}_{kl}e^{-\frac{\pi^{2}}{4}\sqrt{\frac{V_{\mathrm{L}}}{E_{\mathrm{R}}}}}\right),$

which in the "infinitely" deep lattice limit simplifies to Eq. (17). All
cavity-induced tunneling terms are suppressed by the exponential and tend to
be smaller than density-density terms. Also, since $U_{\mathrm{Q}}$ is
maximally on the order of $J$ (at the end of the sweep), cavity-induced
tunneling terms are smaller than $H_{\mathrm{kin}}$. However, also the
density-density terms can be small for example when $h_{m}(x_{i},y_{j})=0$,
which is why we still include cavity-induced tunneling in the simulations.

In the main text we chose uniform $f_{m}=1/M$. To compensate for Gaussian
smoothing one might want to include the exponential as correction

$\tilde{f}_{m}=f_{m}e^{2\left(\frac{k_{m}^{x}}{2k_{\mathrm{L}}}\right)^{2}\sqrt{\frac{E_{\mathrm{R}}}{V_{\mathrm{L}}}}},$
(49)

which leads to even better results ($\tilde{A}$ is closer to $A$ for finite
lattice depths). Note that this correction does only depend on
$V_{\mathrm{L}}$ and not on the problem size or number of modes in our
implementation, since the range of $k_{m}^{x}$ is fixed.

##  Appendix D Shape of the interaction

Here we reformulate the interaction in the infinitely deep lattice limit from
Eq. 22 with Fourier transforms by defining a real envelope function $f(k)$
such that $f(k_{m})=f_{m}$. To get back the discrete wave numbers, this
function is sampled with a Dirac comb at the lines $m\Delta k+k_{s}$ with
$m\in\mathbb{Z}$, where $k_{s}$ is a constant shift and $\Delta k$ is the
spacing between the pumped modes. In the main text we only consider the case
$k_{s}=0$ for simplicity. We define the Fourier transform as
$\mathcal{F}\\{f\\}(r)=\int_{-\infty}^{\infty}\mathrm{d}kf(k)e^{ikr}/\sqrt{2\pi}$
and denote the convolution as
$(f*g)(t)=\int_{-\infty}^{\infty}f(\tau)g(t-\tau)\mathrm{d}\tau$.

For simplicity, we take parallel wave vectors $\boldsymbol{k}_{m}$. Along this
direction $\boldsymbol{r}\parallel\boldsymbol{k}_{m}$ we write

$\displaystyle\tilde{\mathcal{A}}(r)$
$\displaystyle=\sum_{m}f_{m}\cos(k_{m}r)=\int_{-\infty}^{\infty}\mathrm{d}k\,f(k)\sum_{m=-\infty}^{\infty}\delta(k-m\Delta
k-k_{s})\cos(kr)$ (50)
$\displaystyle=\mathrm{Re}\left[\sqrt{2\pi}\mathcal{F}\left\\{f(k)\sum_{m=-\infty}^{\infty}\delta(k-m\Delta
k-k_{s})\right\\}(r)\right]=\mathrm{Re}\left[\sqrt{2\pi}\mathcal{F}\\{f\\}(r)*\sum_{l=-\infty}^{\infty}\delta\left(r-l\frac{2\pi}{\Delta
k}\right)e^{ik_{s}r}\right]$
$\displaystyle=\sum_{l=-\infty}^{\infty}\mathrm{Re}\left[\sqrt{2\pi}\mathcal{F}\\{f\\}\left(r-l\frac{2\pi}{\Delta
k}\right)e^{i2\pi l\frac{k_{s}}{\Delta k}}\right],$

where $r=|\boldsymbol{r}|$. We used the convolution and shift theorem from
Fourier analysis in the second to last line and evaluated the convolution
integrals by pulling out the sum in the last line.

For symmetric envelopes centered around $k_{c}$ we can further simplify using
the shift theorem

$\tilde{\mathcal{A}}(r)=\sum_{l=-\infty}^{\infty}\sqrt{2\pi}\mathcal{F}\\{\tilde{f}\\}\left(r-l\frac{2\pi}{\Delta
k}\right)\cos\left(k_{c}\left(r-l\frac{2\pi}{\Delta k}\right)+2\pi
l\frac{k_{s}}{\Delta k}\right),$ (51)

where $\tilde{f}(k)=f(k+k_{c})$ is the shifted envelope centered around $k=0$,
whose Fourier transform is real.

Let us apply this to our example described in Section 3.2 and find out why the
interaction has the desired property given in equation Eq. (24). There we had
uniform $f_{m}=1/M$ and wave numbers

$k_{m}^{0}=k_{\mathrm{L}}\left(1+\frac{2m+1}{2M}\right)$ (52)

with $m=0,...,M-1$. These have a mode spacing of $\Delta k=k_{\mathrm{L}}/M$
and are centered around $k_{c}=3k_{\mathrm{L}}/2$. One can see that only the
odd modes of the cavity (wave numbers $k_{n}=n\Delta k_{\mathrm{FSR}}$ with
$n$ odd and free spectral range $\Delta k_{\mathrm{FSR}}$) are used.
Therefore, $\Delta k=2\Delta k_{\mathrm{FSR}}$ and $k_{s}=\Delta k/2$, because
the comb has to be shifted to fit the odd modes. Due to the uniform $f_{m}$
the envelope is a rectangular function with width $k_{\mathrm{L}}$ and height
$1/M$ centered at $k_{c}$

$f(k)=\mathrm{rect}((k-k_{c})/k_{\mathrm{L}})/M=\begin{cases}1/M&\mathrm{for}\;k\in[k_{c}-k_{\mathrm{L}}/2,k_{c}+k_{\mathrm{L}}/2]\\\
0&\mathrm{otherwise}.\end{cases}$ (53)

The Fourier transform of a rectangular function centered around zero with unit
width and height is a sinc function $\mathrm{sinc}(x)=\sin(x)/x$. Using the
addition theorem for the cosine and noting that $\cos(\pi l)=(-1)^{l}$ and
$\sin(\pi l)=0$ we obtain an analytical expression for the interaction

$\tilde{\mathcal{A}}(r)=\sum_{l=-\infty}^{\infty}(-1)^{l}\mathrm{sinc}\left(\frac{k_{\mathrm{L}}}{2}\left(r-l\frac{2\pi}{\Delta
k}\right)\right)\cos\left(k_{c}\left(r-l\frac{2\pi}{\Delta k}\right)\right).$
(54)

For $l=0$ and at lattice site spacings $r_{j}=j\pi/k_{\mathrm{L}}$ it takes
the values

$\mathrm{sinc}(\pi j/2)\cos(3\pi j/2)=\begin{cases}1\quad\mathrm{for}\;j=0\\\
0\quad\mathrm{otherwise},\end{cases}$ (55)

as desired. This comes from well known properties of sinc and cosine

$\displaystyle\mathrm{sinc}(\pi j/2)$
$\displaystyle=\begin{cases}1&\mathrm{for}\;j=0\\\
0&\mathrm{for}\;j\;\mathrm{even}\\\ (-1)^{\frac{j-1}{2}}\frac{2}{\pi
j}&\mathrm{for}\;j\;\mathrm{odd}\\\ \end{cases}$ $\displaystyle\cos(3\pi j/2)$
$\displaystyle=\begin{cases}(-1)^{\frac{j}{2}}\quad&\mathrm{for}\;j\;\mathrm{even}\\\
0\quad&\mathrm{for}\;j\;\mathrm{odd}.\end{cases}$

The other summands have the same form, but are shifted by $R=2\pi/\Delta
k=2\pi M/k_{\mathrm{L}}=r_{2M}$ ($2M$ lattice sites) and have alternating
signs. Since $R$ is an integer multiple of the lattice spacing this adds up to
the desired interaction given in Eq. (24) in the main text.

Thus for rectangle envelopes the bandwidth $\Delta k_{\mathrm{BW}}$ determines
the zeros of the interaction. Taking $\Delta k_{\mathrm{BW}}=2k_{\mathrm{L}}$
would lead to zeros at all lattice sites. For the smaller bandwidth $\Delta
k_{\mathrm{BW}}=k_{\mathrm{L}}$ used here, only even sites become zero. This
can be compensated by choosing a central wave number $k_{c}=nk_{\mathrm{L}}/2$
with $n$ odd, which is responsible for the zeros at odd sites. The mode
spacing $\Delta k$ determines the peak distance. Finally, using odd cavity
modes (specifying $k_{s}$) leads to alternating peaks, which does not have an
effect in our implementation, since $-2M<j<2M$.