A Quantum N-Queens Solver

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. The problem is proven NP-complete and for the excluded-diagonals variation the parameter regime to generate hard instances that are intractable with current classical algorithms are known. We propose a special purpose quantum simulator with the aim to solve the excluded diagonals N-queens completion problem using atoms 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.


I. 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]. 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 [10]. Quantum simulation with the aim to solve Hamiltonian systems may serve as a building block of such a special purpose quantum computer [11][12][13]. In particular, adiabatic quantum computing [14][15][16] has been proposed to solve computationally hard problems by finding the ground state of Ising spin glasses [17]. The speedup is expected to be polynomial and not exponential, and the question whether these protocols can show quantum speedup at all is still open. Thus, demonstrating quantum advantage by solving optimization problems using quantum simulation tools is a crucial step towards the development of general programmable quantum optimizers [18,19].
Here, we present a scheme that aims at solving the Nqueens problem, and variations of it, using atoms with cavity-mediated long-range interactions [20][21][22][23][24]. Our proposed setup consists of N ultracold atoms in an optical lattice that represent the queens on the chess board [25]. The non-attacking conditions are enforced by a combination of restricted hopping and interactions between the atoms stemming from collective scattering of pump laser light into a multi-mode cavity [26][27][28][29][30][31][32][33] (see Fig. 1). For the excluded diagonals variation of the Nqueens problem, additional repulsive optical potentials * w.lechner@uibk.ac.at 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). 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 continuous measurement of photons that leave the cavity [34] 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 vi-able candidate to test quantum advantage in near term devices. (a) The problem is proven to be NP-complete and hard instances are known from computer science literature [35], (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. [35]) is available in the lab. As longrange interactions are implemented with infinite-range cavity-mediated forces the intermediate step of embedding the optimization problem in an Ising problem is removed. With methods such as minor embedding [19,36], LHZ [18,37,38] or nested embedding [39] there is at least a quadratic overhead for all-to-all connected models as the N -queens problem. Thus, due to the cavity mediated long range interactions, 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 [35]. Finally, verification is computationally trivial as the final state is classical and no quantum tomography is needed. 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. II 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. III. In Sec. IV, we discuss that light leaking out of the cavity can be used for read-out. Finally, this work contains a numerical comparison between model and implementation in Sec. V and we conclude in Sec. VI.

II. QUANTUM SIMULATION OF THE N -QUEENS PROBLEM
Following the idea of adiabatic quantum computation [14][15][16], we construct a classical problem Hamiltonian H 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 timedependent Hamiltonian from t = 0 to t = τ . Initially at t = 0, the system is prepared in the ground state of H(0) = H 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(τ ) at t = τ . If the lowest energy gap of H pr is much larger than the one of H kin , this state is close to the ground state of H pr and thus the solution of the optimization problem. In the following we construct the problem Hamiltonian H pr and the driver Hamiltonian H kin . The system is modeled as a 2D Bose-Hubbard model with annihilation (creation) operators b ij (b † ij ) 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 xdirection j = l, (ii) the y-direction i = k, the diagonals 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 apriori knowledge that a solution has one queen in a row, which reduces the accessible configuration space size from N 2 N to N N configurations. In this vein, the tunneling Hamiltonian is given by where J is the tunneling amplitude 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 For U 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 D + and D − , respectively, and the coefficients are For U T > 0, the second term favors occupations of certain sites. The sites where queens should be pinned to are pooled in the set T and therefore the coefficients are given by The problem Hamiltonian of the N -queens problem with excluded diagonals is then 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 [25,33], 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 n ij 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 kin is a superposition of excitations along each tube. After evolving for a sufficiently large time Jτ / = 49, the system is in the ground state of H pr and thus assumed the solution of the optimization problem.
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 τ 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.

III. IMPLEMENTATION
In this section we propose a physical implementation of Eq. (1) based on ultracold atoms in an anisotropic two-dimensional optical lattice, where tunneling is suppressed in one dimension [40]. Thus the atoms can move in tubes along the x-direction only, effectively creating an array of N parallel 1D optical lattices, each filled with a single atom. For implementing the queens interactions the optical lattice is placed inside the transverse plane of a multi-mode standing-wave resonator (see Fig.  1(b)). The trapped atoms are illuminated by running wave laser beams with multiple frequencies from three different directions within the transverse plane. Such a multi-frequency laser beam can be created by a frequency comb. Since all light frequencies are well separated compared to the cavity line width, each is scattered into a distinct cavity mode and the cavity modes are not directly coupled [33]. This collective positiondependent scattering into the cavity introduces effective 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. Additionally, light sheets and optical tweezers can be used to make certain diagonals energetically unfavorable (for the excluded diagonal problem) or pin certain queens.
In the first part (Sec. III A), we will derive a Bose-Hubbard-type Hamiltonian for the trapped atoms interacting via collective scattering for general illumination fields. For this we first consider the full system of coupled atoms and cavity modes, and later adiabatically eliminate the cavity fields to obtain a light-mediated atomatom interaction.
In the second part (Sec. III B) we discuss how to implement the non-attacking conditions of the N -queens problem with this light-mediated atom-atom interaction. Specifically, we consider the limit of a deep optical lattice leading to simple analytical expressions. These formulas allow us to find an example of a pump configuration with specific wave numbers and pump strengths leading to the queens interaction Hamiltonian in Eq. (3). However, increasing the lattice depth slows down atomic movement, which leads to huge sweep times. Luckily, it turns out that one can still use the pump configuration obtained in the deep lattice limit also for a moderate lattice depth allowing for tunneling. We show with numerical simulations in Sec. V that while the interaction is altered, its shape does not significantly change and it still well approximates the N -queens interaction.

A. 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 [41].
Single-particle Hamiltonian. For a single particle of mass m A , the motion of the atoms in the x-y-plane is described by [21,42] The first line contains the kinetic term with the momentum operatorsp x andp y . Classical electric fields create optical potentials with depths V x L , V y L and V bias . The first two create the optical lattice with wave number k L and lattice spacing a = π/k L , while V bias is much smaller and only responsible for a bias field on certain sites, for instance for excluding diagonals. Thereby 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 ), 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 η m = g m Ω m /∆ a,m , with the pump laser Rabi frequencies Ω m , the atom-cavity couplings g m and the detunings between pump lasers and atomic resonance frequency ∆ a,m . The effective cavity detunings∆ c,m = ∆ c,m − N U 0,m are given by the detunings between pump laser and cavity mode frequencies ∆ c,m and the dispersive shifts of the cavity resonance due to the presence of the atoms in the cavity N U 0,m [20].
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 TEM 00 modes, the atomcavity 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 atomatom interactions are taken into account by introducing bosonic field operatorsΨ( [21,43]. 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 x L and V y L 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 [42]. In this limit we can expand the bosonic field operators in a localized Wannier basiŝ Ψ(x) = i,j w 2D (x − x ij )b ij with the lowest-band Wannier functions w 2D (x) coming from Bloch wave functions of the lattice [44]. We split the resulting Hamiltonian in three termsH which will be explained in the following. As in the standard Bose-Hubbard model, one obtains a tunneling term H kin as in Eq. (2). Tunneling in ydirection is frozen out by ensuring V y The other terms H cav and H pot originate from the weak cavitypump interference fields and the bias fields introduced above and should resemble H 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 cav and H pot larger) or the lattice depth (make H kin smaller).
The cavity-related terms in Eq. (9) give rise to with the order operator of cavity mode m The structure of the fields enters in the on-site and nearest neighbor atom-mode overlaps 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 y L V x L we can approximate the y-dependence of the Wannier functions by a Dirac delta: The last term H pot describes all extra fields responsible for local energy off-sets at certain sites that stem from the weak classical fields with the distribution F(x, y) The off-sets ought to be calculated from the overlaps of fields and Wannier functions, analogously to v ij m . 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 pot for numerical simulations.
Atom-atom interaction Hamiltonian. The main focus of this work is to show how to create the tailored allto-all particle interactions via collective scattering. We derive this interaction by eliminating the cavity fields introduced in the previous section [42,45,46]. This can be done because the cavity fields decay through the mirrors with the rates κ m and thus end up in a particular steadystate for each atomic configuration. Assuming that the atomic motion is much slower than the cavity field dynamics, i.e. J/ |∆ c,m + iκ m |, this steady-state is a good approximation at all times. The stationary cavity field amplitudes are given by and thus replaced by atomic operators (see Appendix B).
In the coherent regime |∆ c,m | κ m , the atom-light interaction is then described by an effective interaction Hamiltonian for the atoms [24] The collective, state-dependent scattering induces interactions between each pair of sites (i, j) and (k, l): Density-density interactions due to the terms containinĝ n ijnkl and a modified tunneling amplitude due an occupation or a tunneling event somewhere else in the lattice due ton ijBkl ,B ijnkl andB ijBkl .
Since the Wannier functions are localized at the lattice sites, the on-site overlaps v ij m tend to be much larger than the nearest-neighbor overlaps u ij m . Thus density-density interactions are expected to be the dominant contribution to H eff cav . 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 ∝ (V x L ) −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 x L ) when the width tends to zero we get v ij m = h m (x i , y j ) and u ij m = 0 [24], and an interaction Hamiltonian [from Eq. (16)] which only depends on density operators, and hence does not include cavity induced-tunneling.

B. 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 Q [Eq. (3)] containing the non-attacking conditions. Using these pump fields we later show numerically in Sec. V, 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 k respectively, with the wave numbers k 0 m . Therefore the pump fields are given by where x = (x, y) T is the position vector and k m is a wave vector in any of the three directions.
With running wave pump fields, Eq. (17) can be written as This formally corresponds to H Q [Eq. (3)], where the quantities now have a physical meaning: The interaction matrix is given bỹ with lattice site connection vectors x ij − x kl and 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Ã approximates A [Eq. (4)]. The overall strength of the interaction term is captured by the energy U 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 functioñ which returns the interaction matrix when evaluated at lattice site connection vectorsÃ ijkl =Ã(x ij − x kl ). We note that one set of parallel k µ m (µ ∈ {x, +, −}) creates an interactionÃ which is constant and infinite range (only limited by the laser beam waist) in the direction perpendicular to the propagation direction r ⊥ k µ m . Along the propagation direction r k µ m instead, the interaction is shaped according to the sum of cosines, and can be modified by the choice of wave numbers k µ m = |k µ m | and their relative strengths f m .
In the following we will use the example wave numbers 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 Q > 0 and thus∆ c,m > 0), the wave vectors k x m create the non-attacking interaction along the y-direction (Ã ijkl = 1 if i = k and 0 otherwise) as long as N ≤ 2M . This is illustrated in Fig. 4(a) for M = N = 5. Analogously, k ± m cause the non-attacking interactions along the diagonals. In a square lattice the diagonals have the distance r j / √ 2, which is compensated by k ± m = |k ± m | = √ 2k 0 m . Since there are 2N − 1 diagonals, one has to make sure that 2N − 1 ≤ 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 tot = 3M = 3N frequencies in our example.
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∆k with n ∈ Z containing all k µ m . Considering one illumination direction for simplicity, the interaction [Eq. (22)] along r k µ m with r = |r| can be written as 2π is the Fourier transform of f (k) and δ(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π/∆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 σ ∼ 2π/∆k BW .
For the (approximate) non-attacking condition (Ã(ja) ≈ 1 for j = 0 and |Ã(ja)| 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 ≥ N a (R ≥ √ 2N a). Secondly, the width of one peak has to be smaller than the lattice spacing σ a. Combining these conditions to R N σ, we see that the minimum number of modes per direction M scales linearly with N M ≈ ∆k BW /∆k N.
Therefore, with only ∼ 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 ∆k BW k L .

IV. 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 [47,48]. 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, the lattice depth can be in principle increased to some high value in the deep lattice regime.

A. 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 |ψ In the last line we assumed that ζ = 2κ m /∆ 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 = 3N U Q ζ. 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Ã leads to an increase of the photon flux by ∆P = 2U Q ζ. The two atoms create an energy penalty for one another, explaining the factor 2. The relative difference of photon flux due to a state with L attacking pairs and a solution is given by 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 κ m /∆ c,m , photons from different modes have to be distinguished.

B. Field measurement
More information can be gained from a measurement of the output fields (see also Fig. 6). The phase of the scattered light gives insight about the absolute position of the atoms projected onto the pump laser propagation direction. The output field phase can be measured by interference with reference lasers, for example by homodyne detection.
Let us illustrate this by considering only light scattered from the x-direction with incident wave vectors k x m . Neglecting cavity-assisted tunneling, the field quadratures for a phase difference φ are Here we used that for plane wave pumps in x-direction, the atom-field overlaps do not depend on j. Thus the cavity fields are determined by the total occupation on a vertical lineN x i = jn ij . For at least N modes (M ≥ N ) this system of equations can be inverted yielding the line occupations N x i . By measuring the cavity output fields scattered from the diagonal pump light we obtain the occupations of each diagonal N + i and N − i . Inverting the system of equations for diagonals demands at least as many pump modes as diagonals, that is 2N − 1.
A solution of the queens problem has maximally one atom on each diagonal and exactly one atom on each vertical line. Thus a necessary condition is This condition is also sufficient for classical configurations, like occupation number basis states |φ ν . However, some superpositions |ψ = ν c ν |φ ν which are no solutions might also fulfill the above criterion, because summands in the field expectation values a st m = ν |c ν | 2 φ ν |a st m |φ ν can cancel each other. For instance, for U T = 0 the solution from our example in Fig. 2 |ψ sol = |1, 4, 2, 5, 3 scatters the same fields a st m as the superposition |ψ nosol = (|ψ 1 nosol + |ψ 2 nosol )/ √ 2 with |ψ 1 nosol = |1, 3, 2, 5, 4 and |ψ 2 nosol = |1, 4, 5, 2, 3 , both of which are no solution. In this notation the state |i 1 , i 2 , ..., i N has one atom on each site (i j , j). However, these macroscopic superpositions are highly unstable. Even theoretically the measurement back-action [34,49] projects superpositions of states scattering different fields (such as |ψ nosol ) to one of its constituents. The inclusion of measurement back-action and noise due to photon loss might thus lead to intriguing phenomena beyond those presented here and is subject to future work.
For classical states the above measurement determines if we found a solution or not, which is the answer to the combinatorial decision problem. But it does not contain information about the configuration of atoms. This can be measured with single site resolution as demonstrated in several experiments [47,48].

V. NUMERICAL JUSTIFICATION OF ASSUMPTIONS
We compare the ideal model Hamiltonian described in Sec. II [Eq. (1)] to the physically motivated tight-binding Hamiltonian for finite lattice depths introduced in Sec. III [Eq. (10)]. Like for the ideal model in Fig. 2, we consider the time evolution during a slow linear sweep of U Q , U T and U D by numerically integrating the timedependent 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.
In the following we use a realistic lattice depth of V x L = 10E R with the recoil energy E R = 2 k 2 L /(2m A ). For example, for rubidium 87 Rb and λ L = 785.3 nm it is E R / = 23.4 kHz [22]. The chosen lattice depth leads to a tunneling amplitude J ≈ 0.02E 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. III B 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.

A. Time evolution for finite lattice depths
The energy spectrum 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 lat-  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. (31)]. 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).
tices 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τ / = 49 for the same parameters by integrating the time-dependent Schrödinger equation [Eq. (10)]. Snapshots of the site occupations n ij for several times are shown in Fig. 5(a), where we observe a similar behavior as for the ideal case (Fig. 2). This suggests that the system is robust against the errors introduced by the moderate lattice depth.
During the time evolution, the cavity output gives insight about the atomic state. Figure 6(a) shows the cavity field expectation values evolving in the complex plane for the uniform cavity parameters κ m = 400E R / and ∆ c,m = 4κ m . The time evolution of the rate of photons at a detector is shown in Fig. 6(b). They vary around a linear increase, which originates in the linear pump intensity ramp of our protocol. The final values are depicted in Fig. 6(c). We observe that the light scattered from k x m -direction has the same intensity distribution as the light scattered from the k − m -direction. This is because the atomic distribution projected on the respective directions has the same form. The fields and intensities can be used to determine if a state is a solution or not, as discussed in Sec. IV. The phase information contained in the cavity fields additionally gives insight about occupations on vertical lines and diagonals.

B. Classical cavity fields
Here we substitute the field operator a st m by its expectation value, that is, we consider classical cavity fields. Instead of the full quantum Hamiltonian shown in Eq. (16) we obtain the semi-classical Hamiltonian where the expectation values have to be calculated selfconsistently with the current atom state vector. This substitution amounts to considering only first order fluctuations around the mean ofΘ m in Eq. (16).
Consequently, the dynamics are described by a differential equation which is non-linear in the state vector |ψ . 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τ / = 49 is depicted in Fig. 5(b). The discrepancy shows the necessity of quantum effects in our procedure.

VI. 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 testbed to study possible quantum advantage in solving classical combinatorial problems 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-ofthe-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 cannot be neglected.
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 [47].
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 [20,50] 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 [34].   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 Q , the final trapping energy U T , the tunneling amplitude J and the final excluded diagonals penalty U 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 between the state |φ 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 |ψ (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 Q , U D and U T . Figure 7(a) suggests that in order to increase the minimal gap the ratio U Q /U D has to be chosen as small as possible. We vary the ratio by fixing U D and varying U Q . Therewith, Fig. 7(a) indicates that U Q should be as small as possible. However, as it can be seen in Fig. 7(b), a small U 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 D to 5J and U T to 2J we find that for U Q = 1J the overlap is F ≈ 0.93 and the minimal gap is around 0.44J. These values were used for Figs. 2 and 3.
Note that the effective Hamiltonian can also be written in the form which allows for a simple interpretation: For∆ c,m > 0 the lowest energy states tend to minimize the intensity of the cavity fields (a st m ) † a st m .
Appendix C: Harmonic approximation of potential wells In this section we investigate the limit of a deep lattice in more detail. In Section III A 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 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 x m . For V L /E R 1, we obtain v i,j m = h m (x i , y j ), as in the main text. For the off-site overlaps we obtain 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 L /E R 1, leading to u ij m = 0. The order operator is then leading to an interaction Hamiltonian [Eq. (16)] given by 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 Q is maximally on the order of J (at the end of the sweep), cavity-induced tunneling terms are smaller than H 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 correctionf which leads to even better results (Ã is closer to A for finite lattice depths). Note that this correction does only depend on V L and not on the problem size or number of modes in our implementation, since the range of k x m is fixed.
The other summands have the same form, but are shifted by R = 2π/∆k = 2πM/k 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 ∆k BW determines the zeros of the interaction. Taking ∆k BW = 2k L would lead to zeros at all lattice sites. For the smaller bandwidth ∆k BW = k L used here, only even sites become zero. This can be compensated by choosing a central wave number k c = nk L /2 with n odd, which is responsible for the zeros at odd sites. The mode spacing ∆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 .