QUANTUM COMPUTING FOR COMBINATORIAL OPTIMIZATION PROBLEMS USING PROGRAMMABLE ATOM ARRAYS

CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This application claims the benefit of priority to U.S. Provisional Application No. 62/725,874, entitled“QUANTUM OPTIMIZATION FOR MAXIMUM INDEPENDENT SET USING RYDBERG ATOM ARRAYS,” filed on August 31, 2018, the disclosure of which is hereby incorporated by reference in its entirety.

COPYRIGHT NOTICE

[0002] This patent disclosure may contain material that is subject to copyright protection. The copyright owner has no objection to the facsimile reproduction by anyone of the patent document or the patent disclosure as it appears in the U.S. Patent and Trademark Office patent file or records, but otherwise reserves any and all copyright rights.

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR

DEVELOPMENT

[0003] This invention was made with government support under Grant Nos. 1506284, PHY-l 125846, and PHY-1521560 awarded by the National Science Foundation; FA9550-17- 1-0002 awarded by the U.S. Air Force Office of Scientific Research; and N00014-15-1-2846 awarded by the U.S. Department of Defense/Office of Navy Research. The government has certain rights in the invention.

TECHNICAL FIELD

[0004] This patent relates to quantum computing, and more specifically to preparing and evolving an array of atoms for quantum computations. BACKGROUND

[0005] As quantum simulators, fully controlled, coherent many-body quantum systems can provide unique insights into strongly correlated quantum systems and the role of quantum entanglement and enable realizations and studies of new states of matter, even away from equilibrium. These systems also form the basis for the realization of quantum information processors. While basic building blocks of such processors have been demonstrated in systems of a few coupled qubits, increasing the number of coherently coupled qubits to perform tasks that are beyond the reach of modern classical machines is challenging.

Furthermore, existing systems lack coherence and/or quantum nonlinearity for achieving fully quantum dynamics.

[0006] Neutral atoms can serve as building blocks for large-scale quantum systems, as described in more detail in PCT Application No. PCT/US 18/42080, titled“NEUTRAL ATOM QUANTUM INFORMATION PROCESSOR.” They can be well isolated from the environment, enabling long-lived quantum memories. Initialization, control, and read-out of their internal and motional states is accomplished by resonance methods developed over the past four decades. Arrays with a large number of identical atoms can be rapidly assembled while maintaining single-atom optical control. These bottom -up approaches are

complementary to the methods involving optical lattices loaded with ultracold atoms prepared via evaporative cooling, and generally result in atom separations of several micrometers. Controllable interactions between the atoms can be introduced to utilize these arrays for quantum simulation and quantum information processing. This can be achieved by coherent coupling to highly excited Rydberg states, which exhibit strong, long-range interactions. This approach provides a powerful platform for many applications, including fast multi-qubit quantum gates, quantum simulations of Ising-type spin models with up to 250 spins, and the study of collective behavior in mesoscopic ensembles. Short coherence times and relatively low gate fidelities associated with such Rydberg excitations are challenging. This imperfect coherence can limit the quality of quantum simulations and can dim the prospects for neutral atom quantum information processing. The limited coherence becomes apparent even at the level of single isolated atomic qubits.

[0007] PCT/US 18/42080 describes exemplary methods and systems for quantum computing. These systems and methods can involve first trapping individual atoms and arranging them into particular geometric configurations of multiple atoms, for example, using acousto-optic deflectors. This precise placement of individual atoms assists in encoding a quantum computing problem. Next, one or more of the arranged atoms may be excited into a Rydberg state, which can produce interactions between the atoms in the array. After, the system may be evolved under a controlled environment. Finally, the state of the atoms may be read out in order to observe the solution to the encoded problem. Additional examples include providing a high fidelity and coherent control of the assembled array of atoms.

SUMMARY

[0008] In one or more embodiments, a method includes selectively arranging a plurality of qubits into a spatial structure to encode a quantum computing problem, wherein each qubit corresponds to a vertex in the quantum computing problem, and wherein spatial proximity of the qubits represents edges in the quantum computing problem; initializing the plurality of qubits into an initial state; driving the plurality of qubits into a final state by applying a sequence of resonant light pulses with a variable duration and a variable optical phase to at least some of the plurality of qubits, wherein the final state comprises a solution to the quantum computing problem; and measuring at least some of the plurality of qubits in the final state.

[0009] In one or more embodiments, the spatial structure comprises a one-dimensional, two-dimensional or three-dimensional array of qubits. [0010] In one or more embodiments, the encoded quantum computing problem comprises one or more of an unweighted maximum independent set problem, a maximum-weight independent set problem, a maximum clique problem, and a minimum vertex cover problem.

[0011] In one or more embodiments, weights in the maximum-weight independent set problem are encoded by applying light shifts to at least some of the plurality of qubits.

[0012] In one or more embodiments, the final state of the plurality of qubits comprises one or more of a solution to the encoded unweighted maximum independent set problem, a solution to the encoded maximum-weight independent set problem, a solution to the encoded maximum clique problem, and a solution to the encoded minimum vertex cover problem.

[0013] In one or more embodiments, the solution to the quantum computing problem comprises an approximate solution to the quantum computing problem.

[0014] In one or more embodiments, a method includes selectively arranging a plurality of qubits into a spatial structure comprising a plurality of vertex qubits and a plurality of ancillary qubits to encode a quantum computing problem using spatial proximity of the plurality of qubits, wherein each vertex qubit corresponds to a vertex in the quantum computing problem and wherein subsets of the ancillary qubits correspond to edges in the quantum computing problem; initializing the plurality of qubits into an initial state; driving the plurality of qubits into a final state, wherein the final state comprises a solution to the quantum computing problem; and measuring at least some of the plurality of qubits in the final state.

[0015] In one or more embodiments, the driving the plurality of qubits into the final state comprises applying light pulses with a constant or variable Rabi frequency W and a constant or variable detuning D to at least some of the plurality of qubits.

[0016] In one or more embodiments, the applying light pulses to the at least some of the plurality of qubits further includes: applying at least one light pulse with a detuning D ₀ to a vertex qubit comprising a corner vertex or a junction vertex; and applying at least one light pulse with a detuning D _£ to each of i ancillary qubits adjacent to the vertex qubit on an edge connected to the vertex qubit.

[0017] In one or more embodiments, the applying the light pulses to the at least some of the plurality of qubits further comprises applying light shifts to selected qubits of the at least some of the plurality of qubits.

[0018] In one or more embodiments, the driving the plurality of qubits into the final state comprises applying a sequence of resonant light pulses with a variable duration and a variable optical phase to at least some of the plurality of qubits.

[0019] In one or more embodiments, the arranging the plurality of qubits into the plurality of vertex qubits and the plurality of ancillary qubits comprises arranging the plurality of qubits onto a grid.

[0020] In one or more embodiments, the encoded quantum computing problem comprises one or more of an unweighted maximum independent set problem, a maximum-weight independent set problem, a maximum clique problem, and a minimum vertex cover problem.

[0021] In one or more embodiments, weights in the maximum-weight independent set problem are encoded by applying light shifts to a plurality of qubits.

[0022] In one or more embodiments, the final state of the plurality of qubits comprises one or more of a solution to the encoded unweighted maximum independent set problem, a solution to the encoded maximum-weight independent set problem, a solution to the encoded maximum clique problem, and a solution to the encoded minimum vertex cover problem.

[0023] In one or more embodiments, the method further includes renumbering at least two vertices in the quantum computing problem prior to the encoding the quantum computing problem.

[0024] In one or more embodiments, the solution to the quantum computing problem comprises an approximate solution to the quantum computing problem.

[0025] In one or more embodiments, a method includes: selectively arranging a plurality of qubits into a spatial structure to encode a quantum computing problem, wherein each qubit corresponds to a vertex in the quantum computing problem; initializing the plurality of qubits into an initial state; stroboscopically driving the plurality of qubits into a final state, wherein the final state comprises a solution to the quantum computing problem; and measuring at least some of the plurality of qubits in the final state.

[0026] In one or more embodiments, stroboscopically driving the plurality of qubits into a final state comprises applying light pulses sequentially and selectively in an order to subsets of the plurality of qubits, the order of light pulses corresponding to the graph structure of the quantum computing problem.

[0027] In one or more embodiments, the driving the plurality of qubits into the final state comprises applying light pulses with a constant or variable Rabi frequency W and a constant or variable detuning D to at least some of the plurality of qubits.

[0028] In one or more embodiments, the driving the plurality of qubits into the final state comprises applying a sequence of resonant light pulses with a variable duration and a variable optical phase to at least some of the plurality of qubits.

[0029] In one or more embodiments, the encoded quantum computing problem comprises one or more of an unweighted maximum independent set problem, a maximum-weight independent set problem, a maximum clique problem, and a minimum vertex cover problem.

[0030] In one or more embodiments, weights in the maximum-weight independent set problem are encoded by applying light shifts to a plurality of qubits.

[0031] In one or more embodiments, the final state of the plurality of qubits comprises one or more of a solution to the encoded unweighted maximum independent set problem, a solution to the encoded maximum-weight independent set problem, a solution to the encoded maximum clique problem, and a solution to the encoded minimum vertex cover.

[0032] In one or more embodiments, the method further includes renumbering at least two vertices in the quantum computing problem prior to the encoding the quantum computing problem.

[0033] In one or more embodiments, the solution to the quantum computing problem comprises an approximate solution to the quantum computing problem.

[0034] In one or more embodiments, a method includes: arranging a plurality of qubits to encode a quantum computing problem; applying a sequence of q levels of light pulses to the plurality of qubits, wherein the q levels of light pulses comprises at least a first set of q variational parameters and a second set of q variational parameters; measuring the state of one or more of the plurality of qubits; optimizing, based on the measured state of at least some of the one or more of the plurality of qubits, the first set of q variational parameters and the second set of q variational parameters of the q levels of light pulses; optimizing, based at least on the first set of q optimized variational parameters and the second set of q optimized variational parameters of q levels of light pulses, a first set of p variational parameters and a second set of p variational parameters of p levels of light pulses, wherein q<p ; and measuring at least some of the plurality of qubits in a final state.

[0035] In one or more embodiments, optimizing the first set of p variational parameters and the second set of p variational parameters of p levels of light pulses further comprises computing a first set of p variational parameter starting values and a second set of p variational parameter starting values of the p levels of light pulses.

[0036] In one or more embodiments, computing of the first set of p variational parameter starting values of the p levels of light pulses, wherein p> 1, comprises: performing a Fourier transform on the first set of q variational parameters of the q levels of light pulses, into a plurality of k frequency components, each of the k frequency components having an amplitude u _k and computing the first set of p variational parameter starting values of the p levels of light pulses based on the amplitudes u _k ;

[0037] In one or more embodiments, computing of the second set of p variational parameter starting values of the p levels of light pulses, wherein p> 1, comprises: performing a Fourier transform on the second set of q variational parameters of the q levels of light pulses, into a plurality of k frequency components, each of the k frequency components having an amplitude v _k and computing the second set of p variational parameter starting values of the p levels of light pulses based on the amplitudes v _k.

[0038] In one or more embodiments, computing of the first set of p variational parameter starting values and computing of the second set of p variational parameter starting values of the p levels of light pulses, comprises: extrapolating the first set of p variational parameter starting values of the p levels of light pulses based on the first set of q variational parameters of the q levels of light pulses; and extrapolating the second set of p variational parameter starting values of the p levels of light pulses based on the second set of q variational parameters of the q levels of light pulses.

[0039] In one or more embodiments, the method further includes applying a sequence of p levels of light pulses to the plurality of qubits with a first set of p optimized variational parameters and a second set of p optimized variational parameters, wherein the measuring the at least some of the plurality of qubits in the final state comprises measuring the at least some of the plurality of qubits after the applying the sequence of p levels of light pulses to the plurality of qubits.

[0040] In one or more embodiments, the encoded quantum computing problem comprises a MaxCut problem, and wherein the final state of the plurality of qubits comprises a solution to the MaxCut problem.

[0041] In one or more embodiments, the encoded quantum computing problem comprises a maximum independent set problem, and wherein the final state of the plurality of qubits comprises a solution to the maximum independent set problem.

BRIEF DESCRIPTION OF THE FIGURES

[0042] Various objectives, features, and advantages of the disclosed subject matter can be more fully appreciated with reference to the following detailed description of the disclosed subject matter when considered in connection with the following drawings, in which like reference numerals identify like elements.

[0043] FIGS. 1 A-1E illustrate an exemplary scheme for encoding and finding solutions to optimization problems using an array of qubits, according to some embodiments.

[0044] FIGS. 2A-2B illustrate and characterize arrangements of vertex and ancillary qubits, according to some embodiments.

[0045] FIGS. 3 A-3B are graphs illustrating the performance of classical branch and bound algorithms for MIS on random UD graphs, according to some embodiments.

[0046] FIGS. 4A-4D are graphs illustrating the performance of quantum algorithms for MIS on random UD graphs, according to some embodiments.

[0047] FIG. 5 is an example random unit disk graph, according to some embodiments.

[0048] FIG. 6 is a graph showing a method to extract the adiabatic time scale TLZ, according to some embodiments.

[0049] FIGS. 7A-7C are graphs showing aspects of the adiabatic time scale TLZ ,

according to some embodiments.

[0050] FIG. 8 is a phase diagram for a quantum algorithm in terms of an approximation ratio r, according to some embodiments.

[0051] FIG. 9 is an illustration and graph showing behavior of a Rydberg problem, according to some embodiments. [0052] FIGS. 10A-10B are grid representations of planar graphs of maximum degree 3 and a transformation to a unit disk graph, according to some embodiments.

[0053] FIGS. 11 A-l 1C are grid representations of planar graphs of maximum degree 3 taking into account non-nearest-neighbor interactions, according to some embodiments.

[0054] FIGS. 12A-12E are grid representations of transformations to an effective spin model, according to some embodiments.

[0055] FIGS. 13A-13B are illustrations of local non-maximal independent set configurations and domain walls, according to some embodiments.

[0056] FIG. 14A is a schematic of a /i-level QAOA algorithm, according to some embodiments.

[0057] FIG. 14B is an illustration of an example MaxCut problem, according to some embodiments.

[0058] FIGS. 15A-15H are graphs showing QAOA optimal variational parameters, according to some embodiments.

[0059] FIGS. 16A-16B are graphs showing comparisons of embodiments of disclosed heuristic strategies to brute force methods, according to some embodiments.

[0060] FIGS. 16C-16D are graphs showing the average performance of QAOA, according to some embodiments.

[0061] FIGS. 17A-17J are plots showing the performance of QAOA as compared to QA, according to some embodiments.

[0062] FIGS. 18A-18F are graphs and illustrations showing nonadiabatic aspects of QAOA, according to some embodiments.

[0063] FIG. 19 is a graph showing the results of Monte-Carlo simulations of QAOA, according to some embodiments. [0064] FIGS. 20A-20F are graphs and illustrations showing exemplary vertex renumbering techniques and Rydberg implementations, according to some embodiments.

[0065] FIGS. 21A-21F are graphs showing QAOA optimal variational parameters, according to some embodiments.

[0066] FIG. 22 is a schematic of two variants of a FOURIER[c/, A] heuristic strategy, according to some embodiments.

[0067] FIG. 23 is a graph showing a comparison between various heuristic strategies, according to some embodiments.

[0068] FIGS. 24A-24C are graphs showing energy spectra and time to solution for QAOA and QA, according to some embodiments.

[0069] FIGS. 25A-25B are graphs showing exemplary instantaneous eigenstate populations and couplings, according to some embodiments.

DETAILED DESCRIPTION

[0070] Optimization algorithms are used for finding the best solution, given a specified criterion, for a specified problem. Combinatorial optimization involves identifying an optimal solution to a problem given a finite set of solutions. Quantum optimization is a technique for solving combinatorial optimization problems by utilizing controlled dynamics of quantum many -body systems, such as a 2D array of individual atoms, each of which can be referred to as a“qubit” or“spin.” Quantum algorithms can solve combinatorially hard optimization problems by encoding such problems in the classical ground state of a programmable quantum system, such as spin models. Quantum algorithms are then designed to utilize quantum evolution in order to drive the system into this ground state, such that a subsequent measurement reveals the solution. In other words, a problem can be encoded by placing qubits in a desired arrangement with desired interactions that encode constraints set forth by the optimization problem. When properly encoded, the ground state of the many- body system comprises the solution to the optimization problem. The problem can therefore be solved by driving the many -body system through an evolutionary process into its ground state.

[0071] Without being bound by theory, assuming complete control of the interactions between the qubits, it is possible to encode nondeterministic polynomial (“NP”)-complete optimization problems into the ground states of such systems. In most realizations, however, not all interactions are fully programmable. Instead, such interactions are determined by properties of specific physical realizations, such as, but not limited to locality, geometric connectivity, or controllability, which either constrain the class of problems that can be efficiently realized or imply that substantial overhead is required for their realization. Thus, one of the challenges in understanding and assessing quantum optimization algorithms involves designing methods to encode specific and larger classes of combinatorial problems in physical systems in an efficient and natural way.

[0072] In some implementations, quantum optimization can involve: (1) encoding a problem by controlling the positions of individual qubits in a quantum system given a particular type and strength of interaction between pairs of qubits and (2) steering the dynamics of the qubits in the quantum system through an evolutionary process such that their evolved final states provide solutions to optimization problems. The steering of the dynamics of the qubits into the ground state solution to the optimization problem can be achieved via multiple different processes, such as, but not limited to the adiabatic principle in quantum annealing algorithms (QAA), or more general variational approaches, such as, but not limited to quantum approximate optimization algorithms (QAOA). Such algorithms may tackle computationally difficult problems beyond the capabilities of classical computers. However, the heuristic nature of these algorithms poses a challenge to predicting their practical performance and calls for experimental tests. In addition, such systems, in their full generality, are inefficient and difficult to implement owing to practical constraints as described above, and can only be used on a subset of optimization problems.

[0073] Some aspects of the present disclosure relate to systems and methods for arranging qubits in programmable arrays that can encode or approximately encode in an efficient way a broader set of optimization problems. In some embodiments, chains of even numbers of adjacent“ancillary” qubits are used to encode interactions between distant qubits by connecting such distant qubits with chains of ancillary qubits, for example as described in more detail with reference to FIG. 2A. As described in more detail throughout the present disclosure, these chains of“ancillary” qubits can be used to encode interaction between certain“vertex” qubits, but not other vertex qubits and to reduce the strength of long-range interaction between two vertex qubits that are not intended to be connected via an edge. In some embodiments, the effects of long-range interactions can be further reduced by introducing a detuning parameter to a chosen control technique to selectively control interaction between particular qubits. For example, for comer and junction qubits, detuning patterns described in the present disclosure can reduce the effects of long-range interactions such that the ground state of the system is the optimal solution to the encoded problem. The techniques described herein can permit efficient encoding of a larger set of optimization problems beyond simple unit disk graphs.

[0074] Some additional or alternative aspects of the present disclosure relate to systems and methods for coherently manipulating the internal states of qubits, including excitation. In some embodiments, techniques are disclosed that can be used to evolve an encoded problem to find an optimal (or an approximately optimal) solution. For example, embodiments of the present disclosure relate to optimal variational parameters and strategies for performing the Quantum Approximate Optimization Algorithm (“QAOA”), some embodiments of which are described, for example, with reference to FIG. 14 A. For example, embodiments include heuristics for classical feedback loops that can improve the performance of brute-force QAOA implementations. In some embodiments, these strategies perform at least as well if not better than existing algorithms. Some aspects of the present disclosure focus on implementations of using QAOA to solve MaxCut combinatorial problems, but the disclosed techniques are not limited thereto.

Exemplary optimization problems and encodings

[0075] In some embodiments, particular types of optimization problems can be encoded with an arrangement of qubits. For example, FIGS. 1 A-1E show an exemplary scheme for encoding and finding solutions to optimization problems using an array of qubits, according to some embodiments. FIG. 1 A shows aspects of a Rydberg blockade mechanism and maximum independent set on unit disk graphs, according to some embodiments. One exemplary optimization problem that can be solved using the techniques described in the present disclosure is a maximum independent set (“Mis”) problem. Given a graph G with vertices V and edges E, an independent set can be defined as a subset of vertices where no pair is connected by an edge. FIG. 1 A shows an example graph with vertices such as 102, 104. Vertices 102, 104 can be connected via an edge, such as edge 106. The computational task of a Mis problem is to find the largest such set, called the maximum independent set (Mis). AS shown in the graph of FIG. 1 A, the maximum independent set is denoted via black vertices such as 102, none of which are connected. In the example of FIG. 1A, the size of the maximum independent set is 6. Determining whether the size of Mis is larger than a given integer a for an arbitrary graph G is a well-known NP-complete problem. Furthermore, even approximating the size of the Mis is an NP-hard problem. In some embodiments, the MIS problem is also equivalent to the maximum clique problem and the minimum vertex cover problem. Thus, a solution to the MIS problem will constitute a solution to the corresponding maximum clique problem and the minimum vertex cover problem. [0076] Without being bound by theory, the embodiment of FIG. 1 A can be referred to as a unit disk (“UD”) graph. UD graphs are geometric graphs in which vertices are placed in a 2D plane and connected if their pairwise distance is less than a unit length, r. In other words, UD graphs are graphs where any two vertices within a distance r from one another are connected with an edge, such as vertices 102, 104 which are connected via edge 106. Vertex 108 is too far from vertices 102, 104 to be connected therewith with an edge. The Mis problem on UD graphs (UD-Mis) is still NP-complete and can be used to find practical situations ranging from, for example, wireless network design to map labelling in various industry sectors.

[0077] In some embodiments, a Mis problem can be formulated as an energy

minimization problem, by associating a spin- 1/2 with each vertex v E V. Vertices like those shown in FIG. 1 A can be prepared such that after a driving sequence, such as one with a Rabi frequency W and detuning D that vary over time (as shown in FIGS. 1D and 1E), the state |l) of each qubit is energetically favored unless a nearby vertex is also in the state 11), in which case it is energetically favored to have one vertex transition to the state |0). Without being limited by theory, the Hamiltonian (“///”) of such a system can be represented as follows: where a spin-l/2 system is assigned with states |0) and |l) to each vertex, = |l) _u (l|, D is the detuning on the spin, and U is the energy penalty when two spins ( v , w ) connected by an edge (E) are both in state |l). Initially, all vertices can be prepared in the |0) state. Driving causes at least some of the vertices to transition to the 11) state. For D > 0 , Hp favors qubits to be in state |1). However, if //> D, it is energetically unfavorable for two qubits, u and v, to be simultaneously in state |1) if they are connected by an edge u, v E E. Thus, each ground state of Hp represents a configuration where the qubits that correspond to vertices in the maximum independent set are in state |1), and all other qubits are in state |0).

[0078] FIG. 1C is a graph showing interatomic interaction potentials between two adjacent vertices in the limit of weak driving, where W « D and D > 0, according to some embodiments. Changing the detuning and Rabi frequency over time can produce quantum evolution that changes the system from the initial state to the final state, which can include a solution (or one or more approximate solutions) to the encoded problem. As described above, under such conditions, for qubits closer than re (the Rydberg blockade radius) it is energetically favorable for one of the qubits to stay in the |0) state. For example, when II D > 0, the Hamiltonian Hp energetically favors each spin to be in the state 11) unless a pair of spins are connected by an edge (i.e., within the Rydberg radius). Thus, in the ground state of the Hamiltonian Hp , only the vertices in the Mis are in state 11). Such a state can be referred to as a Mis-state, and Hp can be referred to as a Mis-Hamiltonian. In some embodiments, the NP-complete decision problem of Mis becomes deciding whether the ground state energy of Hp is lower than -a A.

[0079] In some embodiments, a quantum annealing algorithm (“QAA”) can be used to evolve the quantum state from the initial state to the final state, which encodes the solution of the optimization problem. For example, a simple QAA can be implemented by adding a transverse field H _D = å _v Ws* with cf = |0) (l| + |l) (0|, that induces quantum tunneling between different spin configurations. FIG. 1B shows a transition between spin states, according to some embodiments. For example, the arrow labelled as Rabi frequency W describes a transition between the spin states |0) and |l) controlled, for example, by a laser beam with the Hamiltonian Ws * for a spin v. When the detuning D < 0, the ground state of the Hamiltonian Hp is the trivial state where all qubits are in the |0) state. When the detuning

D > 0, the ground states of Hp are the MIS states where qubits in state |l) form the MIS. By slowly changing W and D in a time-dependent fashion, for example as shown in FIGS. 1E and 1D, respectively, the trivial ground state with all qubits in |0) (D < 0, W = 0) can be adiabatically connected to a final state encoding the Mis (D > 0, W = 0), resulting in the qubits being left in a solution to the optimization problem. Note that in general, such a procedure can involve transitions across vanishing energy gaps as discussed in more detail throughout the present disclosure. Exemplary non-limiting embodiments of simulations are described in more detail below in the section titled“Quantum annealing for random ETD- Mis.”

[0080] In some embodiments, MIS problems can be implemented using Rydberg interactions between individual atoms. For example, as discussed in more detail in

PCT/ETS18/42080, graphs like that shown in FIG. 1A can be implemented using two-level qubits, the states of which are shown in FIG. 1B, according to some embodiments. ETsing optical tweezers, atoms (qubits) can be individually and deterministically arranged in fully programmable arrays in one, two and even three dimensions. Such a system can use individually trapped homogeneously excited neutral atoms interacting via the so-called Rydberg blockade mechanism. Each atom realizes a qubit, u, with an internal ground state, |0), and a highly excited, long-lived Rydberg state, 11), which can be coherently manipulated by external laser fields. If two atoms are both in this Rydberg state and within the Rydberg blockade radius, they interact via strong van der Waals interactions, which is energetically unfavorable. This makes it possible to encode a ETD-Mis problem using qubits as vertices and the disk size set by the Rydberg radius, as explained in more detail below. FIG. 9 shows an exemplary solution to a ETD-Mis problem, where some qubits 904 are found in the internal ground state, |0), and, other qubits 902 in the long-lived Rydberg state, 11), according to some embodiments. [0081] Without being limited by theory, the Hamiltonian governing the evolution of embodiments of such a system can be represented as follows:

where W _u and D _u are the Rabi frequency and laser detuning at the position x _v of qubit v.

While individual manipulation is feasible, such a system can also be implemented with a homogeneous driving laser, for example, where W _u = W and D _u = D. The operator s = |0) _v (l| + |l) _v (0| can give rise to coherent spin flips of qubit v and n _v = | l) _v (l| counts if the qubit v is in the Rydberg state. In some embodiments, for isotropic Rydberg states, the interatomic interaction strength depends only on the relative atomic distance, x, and is given by Ryd(x) = C/x^, where C is a constant. The strong interactions at short distances energetically prevent two qubits from being simultaneously in state 11) if they are within the

Rydberg blockade radius TB = (c / i ) ¹ ' ⁶ , as shown in FIG. 1C, resulting in the so-called blockade mechanism, according to some embodiments. This can provide a connection between the Rydberg Hamiltonian HRyd and solutions to UD-MlS problems. In other words, the Rydberg blockade causes it to be energetically favorable for adjacent qubits within the Rydberg radius (r/s) to not both be in the Rydberg state |l). Like in UD graphs, in some embodiments of this configuration the qubits cannot be found to be both in state |l) when separated by a distance smaller than the Rydberg blockade radius.

[0082] The MlS-Hamiltonian Hp shares some features with the Rydberg Hamiltonian i¾ _yd in the classical limit, W _n = 0. In some embodiments, the main difference lies in the achievable connectivity of the pairwise interaction, for example, when arbitrary graphs are allowed in Hp. A special, restricted class of graphs can be considered that are most closely related to the Rydberg blockade mechanism. These so-called unit disk (UD) graphs, as discussed above, are constructed when vertices can be assigned coordinates in a plane, and only pairs of vertices that are within a unit distance, r, are connected by an edge. Thus, the unit distance r plays an analogous role to the Rydberg blockade radius TB in i¾ _yd. In other words, spatial proximity of the qubits is used to encode the edges in the UD-MlS problem. Mis is NP-complete even when restricted to such unit disk graphs. While embodiments of the present disclosure discuss 2D problems and 2D arrangement of qubits, a person of skill in the art would understand, based on the present disclosure, that aspects of the problem encoding described herein would be applicable to other spatial structures such as a one- dimensional or three-dimensional structure.

[0083] Without being bound by theory, the maximum-weight independent set problem is a MIS problem where each vertex v has a weight D _n that replaces the homogenous weight D in equation 1. The maximum-weight independent set problem is to find an independent set with the largest total sum of weights. In some embodiment implementation, these weights D _n can be encoded by applying corresponding light shifts to each qubit. In some embodiments, these light shifts can be AC Stark shifts created by off-resonant laser beams or spatial light modulators.

Exemplary qubit arrangements and detunings

[0084] Although Rydberg interactions decay significantly beyond the Rydberg radius, there are still long-range interaction tails between distant qubits, such as 102 and 108, shown in FIG. 1 A. Interaction tails, taken alone or in aggregate, can reduce the likelihood that the system will be found in the solution after driving, such as that shown in FIGS. 1D and 1E. Therefore, implementations such as that shown in FIG. 1 A may not always perfectly encode an NP-complete Mis problem. Furthermore, the range of NP-complete problems that can be encoded using the technique described with reference to FIG. 1 A is limited to those related to unit disks. [0085] In some embodiments, one or more of these problems can be solved by choosing atom positions in two dimensions and laser parameters such that the low energy sector of the Rydberg Hamiltonian HR _yd reduces to the (NP-complete) MlS-Hamiltonian Hp on planar graphs with maximum degree 3. In some aspects of the present disclosure, antiferromagnetic order can be formed in the ground state of (quasi) 1D spin chains of ancillary qubits at positive detuning, due to the Rydberg blockade mechanism. Such a configuration can effectively transport the blockade constraint between distant vertex qubits. In other aspects of the present disclosure, a detuning pattern, {D _n } can be introduced to eliminate the effect of undesired long-range interactions without altering the ground state spin configurations. Some embodiments allow efficient encoding of NP-complete problems in the ground state of arrays of trapped neutral atoms. Without being bound by theory, the ground state energy of

Rydberg interacting atoms in 2D array in such embodiments is NP-hard (and NP-complete when W _n = 0).

Exemplary qubit arrangements with ancillary qubits

[0086] As described with reference to FIGS. 2A and 2B, NP-complete problems can be encoded with“ancillary” qubits or vertices that implement edges between“vertex” qubits with vertex qubits having a maximum degree of 3, according to some embodiments. For example, as shown in FIG. 2 A, a plurality of vertex qubits 202, 204, 212, 214, and 216 can be arranged in a graph. Edges between pairs of vertex qubits, such as vertex qubits 202 and 204, can be implemented using an even number of ancillary qubits 206, each of which is separated by a unit length r. This technique can be used to embed a planar grid graph with maximum degree 3 into a UD graph. For example, a planar graph with maximum degree 3 can be efficiently embedded on a grid (with grid unit g ), Q, and transformed to a UD graph, G, by introducing an even number of ancillary vertices along each edge, as shown in FIG. 2A, according to some embodiments. The UD radius r— can be determined by a parameter k , proportional to the linear density of ancillary vertices along each edge. According to some embodiments, it is desirable to implement approximately the same density of ancillary qubits along any given edge to ensure that the interactions that form each vertex are roughly equal in strength.

[0087] In the example of vertex qubits 202 and 204, vertex qubit 202 can interact with the leftmost ancillary qubit 206 as if it were vertex qubit 204, and vertex qubit 204 can interact with the rightmost ancillary qubit 206 as if it were vertex qubit 202. In this way, edges can be implemented between vertex qubits outside the Rydberg radius, and in ways that cannot be implemented purely as a unit disk graph. Furthermore, whereas in unit disk graphs like that discussed with reference to FIG. 1 A, it would not be possible to have some vertices separated by the same distance connected by an edge (e.g., vertex qubits 202, 204) while other vertices separated by the same distance are not connected by an edge (e.g., vertex qubits 212, 214), such configurations can be realized using the ancillary vertices described in the present disclosure.

Exemplary detuning patterns

[0088] In some embodiments, when a Mis graph is implemented as shown in FIG. 2A, the tail of the Rydberg interactions does not affect interactions between qubits except in the vicinity of comers and junctions, where arrays of ancillary qubits meet at an angle. Vertex qubits 202 is an exemplary corner, while vertex qubit 216 is an exemplary junction. As described in more detail with reference to FIG. 2B, the detuning pattern during driving of the qubits can be adjusted around these structures to compensate for the effect of interaction tails.

[0089] FIG. 2B shows an example of such a detuning pattern, according to some embodiments. The inset of FIG. 2B shows a vertex qubit 216 with degree of 3, such as the vertex qubit 216 in FIG. 2A, which is adjacent ancillary qubits 222, 224 in a first direction, 232, 234 in a second direction, and 242, 244 in a third direction. FIG. 2B plots detunings A _jpQi) of qubits at distances j from the vertex qubit, with detuning Ac, along the vertical (horizontal) direction. As shown in FIG. 2B, the detunings A _j can individually address each ancillary qubit (for example, using individual addressing techniques described in PCT/US18/42080), and set in a way that eliminates the effects of interaction tails. For example, detunings A _v can be selected in a range 0.5lC/r ⁶ < A _v < Amax º C/r ⁶ , and with

sufficiently large k. In some embodiments, detunings can be selected arbitrarily within this range. Such detunings can compensate for undesired effects of long-range interactions.

Under such conditions, in the ground state of i¾ _yd (at W _n = 0), each array of ancillary qubits is in an ordered configuration, alternating between |0) and |l) (with at most one domain wall). This ensures that the ground state of i¾ _yd (2) for W _n = 0, coincides with the ground state of the Mis Hamiltonian (1) on G. This perturbation of the detunings for individual ancillary qubits can be performed with multiple different types of global driving patterns (such as QAA and QAOA described in more detail throughout the present disclosure) and without prior knowledge of what the solution is to the Mis. Furthermore, since the size of the Mis on G is identical to that of the Mis on G , up to half of the number of ancillary vertices, it is still NP-hard to find the ground state of the Rydberg Hamiltonian (2) when a Mis problem is implemented using ancillary vertices.

[0090] In some embodiments, implementations can be expanded beyond unit disk graphs to more general graphs. Such implementations can be referred to as“stroboscopic” implementations, which can involve an arbitrary graph implemented without ancillary qubits described with reference to FIG. 2 A. For example, in some embodiments, various optical techniques can be employed to enlarge the class of problems that can be addressed in these systems. One exemplary approach involves qubits encoded in hyperfme ground states, rather than ground-Rydberg qubit encoding, and selective excitation (for example, with individual addressing of qubits) into various kinds of Rydberg states, such as Rydberg S and P states. In this example, strong and long-range dipolar interactions between Rydberg states with different parities S and P can be used to efficiently realize rotations of a single qubit, controlled by its multiple neighbors. For example, Rydberg S and P states can be used to realize multi-qubit-controlled rotations. Rotations (or qubit flips) can be used to evolve the quantum state from the initial state to the final MIS-state. By controlling these interactions, it is possible to ensure that evolution of a system does not violate the independent set constraint, for example, by ensuring that two qubits do not evolve to both be in the Rydberg state if they are connected by an edge. In some examples, all neighbors of a given (central) qubit can be selectively excited from the state 11) to an Estate. Subsequently, a two-step transition between the two hyperfme states of the central qubit, with an intermediary step in the Rydberg -state, can implement a multi-qubit-controlled rotation. Combined with Trotterized time evolution (e.g., splitting the total Hamiltonian evolution into small discrete steps), these techniques provide means to go beyond UD graphs, and implement various quantum optimization algorithms for Mis problems on arbitrary graphs. In some

embodiments, this technique can implement arbitrary graphs, which can include a broader range of graphs than what is depicted in FIG. 2A. In some embodiments, the“stroboscopic implementation” (selective and sequential excitation of qubits into Rydberg states) allows for the implementation of multi-qubit-controlled rotations. This selective excitation can facilitate controlled rotations according to the graph structure, which is not necessarily limited to the unit disk geometry.

Comparison of exemplary classical and quantum algorithms

[0091] According to some embodiments, it is desirable to identify the types of UD-MlS problems for which finding solutions with quantum computing would show large

improvements over standard computing approaches. Numerical simulations of both classical and quantum algorithms can be used to identify regimes and system sizes where quantum algorithms are well suited for UD-MlS. In one embodiment, Mis on random UD graphs constructed by placing N vertices randomly with a density p in a 2D box of size L x L, where

L = JN/P and UD radius r = 1 can be considered. The hardness of this problem can depend on the vertex density: if the density is below the percolation threshold p _c ~ 1.436, the graph decomposes into disconnected, finite clusters, allowing for efficient polynomial time algorithms. In the opposite limit, for example if the density is very large, p ® ¥ (and N oc p), the problem becomes essentially the closest packing of disks in the (continuous) 2D plane, with the known best packing ratio of n/(2 [3).

[0092] FIGS. 3 A-3B show the performance of an exemplary classical optimization algorithm, according to some embodiments. In particular, it shows that in some density regions it is difficult for classical algorithms to solve optimization problems, requiring exponential time to solve exactly. In these regimes, quantum advantage is more beneficial than other regimes. FIGS. 4A-4D show the corresponding performance of quantum algorithms (both QAA and QAOA), according to some embodiments. Though implemented with small-size simulations, FIGS. 4A-4D indicate that QAOA can solve problems fast and beat QAA. FIGS. 4A-4D show quantum advantage for certain regimes, according to some embodiments. More particularly, FIG. 3 A is a graph showing the performance of classical optimization algorithms, focusing on a branch and bound algorithm, according to some embodiments. The median runtime 7iun to find the Mis in CPU time is shown for bins at given numbers (vertical axis) and densities (horizontal axis). The exemplary statistics shown in FIG. 3 A were obtained from 50 graphs per data point. The dashed line indicates the percolation threshold described above, which in this case is shown to hep = p _{c «} 1.436. As can be seen in FIG. 3 A, in certain regimes of the parameter space of N and >, problems are difficult for classical optimization algorithms to solve because, for example, the run time is large. [0093] FIG. 3B is a graph showing the run time Tnm as a function of N, according to some embodiments. While there is a polynomial scaling for p < p _{c «} 1.436, at an intermediate density p ~ l, the runtime exhibits an exponential dependence with system sizes exceeding N - 150 vertices. This exponential dependence can be seen in both worst cases and typical instances for any given number N. Note that at p ~ 7, the classical algorithm often could not find a solution for N ³ 440 within 24 hours of CPU time on a single node of the Odyssey computing cluster (managed by Harvard University). For more general, non-UD graphs, this limitation of the classical algorithm is observed aXN ³ 320. For comparison, the

performance on Erdos-Renyi graphs with edge probability p = 0.05 are also shown, demonstrating a similar exponential runtime. Error-bars are 5th and 95th percentile values among 100 graphs. FIG. 3B shows that it can take a long time (for example, 24 hours) for classical algorithm to solve problems of a certain size (for example, -N = 400). This informs the problem size to test quantum algorithms such that the problem is nontrivial to solve classically.

[0094] FIGS. 4A and 4B are graphs showing performance of quantum algorithms for Mis on random UD graphs, according to some embodiments. FIG. 4A is a hardness diagram for a quantum adiabatic algorithm in terms of adiabatic time scale TLz, according to some embodiments, plotted against N on the vertical axis and density p on the horizontal axis. For example, FIG. 4A shows the results of numerical simulations of the adiabatic QAA, with the Mis annealing Hamiltonian HD + Hp described above, while sweeping the parameters as:

Equation 3: ri(t) = A ₀ (2t/T— 1), Equation 4: W(ί ) = W ₀ sin ² (nt/T )

[0095] The limit U /W ₀ ® ¥ can be observed, where the dynamics are restricted to the independent set subspace, allowing numerical simulation of system sizes up to A— 40 qubits. Sloped dashed lines parallel to the stripe pattern correspond to optimal disk packing. The fully connected region has trivially |MIS| = 1. The Landau-Zener time scale, TLz, required for adiabaticity can be extracted by fitting numerical results to the expected long-time behavior of the ground state probability MIS = 1— e ^{a_T//T z} . For example, it is possible to fit MIS for large T to extract 7LZ . FIG. 4 A shows a clear transition at the percolation threshold p = 1.436 from the regime in which the problem is easy to solve to that where it is not. At high density the size of the Mis can be determined by L = N / p owing to the relation to the optimum packing problem. This is also reflected in TLz as visible by the stripe pattern in FIG. 4A. The exponential Hilbert space size limits numerical simulations, making a conclusive comparison with the scaling of the classical algorithm above difficult. In some implementations, finite coherence times can further limit the performance of the adiabatic algorithms with exponential scaling of the minimum energy gap, e _gap corresponding to the energy gap between the ground state and first excited state along a quantum adiabatic path.

[0096] Several approaches can be implemented to try to overcome these potential limitations. Such approaches include heuristics to open up the gap, the use of diabatic (non- adiabatic) transitions in QAA, and variational quantum algorithms such as QAOA studied below. FIG. 4B is a graph showing the probability to find the Mis, MIS for simulated non- adiabatic QAA, according to some embodiments. In some embodiments, 7LZ and /Vis can convey similar information, and is related by .PMIS = 1— e ^{a_T//T z} for large T (in the adiabatic regime). FIG. 4A shows an exemplary adiabatic regime (large T) while FIG. 4B shows an exemplary nonadiabatic regime (short T). For example, FIG. 4B shows a single sweep (rather than sweeping multiple times for total time T, as shown in FIG. 4A) for a single T, which is in an exemplary short-time limit ( e.g ., nonadiabatic) with profile in Eq. (3) with a short annealing time T= 10/Wo « TLz. FIG. 4B shows that substantial overlap with the Mis state can be obtained with evolution time T « TLz, while the success probability resembles qualitatively the adiabatic hardness diagram. These results are consistent with observations that non-adiabatic annealing can be advantageous, especially when the minimum spectral gap € _gaP is small. In some embodiments, if the gap is small (compared to other energy scales), it can require a long coherence time to run the quantum adiabatic algorithm (for example, T »

1/ <ºgap)·

[0097] As described throughout the present disclosure, it is possible to take advantage of a direct connection between the many-body problem of spins interacting via van der Waals interactions and computational complexity theory. Individual control over the positions of such spins allows for NP-complete optimization problems to be directly encoded into such quantum systems. This result can be obtained from a reduction from MIS on planar graphs with maximum degree of 3. Quantum optimizers based on the techniques described in the present disclosure, in combination with techniques to trap and manipulate neutral atoms, can address NP-hard optimization problems as an improvement over traditional computing techniques.

Exemplary quantum algorithms

[0098] As discussed above, after encoding a combinatorial problem using the position of a quantum system, whether using the“ancillary” qubit technique described above or not, the next step is to evolve the system in a way that produces a ground state that is a solution to the combinatorial problem. Some examples include QAA and QAOA.

QAOA

[0099] In some embodiments, QAOA can be applied to quantum optimization problems like those described in the present disclosure. For example, a QAOA (of level p) can consist of applying a sequence of p resonant pulses to all qubits (or with some detuning and energy adjustments for specific qubits as described in more detail in the present disclosure) of varying duration, fe, and optical phase, pk on an initial prepared state. This can generate a variational wavefunction i p({t _k , _k })) = n _k=1 U _k \0)® ^N , where U _k = e ^~lHktk and H _k = Without being bound by theory, the present disclosure describes some additional or alternative optimization parameters (g, b), which can in some embodiments be related to the duration, fe, and optical phase, p as described in more detail in the section below titled QAOA for MIS.

[0100] The resulting quantum state can be measured in the computational basis and then used for optimization. For example, the variational parameters t and can be optimized in a classical feedback loop from measurements of ( H _p ) in this quantum state. Examples of this optimization are described in more detail below.

[0101] The performance of QAOA depends in part on the chosen classical optimization routine. Before explaining exemplary implementations of such techniques, the performance of QAOA can be shown, which marks an improvement upon classical computation techniques in some embodiments. For example, FIG. 4C compares optimized QAOA with quantum annealing, on an exemplary graph instance, Go, (with N = 32 and independent set space dimension dimis = 17734, p = 2.4 and |MIS| = 8, e _{gap »} 0.0012Wo) for which the adiabatic timescale can be observed to be approximately ~10 ⁶ /Wo, according to some embodiments. As shown in FIG. 4C, QAOA achieves 0(1) success probability at depth p < 40 (i.e., less than 40 resonant pulses) thereby displaying significantly shorter evolution times, T = å _k t _k ~10/n ₀ « T _Lz and showing exceptionally good variational states, which can indicate that the algorithm is effective at solving the encoded problem. And discussed in more detail below, this is also reflected in implementations of QAOA on large systems. For exemplary system sizes that can be assessed using common numerical analytics, however, the fluctuations associated with measurement projection noise can imply that the MIS size can be easily found during the classical optimization routine, even at smaller p. In some

embodiments, this“measurement projection noise” is not noise in the measurement. For example, each measurement can produce a configuration space. To estimate the expectation value of an observable and perform QAOA optimization, the measurement has to be repeated multiple times. For small system sizes, the measurement result may hit the MIS before much optimization since each measurement has a high probability to produce the MIS.

[0102] FIG. 4D shows the simulated performance of various quantum algorithms on the graph Go, according to some embodiments. For example, FIG. 4D shows the average difference between the true MIS and the largest independent set found after M measurements for simulated QAOA and QAA on Go. As shown in FIG. 4D, both QAOA with heuristic ansatz as well as random guess methods are shown. To obtain such results, simulations for QAOA were performed with Monte Carlo simulations techniques, including sampling from projective measurements. The size of the largest measured independent set is plotted as a function of the number of measurements (averaged over Monte Carlo trajectories). With a properly chosen heuristic ansatz for tk and pk, QAOA (at p = 3) finds the MIS already after ~ 10 ² -10 ³ measurements. This is comparable to the number of measurements required for accurately estimating (HP), as part of the classical optimization routine. This performance can be compared to repeated runs of QAA with sweeps described by equation 3 above (and Do/Ώo = 6) for non-adiabatic annealing times T= 1/Wo, 10/Wo, 100/Wo. For this small system size, QAOA performs similarly as non-adiabatic QAA at the annealing time T ~ 10/Wo.

[0103] As explained above, QAOA is a quantum processor that prepares a quantum state according to a set of variational parameters. Using measurement outputs (such as the measured qubit state (|0) or 11)) of each qubit), the parameters can then be optimized, for example via a classical computer and fed back to the quantum machine in a closed loop. In QAOA, the state is prepared by a /i-level circuit specified by 2 p variational parameters. In other words, each pulse has two parameters ( t _k , <p _k ), for k from 1 to p. Even at the lowest circuit depth (p = 1), QAOA has non-trivial provable performance benefits (for example, better performance than a simple classical algorithm) but cannot be efficiently simulated by classical computers.

[0104] However, very little is known about QAOA with p > 1. One hurdle lies in the difficulty to efficiently optimize in the non-convex, high-dimensional parameter landscape. Without constructive approaches to perform the parameter optimization, any potential advantages of the hybrid algorithms could be lost. Furthermore, although QAOA can be shown to succeed in the/? ¥ limit due to its ability to approximate adiabatic quantum annealing (i.e., the adiabatic algorithm), its performance when 1 <p « oo is largely unexplored.

Exemplary parameter optimization for QAOA

[0105] Aspects of the present disclosure detail techniques to efficiently optimize QAOA variational parameters. In some examples, given a set of qubits in a particular arrangement, QAOA proceeds by applying a series of operations to the qubits, each operation having at least two variational parameters. The evolved state of the qubits is measured, which is fed back into an optimization routine (such as a classical algorithm) to adjust the variational parameters. The process then repeats on the qubits until it is determined that the measured state of the qubits is a solution to the encoded problem or an approximation thereof. In some embodiments, techniques for classical optimization routines are disclosed. These techniques are quasi-optimal in the sense that they typically produce known global optima, and generically require 2 ^oip) brute-force optimization runs to surpass performance. Aspects of the present disclosure also disclose implementations of QAOA with such optimization techniques, such as an example involving a 2D physical array of a few hundred Rydberg- interacting atoms that presents potential advantages over classical algorithms for solving MaxCut problems. [0106] As discussed above, many real-world problems can be framed as combinatorial optimization problems. In some embodiments, these are problems that can be generally defined on iV-bit binary strings z = z _x ··· z _N , where the goal is to determine a string that maximizes a given classical objective function C(z ) : {+1, - 1 } ^v ·® _³0 . An approximate optimization algorithm aims to find a string z that achieves a desired approximation ratio

where Cmax = max _z C(z).

[0107] The Quantum Approximate Optimization Algorithm (QAOA) is a quantum algorithm that can tackle these combinatorial optimization problems. To encode the problem, the classical objective function can be converted into a quantum problem Hamiltonian by promoting each binary variable zt into a quantum spin af

[0108] FIG. 14A is a flow diagram of an example QAOA algorithm, according to some embodiments. As shown in FIG. 14A, the quantum circuit takes inputs |+) ^®n 1446 and alternately applies sets of p levels (1440, . . .,1450) of e ^~lYiHc (1442, . . ., 1452) and C _b. = 1454). In some embodiments, e ^~lYiHc and ¾ represent 2 different types of evolution of the system, which can be induced physically by laser pulses. The parameters g and b can quantify for how long each evolution occurs. Each of these levels applies variational parameters 1480, which can differ between each of the set of p levels. The outputs are then measured at 1460 with respect to a function ( H _c ). In some embodiments, the outputs are a superposition of spin states. When measured the system collapses to one of these spin states, so the measurement result can be a single spin configuration. In some embodiments, the expectation value ( H _c ), which is a function of the spin configuration, is the number arrived at when evaluating the Hamiltonian He in this configuration, and averaged for multiple experimental runs, which is then fed through an optimizer 1470 (which can be a classical optimization function), which calculates the parameters (g, /?) for each level 1440, . . ., 1450 in the set that maximizes ( H _c ). These parameters ( gb ) are then provided as variational parameters 1480 for the next iteration (for example, as a feedback loop). Note that the variational parameters 1480 may be different at each level 1440, . . .,1450 in the set of p levels. The inputs |+) ^®n 1446 are then provided again and the process repeats. In some embodiments, after a number of iterations, the measurement 1460 can be considered to be the solution to the encoded problem.

[0109] For the />-level QAOA described in FIG. 14 A, the quantum processor can be initialized in the state |+) ^®w (1446). The problem Hamiltonian Hc (1442, 1452) and a mixing Hamiltonian H _B = 1454) are applied alternately (p times) with controlled durations to generate a variational wavefunction:

which is parameterized by 2 p parameters g; and bί (i = 1, 2, p) (one for each level in the set of p level). The expectation value He in this variational state can be determined as follows: which can be identified by repeated measurements 1460 of the quantum system in the computational basis and taking the average (e.g., with two orthogonal spin configurations |0) and 11), which can be measured using existing measurement techniques).

[0110] Once the measurements 1460 are performed, the results (e.g., the calculated F _p determined by taking the average over many Hc) can be fed back to an optimizer 1470, such as a classical computer, to search for the optimal parameters (g ^* b ^* ) so as to maximize

Equation !·. (y ^* ,/? ^* )= arg max _p (yj ?).

Y. b

[0111] The performance of QAOA can, in some embodiments, be benchmarked based on the approximation ratio:

[0112] In some embodiments, r characterizes how good the solution provided by QAOA is. The higher the r value, the better the solution.

[0113] In some embodiments, without being bound by theory, this QAOA framework can be applied to general combinatorial optimization problems. In one example, an archetypical problem called MaxCut can be considered.

[0114] FIG. 14B shows an example of a MaxCut problem, according to some

embodiments. A MaxCut problem can be described for an input graph G = (V, A). Here, V = { 1, 2, . . . , A) denotes the set of vertices (shown as up and down spins 1492A, 1492B,

1492C, 1494A, 1494B) and E = {((iQ), w _i;· )} is the set of edges (shown with solid and dotted lines), where w _Lj e M _³0 is the weight of the edge (i,j) connecting vertices i and j. For example, the edge connecting the two up spins 1492A, 1492B has a weight of wis. The goal of MaxCut is to maximize the following objective function:

where an edge (i,j) contributes with weight Wij if and only if spins of and o are anti aligned. For example, the curved dashed line shows a cut through all edges of spins that are anti-aligned, excluding edges wn and W25, which are aligned. Exemplary combinatorial problems

[0115] While embodiments of the present disclosure discuss MaxCut on 6/-regular graphs, where every vertex is connected to exactly d other vertices, based on the present disclosure a person of skill in the art would understand that the aspects of QAOA described herein would be applicable to other types of combinatorial problems such as, but not limited to MaxCut problems on other types of graphs, Maximum Independent Set problems, and others. Two types of d- regular MaxCut graphs are considered: (1) unweighted 6/-regular graphs (ur/R), where all edges have equal weight w _Lj = 1; and (2) weighted _< i-regular graphs (wr/R), where the weights w _Lj are chosen uniformly at random from [0,1] (though other weights other than the interval [0,1] can be selected in other embodiments).

[0116] It is NP-hard to design an algorithm that guarantees a minimum approximation ratio of r * > 16/17 for MaxCut on all graphs, or r* > 331/332 when restricted to unweighted 3 regular graphs (“u3R”) discussed above.

[0117] According to some embodiments, QAOA presents several benefits. For certain cases, it achieves a guaranteed minimum approximation ratio when p = 1. Additionally, under some reasonable complexity -theoretic assumptions, QAOA may not be efficiently simulated by any classical computer even when p = 1, making a candidate algorithm for “quantum supremacy,” or the ability of a quantum computer to perform calculations that a traditional computer cannot. The square-pulse (“bang-bang”) ansatz of dynamical evolution, of which QAOA can be one example, can be optimal given a fixed quantum computation time. In general, the performance of QAOA can improve with increasing p , achieving r 1 as p ¥ since it can approximate adiabatic quantum annealing via Trotterization. This monotonicity makes it more attractive than quantum annealing, whose performance may decrease with increased run time. [0118] While some embodiments of QAOA have a simple description, not much is currently understood beyond p = 1. For the example problem of MaxCut on u2R graphs, (such as 1D antiferromagnetic rings,) QAOA may yield r > (2p+ 1 )/(2p+2) as determined by numerical evidence. In another example, such as Grover’s unstructured search problem among n items, QAOA can find the target state with p achieving the optimal query complexity within a constant factor. In some embodiments, for more general problems, a simple brute-force approach can be used by discretizing each parameter into 0(poly(/V)) grid points. However, this technique requires examining N ^oip) possibilities at level p, which can become impractical to calculate using typical computers as p grows. Embodiments of the present disclosure therefore address efficient optimization of QAOA parameters and understanding of the algorithm for 1 « p < ¥.

Exemplary heuristics for parameter optimization

[0119] Some embodiments of the present disclosure relate to techniques for optimizing variational parameters. As described in more detail below, patterns in optimal parameters can be exploited to develop a heuristic optimization strategy for more quickly identifying the optimal variational parameters. In some examples described below, parameters identified for level -p QAOA can be used to more quickly optimize parameters for level -(p + 1) QAOA, thereby producing a good starting point for optimization. These techniques provide improvements over brute-force techniques. In some embodiments, parameters identified for level- QAOA, for any q<p , can be used to more quickly optimize parameters for level-p QAOA. Further, while some examples discuss randomly generated instances of u3R and w3R, similar results can be found when applying these techniques to u4R and w4R graphs, as well as complete graphs with random weights (or the Sherington-Kirkpatrick spin glass problem. These other exemplary u3R, w3R, u4R, and w4R graphs are regular graphs, meaning, for example, that each vertex has the same number of neighbors (3 or 4 respectively). The letters u and the w can refer to whether one considers unweighted or weighted graphs, respectively. In some embodiments, these graphs are useful as test samples. The patterns in the optimal parameters identified herein are used to develop example heuristic strategies that can efficiently find quasi-optimal solutions in 0(poly(p)) time.

[0120] In some embodiments it is possible to eliminate degeneracies in the parameter space due to symmetries. For example, generally, QAOA can have a time-reversal symmetry, F _p (j b ) = _p (-y,— /?), since both HB and He are real-valued. For QAOA applied to MaxCut, there can be an additional Z 2 symmetry, as _{e -} ^l O/2

commutes through the circuit. Furthermore, the structure of the MaxCut problem on u<iR graphs can create redundancy since e ^mHc = 1 if d is even, and (<r ^z ) ^®w if d is odd. These symmetries make it possible to restrict bi E in general, and gi e for ut/R graphs.

[0121] In some embodiments, it is possible to numerically investigate the optimal QAOA parameters for MaxCut on random u3R and w3R graphs with vertex number 8 < N< 22, using a brute-force approach. For each graph instance and a given level p, a random starting point (seed) in the parameter space can be chosen and a gradient-based optimization algorithm such as Broyden-Fletcher-Goldfarb-Shanno (“BFGS”) can be used to find a local optimum ( g ⁱ , b ⁱ ) starting with this seed. In some embodiments, a local optimum can refer to a case where for a given choice of parameters b and y, the result always gets worse if the values of the parameters b and y are changed slightly. However, for such local optima, the result might get better if the values of these parameters are instead changed drastically. Thus the optimum can be referred to as local optimum, not global optimum. In some

embodiments, for each graph, it is possible to optimize the variational parameters to cause the solution to go to a local optimum by a local optimization method, such as one where the optimization only searches parameters close to the initial starting parameters. This local optimization can be repeated with sufficiently many different seeds to find the global optimum (y ^* ,/? ^* ). In some embodiments, a global optimum may refer to the case where there is not a better choice of parameters. The global optimum can change from graph to graph. The degeneracies of the optimal parameters (y ^* ,/? ^* ) can be reduced using the symmetries discussed above (e.g., by finding some (distinct) values of y and b that are equivalent because they lead to exactly the same result and thus do not need to be considered individually). In some illustrative examples, the global optimum can be non-degenerate up to these symmetries.

[0122] In some embodiments, the process of identifying optimal parameters (y ^* ,/? ^* ) can be repeated for additional random graphs, such as 100 u3R and w3R graphs with various vertex numbers N , which can draw out one or more patterns in the optimal parameters (g , b ). In one example, the optimal gΐ can tend to increase smoothly with each iteration i = 1, 2, · , >, while b can tend to decrease smoothly. For example, FIGS. 15A-15B show an example instance of optimal QAOA parameters (y ^* , /? ^* ) obtained for example instances of MaxCut on a 16-vertex unweighted 3 -regular (u3R) graph at level p 7, according to some embodiments. As shown in FIGS. 15A-15B, the parameters for level p 7 transition smoothly for each iteration i.

[0123] FIGS. 15C-15H shows a parameter pattern visualized by plotting the optimal parameters of 40 instances of 16-vertex u3R graphs (such as the example used to plot FIGS. 15A-15B), for 3 <p < 5 as a function of iteration according to some embodiments. Each dashed line connects parameters for one particular graph instance. For each instance and each p , the classical BFGS optimization algorithm was used from 10 ⁴ random seeds (starting points), with the best parameters being kept. Similar patterns can be found for w3R graphs, which are discussed in more detail below. Furthermore, in some examples, for a given class of graphs, the optimal parameters can be observed to roughly occupy the same range of values as p is varied. This demonstrates a pattern in the optimal QAOA parameters that can be exploited in the optimization, as discussed in more detail below. Similar patterns are found for parameters up to p 5 15, if the number of random seeds is increased accordingly.

[0124] FIGS. 21A-21F show graphs demonstrating optimal QAOA parameters (g ^* ,b ^* ) for 40 instances of l6-vertex weighted 3-regular (w3R) graphs, for 3 <p < 5, according to some embodiments. Each dashed line connects parameters for one particular graph instance. For each exemplary instance and each p , the BFGS algorithm is used to optimize from 10 ⁴ random starting points, and the best parameters are kept as (y ^* ,/T). FIGS. 21A-21F are analogous to the results of unweighted graphs in FIGS. 15C-15H. The spread of f ^* for weighted graphs is wider than that for unweighted graphs shown in FIGS. 15C-15H. In some embodiments, this is because the random weights effectively increase the number of subgraph types. Moreover, the larger value for g ^* for weighted graphs compared to unweighted graphs can be understood via the effective mean-field strength that each qubit experiences.

[0125] Notably, even at small depth, this parameter pattern can be reminiscent of adiabatic quantum annealing where He is gradually turned on while HB is gradually turned off, in some embodiments. However, the mechanism of QAOA can be shown to go beyond the adiabatic principle, as discussed in more detail below. In addition, in some embodiments, the optimal parameters can have a small spread over many different instances. This can be because the objective function F _P (y, /?) can be a sum of terms corresponding to subgraphs involving vertices that are a distance < p away from every edge. At small p, there are only a few relevant subgraph types that enter into F _P and can effectively determine the optimal parameters. As N ® ¥ and at a fixed finite p , the probability of a relevant subgraph type appearing in a random graph can approach a fixed fraction. This implies that the distribution of optimal parameters (y ^* , b can converge to a fixed set of values in this limit.

[0126] In some embodiments, the optimal parameter patterns observed above can indicate that generically, there is a slowly varying continuous curve that underlies the parameters y _t ^* and b[ . In some embodiments, this curve changes only slightly from each level p to p+ 1. Based on these observations, a new parameterization of QAOA can be used, as well as a heuristic optimization strategy that can, without limitation, be referred to as“FOURIER.” In some embodiments, the heuristic strategy uses information from the optimal parameters at level p to help optimization at level p + 1 by producing good starting points. While this heuristic does not necessarily find the global optimum of QAOA parameters, it can produce, in 0(poly(p)) time, quasi-optima that can only be surpassed with 2 ^oip) number of brute-force runs. In some beneficial embodiments, this facilitates evaluation of the performance and mechanism of QAOA beyond p = 1.

[0127] In some embodiments, QAOA can be re-parameterized based on the observation that the optimal QAOA parameters y _t ^* and b[ appear to be smooth functions of their index i. An alternative representation of the QAOA parameters can be implemented using Fourier- related transforms of real numbers, which permits reductions in the number of necessary parameters. Instead of using the 2p parameters (y ^* ,/? ^* ) e M ^2p , FOURIER can instead use 2 q parameters (u, v) E M ^2q , where the individual elements Yi and bi are written as functions of (u, v) through the following transformation:

Equation 10:

Equation 11:

[0128] In some embodiments, these transformations can be referred to as Discrete Sine/Cosine Transforms, where m and vk can be interpreted as the amplitude of k- th frequency component for g and /?, respectively. When q ³p, this new parametrization can describe all possible QAOA protocols at level p.

[0129] Embodiments of the FOURIER strategy work by starting with level p = 1, optimizing using an optimization function such as brute force for level p 1, and then using the optimum at level p to determine a starting point for level p + 1. The starting points can be generated by re-using the optimized amplitudes ( u ^* , v ^* ) of frequency components from level p extrapolated from the optimized parameters (g ^* ,b ^* ) to identify the parameters (g ^* ,b ^* ) for the level p + 1. This can be repeated for increasing p.

[0130] Some embodiments include several variants of this strategy, examples of which are referred to as FOURIER[c/, R] and INTERP, for optimizing -level QAOA. Without limitation, embodiments of one of variants can be referred to as FOURIER[c/, R],

characterized by two integer parameters q and R. The first integer q can refer to the maximum frequency component allowed in the parameters (u, v), which can be the maximum value of k. If q =p, the full power of /-level QAOA can be utilized. However, since the smoothness of the optimal parameters (j, ^* /?) implies that only the low-frequency components are important, it is also possible to consider the case where q is a fixed constant independent of, but smaller than p , so the number of parameters is bounded even as the QAOA circuit depth increases.

[0131] In some embodiments, the second integer R can refer to the number of controlled random perturbations added to the parameters to escape a local optimum towards a better one. For example, where the optimization parameters (g ^* , /?) were identified at a local but not global optimum for the initial value p = 1, perturbations can be introduced to avoid focusing only on that local optimum for increased values of p. Exemplary results discussed in the present disclosure implement the FOURIER[c/, R] strategy with q =p and R = 10 unless otherwise stated, but such as selection is not limiting.

[0132] In embodiments where q is chosen such that q =p, the strategy can be denoted as FOURIER[oo, //], since q grows unbounded with p. In embodiments of the FOURIER[oo, 0] variant of this strategy, a starting point is generated for level p + 1 by adding a higher frequency component, initialized at zero amplitude, to the optimum at level p. For example, as shown in FIG. 22, according to some embodiments, from the optimized variational parameters 2212 (¾_ _!) , i¾-i ₎ ) ^e D¾ ^2/,_i at level p - 1, the starting point 2222A (W _(p) , ^ _(p)) G R ^2(p) is generated according to:

Equation 12: u° _p) = (¾_ _1} , 0), ¾ = (¾_ _1} , 0).

Eising (ϊί _(r) , 2222A as a starting point, a BFGS optimization routine can be performed to obtain a local optimum 2224B (¾ ₎ , ^ _(P) ) for the level p. This is output to the next level of p, as the process continues.

[0133] In some embodiments, improvements over this technique can be gained with the strategy, FOURIER[oo, R > 0], which is also shown in FIG. 22. In some embodiments, the strategy FOURIERfoo, 0] can become stuck at a sub-optimal local optimum. Accordingly, perturbing its starting point as is performed in FOURIER[oo, R > 0] can greatly improve the performance of QAOA already seen with FOURIERfoo, 0]

[0134] As shown in FIG. 22, the technique begins at level p 1, according to some embodiments. In some embodiments of FOURIERfoo, R > 0], in addition to optimizing according to the strategy FOURIER[oo, 0], the / level QAOA can be optimized from R+ 1 extra starting points, R of which are generated by adding random perturbations to the best of all local optima ( u _p _ _!) , ¾-i ₎ ) found at level p-1. Specifically, for each instance at />-level QAOA, and for r = 0, 1, . . . , R, optimization can start from contain random numbers drawn from normal distributions with mean 0 and variance given

Equation 15: Normal

Equation 16: Normal

[0135] In such embodiments, there is a free parameter a corresponding to the strength of the perturbation. In some non-limiting examples derived from trial and error, setting a = 0.6 can yield good results. This choice of a is assumed in the present disclosure unless otherwise stated. However, a person of skill in the art would understand that other values of a can be used, including dynamic values of a as needed depending on particular implementations.

[0136] In some embodiments, additional strategies can be used to take advantages of the parameter pattern indicated above. One exemplary strategy can use linear interpolation of optimal parameters at lower level QAOA to generate starting points for higher levels, which can be referred to without limitation as“INTERP.” Both INTERP and FOURIER strategies are effective for the instances discussed throughout the present disclosure and are applicable to others as well. While FOURIER has demonstrated a slight edge in its performance in finding better optima when random perturbations are introduced, a person of skill in the art would understand from the present disclosure that INTERP is also an efficient way of improving QAOA and can present additional benefits. Embodiments of FOURIER[q, R] and INTERP are described below in more detail. However, additional techniques are

contemplated by the present disclosure, such as the use of machine learning. Furthermore, although in aspects of the present disclosure, the heuristic strategies makes use of optimal variational parameters found at level -(/ -/) QAOA to find initial variational parameters at level-/; QAOA, a person of skill in the art would understand that optimal variational parameters found at level-/??, for any m<p, can be used to design initial variational parameters at level-p QAOA.

[0137] In some variants of the FOURIER strategy, the number of frequency components q is fixed. These variants can be treated similarly to the strategies where q =p discussed above, except all u and v parameters can be truncated at the first q components. For example, when optimizing QAOA at level p ³ q with the FOURIER[c/,0] strategy, no further higher frequency components are added, and the starting point begins at G

[0138] In some embodiments of the optimization strategy referred to as INTERP, linear interpolation can be used to produce a starting point for optimizing QAOA and an

optimization routine can iteratively increase the level p. However, for purposes of the present discussion, p should be considered the same as 1 in the discussion for the FOURIER strategy, as this is simply a matter of semantics for where to begin the algorithm. In some embodiments, this is based on the observation that the shape of parameters (7 _(p+1) , b _{(r+1) )} closely resembles that of (7 _(p) , /¾, ₎ ) _· For a given instance, QAOA is iteratively optimized by starting from p = 1, obtaining local optimum parameters (7 _(p) , ¾ ₎ ), and incrementing p top + 1. To optimize parameters for level p+l, optimized parameters from level p are used to produce starting points (7 ₍ ° _{p+1) 1)} ) according to:

Equation 17:

for i = 1, 2, ·,r + 1, where i denotes the ?-th element of the vector. Here, [7] _{/ º} 7, denotes the ?-th element of the parameter vector 7, and [7 _(P) ] _Q = [7 _(p) ] _p+1 ^º 0 _· The expression for b(r+i) ^can the same as above after swapping y <-> b. Starting from (K ₍ ° _r+1) , ? ₍ ° _p+i) X the BFGS optimization routine (or any other optimization routine) can be performed to obtain a local optimum (K _(r+1) , /¾ _>+1) ) for the (p + l)-level QAOA. Finally, » can be incremented by one and the same process can be repeated until a target level is reached.

[0139] The INTERP strategy can also get stuck in a local optimum in some embodiments. Adding perturbations to INTERP can help but may not be as effective in some embodiments as with FOEIRTER. This may occur because the optimal parameters are smooth, and adding perturbations in the (u, u)-space modify (y, b) in a correlated way, which can enable the optimization to escape local optima more easily. However, a similar perturbation routine is contemplated.

[0140] As discussed above, the heuristic approaches described in the present disclosure constitute a significant improvement over brute force QAOA techniques. Non-limiting comparisons of example implementations are discussed below in the sections below.

[0141] Based on the present disclosure, a person of skill in the art would understand that the disclosed heuristic strategies could be implemented on a number of technical platforms.

In the section below titled“Example QAOA Implementations” the MaxCut problem is considered as an example, although it can also be applied to solve other interesting problems.

Example implementations with quantum systems

[0142] In some embodiments, large-size problems are suitable for implementation on quantum systems. Two aspects of such implementations (reducing the interaction range and examples with Rydberg atoms) are discussed in more detail below.

[0143] First, with regard to reducing the interaction range, in some quantum

implementations, as discussed above, each vertex can be represented by a qubit. For a large problem size, a major challenge to encode general graphs is the necessary range and versatility of the interaction patterns (between qubits). The embedding of a random graph into a physical implementation with a 1D or 2D geometry may require very long-range interactions. By re-labelling the graph vertices, it is possible reduce the required range of interactions. Without being bound by theory, this can be formulated as the graph bandwidth problem. Given a graph G = (V, E) with N vertices, a vertex numbering is a bijective map from vertices to distinct integers, f: V ® {1,2, ··· , N}. The bandwidth of a vertex numbering / is, Bf(G ) = max{\f(u )— f(v) |: ( u , v ) e E(G)} which can be understood as the length of the longest edge (in 1D). The graph bandwidth problem is then to find the minimum bandwidth among all vertex numberings, i.e., B(G ) = min/B/(G); namely, it is to minimize the length of the longest edge by vertex renumbering.

[0144] In general, finding the minimum graph bandwidth is NP-hard, but good heuristic algorithms have been developed to reduce the graph bandwidth. FIGS. 20A-20E show a simple example of bandwidth reduction, according to some embodiments. FIGS. 20A, 20B illustrate the vertex renumbering with a 5-vertex graph. FIG. 20E shows the histogram of graph band-widths for 1000 random 3 -regular graphs of N= 400 each. ETsing the Cuthill- McKee algorithm, the graph bandwidth can be reduced to around B ~ 100. While this still may require quite a long interaction range in 1D, the bandwidth problem can also be generalized to higher dimensions. In 2D arrangements, the radius of interaction can be considered to be B _2D ~ 100 for N = 400 3-regular graphs, which is possible for quantum devices. FIGS. 20C and 20D show sparsity patterns of the adjacency matrix before and after vertex renumbering for one exemplary graph. While, in this illustrative example, the vertex renumbering technique is applied to the MaxCut problem, a person of skill in the art would understand based on the present disclosure that the same technique can be used for other combinatorial optimization problems such as maximum independent set problems. [0145] In some embodiments, a general construction can be used to encode any long- range interactions to local fields by including additional physical qubits and gauge constraints. It is also possible to restrict to special graphs that exhibit some geometric structures. For example, unit disk graphs are geometric graphs in the 2D plane, where vertices are connected by an edge only if they are within a unit distance. These graphs can be encoded into 2D physical implementations, and the MaxCut problem is still NP-hard on unit disk graphs.

[0146] In some embodiments, the above discussion of QAOA has been platform independent, and is applicable to any state-of-the-art platforms. Exemplary platforms include neutral Rydberg atoms, trapped ions, and superconducting qubits. Although the following discussion focuses on an implementation of QAOA with neutral atoms interacting via Rydberg excitations, where high-fidelity entanglement has been recently demonstrated, other implementations are contemplated.

[0147] In some embodiments, it is possible to implement control over the interaction between individual atoms according to a given graph. In some embodiments, in this way it is possible to control which types of problems are to be solved. Since the interactions can specify the problem, controlling the interactions is one way to control the problem. As discussed in more detail above, in an exemplary Rydberg implementation, the hyperfme ground states in each atom can be used to encode the qubit states |0) and |l), and the 11 ) state can be excited to the Rydberg level | r) to induce interactions. The qubit rotating term, exp can be implemented by a global driving beam with tunable durations. The

interaction terms a can be implemented stroboscopically for general graphs; this can be realized by a Rydberg-blockade controlled gate, as illustrated in FIG. 20F. For example, FIG. 20F shows a protocol to use Rydberg-blockade controlled gate to implement the interaction term exp For example, FIG. 20F shows how to realize a specific component of the interaction term by first applying a pi-pulse to a control qubit, second applying a pulse with coupling strength W and detuning D to a target qubit, and following this pulse with another pi-pulse on the control qubit. By choosing a proper gate time for the second step, such as (t = 2p / Vl2 ² + D ² ), the population on the Rydberg level |r) may not be accumulated. With tunable coupling strength W and detuning D, it is possible to control the interaction time y .

[0148] By controlling the coupling strength W, detuning D, and the gate time, together with single-qubit rotations, it is possible to implement exp(— iyafa ), which can then be repeated for each interacting pair. In this way, it is possible to choose which pairs should interact, and thus control which problem to solve. In some embodiments, an additional advantage of the Rydberg-blockade mechanism is its ability to perform multi-qubit collective gates in parallel. This can reduce the number of two-qubit operation steps from the number of edges to the number of vertices, /V, which means a factor of N reduction for dense graphs with ~N ² edges. While the falloff of Rydberg interactions may limit the distance two qubits can interact, MaxCut problems of interesting sizes can still be implemented by vertex renumbering or focusing on unit disk graphs, as discussed above. Furthermore,

implementing ancillary vertices discussed in the present disclosure can be used to increase the length of interactions.

[0149] Without being bound by theory, in some embodiments, for generic problems of 400-vertex regular graphs, the interaction range can be on the order of 5 atoms in 2D. This can be determined by assuming a minimum inter-atom separation of 2 «m, which means an interaction radius of 10 «m, which is realizable with high Rydberg levels. Given examples of coupling strength W ~ 2p X 10-100 MHz and single-qubit coherence time t ~ 200 u (limited by Rydberg level lifetime), with high-fidelity control, the error per two-qubit gate can be made roughly (Wt) ^-1 ~ 10 ^-3 - 10 ^-4 . For 400-vertex 3-regular graphs, QAOA of level p— Wt/N~ 25 can be implemented with a 2D array of neutral atoms. Advanced control techniques such as pulse-shaping would increase the capabilities of QAOA in such systems. Furthermore, QAOA may not be sensitive to some of the imperfections present in existing implementations with Rydberg atoms.

[0150] The following sections explore additional examples and embodiments of the present disclosure. The present disclosure is not limited by the theory described herein, which is merely meant for illustration of some aspects of operational principles that underly some embodiments of the present disclosure.

Quantum annealing for random UD-MlS

[0151] In some embodiments, quantum optimization algorithms such as quantum annealing algorithm (QAA) can be used for random UD-Mis problems. As discussed above, the maximum independent set problem on random unit disk (UD) graphs is only one type of problem that is contemplated by the present disclosure. For QAA with UD graphs, a random UD graphs can be parameterized by two parameters: the number of vertices N and the 2D vertex density p. As shown in FIG. 5, the unit distance can be taken to be r = 1, and the vertices are put into a box of L x L, where L = ^N/p, in some embodiments. In particular, FIG. 5 shows an example of a random unit disk graph with N= 40, p = 1.5, |MIS| = 14, and 93 edges. Without being bound by theory, for a random UD graph with density p, the average vertex degree is approximately pr. To minimize the finite-size boundary effect, periodic boundary conditions are used for UD graphs in numerical simulations discussed in the present disclosure. However, such conditions are not necessary.

[0152] As discussed above, a QAA for MIS can be performed using the following

Hamiltonian: Equation 18:

[0153] The QAA can be designed by first initializing all qubits in |0) at time / = 0, which is the ground state of F/QA(/ = 0) when D(/ = 0) < 0 and Q(t = 0) = 0 (with U > 0). The parameters can then be changed, for example by first turning on W(I) to a non-zero value, sweeping D(/) to a positive value, and finally turning off W(ή again. Without being bound by theory, the exemplary annealing protocol discussed throughout the present disclosure can be specified by

Equation 19: s = t/T.

[0154] If the time evolution is sufficiently slow, then by the adiabatic theorem, the system can follow the instantaneous ground state, ending up in the solution to the MIS problem. Wo = 1 can be considered the unit of energy, and it is possible to fix Do/Wo = 6, which in non-limiting examples has been identified as a good ratio to minimize nonadiabatic transitions.

[0155] In some embodiments, quantum annealing can be explored on random unit disk graphs, with N vertices and density p. In some embodiments, in the limit of Do, Wo « U, the non-independent sets are pushed away by large energy penalties and can be neglected. In some implementations, this can correspond to the limit where the Rydberg interaction energy is much stronger than other energy scales. Without being bound by theory, in some examples in this limit, the wavefunction can be restricted to the subspace of all independent sets, such as:

Equation 20: H _IS = {\ip): n _v n _w \ip) = 0 for any (v, w) E E} in the exemplary numerical simulation discussed herein, which allows for access to a much bigger system size up to N~ 50 since dim (H _IS ) « 2 ^v . First, in an example, the subspace of all independent sets can be found by a classical algorithm, the Bron-Kerbosch algorithm, and the Hamiltonian in equation 18 can then be projected into the subspace of all independent sets. The dynamics with the time-dependent Hamiltonian can be simulated by dividing the total simulation time t into sufficiently small discrete time steps t and at each small-time step, a scaling and squaring method with a truncated Taylor series approximation can be used to perform the time evolution without forming the full evolution operators.

[0156] In some embodiments, exemplary time scales for adiabatic quantum annealing to perform well can be explored. In some examples, this time scale can be governed by the minimum spectral gap, £¾ap, where the runtime required can be T = 0 However, the minimum spectral gap can be considered to be ambiguous when the final ground state is highly degenerate, since it is possible for the state to couple to an instantaneous excited state as long as it comes down to the ground state in the end. For an example generic graph, there can be many distinct maximum independent sets (the ground state of Hp can be highly degenerate). So instead of finding the minimum gap, a different approach can be taken to extract the adiabatic time scale.

[0157] In some embodiments, in the adiabatic limit, the final ground state population (including degeneracy) can take the form of the Landau-Zener formula PMIS ~ 1 - e ^a~T ^ ^TLZ , where a is a constant and TLZ is the adiabatic time scale. In embodiments of the

nondegenerate case, it is typical to have TLZ = 0(l /€ _gap ). In the more general case, the time scale TLZ can be extracted by fitting to this expression. However, in some examples, the simple exponential form holds only in the adiabatic limit, where T > TLZ. Hence, for each graph instance, it is possible to search for the minimum T such that the equation holds. For example, T can be adaptively doubled iteratively (from 5) until the minimum T* is

found such that MIS > 0.9, at which it is possible to assume, in some embodiments, that the time evolution lies in the Landau-Zener regime. The dynamics can then be simulated for another three time points 1.57\ 2/\ and 2.57\ before finally fitting to the equation from T* to 2.5 T* to extract the time scale TLZ.

[0158] The fitting has been shown in some examples to be effective for most instances. For example, FIG. 6 shows the Landau-Zener fitting to 1 -/Vis = e ^a~T R ^Lz to extract the adiabatic time scale TLZ. Here, four random unit disk graphs with N= 10, 20, 30, 40 are plotted. For each exemplary instance, the first T is found iteratively such that / tis( 7) > 0.9, denoted as T*. The fitting is then performed on four points T*, 1.5 T*, 2 T*, 2.5 /-' to extract the time scale TLZ. It is possible to drop the few graphs where the goodness-of-fit R ² < 0.99. In some embodiments, this procedure can be performed to extract TLZ for up to 30 graph instances at each N and p. From these, the median can be taken, which can be used to produce the full phase diagram in terms of TLZ as plotted in FIG. 4A described above.

[0159] FIGS. 7A-7C show the adiabatic time scale TLZ at some fixed densities. More particularly, FIGS. 7A-7C show 200 random unit disk graphs that are simulated for QAA for system sizes up to N = 46. FIG. 7A plots the median TLZ, while FIGS. 7B and 7C plot TLZ for individual instances for p = 0.8 and p = 3. FIGS. 7A-7C show the scaling of TLZ with N at some fixed densities p = 0.8 (below the percolation threshold) and p = 3 (above the percolation threshold). FIG. 7A shows a clear separation between p = 0.8 and p = 3, according to some embodiments. However, the scaling of N is unclear, which can be due to a finite-size effect: N ³ 100 may be a better condition to show scaling. FIGS. 7B and 7C also show the spread of TLZ for each instance. It can be seen that some exemplary hard instances require significantly longer TLZ than the typical instance, though even on average TLZ > 10 for p = 3, N ³ 20.

[0160] While some embodiments discussed above focused mainly on the capacity of exemplary algorithms to solve MIS exactly, it is also possible to identify whether the algorithms can solve MIS approximately, for example in the sense of finding an independent set as large as possible. For some exemplary quantum algorithms, without being bound by theory, the approximation ratio r can be used to gauge performance for approximations. For example, for a quantum algorithm (such as a quantum annealer) that outputs a state \ipf), it is possible to write r = where |MIS| is the size of the MIS. In other words, r can quantify the ratio of the average independent-set size found by measuring the output quantum state to the maximum independent-set size. FIG. 8 shows an exemplary analogous phase diagram in terms of the approximation ratio r by simulating quantum annealing at a fixed time T= 10/Wo, averaged over 30 graph instances per N and p (the same instances as in FIG. 4B discussed above), according to some embodiments. The vertical dashed line corresponds to the percolation ratio at = p _{c «} 1.436, while the sloped dashed lines correspond to optimal disk packing. FIG 8 shows qualitatively the same features as the ground state population discussed above with reference to FIG. 4B, but the finite-size effect is stronger due to the small discrete |MIS| values at large densities.

Generalizations to arbitrary graph structure beyond UD graphs

Stroboscopic evolution

[0161] As discussed above, it is possible to generalize the exemplary implementation discussed above to address MIS problems on graphs, G = (V, E) that are beyond the UD paradigm.

[0162] In some embodiments, the quantum algorithms discussed above with reference to the UD paradigm involve evolution with a Hamiltonian H(t) = å _n W _n (ΐ s^— D(ΐ)h _n + å _{(u V)EE} Un _u n _{v .} In exemplary embodiments where U »|W|,|D|, the dynamics can be effectively restricted to the independent set space His. Without being bound by theory, in some embodiments to generate such evolution with a Hamiltonian corresponding to a general graph structure, a Trotterized version of the time evolution operator can be considered: Equation 21:

where the time interval [0,7] is sliced defining times tj such that å _j t _j = T and t _j+1 — t _j « ^]D _max n(t _j ), \A(t _j )\. Here /) max can denote the maximum degree of the graph. Each U(/y) can be further Trotterized as follows

[0163] In other words, it can be split into a product of terms U _v that each are associated with the evolution of one spin, v. Here N (v) can denote the neighbors of v on the graph.

Note that in the U ¥ limit this can be written as

Equation 23: ] ^~ [ |0) _u (0|

UEN( _p )

[0164] This is a simple single qubit rotation of atom v , conditioned in some embodiments on the state of the atoms corresponding to neighbors of v being |0). If at least one of the neighbors is in state 11), atom v may not evolve.

Implementation using qubit hyperfme encoding

[0165] In some embodiments it is desirable to add further control of a quantum system in order to specify interactions between particular pairs or subsets of qubits. One exemplary approach to implement the corresponding dynamics with individually controlled neutral atoms can take advantage of qubit states |0) and |1) encoded in two (non-interacting) hyperfme states, in the internal atomic ground state manifold. In other words, when implementing other forms of graphs, it is possible to iterate through desired interactions between individual pairs (or sets) of qubits.

[0166] In some embodiments, the atoms can be positioned on the points of a 2D square lattice with distance g. To realize a single step, U _v (t _j ), the atoms, w, can be excited that correspond to neighbors of v on the graph (u e K (v)), selectively from the state 11) to a Rydberg Estate \ V). In some embodiments, a grid length g « r _B can be chosen such that none of the atoms u £ ( v ) interact during this process. Atom v can then be driven to realize the single qubit rotation in the hyperfme manifold, for example with a unitary, such as, but not limited to, a physical operation that can be performed deterministically such as driving the atom with a laser to cause a quantum state evolution corresponding to an evolution with W _n (ΐ^s^— D _n (ί _/· )h _n , where s couples the two hyperfme qubit states of atom v, and n _v counts if atom v is in hyperfme state |1). In some embodiments, to realize this rotation, an individually addressed two-step excitation can be used that can couple the two hyperfme states |0) and 11) of atom v via a transition through a Rydberg / ^J - state (for example, by first exciting the atom from the state |0) to a Rydberg P state with a first laser pulse, and then changing its state from the Rydberg P state to the state 11) with a second laser pulse). In some embodiments, if all atoms u are in the state |0), then this process may not be disturbed, but if at least one of the neighbors is the Rydberg S state, the strong S I* interaction can give rise to a blockade mechanism that prevents the rotation of the qubit v, thus realizing exactly equation 23. Note that in some embodiments this may require a different scale of interaction length of the blockade radius for two atoms in the Rydberg S- states on one hand, and the S-P blockade radius on the other hand because these two interactions scale differently with separation of the atoms: the S-P interactions decay as l/x ³ (much slower than the S-S interactions that scale like l/x ⁶ ), which makes it possible to implement collective gate with high fidelity. Accordingly, it is possible to apply quantum optimization to MIS problems on graphs that are more general than the unit disk graphs.

Maximum Independent Sets for Unit Disk Graphs

[0167] Without being bound by theory, this section addresses some aspects of Mis problems on unit disk graphs that can be solved according to the techniques described in the present disclosure. FIG. 5 shows a unit disk graph with vertices 502, 504, some of which are connected by edges 506 when they fall within a radius r of one another, as represented by disks 508, according to some embodiments. The problem can be formulated as minimizing the energy of:

Equation 24:

where:

Equation 25:

which is a subclass of Hp (2). This problem can be proved to be NP-complete by reducing it from Mis on planar graphs of maximum degree 3. Since, in some embodiments, analysis of the computational complexity associated with the Rydberg Hamiltonian (Equation 2) discussed herein is based on a similar reduction, it can be instructive to review the following theorem and its proof: Mis on unit disk graphs is NP-complete .

[0168] Without being bound by theory, this theorem can be proven as follows: (1) Mis on planar graphs with maximum degree 3 is NP-complete. (2) A planar graph Cj = (V, e ), with vertices V and edges 8, with maximum degree 3 can be embedded in the plane using 0( \V\) area in such a way that its vertices are at integer coordinates, and its edges are drawn by joining line segments on the grid lines x = i x g or y = j x g, for integers i and j, and grid spacing g. (3) One can replace each edge {u, v} of Q by a path having an even number, 2ku,v , of ancillary vertices, in such a way that a UD graph, G = ( V, E) with vertices V and edges E, can be constructed. It is straightforward to verify that Q has an independent set S i V such that \S\ < a if and only if G has an independent set S a V such that |S | < a' = a +

[0169] In some embodiments, this theorem shows that it is NP-complete to decide whether the ground state energy of HUD is lower than -a’A. In some embodiments, the transformation in the proof of this theorem does not fully determine the actual positions of the ancillary vertices in the 2D plane. In some embodiments, a particular arrangement can be specified consistent with the requirements of this transformation. Once Rydberg interactions are considered, the interaction strength between each pair of qubits can be fixed in a way that takes into account the distance of the atoms.

[0170] Given an edge {u, v} of the graph Cj embedded on a grid, the length of this edge can be denoted by g x £ _{u v} with £ _{u v} an integer, and g the grid unit length, in some embodiments. First, an ancillary vertex can be placed on the £ _{u v} — 1 grid point along the edge, separating the edge in £ _{u v} segments of length g. An integer k > 3 can be selected and equally spaced 2k ancillary vertices can be placed along each segment, dividing it into 2k + 1 pieces of equal length d = g l{2k + 1). In some embodiments, irregular vertices can be used to ensure that the number of ancillary spins on each edge is even. In some embodiments, however, the total number of spins on each edge might be different. If £ _{u v} is even, one segment can be chosen and the 40 + 2 ancillary vertices close to the center of this segment can be replaced with 40 + 1 ancillary vertices, that are equally spaced by a distance D = d + for some integer ø to be determined. Such segments can be referred to as“irregular segments”, and to the vertex at the center of the irregular segment can be referred to as an irregular vertex. In some embodiments, these exceptions are made to ensure that the total number of ancillary vertices along each edge {u, v} e £, 2 k _u ,v, is even in order to ensure that the ancillary vertices transfer the independent set constraint without changing the nature of the problem to be solved. Following this arrangement, the nearest-neighbor distance of the ancillary vertices can be either d or D. Setting the unit disk radius to r = I) + CT can produce the unit disk graph G. The positions of the vertices can be labelled by x _{v .} In some embodiments, arrangement depends on the freely chosen parameters k and f. Accordingly, as described throughout the present disclosure, a hard problem can be transformed into an MIS problem on an arrangement of vertices that form a unit disk graph.

[0171] Without being bound by theory, a few properties of the maximum independent set on the types of unit disk graphs constructed in this way can be noted, according to some embodiments. First, in some embodiments, the maximum independent set on G is in general degenerate, even if the maximum independent set on Cj is unique. The properties of a particular set of MIS-states on G can be considered, such as MIS-states on G, | ip _G ), that coincide with MIS-states on Cj , \ipg), on the vertices V. Without being bound by theory, to see that such states can always exist, \tp _G ) can be constructed from | ipg) as follows: given the state of two qubits corresponding to two vertices, w, v e V, if u is in the state |0) and v is in state |1) (or vice versa) then the 2 k _u ,v ancillary vertices can be placed on the edge connecting u and v in an (antiferromagnetically) ordered state, where the qubits are alternating in states |1) and |0). In the other case, if u and v are both in states |0), the ancillary vertices can be placed in an analogously ordered state. In embodiments of this latter case, a domain wall, such as an instance where two neighboring qubits are both in the state |0) can be introduced. The position of this domain wall along the edge is irrelevant. In both cases half of the 2 k _u ,v ancillary vertices along the edge are in state |1). Therefore, by the above theorem, the state constructed from | i g) by applying this process on all edges of Q, is a MIS-state on G.

[0172] Without being bound by theory, the structure of these MIS-states around points where edges meet under a 90° angle can be further specified, according to some

embodiments. Such points can be either junctions, where 3 edges meet at a vertex, such as point 216 in FIG. 2A, or corners, such as 202 in FIG. 2A. There is a MIS-state, \ip _G ), such that the qubits close to each corner or junction are ordered (for example, in a state where every other vertex is in state |1) and the others are in state |0)). For example, there exists a maximum independent set, such that for every corner and junction, all vertices within a distance d < g/4 are in one of the two possible ordered configurations with no domain wall. Without being bound by theory, this follows from the above discussion by noting that the position of any domain wall along an edge £ is irrelevant, in some embodiments. For example, any domain wall can be moved along an edge £ such that its distance to any vertex on a grid point is larger than g/4. The possibility to move domain walls is also exploited in the discussion of the full Rydberg problem discussed, for example below in the section titled NP-Completeness of the Rydberg Problem.

Model Detunings for Corners and Junctions

[0173] Without being bound by theory, to illustrate applications of the above reduction to implementations that employ Rydberg interactions, a simple model can be implemented to explain aspects of some implementations. This model can be used to show some aspects and benefits of embodiments of the present disclosure, including, but not limited to treatment of special vertices. For example, a Hamiltonian similar to HUD, the MIS-Hamiltonian for UD graphs, can be considered with the introduction of interactions beyond the unit disk radius. Similar to the situation in the Rydberg system, these additional interactions can cause energy shifts that can result in a change of the ground state, thus invalidating the encoding of the Mis. In other words, such additional interaction can cause the ground state of an encoded Mis problem to be mismatched with the solution to the Mis problem. To resolve this issue, local detunings can be used to compensate for the additional interactions.

[0174] Without being bound by theory, embodiments of the model can be expressed as: Equation 26:

with interactions given by:

where W < U. For W = 0 and D _n = D, //model can reduce to the Hamiltonian HUD described in Equation 25. For W> 0 it includes interactions beyond the unit disk radius, r, and can thus be considered as a first approximation to the Rydberg Hamiltonian.

[0175] In some embodiments, in the case of a qubit arrangement described in the section titled“Maximum Independent Sets for Unit Disk Graphs” corresponding to a unit disk graph, G, and the case 2r < R < 2r, most qubits interact only with their neighbors on G, with an exception being qubits that are close to corners (such as qubit 202 in FIG. 2A) or junctions (such as qubit 216 in FIG. 2A). In such geometric arrangements, some qubits that are not neighboring on G interact with strength W> 0. These interactions can, in some embodiment, change the ground state, as compared to the MIS-Hamiltonian, Hr.

[0176] FIGS. 11 A-l 1C show examples of vertex arrangements corresponding to a unit disk graph, where the ground state of the model //model (Equation 26) does not correspond to the maximum independent set, according to some embodiments. The maximum independent set on this graph in FIG. 11 A is indicated by the black vertices, such as vertices 1106 and 1156. The white vertices, such as vertices 1102, 1104, 1152 indicate another independent set whose size is smaller than the size of the maximum independent set by one vertex. If W is too large, the additional interactions on the diagonals 1130, 1160 close to comers 1102 and junctions 1152 disfavor the MIS-state; instead, the state where the black qubits are in state 11), becomes the ground state. As described throughout the present disclosure, this problem can be resolved by changing the detuning locally at the problematic structures as indicated FIGS. 11B and 11C for corners and junctions, respectively. [0177] In some embodiments, it is desirable to find a detuning pattern, D _n , such that the MIS-state of HVO remains the ground state of if model, even at finite W. In other words, it is desirable to find a detuning pattern that renders the MIS solution of the graph more energetically favorable (i.e., coupling it to the ground state of the system) than other solution of the graph. For the relevant qubit arrangements, the interactions of the H\ u and if model differ only around comers and junctions, in some embodiments. Thus, D _n can be set to D _n =

D everywhere (e.g., for all qubits) except at these structures, which can be considered individually and separately.

[0178] In some embodiments, with reference to corners as shown in FIG 11B, as discussed above in the section titled“Maximum Independent Sets for Unit Disk Graphs,” there exists at least one MIS-state, |i/> _G ), such that the comer qubit 1102 and its two neighbors 1106 are in one of the two ordered configurations: either the corner qubit 1102 is in state |0) and its two neighbors 1106 are in state |1) (|011)), or the corner qubit 1102 is in state |1) and its two neighbors 1106 are in state |0) (| 100)), where 1100) is a state of three vertices where the first vertex is in the state |1) , and the other two vertices are in the state |0). Detunings can be chosen for these three qubits in such a way that only these two configurations are relevant for the ground state of ffmodei. That is, detunings can be specified such that the energy of all possible qubit configurations can always be lowered by arranging the qubits on the comers in one of the two ordered states, such as, but not limited to 01010 or 10101. In some embodiments, this can be achieved if the detuning of the corner qubit 1102 is selected to be Ac = D + W+ 2e, and the detuning of its two neighbors 1106 is selected to be AN = D + W+ e. Here e > 0 can be chosen freely (up to the trivial constraint D _n < U, where U is the strength of interaction for two vertices that are very close). In some embodiments, there might be some advantages for choosing different e, and the values described herein are discussed as examples. In some embodiments, this choice of the detuning restores the energy difference between the two ordered configurations on the three qubits to D, corresponding to the additional vertex qubit in state |1). Without being bound by theory, up to a trivial constant (and contributions from junctions discussed below), //model and //UD are identical on all states where each comer is in one of the above ordered configurations. This includes in particular one ground state of UD, which corresponds to a MIS-state. All states that are not of this type have higher energy with respect to //model by construction, according to some embodiments.

[0179] In some embodiments, with reference to FIG. 11C, junctions can be treated with a similar detuning technique. Without being bound by theory, there exists at least one MIS such that the central qubit 1152 is in the state |0) and its three neighbors 1156 are in the state |1), or the central qubit 1152 is in the state |1) and its three neighbors 1156 are in the state |0). A detuning of the degree-3 junction qubit 1152 (denoted A/) can be chosen as A/ = D + W+ 3e, and the detuning of its three neighbors 1156 (denoted by Dn ) can be chosen as A v = A + W+ e. In some embodiments, this choice renders one of the two ordered states on the junction energetically more favorable than any other state and restores their energy difference to exactly 2D. It should be appreciated, based on the present disclosure, that other values of e can be selected.

[0180] In some embodiments, by using the detuning patterns described above, the actions of //UD and //model are, for non-limiting theoretical purposes, identical for at least one MIS- state. In addition, some embodiments ensure that the chosen detunings do not lower the energy of any other configurations. Therefore, a ground state of //model is a ground state of //UD, encoding an MIS problem on the corresponding unit disk graph such that the ground state is a solution to the MIS problem. NP-Completeness of the Rydberg Problem

[0181] Without being bound by theory, in some embodiments the detuning model described above can be applied to the case of the Rydberg Hamiltonian, thereby showing that it is NP-complete to decide whether the ground state energy of HR _yd is below a given threshold, when W _n = 0, where the atoms can be positioned arbitrarily in at least two dimensions. While the main idea is similar, the infinite ranged Rydberg interactions warrant more explanation.

[0182] As described above, it is possible to apply a detuning pattern that separates length and interaction scales, in some embodiments. This is possible in part because the Rydberg interactions decay fast, such that interactions between qubits that are close can be separated from the interactions between qubits that are far apart on the graph G. For example, the interactions between distant qubits can be neglected, if k is chosen large enough. In such embodiments, individual substructures of the system can be treated separately for purposes of theory.

[0183] In some embodiments, the low energy sector of the Rydberg Hamiltonian can be mapped to a much simpler effective spin model. For example, clusters of qubits can be addressed with specific detuning patterns such that only two configurations are relevant for each cluster. In some embodiments, the resulting, effective pseudo-spins are described by a MIS Hamiltonian. This makes it possible to encode Mis problems on planar graphs with maximum degree 3 such that the ground state of the Rydberg Hamiltonian is the Mis solution. The discussion in this section proves the details of the mapping to the effective model.

[0184] First, with regard to separation of length scales, the Rydberg interactions can decay as a power law with distance, and thus do not define a length scale where the interactions cease to exist. Nevertheless, in some embodiments, an exemplary qubit arrangement introduces two length scales: (1) the closest distance between two qubits, d , and (2) the grid length, g. These length scales are separated by a factor, g = d (2k + 1), that can be chosen arbitrarily in the transformation of the planar graph, Q, to the unit disk graph, G. In some embodiments, with such implementations, interactions between qubits that are“close” can be separated from interactions between qubits that are“distant”, for example when k is selected to be sufficiently large, as described in the present disclosure.

[0185] The closest qubits in the system are separated by a distance d and thus interact with a strength of EL = C/d ⁶ . This defines a convenient unit of energy, 1 L, which depends on the choice of the parameter k specifying the transformation from Cj to G.

[0186] Without being bound by theory, two qubits can be described as“distant” if their x or y coordinate differs by at least g (the grid length). The interaction energy of a single qubit with all distant qubits can be upper bounded by:

Equation 28:

[0187] In some embodiments, this bound can be derived by considering a system where a qubit can be placed in state |1) on each possible qubit position in the 2D plane (i.e. along all grid lines). Edist is the interaction energy of a single qubit in such a system with all other qubits that are distant in the above sense. This is an upper bound to the maximum interaction energy of a single qubit with all distant qubits on arbitrary graphs, that itself can be upper bounded by

Equation 29:

where z(c) is the Riemann zeta function. Since in some embodiments, dig = 1/(2 k + 1), the interaction energy of a spin with all distant spins can scale as Edist = 0 ( Wk ⁵ ). The total number of spins can scale as \ V\ ~ 0 (k\V\ , as the grid representation of the original planar graph Cj in the plane can be done with 0( \V\) area. Thus, in some embodiments, the contribution of interactions between all pairs of spins that are distant with respect to each other can be bounded by Edia\V\ ~ 0(lL\V\/k ⁴ ). In some embodiments, k can be chosen large enough such that these contributions can be neglected. In such embodiments, the ground state does not change due to the long range interactions. Note that the required k can scale polynomially with the problem size | |.

Effective Spin Model

[0188] As described in more detail below, without being bound by theory, an effective spin model can be developed in order to show that for each planar graph Cj = ( , £) of maximum degree 3, an arrangement of 0{k\V\) = 0( | | ^5/4 ) atoms can be identified, and detuning patterns applied such that the ground state of the corresponding Rydberg

Hamiltonian directly reveals a maximum independent set of Q. Furthermore, without being bound by theory, it can be seen that it is NP-complete to decide whether the ground state energy of H _Ryd (D. _v = 0) is lower than some threshold. In addition, without being bound by theory, it is possible to treat certain arrangements as pseudo-spins to simplify nonlimiting discussion of aspects of some embodiments.

[0189] In some embodiments, similar to the model discussed above, interactions beyond the unit disk radius (e.g., interaction tails) can be problematic. For example, longer range interactions can cause the ground state of systems with encoded MIS problems described above to differ from MIS solutions. These interactions can be more prominent close to vertices of degree 3, such as at vertex qubit 216 in FIG. 2 A or corners, such as vertex qubit 202 in FIG. 2A.

[0190] In some embodiments, as discussed throughout the present disclosure, to address these complications length and interaction scales can be separated to isolate these structures and control them individually. For example, for some graphs discussed in the present disclosure, various“special” vertices can be defined. These are special vertices can correspond to those that that either have integer grid coordinates or are irregular vertices. Special vertices thus include all vertices V of Q, but also a subset of vertices introduced to transform Cj to the UD graph G. FIG. 12A shows an example of a planar graph with maximum degree 3, Cj embedded on a grid, according to some embodiments. The planar graph of FIG. 12A includes vertex qubits 1202 A, 1202B that can be connected with edge qubits as discussed above.

[0191] FIG. 12B shows a spin arrangement implementing the graph of FIG. 12A with k = 7 and 0 = 1, according to some embodiments. For each special vertex, denoted as a set of vertices, Ai, can be defined that includes i and its 2q neighbors on each leg. The value of q can be chosen such that 2p < 2q < k!A.

[0192] FIG. 12B shows the regions At, corresponding to special vertices and their 2 q neighbors on each leg with a value of q = 1. The regions Ai, A ₄ , A ₆ and A9 can be referred to as corners, and regions A3, A7, and A11 can be referred to as junctions. The region A12 can be referred to as an open leg, while regions A2, As, As, A10, A13 and An can be referred to as special vertices on straight segments. Regions A10 and An can be referred to as irregular regions (e.g., those in which the spacing between vertices is different than in the other regions). In a spin model, each region A can be represented by a pseudo-spin A that interacts with neighboring pseudo-spins. In other words, the spins in a certain region can be aggregated as a pseudo-spin, which, as explained in more detail below, is limited to a specified number of states in the ground state. The two pseudo-spin states correspond to the two ordered configurations in A,. These two configurations, denoted by A = 0 and = 1, correspond to states where the special spin, is in state |S _j ) and all other spins in A, are in the corresponding perfectly ordered state. For each region At a pseudo-spin A can be defined and E _{a .} can be calculated for both pseudo-spin states (denoted E _Af ). Note that the number of spins in Ai that are in state 11) can be directly determined by the pseudo-spin state A, as it is given by A +m,q (where m denotes the degree of the special vertex i). As shown in FIG. 12B, sets Bi _j can be defined as the vertices along the path connecting A, and A _j . In such

embodiments, by construction Bij contains an even number of vertices. The ground states of different types of pseudo-spins A and B are described in more detail with reference to the section titled“Behavior at Low Energy.” For example, FIG. 12D shows the 3 possible configurations of vertices along a straight segment of ancillary spins that can appear in a ground state configuration. FIG. 12E shows the effective spin model corresponding to the graph Cj ^eff , according to some embodiments.

[0193] Without being bound by theory, due to the separation of length and interaction scales described throughout the present disclosure, interactions between two spins, u and v, are only relevant if they belong to the same set, u, v e At (or u,v e B _; y), or if they belong to adjacent sets, u e At and v e Bij. All other interaction terms will give contributions that vanish like 0(U\V\/k ⁴ ). Without being bound by theory, the total energy can be written as

Equation 30:

Where, E _c = gives the energy of the system when only spins inXare taken into account, and E _{X Y} =

quantifies the interaction energy between spins in two regions X and Y. The sum runs over i and j, such that A, is connected to Aj, by a non-empty set Bij, and the factor 1/2 compensates for double counting.

[0194] It is possible to apply an appropriate choice of the detunings, D _n , such that only a few configurations of spins in regions At and Btj are relevant to construct the ground state of the entire system, in some embodiments. In some such configurations it is possible to consider only a few configurations as candidates for the ground state, rather than

exponentially many. For example, as discussed in more detail below, in such embodiments only two configurations of spins in the regions At shown in FIG. 12C are relevant.

Furthermore, in such embodiments, only four configurations of spins in each region Bij are relevant, as shown in FIG. 12D, according to some embodiments. Moreover, these configurations can be completely determined by the pseudo-spin state & and .sy. If s, = 1 and .sy = 0 (or vice versa) the total energy is minimized if Bij is in the perfectly ordered state. If & = Sj = 0 or Si = j = 1, the lowest energy state is also ordered but with a single domain wall in the center of the region Bi _j. Without being bound by theory, the energy E _{B 1.} J . for these

S S

configurations can be denoted by by åJ ^{’ 1} , which yields:

1 _> J

Equation 31:

where note E _g = E _g The number of spins in state |1) in Bi _j can be written as b _J— .sv.sy

(where 2 bij gives the number of spins in Bij). The interaction energy E _{A .B . .} can be constant for all of the above combinations of relevant states [up to 0(11 / k ⁵ )].

[0195] In some embodiments, the total energy in the relevant configuration sector containing the ground state of i¾ _yd can be given by an effective spin model for the pseudo- spins:

Equation 32:

where the sum runs over neighboring pairs of pseudospins (without double-counting). Effective detunings A ^ff (e.g., detunings that have a similar effect on pseudospins as detunings for normal spins) can be introduced and pseudo-spin interactions U can be given by:

Equation 33: Af* = £¾ - ¾ - å (¾ - ¾) ,

jEM(i)

Equation 34: Equation 35:

[0196] In some embodiments, detuning patterns, D _n , can be chosen such that the effective pseudo-spin detunings and interactions are homogeneous, A ^ff = A ^eff and U = t/ ^eff , and satisfy 0 < A ^eff U ^eff . In some embodiments, when these parameters are chosen such that the above equations are satisfied the problem of the Rydberg Hamiltonian can be treated as equivalent to the MIS problem. To ensure that the long-range correction is bounded by some small constant n « A ^eff , k can be selected such that k ³ 0((U\V\/q) ^l/4) ) x can be computed based on the chosen detuning pattern.

[0197] This effective spin model in equation 32 can be related back to the original MIS problem on Q. Without being bound by theory, embodiments of this effective spin model corresponds exactly to a MIS problem on a graph ^eff obtained from = (¥, £) by replacing each edge {w,v} e E by a string of an even number 2 K _{u v} of vertices as discussed in more detail with reference to FIGS. 12A-12E. These additional vertices transport the independent set constraint giving an exact correspondence between the maximum independent sets of Q and ^eff , in complete analogy to the discussion above in the section titled“Maximum

Independent Sets for Unit Disk Graphs.”

[0198] Therefore, for each planar graph Q = (V, E) of maximum degree 3, an arrangement of 0{k\V\) = 0( |F| ^5/4 ) atoms can be identified, and detuning patterns applied such that the ground state of the corresponding Rydberg Hamiltonian directly reveals a maximum independent set of Q. Thus, without being bound by theory, it can be seen that it is NP- complete to decide whether the ground state energy of H _Ryd (D. _v = 0) is lower than some threshold. First, this question is in NP since the ground state of H _Ryd (D. _v = 0) has a classical description. In addition, for reduction from the NP-complete problem of deciding whether the size of MIS on Q, a planar graph of maximum degree 3, is > a. Let a' = a + the size of MIS of Q is > a, then the ground state energy of i¾ _yd is <

— a'A ^eff + x + h. If the size of MIS of Q is < a—

1, then the ground state energy of H _Ryd is >— ( a '— l)A ^eff + x— h. In some

embodiments, by choosing h < A ^eff /2, the MIS of Q has size > a if and only if the ground state energy of i¾ _yd is <— (a'— l/2)A ^eff + x.

Behavior at Low Energy

[0199] Without being bound by theory, embodiments of aspects of the ground state of the Rydberg Hamiltonian are described below. Example aspects of the Rydberg Hamiltonian show that at lower energy, the ground state corresponds to an MIS solution to an encoded problem. Thus, in example implementations, if a problem is properly encoded, transitioning the corresponding qubits into lower energy states can result in the system in the ground state which corresponds to the solution to the encoded problem. Without being bound by theory, discussion of such transitions can be aided by referencing the“pseudo-spins” discussed in the previous section.

[0200] In some embodiments, the ground state of i¾ _yd corresponds to a maximal independent set on the associated UD graph if the detuning of each spin is selected according to the techniques described herein. Maximal independent sets refer to independent sets where no vertex can be added without violating the independence property (e.g., that no two neighboring vertices can both be in the state |1)). The largest maximal independent set is the maximum independent set. For configurations corresponding to such sets, no two

neighboring spins can be in state |1) (independence), and no spin can be in state |0) if all its neighbors are in state |0) (maximality). [0201] In some embodiments of the Rydberg Hamiltonian, no two neighboring vertices on G can be simultaneously in state |1) if the system is in the ground state of i¾ _yd as long as the detuning D _n on all vertices obeys:

Equation 36:

where D is the maximal Euclidean distance between two vertices that are neighboring on G.

[0202] FIG. 13A shows a portion of a graph of qubits at a vertex of degree-3, according to some embodiments. The graph includes qubits 1302A, 1302B, 1304, 1306, and 1310.

FIG. 13B shows a similar graph with two vertices of degree-3, according to some

embodiments. FIG. 13B shows vertices of degree three 1314, 1324, as well as qubits 1312A, 1312B, 1316, 1318, 1320, 1322A, 1322B, 1326, 1328, 1330. While for the MIS Hamiltonian D > 0 is sufficient to energetically favor spins to be in the state |1) and thus maximality of the ground state, this is not necessarily the case for the Rydberg Hamiltonian. This is in part due to long range interactions discussed throughout the present disclosure. In this case, the energy gain D _n has to exceed the energy cost associated with the interaction energy of a spin v, with all other spins that are in state 11). In order to determine the required detuning that guarantees maximality, the interaction energy of a spin, v (such as qubit 1306 in FIG. 13A or qubits 1316 and 1326 in FIG. 13B) whose neighbors (on G) are all in state |0) can be bounded with the rest of the system. In order to obtain a bound, the energy to a contribution can be split from distant spins (bounded by E _dist = 0(U/k ⁵ ), see equation 29) and a contribution from the remaining“close” spins, B . The latter contribution from close spins is maximal if the spin v (1306, 1316, 1326) is directly neighboring to a vertex 1304 of degree 3. Since no two neighboring spins can be simultaneously in the state |1) as long as D _n < Amax, it can be upper bounded by B _{i)2 + 1 3} ⁼ 0.268031 X 'll. The first term corresponds to the maximum interaction of spin v with spins on the same grid line (such as qubits 1304, 1310 in FIG. 13 A, and qubits 1314, 1320 or 1324, 1330 in FIG. 13B), while the second term bounds the interaction with spins on the two perpendicular gridlines (such as 1302 A, 1302B in FIG. 13 A, or 1312A, 1312B, 1322A, 1322B in FIG. 13B).

[0203] The configuration shown in FIG. 13 A corresponds to a non-maximal independent set configuration where three qubits (1304, 1306) are in the same state |0), according to some embodiments. It is energetically favorable to flip the spin 1306 if the detuning is D _n > 0.268031 X T L. Q embedded on a grid. The configuration shown in FIG. 13B is another non- maximal independent set configuration with two domain walls 1318, 1328 along a straight segment. As discussed in more detail below, by merging the two domain walls (i.e., by flipping the spins along the straight segment such that only one domain wall exists), it is possible to obtain a defect-free configuration with one more spin in state |1). The increase in interaction energy between these two configurations can be bounded by 0.504748 X T L, where the state with one domain wall is the lowest energy configuration. Thus, for any given segment connecting two special vertices, such as vertices 1314, 1324, it is energetically favorable to include at most one domain wall.

[0204] Therefore, in some embodiments, the configuration of spins in the ground state of the Rydberg Hamiltonian corresponds to a maximal independent set on the associated UD graph (edge between nearest neighbor) if, C/D ⁶ > A _v > 0.268031 x C/d ⁶ .

Straight Segments

[0205] FIG. 13B shows an exemplary straight segment pseudo-spin, according to some embodiments. As explained below, such segments can be found in at most three states, shown in FIG. 12D.

[0206] The spins on each edge £ of the graph Cj are in an ordered configuration, such as a configuration where spins are alternating in state |0) and | 1), up to so called domain walls, where two (but not more) neighboring spins are in state |0). A MIS configuration on G has at most one such domain wall on any array of spins connecting two special vertices i and j such as 1314, 1324 in FIG. 13B, according to some embodiments. In order to ensure that the ground state of the Rydberg system respects this property, the detuning can be chosen to be large enough. For example, assuming such a straight segment is in an ordered state with n domain walls, when two such domain walls (1318, 1328) are combined an additional spin, v, can be flipped from state |0) to state | 1), lowering the single-particle contribution to the total energy by D _n. The corresponding difference in interaction energy grows with the distance of the domain walls, for example, it is maximal if the two domain walls are initially on opposite ends of a segment connecting two junctions as shown in FIG. 13B. It can be bounded by Ei + f iist, where

Equat ⁱ on 37:

£ 0.490084.

[0207] This bound can be understood with reference to FIG. 13B, where the zero- domain-wall configuration can be obtained from one with two domain walls by shifting the ordered spins on the middle bar 1340 (i.e., between the domain walls) one atom to the left in, and then flipping a spin next to the right junction to the state |1), according to some embodiments. Thus, the first sum corresponds to the interactions between the middle bar of the H-shape with the sides, and the last sum is the extra interaction from a new flipped spin. Thus (for large enough k where £dist is negligible) it is favorable to merge two neighboring domain walls along a straight segment if

Equation 38: A _v ³ A _min º 0.490084 X 11.

[0208] Therefore, for Amm < Dn < Amax, the ground state of the Rydberg Hamiltonian is ordered on all segments connecting two special vertices, with at most one domain wall per segment. This justifies the assumption made with reference to FIG. 12D, that it is possible to restrict the analysis of the low energy sector to only 3 configurations of the structures Bij. In the first configuration, all spins alternate with different states on either end of the structure Bij. In the second configuration, the spin states on either end of the structure Bij are equal to 0, with one domain wall in between. In the third configuration, the spin states on either end of the structure B are equal to 1, with one domain wall in between. The detuning for all those spins can be set to be equal and in the above range, which can be denoted by As.

[0209] Below, embodiments of special structures L consisting of a vertex i and its 2 q neighbors on each leg are described in more detail. The discussion can be restricted to states that are ordered in At with at most one domain wall per leg. Detuning patterns can be selected for the spins in At such that any such domain wall is energetically pushed away from spin, out of the structure A,. Without being bound by theory, with such a detuning only the two ordered states on At have to be considered for the ground state.

Corners

[0210] A comer vertex and the spins on the two legs, which can be collectively treated as a pseudo-spin, are considered first, such as those shown in embodiments of FIG. 12C. A detuning pattern consistent with Amax > A _v > Amin can be applied such that only two configurations are relevant in the low energy sector (e.g., for the ground state). These are the two ordered states across the comer which correspond to the state of the corner spin, |0) or 11), which can be denoted by 5 = 0 and 5 = 1 respectively.

[0211] In some embodiments, if the corner spin is in state |0), a domain wall can be formed where on one of the legs the two spins at distance 2p - 2 and 2p - 1 from the corner are in state |0). The domain wall can be moved by one unit (i.e. p p +1) by flipping the state of the spins 2p - \ and 2 p. The interaction energy increases in this process. Its amount can be bounded by T LF _Q (p) + 0(U/k ⁵ ) with [0212] Without being bound by theory, this bound can be understood as follows.

Effectively, the one spin excitation is moved by one site towards the corner. While the interaction energy of this excitation with all spins on the same leg is reduced (such as where 2 P < k ), and thus is upper bounded by zero, the interaction energy with all spins on the other leg increases. Since any defect on this leg would decrease this interaction energy, the energy can be maximized if all spins on the other leg are in the perfectly ordered state, which gives the bound in equation 39.

[0213] While the interaction energy can increase in this process, the contribution from the single particle term to the total energy changes by Aip-i - Ai _P . Depending on the choice of the detuning, this can lead to an energy gain such that it is energetically favorable to move the domain wall by one unit away from the comer spin. Thus, the energy is minimized if the first 2 q spins on each of the legs are in the perfectly ordered state if the corner spin is in state |0) and the detuning for spins satisfy for/ _? = 1, 2 . . ., q:

Equation 40:

[0214] Similarly, if the comer spin is in state 11), the energy is minimized if the first 2 q +1 spins on each of the leg are in the perfectly ordered state if detunings satisfy:

Equation 41: with

Equation 42:

[0215] The conditions in equations 40, 41 can be satisfied by the choice (for p = 1, . . q):

Equation 43:

Equation 44:

[0216] In some embodiments, these sums are convergent and can be efficiently evaluated numerically. D _U is monotonically decreasing with x and for large x approaches D _¥ +

0(U / x ⁵ ). Moreover, the maximum value of D _U in this sequence can be evaluated to

D ₁ = D _¥ + 0.134682 x Ίί. It is therefore possible to choose D such that the detuning along the legs of the comer are within the required range.

[0217] With the above detuning pattern and the choice D > AB , the ground state configuration is necessarily such that the 4q + 1 spins forming a corner structure are in one of the two ordered states. These states can be labeled as 5 = 0 and 5 = 1 according to the corresponding state of the comer spin (|0) and |1)).

[0218] Note that the choice in equations43, 44 fixes the detuning on all spins except the detuning of the corner spin, denoted Ac. This makes it possible to tune the relative energy between the two relevant configurations, E .— E . for the ground state. Without being bound by theory, to calculate this difference, the quantity Ic can be defined as the difference in interaction energies between the two spin configurations,

[0219] Similarly, without being bound by theory, the difference in energy ( D _c ) of the two configurations due to the single particle term can be calculated as:

Equation 46: 0.237094 + 0(ZL/q ⁵ )

[0220] For the corner structures, f¾.— E . =—D _c — I _c , this can be evaluated to Equation 47: E — E . =— A _c + 0.143797 X Ί1 + 0(ZL/q ⁵ ) which can be fully tuned by the detuning of the corner spin. Thus, in some embodiments, corner structures can be treated as having only two configurations.

Junctions

[0221] In some embodiments, junctions, such as that shown in FIG. 13 A can be treated in a similar way as corners. For example, the spin at the junction and adjacent spins can be treated as a pseudo-spin. A detuning pattern on the legs of a junction can be chosen such that it is energetically favorable to push domain walls away from the junction, such that in the ground state, a junction can only be in one of the two ordered configurations, according to some embodiments. In this way, junctions also can be described by at most two states defined by the state of the 3-way junction vertex.

[0222] First, with reference to central spin, such as the degree-3 vertex 1304 in FIG. 13 A, which can be in state |0), the three legs can be referred to by X Y and Z, where for concreteness X and Z are on the same gridline (e.g., corresponding to spins continuing from vertices 1302A, 1302B) with F being perpendicular thereto (e.g., corresponding to spins continuing from vertex 1306), according to some embodiments. In a situation where two vertices on leg X at a distance 2p - 2 and 2p - 1 from the central vertex are in state |0) and thus form a domain wall, in order to push this domain wall one unit away from the vertex, the state of vertices 2p - \ and 2 p has to be flipped. The interaction energy required to do so is bounded, in some embodiments, by T LF ₀ (p). Note that this bound can hold true regardless of whether (or where) on legs Y and Z a domain wall is present. Similarly, the interaction energy required to push a domain wall on leg 7 by one unit from the corner can be bounded by 2 F ₀ (pl

[0223] In the other case, i.e. if the central vertex is in state |1), similar bounds can be found. If the two vertices on leg X at a distance 2p - l and 2 p from the central vertex are in state |0), the interaction energy will be increased when the states of vertices 2 p and 2 p +1 are flipped, i.e. when the domain wall is pushed away from the center by one unit. The interaction energy required to push the domain wall on leg A" by one unit can therefore be bounded by 'lLF ₁ (p). This bound can be valid independently of the configuration of the spins on legs Y and Z. The interaction energy required to push a domain wall on leg Y by one unit from the corner can be bounded by 21 LF^p).

[0224] Therefore, in some embodiments, the state that minimizes the energy may not contain any domain wall on the first 2q + 1 spins on each leg if the detuning pattern satisfies

for p = 1, 2, . . . q, and D^= D® . Here D^ denotes the detuning of the v-th spin on leg s. Without being bound by theory, this can be achieved by the following:

Where D^= D® . For the choice above, the maximum value of in this sequence can

( y J

evaluate to D^ = D _¥ + 0.269364 x 'll. Thus, all the detunings on the legs of a junction At are within the range [D _¥ , D _¥ + 0.269364 X Ti] It is therefore possible to choose D such that all of them are in the allowed range, Amin < D _n < Amax. With the additional choice D > AB, it can be seen that, in some embodiments, only the two ordered configurations are relevant for the ground state.

[0225] Analogous to the situation of corner structures, the detuning of the junction spin, denoted A / can be used to tune the relative energy between the two relevant configurations of the junction. Without being bound by theory, the difference in interaction energies of the two spin configurations in this structure can be given by

which can be evaluated to I _j = 0.218387 X Ί1 + 0(U/q ⁵ ). Similarly, the difference in energy of the two configurations due to the single particle terms can be calculated as

Equation 57:

which can be evaluated to D _j = Ay— 0.474188 X ZL + Q(ZL/q ⁵ ). Without being bound by theory, the quantity — E {. =—D _j — l _j can therefore be written as

Equation 58: 0.255801 X V. + 0(11/ q ⁵ ).

Other special vertices

[0226] Without being bound by theory, some embodiments include other special vertices described above in addition to vertices at corners and junctions. These are vertices of degree 1 (open ends, such as An in FIG. 12B), vertices on a grid point with two legs on the same grid line (straight structures, such as An in FIG. 12B), and irregular vertices (such as AM in FIG. 12B). These special vertices can be treated as pseudo-spins, as described in more detail below. The analysis in the previous two sections can be repeated for all such vertices, thereby showing that these two can be described as being in one of two states.

[0227] Open ends. In some embodiments, energy can be lowered if a domain wall is moved away from a spin at the end of an open leg if the detuning is constant for the spins on that leg. Therefore, the pseudo-spin states can be restricted to the two ordered configurations on the 2 q spins adjacent to the spin at the end of an open leg. Without being bound by theory, the energy difference between the two states of such an open leg can be denoted by:

Equation 59: 0.015896 X ΊI + 0(11/ q ⁵ )

where Do is the detuning for the spin corresponding to the vertex of degree 1. The homogenous detuning on the 2 q spins adjacent to the latter can be chosen to be equal to D .

[0228] Straight structures. In some embodiments, for a regular special vertex, the relevant configurations for the ground state can be restricted to two ordered states by choosing a detuning for all 4 q spins on the leg as D > AB. With this choice it would be energetically favorable to move a potential domain wall into the adjacent neighboring regions. Without being bound by theory, the energy difference can be denoted by:

Equation 60: 0.015896 ΊI + 0(11/ q ⁵ )

where As denotes the detuning of the special vertex.

[0229] Irregular structures. In some embodiments, irregular vertices can be treated identically to straight structures. Since the spacing of the spins close to the irregular vertex is slightly larger than elsewhere (e.g., because irregular structures can be defined as structures where the spacing between vertices is larger than in ordinary structures to ensure that the number of ancillary vertices on each edge is even), any domain wall will be pushed away from the irregular structure naturally, if the detuning in the irregular structure is larger than AB. Thus, only the two ordered configurations are relevant for the ground state. The corresponding energy difference can be numerically evaluated for every choice of f (which can indicate how ancillary vertices are positioned). In the large f (and q) limit, it is possible to obtain the same analytic expression as in equation 60.

Effective energy

[0230] In some embodiments, it is possible to determine the effective detuning A of the pseudo-spins (such as the straight segments, corners, junctions, and other special vertices described above that can be described as having a limited set of states in the ground state) and their effective interaction energies Uf^ . In some embodiments, knowing the effective detunings and effective interactions of pseudospins allows for encoding of the problem to be effectively solved.

Effective interactions

[0231] In some embodiments and without being bound by theory, the discussion above with reference to straight segments, corners, junctions, and other special vertices can rely on the concept that only four spin configurations are relevant in each region Bij to describe the ground state. These correspond to the four possible configurations of the spins in the adjacent regions At and Aj. For example, if the detuning in Bij is set to be homogeneous, AB, if i = 1 and s, = 0 (or vice versa), the lowest energy configuration should correspond to the perfectly ordered state on Bij (with bij spins in state 11), with b vertices in the set B ), if Amax > AB > Amin. Any other configuration would require at least two domain walls. Second, if Si =

Sj =0 then energy can be lowered by arranging the spins in Bi _j in an ordered configuration with bij spins in state |1), and one domain wall. The position of this domain wall does not change the energy up to 0(11 / k ⁵ ), such that there is no need to distinguish between the different domain wall configurations. Finally, if Si = Sj = 1 the lowest energy configuration is similarly achieved by an ordered configuration with one domain wall, if Amax > As > Amin. While the position of the domain wall can be irrelevant, only bij - 1 spins are in state 11) in this state.

[0232] Without being bound by theory, in some embodiments, the relevant energy differences between these different relevant configurations can be readily calculated, as

[0233] If Ai and Aj are not connected, this can be zero. Else, this term becomes independent of i and j in the large q and k limit, since = 0(k ), and evaluates to

Equation 62: ¾'° - £°'° = 0.0146637 X i + 0{U/k ^s )

[0234] Similarly, E _g ¹ — E _B ° can be calculated as

[0235] Thus, in some embodiments the effective interaction between pseudo-spins U = i/ _eff = E ¹ + E _b ° can be written as

Equation 64: t/ ^eff = A _B - 0.0134313 x U + o(u/k ^s ).

[0236] In such embodiments, the effective interaction depends only on the choice of the detuning in the connecting structures, As, thereby allowing for control of the effective interactions for specifying problems.

Effective detuning [0237] In some embodiments, the effective detuning for a pseudo spin can be given by where , is the degree of the special vertex, i (i.e. , = 2 for corner vertices and straight structures, rm = 3 for junctions, and m, = 1 for open legs). The value of E — E . can be fully tuned by the choice of the detuning of spin i. This can be chosen such that a homogeneous effective detuning, A ^ff = A ^eff , can be used for all four types of pseudo-spins. This can be achieved by

Equation 65: _c = A ^eff + 0.173124 X 'll,

Equation 66: _j = A ^eff + 0.299792 X 'll,

Equation 67: D ₀ = A ^eff + 0.0305597 X 'll,

Equation 68: D ₅ = A ^eff + 0.0452234 X 'll.

[0238] This is compatible with Amax > Dn > Amm and the realization of an effective spin model with 0 < A ^eff < l / ^ciT . In some embodiments, where this inequality is satisfied, the solution of the effective model can coincide with the solution of the hard problem to be solved.

Example QAOA Implementations

[0239] In some embodiments, additional considerations are relevant to implementing the QAOA techniques described in the present disclosure. The framework of QAOA is general, however, and can be applied to various technical platforms to solve combinatorial optimization problems.

Finite measurement samples

[0240] In some embodiments, quantum fluctuations (such as projection noise) can result in finite precision since the precision can be obtained via averaging over finitely many measurement outcomes that can only take on discrete values. Hence, there can be a trade-off between measurement cost and optimization quality: finding a good optimum can require good precision at the cost of a large number of measurements. Additionally, large variance in the objective function value can demand more measurements but may help improve the chances of finding near-optimal MaxCut configurations.

[0241] Without being bound by theory, as an example, the effect of measurement projection noise can be demonstrated with a full Monte-Carlo simulation of QAOA on some example graphs, where an objective function is evaluated by repeated projective

measurements until its error is below a threshold. Exemplary implementation details of this numerical simulation are discussed in the section titled“Exemplary simulation with measurement projection noise.”

[0242] FIG. 19 shows Monte-Carlo simulation of QAOA accounting for measurement projection noise, on the example instance studied in FIG. 18B, according to some

embodiments. Simulated QAOA is optimized with embodiments of the disclosed FOEIRIER heuristic starting at level either p = 1 or p = 5 using an educated choice of initial parameters based on known optima found for small size instances. The approximation ratio n = |Cuh| / |MaxCut| is tracked from the z-th measurement, and the minimum fractional error 1 - n found after M measurements is plotted, averaged over many Monte-Carlo realizations. The solid and dashed lines correspond to QAOA optimized with the FOEIRIER strategy starting with an educated guess of (y,/?) atp = 1 and p = 5, respectively, and increasing p by 1 after a local optimum is found. A local optimum can be presumed to be found if the new candidate parameter vector norm changes by less than 0.01, or when the average Cut changes by less than 0.1. The dash-dot line corresponds to QAOA optimized starting with randomly guessed (y, /?) at p = 1 , and increasing p by 1 after a local optimum is found. The three labelled dotted lines are results from repeatedly running QA with total time T= 10, 10 ² , 10 ³ , and performing measurements on the output state. [0243] As shown in FIG. 19, QAOA can find the MaxCut solution in -10-10 ² measurements, and starting embodiments of the disclosed optimization at intermediate level (p = 5) is better than starting at the lowest level (p = 1). In comparison, some embodiments with random choices of initial parameters starting at/> = 1 perform worse, which fail to find the MaxCut solution until 10 ³ -10 ⁴ measurements have been made. Moreover, if QAOA is compared to QA with various annealing time, it appears that the choice of annealing time T = 100 can perform just as well as QAOA on this instance. Nevertheless, running QAOA at level p = 5 is still more advantageous than QA at T= 100, in some embodiments, when coherence time is limited.

[0244] While this exemplary simulation is limited to small-size instances, the good performance of QAOA and QA can be complicated by the small but significant ground state population from generic annealing schedules. Since it can take 10 ² measurements to obtain a sufficiently precise estimate of the objective function, a ground state probability of ³ 10 ^-2 would mean that one can find the ground state without much parameter optimization. In some embodiments, for larger problem sizes with 10 ² - 10 ³ qubits, the ground state probability from generic QAOA/QA protocol without optimization can decrease

exponentially with system size, whereas the measurement cost of optimization grows merely polynomially with the problem size. The results here indicate that the parameter pattern and the disclosed heuristic strategies are practically useful guidelines in realistic implementation of QAOA that leverages optimization to significantly increase the probability of finding the ground state.

Implementation for large problem sizes

[0245] Implementing a solution to the MaxCut problems with quantum machines can be limited by quantum coherence time and graph connectivity. In some embodiments, in terms of coherence time, QAOA is highly advantageous: the hybrid nature of QAOA as well as its short- and intermediate-depth circuit parametrization makes it useful for quantum devices. In addition, QAOA is not generally limited by the small spectral gaps, which demonstrates that interesting problems can be solved (or at least approximated) within the coherence time.

Performance of Heuristically Optimized QAOA

Comparison between heuristics embodiments and brute-force

[0246] According to some embodiments, non-limiting comparisons of exemplary heuristic strategies to exemplary brute-force approaches can be evaluated. For example, FIG. 16A shows a comparison between an exemplary FOURIER heuristic and the brute-force (BF) approach for optimizing QAOA, on an example instance of 16-vertex w4R graph, according to some embodiments. The value of 1— r, where r is the approximation ratio, is plotted as a function of QAOA level p on a log-linear scale. The brute-force points are obtained by optimizing from 1000 randomly generated starting points, averaged over 10 such realizations. At low p, the exemplary FOURIER heuristic strategies perform just as well as the best out of 1000 brute-force runs— both are able to find the global optimum. The average performance of the brute-force strategy, on the other hand, is much worse than the disclosed heuristics. This indicates that the QAOA parameter landscape can be highly non-convex and filled with low-quality, non-degenerate local optima. When p ³ 5, the exemplary FOURIER heuristic outperforms the best brute-force run.

[0247] To estimate the number of brute-force runs needed to find an optimum with the same or better approximation ratio as the exemplary FOURIER heuristics, brute-force optimization was performed on 40 instances of 16-vertex u3R and w3R graphs with up to 40000 runs, as shown in FIG. 16B, according to some embodiments. In particular, FIG. 16B shows the median number of brute-force runs needed to obtain an approximation ratio as good as an exemplary FOURIER heuristic, for 40 instances of l6-vertex u3R and w3R graphs. A log-linear scale is used, and exponential curves are fitted to the results. Error bars are 5th and 95th percentiles. In FIG. 16B, the median number of brute-force runs needed to beat the exemplary FOURIER heuristic is shown to scale exponentially with p. Hence, the embodiments of the disclosed heuristics offer a dramatic improvement in the resource required to find good parameters for QAOA. As verified with an excessive number of brute- force runs, the heuristics often find the global optima.

[0248] Although most examples discussed in the present disclosure use gradient-based methods (BFGS) in numerical simulations, non-gradient based approaches, such as the Nelder-Mead method, can be used with the disclosed heuristic strategies. The choice to use gradient-based optimization can be motivated by the simulation speed, which in some implementations is faster with gradient-based optimization. In other embodiments, other procedures can be used.

Performance of QAOA on typical instances

[0249] Using embodiments of the disclosed heuristic optimization strategies in hand, the performance of intermediate /i-level QAOA on many graph instances can be examined. For example, it is possible to consider many randomly generated instances u3R and w3R graphs with vertex number 8 < N< 22 and use embodiments of the disclosed FOURIER strategy to find optimal QAOA parameters for level p < 20. In the following discussion, the fractional error 1— is used to assess the performance of QAOA.

[0250] In one example, the case of unweighted graphs, specifically u3R graphs can be considered. For example, FIGS. 16C and 16D plot the fractional error 1 - r as a function of QAOA’s level p. The exemplary results shown in FIGS. 16C and 16D were obtained by applying embodiments of the disclosed heuristic FOURIER optimization strategies to up to 100 random instances of u3R graphs (FIG. 16C) and w3R graphs (FIG. 16D). Lines with different shapes correspond to fitted lines to the average for different system size A, where the model function is 1 r oc e ^P//p ° for unweighted graphs and 1— r oc e ^ ^P//p ° for weighted graphs. Insets show the dependence of the fit parameters po on the system size N.

[0251] As shown in FIG. 16C, on average, 1— r oc e ^p / ^p ° appears to decay

exponentially with p in some embodiments. Note that since the instances studied were u3R graphs with system size N < 22, where Cmax < A| < 33, such embodiments prepared the MaxCut state in virtually all scenarios where the fractional error goes below ~l0 ^-2 . Without being bound by theory, this good performance can be understood by interpreting QAOA as Trotterized quantum annealing, especially when the optimized parameters are of the pattern seen in FIGS. 15C-15H, where initialization is in the ground state of—HB and the system then evolves with HB (and He) with smoothly decreasing (and increasing) durations. The equivalent total annealing time T can be approximately proportional to the level p, since the individual parameters gi, bi = 0(1) and can correspond to the evolution times under He and HB. If Tis much longer than 1/D^ _ίh , where Amin is the minimum spectral gap (which can govern the application of quantum annealing algorithm (QAA), for example, by requiring a long time T to run where the spectral gap is small), quantum annealing can find the exact solution to MaxCut (ground state oi—Hc) by the adiabatic theorem, and achieve

exponentially small fractional error as predicted by the Landau-Zener formula. Numerically, the minimum gaps of these u3R instances when running quantum annealing can be determined to be on the order of A _min > 0.2 in some embodiments. It is thus not surprising that QAOA can achieve exponentially small fractional error on average of exemplary embodiments, since it can prepare the ground state of— He through the adiabatic mechanism for these large-gap instances. Nevertheless, this exponential behavior can break down for some hard instances, where the gap is small.

[0252] As shown in FIG. 16D, in the case of w3R graphs, the fractional error appears to scale as to some embodiments. In some embodiments, the stretched-exponential scaling is true in the average sense, while individual instances have very different behavior. For easy instances whose corresponding minimum gaps Amin are large, exponential scaling of the fractional error can be found. For more difficult instances whose minimum gaps are small, fractional errors reach plateaus at intermediate p, before decreasing further when p is increased. These hard instances are discussed in more detail in the section below titled“Adiabatic Mechanisms, Quantum Annealing, and QAOA.” Notably, when averaged over randomly generated instances, the fractional error is fitted remarkably well by a stretched-exponential function.

[0253] These results of average performance of embodiments of QAOA are notable despite considerations of finite-size effect. While the decay constant po does appear to depend on the system size N as shown in the insets of FIG. 16C, 16D, the exemplary finite- size simulations shown in the present disclosure may not conclusively determine the exact scaling. In some embodiments, if the trend continues to system size of N= 100 ~ 1000, then QAOA will be practically useful in solving interesting MaxCut instances, better or on par with other known algorithms, in a regime where finding the exact solution will be infeasible (i.e., it will be an effective quantum approximator). Even if QAOA fails for some worst-case graphs, it can still be practically useful, if it performs well on a randomly chosen graph of large system size.

Comparison between heuristics

[0254] The difference in the performance among embodiments of the various heuristic strategies proposed in the present disclosure can be examined on an example instance of 14- vertex w3R graph. As shown in FIG. 23, exemplary embodiments of INTERP and

FOURIER[oo, 0] strategies disclosed herein have essentially identical performance (except for small variations barely visible at large p ), according to some embodiments. This can be explained by considering that both strategies generate starting points for optimizing level p +1 based on smooth deformation of the optimum at level p , in some embodiments.

Furthermore, the FOURIERfoo, 10] strategy outperforms INTERP and FOETRIERfoo, 0] beginning at level p ³ 20. Notably, even when restricting the number of QAOA parameters to 2 q = 10, the FOURIER[c/=5,/i= 10] strategy closely matches the performance of other heuristics at low p and beats the R = 0 heuristic at large p. This suggests that the optimal QAOA parameters may in fact exist in a small-dimensional manifold. Therefore, optimization for intermediate p-level QAOA can be dramatically simplified by using the disclosed parameterization (u, v).

Adiabatic mechanism, quantum annealing, and QAOA

[0255] The previous section discussed the performance of embodiments of QAOA for MaxCut on random graph instances in terms of the approximation ratio r. Although useful for approximate optimization, QAOA is often able to find the MaxCut configuration— the global optimum of the problem— with a high probability as level p increases. In this section, the efficiencies of example embodiments of the disclosed algorithms are shown to find the MaxCut configuration and compare it with quantum annealing. In some embodiments, QAOA is not necessarily limited by the minimum gap in the quantum annealing and explain a mechanism at work that allows it to overcome the adiabatic limitations.

Comparing QAOA with quantum annealing

[0256] A predecessor of QAOA, quantum annealing (QA) can be used for solving combinatorial optimization problems. Without being bound by theory, to find the MaxCut configuration that maximizes ( H _c ), the following simple QA protocol can be considered:

Equation 69: H _QA (s ) =—[sH _c + (1— s)H _B ], s = t/T, where / e [0, 7] and / ^' is the total annealing time. The initial state can be prepared to be the ground state of HQA(S = 0), i.e., \y(0)) = |+) ^®w . The ground state of the final Hamiltonian, HQA(S = 1), can correspond to the solution of the MaxCut problem encoded in He. In adiabatic QA, the algorithm can rely on the adiabatic theorem to remain in the instantaneous ground state along the annealing path and solves the computational problem by finding the ground state at the end. To ensure a high probability of success, the ran time of the algorithm typically scales as T= 0(1/D^ _ίh ), where Amin is the minimum spectral gap. Consequently, adiabatic QA becomes inefficient for instances where extremely small. These graph

instances can be referred to as hard instances for adiabatic QA.

[0257] In some embodiments beyond the completely adiabatic regime, there is often a tradeoff between the success probability (ground state population pcsil)) and the run time (annealing time 7): either algorithm can be run with a long annealing time to obtain a high success probability, or it can be run multiple times at a shorter time to find the solution at least once. One metric that can be used to determine the best balance can be referred to as the time-to-solution (TTS):

Equation 70:

[0258] TTSQA(T) can measure the time required to find the ground state at least once with the target probability pd (taken to be 99% in the present disclosure), neglecting non- algorithmic overheads such as state-preparation and measurement time. The adiabatic regime where ln[l— per Landau-Zener formula can yield TTSQ _A OC 1/D^ _ίh which is independent of T. In some cases, it can be better to run QA non-adiabatically to obtain a shorter TTS. By choosing the best annealing time T, regardless of adiabaticity, TTS _QA can be determined as the minimum algorithmic run time of QA. Without being bound by theory, a similar non-limiting metric can be defined for QAOA for purposes of benchmarking. The variational parameters g I and b[ can be regarded as the time evolved under the Hamiltonians He and HB, respectively. The sum of the variational parameters can be interpreted to be the total“annealing” time such that T _p = l + \b[ I), and

TTS _QA0A (p) and TTS _QA0A can be defined as:

Equation 72:

Equation 73: TTS _QA0A = min TTS _QA0A (P). where pcs(p) is the ground state population after the optimal p-level QAOA protocol, in some embodiments. This quantity need not take into account the overhead in finding the optimal parameters. TTSQAOA(P) can be used to benchmark the algorithm but should not be taken directly to be the actual experimental ran time.

[0259] To compare examples of the algorithms, TTS^and TTS _QA0A can be computed for many random graph instances. For each even vertex number from N = 10 to N = 18, 1000 instances of w3R graphs are generated. FIGS. 17A-17E are plots of the relationship between TTS _QA0A and the minimum gap Amin and vertex size N in quantum annealing for each exemplary instance, according to some embodiments. The maximum cutoff p is taken to be /½ax = 50, 50, 40, 35, 30 for N= 10, 12, 14, 16, 18, respectively, for FIGS. 17A-17E. The dashed line corresponds to adiabatic QA run time of 1/A^ _in predicted by Landau-Zener. Most of the exemplary random graphs have large gaps (A _min > 10 ^-2 ). The optimal TTS can be observed to follow the Landau-Zener prediction of 1/A^ _in min for these graphs. This indicates that a quasi-adiabatic parametrization of QAOA can be the best when

reasonably large, in some embodiments. Many exemplary graphs, however, exhibit very small gaps (A _min < 10 ^-3 ), and thus require exceedingly long run time for adiabatic evolution. For some graphs, Amin is as small as 10 ^-8 , which can imply that an adiabatic evolution requires a run time T > 10 ¹⁶ . Nevertheless, QAOA can find the solution for these hard instances faster than adiabatic QA. Notably, TTS _QA0A appears to be independent of the gap for many graphs that have extremely small gaps and beats the adiabatic TTS (Landau- Zener line) by many orders of magnitude, in some embodiments. Thus, an exponential improvement of TTS is possible with non-adiabatic mechanisms when adiabatic QA is limited by exponentially small gaps.

[0260] Similarly, for QA, the optimal annealing time T need not be in the adiabatic limit for small-gap graphs. FIGS. 17F-17J show TTS _QA0A versus TTS _QA for each random graph instance, according to some embodiments. Embodiments of QAOA outperform QA for instances below the dashed line. The (Pearson’s) correlation coefficient between QAOA TTS and QA TTS is p (ln 0.91. As shown in FIGS. 17F-17J , there can be a strong correlation between TTS _QA0A and TTS _QA for each graph instance. Without being bound by theory, this suggests that QAOA is making use of the optimal annealing schedule, regardless of whether a slow adiabatic evolution or a fast diabatic evolution is better. Without being bound by theory, the following section titled“Beyond the adiabatic mechanism: a case study” explores a non-limiting example to explain aspects of the results observed in FIGS. 17A-17H and a mechanism of QAOA that go beyond the adiabatic paradigm.

Beyond the adiabatic mechanism: a case study

[0261] To understand the behavior of QAOA, graph instances that are hard for adiabatic QA can be addressed in more detail. For example, a representative instance is used to explain how embodiments of QAOA as well as diabatic QA can overcome the adiabatic limitations. As illustrated in 18 A, QAOA can learn to utilize diabatic transitions at anti-crossings to circumvent difficulties caused by very small gaps in QA. FIG. 18A shows a schematic of how QAOA and the interpolated annealing path can overcome the small minimum gap limitations via diabatic transitions (small dashed line). Naive adiabatic quantum annealing path leads to excited states passing through the anti-crossing point (large dashed line).

[0262] FIG. 18B shows a representative instance of weighted 3-regular graph with N =

14, which has a small minimum spectral gap Amin < 10 ³ , according to some embodiments. For example, FIG. 18B shows an instance of weighted 3-regular graph that has a small minimum spectral gap along the quantum annealing path given by equation 69. As an example, the optimal MaxCut configuration is shown with two vertex types (circles and squares), and the solid (dashed) lines are the cut (uncut) edges. For this example, hard instance, the quantum annealing process can be simulated. Without being bound by theory, FIG. 18C shows populations in the ground state and low excited states at the end of the process for different annealing time T, according to some embodiments. The solid line follows the Landau-Zener formula for the ground state population, pcs = 1— exp

where c is a fitting parameter. Since the minimum gap can be very small, the adiabatic condition requires T > 1 / D^ _ίh ~ 10 ⁶ . The Landau-Zener formula for the ground state population pcs = 1— exp fits well with the exact numerical simulation discussed herein, where c is a fitting parameter. However, there is a“bump” in the ground state population at annealing time T » 40. At such a time scale, the dynamics can be considered non-adiabatic; which can be referred to as a“diabatic bump.” This phenomenon has been observed in quantum annealing of other optimization problems as well.

[0263] In some embodiments, it is also possible to simulate QAOA on this hard instance. As mentioned above, although QAOA can optimizes energy instead of ground state overlap, substantial ground state population can still be obtained even for many hard graphs. Using an exemplary embodiment of the disclosed FOURIER heuristic strategy, various low-energy state populations of the output state are shown for different levels p shown in FIG. 18D, according to some embodiments. More particularly, FIG. 18D shows populations in low excited states using QAOA at different level p. The exemplary FOURIER heuristic strategy is used in the optimization shown in FIG. 18D. In such embodiments, QAOA can achieve similar ground state population as the diabatic bump at small p, with substantial

enhancements occurring after p ³ 24.

[0264] Without being bound by theory, to better understand the mechanism of

embodiments of QAOA and make a comparison with QA, the QAOA parameters can be interpreted as a smooth annealing path. The sum of the variational parameters can be interpreted to be the total annealing time, i.e., T _P = as discussed above. A smooth annealing path can be constructed from QAOA optimal parameters as

Equation 74:

Equation 75:

where A can be chosen to be at the mid-point of each time interval (g[ + b[). With the boundary conditions f(t = 0) = 0, f(t = T _p ) = 1 and linear interpolation at other intermediate time /, QAOA parameters can be converted to a well-defined annealing path. This conversion can be applied to the QAOA optimal parameters at p = 40, as shown in FIG. 18E, according to some embodiments. The annealing time parameter can be written as s = t / T _p , where T _p = \) _· The flat upper dotted line labels the location of anti crossing where the gap is at its minimum, at which point / (s) ~ 0.92. The inset shows the original QAOA optimal parameters y _t ^* and b[ for p = 40 With this annealing path, it is possible to follow the instantaneous eigenstate population throughout the quantum annealing process, as shown in FIG. 18F, according to some embodiments. In particular, FIG. 18F shows instantaneous eigenstate populations under the annealing path given in FIG. 18E.

Note that the instantaneous ground state and first excited state swap at the anti-crossing point. In contrast to adiabatic QA, the state population leaks out of the ground state and accumulates in the first excited state before the anticrossing point, where the gap is at its minimum. Using a diabatic transition at the anti-crossing, the two states can swap populations, and a large ground state population is obtained in the end. The final state population from the constructed annealing path can differ slightly from those of QAOA, due to, for example, Trotterization and interpolation, but without being bound by theory, the underlying mechanism can be interpreted as the same, which can also be considered responsible for the diabatic bump seen in FIG. 18C. In addition to the conversion used in the equation above, other prescriptions can be used to construct an annealing path from QAOA parameters, and qualitative features do not seem to change.

[0265] Without being bound by theory, these results indicate that QAOA is closely related to a cleverly optimized diabatic QA path that can overcome limitations set by the adiabatic theorem. Through optimization, QAOA can find a good annealing path and exploit diabatic transitions to enhance ground state population. This explains the observation in FIGS. 17A-17E that TTS _QA0A can be significantly shorter than the time required by the adiabatic algorithm. On the other hand, as seen in FIGS. 17H-17J, TTS _QA0A can be strongly correlated with TTS _QA : QAOA can find a good annealing path, which could be adiabatic or not, depending on what is the best route for the specific problem instance.

[0266] The effective dynamics of QAOA for these exemplary specific instances, as shown in FIG. 18F, can be understood, without limitation, by an effective two-level system discussed in more detail in the following section. In general, the energy spectrum can be more complex, and the dynamics may involve many excited states, which may not be explainable by the simple schematic in FIG. 18 A. Nonetheless, QAOA can find a suitable path via the disclosed heuristic optimization strategies even in more complicated cases. Effective few-level understanding of the diabatic bump

[0267] Without being bound by theory, this section elucidates the mechanism of the diabatic bump discussed above with reference to FIG. 18C via effective few-level dynamics. The intermediate dynamics during quantum annealing can be understood using the basis of instantaneous eigenstates | E _t (t)), where

Equation 76:

[0268] Expanding the time-evolved state in this basis as \ip(t)) = åi cq (t) | e _¾ (t)), the Schrodinger equation can be written as

where h = 1 and the time dependence in the notations is dropped for convenience.

Multiplying the equation by ( E _k \ , the Schrodinger equation becomes

Equation 78:

where ( E _k \ E _k ) = 0 is taken by absorbing the phase into the eigenvector \ E _k ). Written in a matrix form, it can be seen that

where the ground state energy is e ₀ = 0 (by absorbing it into the phase of the coefficients) and Aio=Ei— e ₀ is the instantaneous energy gap from the /th excited state to the ground state. The time evolution starts from the initial ground state with ao = 1 and a, = 0 for i ¹ 0, and the adiabatic condition to prevent coupling to excited states is

Equation 80:

[0269] The first equality can be derived from equation 76. This can produce the standard adiabatic condition T = 0(1/A^ _in ) As discussed above, the minimum gap for some graphs can be exceedingly small, so the adiabatic limit may not be practical. However, is possible to choose an appropriate run time T, which breaks adiabaticity, but is long enough such that only few excited states are effectively involved in the dynamics. This is the regime where the diabatic bump operates and one can understand the dynamics by truncating equation 79 to the first few basis states.

[0270] As an example, FIG. 25A plots the instantaneous eigenstate populations of the first few states, according to some embodiments. In particular, FIG. 25 A plots instantaneous eigenstate populations along the linear-ramp quantum annealing path given by equation 69 for the example graph in FIG. 18B. The annealing time can be chosen to be T= T* = 40, which can correspond to the time where the diabatic bump occurs in FIG. 18C. FIG. 25 A is simulated with the full Hilbert space, but effectively the same dynamics will be generated if the simulation is restricted to the first few basis states in equation 79.

[0271] FIG. 25B shows the strength of the couplings between the instantaneous ground state and the low excited states, according to some embodiments. The plotted quantities measure the degree of adiabaticity (as explained in equation 80). By comparing FIGS. 25A and 25B it is possible to see that T= 7^=40 allows the time evolution to break the adiabatic condition before the anticrossing: population leaks to the first excited state, which becomes the ground state after the anticrossing. Thus, the time scale of T* for the diabatic bump represents a delicate balance between allowing population to leak out of the ground state and suppressing excessive population leakage, which without being bound by theory explains why it happens at a certain range of time scale.

Example graph instance and time-to-solution [0272] In the previous sections, an example representative graph instance where the adiabatic minimum gap is small was considered with reference to FIG. 18B. The low energy spectrum for the graph along the QA path can be seen in FIG. 24A, according to some embodiments. More particularly, FIG. 24A shows an exemplary energy spectrum of low excited states (measured from the ground state energy EGS) along the annealing path for the graph instance given in FIG. 18B. Only states that can, in this example, couple to the ground state are shown, such as states that are invariant under the parity operator P = Y\ ₌₁ of. The inset shows the energy gap from the ground state in logarithmic scale. As shown in the embodiment of FIG. 24 A, only eigenstates that are invariant under the parity operator P = Y\ ₌₁ of are shown, since the Hamiltonian QA(S) commutes with P and the initial state is P- invariant: P |i/>(0)) = \y(b ) _·

[0273] FIGS. 24B and 24C illustrate TTSQA and TTSQAOA for the same graph. In particular, FIG. 24B shows the time-to-solution for the linear-ramp quantum annealing protocol, TTSQA, for the same graph instance. In this embodiment, TTS in the long-time limit follows a line predicted by the Landau- Zener formula, which is independent of the annealing time T. FIG. 24C shows TTS for QAOA at each iteration depth p. For QA, it can be seen that non-adiabatic evolution with T~ 20 yields orders of magnitude shorter TTS than the adiabatic evolution. In addition, TTS in the adiabatic limit is independent of the annealing time T, following the Landau-Zener formula pcs = 1 - exp For QAOA, an exemplary embodiment of the FOURIER[oo, 10] heuristic strategy is used to perform the numerical simulation up to />max = 50, and compute TTSQAOA(P) and TTSQ _A0A (up top £ pmax). For this particular figure, although TTSQ _A0A occurs at p = 49, in some embodiments it may be better to run QAOA at/> = 4 or p = 5 due to optimization overhead and error accumulation at deeper circuit depths. The apparent discontinuous jump in TTSQAOA shown in FIG. 24C is due to the corresponding jump in pGsip), which can be explained by two reasons: first, embodiments of the disclosed heuristic strategy is not guaranteed to find the global optimum, and random perturbations may help the algorithm escape a local optimum, resulting in a jump in ground state population; second, even when the global optimum is found for all level p, there can still be discontinuities in / cs, since the objective function of QAOA is energy instead of ground state population.

Simulation techniques

Details of Simulation with Measurement Projection Noise

[0274] When running QAOA on actual quantum devices, the objective function can be evaluated by averaging over many measurement outcomes, and consequently its precision can be limited by the so-called measurement projection noise from quantum fluctuations, in some embodiments. This effect can be accounted for by performing full Monte-Carlo simulations of actual measurements, where the simulated quantum processor only outputs approximate values of the objective function obtained by averaging M measurements:

Equation 81:

where z _P ,i is a random variable corresponding to the /th measurement outcome obtained by measuring \y _r (g,b )) in the computational basis, and C(z) is the classical objective function. Note when M ¥, F _{p M} ® F _p = (y _r (g, b)\H _€ \y _r (y, /?)). In the exemplary simulation, finite precision _{r M} — F _p \ ^" x can be achieved by sampling measurements until the cumulative standard error of the mean falls below the target precision level x. In other words, for each evaluation of F _P requested by the classical optimizer, the number of measurements performed is set by the following criterion:

Equation 82:

[0275] In some embodiments, it can be expected that M » Var(F _p )/^ ² . To address issues that can appear with finite sample sizes, at least 10 measurements are performed (M > 10) for each objective function evaluation.

[0276] Regarding the classical optimization algorithm used to optimize QAOA parameters, generally, classical optimization algorithms iteratively use information from some given parameter point (y, /?) to find a new parameter point (y\ /?') that produces a larger value of the objective function F _p (y' , b'^ ³ F _p (y, b^. In order for the algorithm to terminate, some stopping criteria can be set. In some embodiments, up to two can be used: First, an objective function tolerance e can be set, such that if the change in objective function <E, the algorithm terminates. Second, a step-tolerance d can also be set so that the algorithm terminates if the new parameter point is very close to the previous one |y'— y| ² + \b'— b\ < d ² . For gradient-based optimization algorithms such as BFGS, d can also be used as the increment size for estimating gradients via the finite- difference method: 5 _p /dyi— [F _{p M>} ( i + d)— F _{p M} (ji) /S. In the simulations discussed herein, the BFGS algorithm is implemented as fminunc in the standard library of MATLAB R20l7b.

[0277] Using the approach described above, it is possible to simulate experiments of optimizing QAOA with measurement projection noise for a few example instances, with various choices of precision parameters (e, x, d) and starting points. For the example representative instance studied in FIG. 19, e= 0.1, x = 0.05, and d = 0.01. Each run of the simulated experiment begins with QAOA of level either p = 1 or p = 5 and continues with optimizing increasing levels of QAOA using embodiments of the disclosed FOURIER[oo, 0] heuristic strategy. The starting point of QAOA optimization can either be randomly selected (when starting at p = 1) or chosen based on an educated guess using optimal parameters from small-sized instances (at p =1 and p = 5). Specifically, the starting points obtained from educated guesses for the example studied in FIG. 19 are ( u° , v°) = (1.4849,0.5409) at level p = 1, and

Equation 83: u° = (1.9212, 0.2891, 0.1601, 0.0564, 0.0292)

Equation 84: v° = (0.6055, -0.0178, 0.0431, -0.0061, 0.0141) at level p = 5. For each such run, the history of all the measurements can be tracked so that the largest cut Cut found after the z-th measurement can be calculated. Each experiment is repeated 500 times with different pseudo-random number generation seeds, and an average over their histories is taken.

Techniques to speed up numerical simulation

[0278] In some embodiments, a number of techniques can be exploited to speed up the numerical simulation for both QAOA and QA.

[0279] For example, first, the symmetries present in the Hamiltonian can be used. For MaxCut on general graphs, the only symmetry operator that commutes with both He and HB is the parity operator P = [H _c , P] = [H _B , P] = 0, and so does [H _QA (S), P] = 0, where HQA(S) is the quantum annealing Hamiltonian equation 69. The parity operator can have two eigenvalues, +1 and -1, each with half of the entire Hilbert space. The initial state for both QAOA and QA are in the positive sector, i.e., P |+) ^®w =|+) ^®w . Thus, any dynamics must remain in the positive parity sector. Hc and HB can be rewritten in the basis of the eigenvectors of P , and the Hilbert space reduced from 2 ^v to 2 ^v - 1 by working in the positive parity sector.

[0280] For QA, dynamics involving the time-dependent Hamiltonian can be simulated by dividing the total simulation time / into sufficiently small discrete time T and implementing each time step sequentially. At each small step, it is possible to evolve the state without forming the full evolution operator, either using the Krylov subspace projection method or a truncated Taylor series approximation. In the simulations discussed herein, a scaling and squaring method is used with a truncated Taylor series approximation as it appears to run slightly faster than the Krylov subspace method for small time steps.

[0281] For QAOA, the dynamics can be implemented in a more efficient way due to the special form of the operators He and i¾, in some embodiments. Work can be performed in the standard s _z basis. Thus, H _c = can be written as a diagonal matrix and the action of e ^~lyHc can be implemented as vector operations. For i¾, the time evolution operator can be simplified as

Equation 85:

[0282] Therefore, the action of e ¹ ^ ^HB can also be implemented as N sequential vector operations without explicitly forming the sparse matrix HB, which both improves simulation speed and saves memory. In addition, in the optimization of variational parameters, the gradient can be calculated analytically, instead of using finite- difference methods.

Techniques similar to the gradient ascent pulse engineering (GRAPE) method are also used, which reduces the cost of computing the gradient from 0(p ² ) to O(p), for a /7-level QAOA. Lastly, in example simulations of embodiments of the disclosed FOURIER strategy, the gradient of the objective function can be calculated with respect to the new parameters (u, v). Since g = A _s u and b = A _c v and for some matrices As and Ac, their gradients can also be related via

QAOA for MIS

[0283] This section discusses exemplary techniques to simulate the Quantum

Approximate Optimization Algorithm to solve Maximum Independent Set Problems, according to some embodiments. Quantum approximate optimization algorithm

[0284] Without being bound by theory, in some embodiments, given a problem to find MIS on a given a graph G = (V, E ), the /i-level QAOA for MIS can be a variational algorithm consisting of the following steps:

[0285] (i) Initialization of the quantum state in \y ₀ ) = |0)® ^w .

[0286] (ii) Preparation of variational wavefunction:

parameters g E M ^p_1 and b E M ^p , specify the variational state.

[0287] (iii) Measurement of Hp .

[0288] The three steps (i)-(iii) can be iterated and combined with a classical optimization of the variational parameters in order to minimize {y _r (9, b)\H _R \y _r (g, b)).

Alternative formulation

[0289] In some embodiments, where U » |W|, |D|, the variational search can be restricted to the sub space Jfis spanned by independent sets, such that the algorithm does not need to explore states that can be directly excluded as MIS candidates. In this limit, it is possible to write:

Equation 87: Equation 88:

where Pis is a projector onto the independent set subspace H is. Evolution with Hp can reduce to simple rotation of individual spins around the z axis. Since

Equation 89: exp(— i gH _R ) exp(— ίbHo) it is possible to commute all the unitaries generated by Hp in Equation 86 to the rightmost side until they act (trivially) on the initial state. Thus, it is possible to rewrite the state

where the following can be identified

Equation 91: Equation 92: t _k = /? _fc

[0290] Thus, some embodiments of the formulation of QAOA given above can be equivalent to equation 90 for U » W.

[0291] Aspects of the present disclosure show that quantum algorithms can be implemented for solving computationally hard problems with coherent quantum optimizers with minimal resources and implementation overhead, which can include, but is not limited to the number of ancillary qubits needed, additional depth of the quantum circuits needed, etc.

[0292] Aspects of the present disclosure show that NP-complete combinatorial optimization problems can be encoded exactly in quantum systems even considering the unwanted interactions between the qubits.

[0293] Aspects of the present disclosure show that quantum algorithms can be implemented by applying light pulses with a variable duration and a variable optical phase to at least some of the plurality of qubits.

[0294] Aspects of the present disclosure show that heuristic optimization strategies can find quasi-optimal variational parameters in variational quantum algorithms in 0(poly(p)) time without the 2 ⁰⁽ · ^r) resources required by brute force approaches using many initial guesses of the parameters. [0295] Aspects of the present disclosure show that quantum approximate optimization algorithms utilize non-adiabatic mechanisms to overcome the challenges associated with vanishing spectral gaps.

[0296] Aspects of the present disclosure show that vertex renumbering of the combinatorial optimization problem permits quantum systems to implement a broader class of problem instances.