Solving optimization problems with local light shift encoding on Rydberg quantum annealers

We provide a non-unit disk framework to solve combinatorial optimization problems such as Maximum Cut (Max-Cut) and Maximum Independent Set (MIS) on a Rydberg quantum annealer. Our setup consists of a many-body interacting Rydberg system where locally controllable light shifts are applied to individual qubits in order to map the graph problem onto the Ising spin model. Exploiting the flexibility that optical tweezers offer in terms of spatial arrangement, our numerical simulations implement the local-detuning protocol while globally driving the Rydberg annealer to the desired many-body ground state, which is also the solution to the optimization problem. Using optimal control methods, these solutions are obtained for prototype graphs with varying sizes at time scales well within the system lifetime and with approximation ratios close to one. The non-blockade approach facilitates the encoding of graph problems with specific topologies that can be realized in two-dimensional Rydberg configurations and is applicable to both unweighted as well as weighted graphs. A comparative analysis with fast simulated annealing is provided which highlights the advantages of our scheme in terms of system size, hardness of the graph, and the number of iterations required to converge to the solution.


I. INTRODUCTION
Using quantum computing for applications requires a faulttolerant scalable architecture [1][2][3][4][5].The current so-called Noisy Intermediate Scale Quantum (NISQ) era is dictated by noisy qubits and gates with comparatively low fidelities [6,7].Despite these limitations, there is ongoing research in trying to find suitable problems for which the quantum approach is superior [7][8][9][10].In particular, there is a class of problems whose exact solutions are hard to obtain in polynomial time, and at best, in certain cases, there exist only approximate solutions to them.They fall under the category of NP-hard problems [11,12] and are regularly investigated beyond the realm of computational complexity theory in order to get better insight into the performance of quantum algorithms [13][14][15].In this work, we consider two such NP-hard problems, namely Max-Cut and MIS, which are combinatorial optimization problems that have many practical applications [16][17][18][19][20].
When solving a combinatorial optimization problem with quantum systems, there are two pertinent challenges that arise: (i) the choice of quantum hardware and (ii) the choice of quantum algorithm to be implemented.An ideal quantum hardware should allow scalability with the number of qubits.Qubits of superconducting devices still suffer from noise and rely on a fixed architecture with localized connectivity [21], while trapped ions are still too sensitive to external fields to allow scalability for more than 50 qubits [22,23].Motivated by recent experiments where hundreds of neutral Rydberg atoms can be trapped and whose interactions can be manipulated [24][25][26], there is a tendency and perspective to use Rydberg platforms to solve optimization problems, which is also the focus of this work.
Very recent developments in solving optimization problems on Rydberg platforms include solving the weighted MIS graph problem on a Rydberg simulator [27] and the Max-Cut problem using Rydberg gates [15] whose approximation ratio suffers due to noisy gates.In both works, QAOA (Quantum Approximate Optimization Algorithm) approach is used to approximate the optimal solution [28,29].In general, variational methods require finding optimal values for a large number of variational parameters (initial state ansatz, choice of the mixer, number of layers), which by itself is found to be an NPhard problem [30,31].This can make the implementation of variational algorithms to solve a particular optimization problem fairly cumbersome.For this purpose, we propose an optimized quantum annealing protocol that is implemented on the Rydberg platform.
Another issue is the choice of encoding which is a scheme by which a real-world optimization problem is mapped onto a system of interacting qubits.The encoding scheme depends on the choice of both, hardware and algorithm.Max-Cut and MIS are optimization problems that can be represented in a QUBO form.In this representation, the objective function along with its constraints has a quadratic form with variables that take binary values [32,33].The advantage of formulating the problem in QUBO form is that it can naturally map to interacting spin models [34].However, it should be mentioned that depending on the type of optimization problem, this mapping is neither unique nor trivial to be physically realized [27,[35][36][37][38].For example, in [27], the MIS problem was encoded on an interacting spin system using the unit-disk encoding [39] that relies on the perfect implementation of the Rydberg blockade effect [40] between the atoms.Based on the unit-disk encoding, a problem graph with N vertices will require ∼ 4N 2 physical spins which can be a substantial over-head of resources [41].
In this work, by taking advantage of localized light shifts on individual atoms, we encode and solve the Max-Cut/MIS problems using a non-unit disk mapping to spin models.The proposed scheme has the benefit of tackling weighted as well as unweighted graphs within the same framework and can in principle be generalized to other QUBO problems.By combining gradient and non-gradient optimal control methods, we implement the temporal evolution of the laser parameters that make the annealing process on the Rydberg platform efficient to solve optimization problems.Thus we obtain optimal solutions for different graphs (with sizes 5 − 15, up to degree 5) with varying hardness within 3 − 15µs with an approximation ratio close to one.The hardness parameter measures the complexity of the problem graph and serves as a convenient metric for benchmarking the performance of our optimal quantum annealing method.For graphs with similar size and complexity, we find that our protocol outperforms the simulated annealing algorithm [42] in terms of accuracy.Our encoding scheme should be applicable on any quantum platform that allows local qubit control along with global driving of the many-body system to its ground state.
We proceed as follows.In Section II we outline the theoretical description for the Rydberg simulator and introduce the Max-Cut/MIS problems.We propose a three-step scheme for non-unit disk encoding and solving graph problems using the Rydberg annealer architecture in Section III.In Section IV we discuss the tools to characterize the complexity of the problems and the quality of the solutions.Finally, in Section V, we analyze our results, and detail the physical implementation of our scheme in Section VI.Section VII contains our conclusions and outlook.

II. THEORY
In Sec.II A, we provide a brief introduction to the manybody Rydberg Hamiltonian which provides the setup for quantum annealing while Sec.II B defines the Max-Cut and MIS problems.

A. Hamiltonian of Rydberg Simulator
We consider ultracold atoms trapped in optical tweezers [43,44].They are described as two level systems consisting of a ground state (|g⟩), and a Rydberg excited state (|e⟩) that are optically coupled using a laser [45,46].The Hamiltonian expressed in the atomic basis reads as follows, where V k j = C 6 /|r j − r k | 6 is the van der Waals (C 6 > 0) interaction between the pair of Rydberg atoms |e⟩ k and |e⟩ j where C 6 is the associated dispersion coefficient.Here r j and r k are the positions of the two atoms labeled as j and k.
The detuning ∆ j is the site-dependent laser parameter describing the difference between the frequency of the applied field and the natural frequency associated with the atomic transition.The Rabi frequency Ω is the global laser parameter that couples the two states and is proportional to the intensity of the driving field and the dipole moment associated with the atomic transition.Using the Pauli spin operator representation and setting Ω = 0, the above equation reduces to the Ising Hamiltonian with a longitudinal field (up to a constant where the negative sign is absorbed in ∆ j .The first term is the local longitudinal field which depends on the detuning as well as N k=1, j k V jk , while the second term corresponds to the interaction between spins.The above Hamiltonian is useful for the connection to optimization problems that can be formulated in the QUBO framework.

B. Max-Cut and MIS Problem Definitions
QUBO problems are often represented by a graph G(V, E) with vertices V and edges E which take different representations depending on the problem as shown in Fig. 1.Below we define the classical cost functions for the Max-Cut and MIS problems respectively.
Max-Cut: Given a graph G with weights w jk associated with each edge ( j, k), a subset of vertices from the graph is referred to as a cut.This cut becomes a Max-Cut when the vertices are split into two sets in a way that maximizes the total weight (Max-Cut value) of edges between them.The problem graph is illustrated in the upper panels of Fig. 1(a-b).Panel (a) displays a weighted graph of size 4 and (b) shows a dashed curve representing the cut that divides the graph into two sets, blue and red.Both sets are interchangeable as the solution is degenerate.The Max-Cut cost function, which is to be maximized, is expressed as, where the sum is over all the edges in the graph G.The variable X j ∈ {0, 1} represents the two sets in a cut.

MIS:
Given a graph, G with weights w j associated with each vertex j, a subset of mutually non-adjacent vertices is said to be a maximum independent set if it has the largest possible sum of weights w j over all vertices in the subset.Figure 1 (de) depicts the MIS problem.MIS can be found by maximizing the following cost function, where the sum is over all the vertices V in the first term and in the second term, the sum is over all the edges E. Here, X j ∈ {0, 1} where X j = 1 indicates that a vertex j is in the independent set.The first term in the cost function contributes only when X j = 1, and the second term ensures that no two vertices with an edge between them belong to the independent set.It is done by suppressing the occurrence of simultaneous X j = 1 and X k = 1 along an edge.This section discusses the protocol adopted in this work to solve the QUBO problems.The first step in solving the Max-Cut/MIS problems using Rydberg annealers is to encode the problem graph onto the physical spins (qubits) of the Rydberg setup.In principle, this encoding scheme is not unique and the most widely used approach, in the context of Rydberg simulators, is the unit disk (UD) encoding [39,47].Although it is the natural choice of encoding owing to the fact that the Rydberg blockade effect successfully implements unit disks in a straight-forward manner [27,35,41], it has its drawbacks.In order to solve a particular QUBO problem (in this case Max-Cut/MIS) graph with N vertices, it requires solving a graph with ∼ 4N 2 unit disks.This significant overhead along with the issue of unwanted interactions calls for alternative encoding schemes.
One of the highlights of this work is to explore an alternative, non-blockade-based encoding scheme that allows us to map QUBO to Rydberg annealers both for weighted and unweighted graphs using a single framework and it possesses a linear scaling with respect to N. In Steps 1 and 2, we outline our scheme of encoding the Max-Cut and MIS problem graphs onto the Rydberg spin model, while in Step 3, we discuss the implementation of optimal quantum annealing using Rydberg atoms.
Step 1: Mapping of cost functions using local detuning The cost functions C Max−Cut and C MIS defined in the earlier section can be directly mapped to the spin Hamiltonians using the standard definitions of the Pauli operators.The details are provided in the Supplemental material S-I and the corresponding spin Hamiltonians with N spins are given in the following, where is the spin Hamiltonian for the Max-cut cost function Eq.( 3).Comparing Eqs. ( 2) and (5), we define the detunings ∆ j in a manner such that they cancel the effective longitudinal field in ĤIsing .Thus, for an atom at j th position, the local light shift is chosen to be ∆ j = −1/2 N k j V jk .The edge weights w jk from the graph are identified as the interaction strength V jk between the atoms labeled as j and k as V jk = 4w jk .This implies that the value of the detunings is related to the weights according to Similarly for the MIS problem, we get the Hamiltonian for interacting spins to be e jk w j w k σz In the case of the MIS problem, e jk = 1, if and only if there is an edge between vertices j and k, otherwise e jk = 0. On comparing Eq. ( 2) with (7), the vertex weights w j from the graph are related to the interaction strength V jk between the atoms labeled as j and k as V jk = e j,k w j w k .The detunings ∆ j get chosen to represent the effective longitudinal field for an atom at j th position.This gives the following relationship, Thus, the mathematical problem of maximizing Eqs. ( 3) and ( 4) is reduced to a physical problem of finding the manybody ground state (MBGS) of the spin Hamiltonians given by Eqs. ( 5) and ( 7) respectively.The scheme involves temporally varying the laser parameters till the specific choice of detuning values are obtained at the end of the protocol, all the while minimizing the energy of the full system.This will be elaborated in Step 3. The mapping of classical cost function to the Rydberg Hamiltonian described here is a direct encoding type which results in the linear scaling of the number of atoms needed to represent the vertices of the graph.Consider a problem graph G(V, E; W) where n(V) and n(E) are the number of vertices and edges respectively.Thus for an all-to-all connected graph, n(E) takes the maximum value which is However, the number of edges for a general graph can vary from 1 to n(E All ) depending on the problem that needs to be solved.This implies that the dependence of n(E) on n(V) can be expressed as where α ∈ [0, 2].In some sense, the degrees of freedom for a graph can be understood in terms of n(E).Hence the number of independent degrees of freedom in w jk for generic graphs does not always scale quadratically with n(V) but can have α < 2. For a direct encoding, the number of atoms required are same as the number of vertices in the problem graph.
If the direct encoding of the graph with N vertices were to be implemented using only interactions between atoms, then unwanted interactions cannot be ignored.This leads to 1 2 43 unavoidably creating unwanted edges that did not exist in the original problem graph.
Step 2: Spatial arrangement of atoms As seen in the previous step, the distance-dependent interactions encode the information about the weights of the graph.Therefore, before any evolution of the many-body Hamiltonian, the initial arrangement of the atoms is crucial to the encoding scheme as is represented schematically in Fig 1 .For example, in the case of a weighted Max-Cut graph with weights ω 34 > ω 14 > ω 24 > ω 12 will corresponds to atoms (purple balls) arranged in a configuration such that r 34 < r 14 < r 24 < r 12 which are trapped in tweezers (shown in green).Similarly, for the weighted MIS problem with weight relationship given by ω 3 > ω 1 > ω 4 > ω 2 is associated with a Rydberg configuration where r 34 < r 14 < r 12 < r 24 .Since the weights ω i j are related to the interactions V(r i j ), which in turn is reflected in the choice of the r i j , the atoms are arranged in a manner that provides a true representation of the graph.However, it is possible to have interactions between atoms that do not share an edge in the original graph problem.Such unwanted interactions become a serious issue for graphs with higher degrees.For this purpose, in 2D geometry, the maximum degree of a single node in the graph is limited to five for which we numerically checked that the unwanted interactions do not play a significant role.Apart from the unwanted interactions, the arrangement of the atoms which do have an edge in the problem graph is also limited.This case appears when the weights in the graph are such that the resulting relative distances between the atoms cannot be realized on a 2D plane.Although the graphs that can be represented by the geometrical arrangement of atoms are limited, we emphasize that our scheme is applicable to a large class of graphs out of which only a few are demonstrated in this work.One of the ways to generate the family of accessible graphs is by considering geometrically implementable graphs as the basis and connecting them through an additional node between them.However, adding a third dimension can also increase the value of the allowed degree and the flexibility of the atom arrangement.
Step 3: Optimal Quantum Annealing with Rydberg atoms The previous two steps outlined the encoding of the graph problem to the Rydberg spin Hamiltonian.The goal is to find accurate solutions to a problem graph as efficiently as possible with respect to the number of iterations and run time.This is achieved by numerically solving the spin dynamics using the following Rydberg Hamiltonian for a specific arrangement of atoms ĤRyd = Ω(t) where σx j = |e⟩ j ⟨g| + |g⟩ j ⟨e| and excitation number operator ne j = 1 2 σz j + I as |e⟩ j ⟨e|.The objective is to reach the target Hamiltonian ĤMax−Cut,MIS (Eqs.5 and 7) while minimizing the energy of the system to obtain the instantaneous ground state whose configuration provides the solution to the optimization problem.Initially, all atoms are in the ground state |gg...g⟩ which corresponds to a large non-zero detuning value ∆ j (t = 0) 0. During the protocol, the detuning value on each atom ∆ j (t) varies but attains a specific value at the end of the protocol as defined by Eqs. ( 6) and (8).The initial and final values of the Rabi frequency are set to zero.At intermediate times, non-zero values of the Rabi frequency Ω(t) provide a transverse field for the above Ising Hamiltonian.This generates quantum fluctuations and causes different many-body states to couple with each other thereby accessing a larger part of the Hilbert space.This is key to the quantum annealing process as it explores different configurations through the application of the non-zero transverse field while traversing the energy landscape to get to the desired many-body ground state configuration (i.e. the solution to the problem graph).The idea is to use well-established techniques of optimal control theory applied to quantum systems [48][49][50][51][52][53][54] in order to reach this state in an efficient manner well within the system lifetime.In particular, a combination of gradient and nongradient-based methods (see Supplemental Material S-III for details) are used to shape the pulses in time which allows us to optimally steer the many-body state towards the true solution.The initial guesses for detuning and Rabi frequency are inspired by the adiabatic evolution and are then optimized in time.The objective function that needs to be minimized during the optimal control is the expectation value of the target Hamiltonian ĤMax−Cut,MIS with respect to the instantaneous many-body state |ψ inst (t)⟩ which is given as E is used as the cost function for the optimization but is not the typical energy of the instantaneous many-body state.Apart from E, we also define the overlap between the instantaneous many-body state with the many-body ground state ψ g of the problem Hamiltonian ĤMax−Cut,MIS during the protocol given as The fidelity F is calculated over all the degenerate ground states for an aposteriori analysis of the quality of the solution obtained.More details about the control methods are provided in Supplemental Material S-III.

IV. CHARACTERIZATION OF MAX-CUT/MIS PROBLEMS
In order to benchmark our protocol (algorithm), we need to ascertain how close is the solution provided by our algorithm to the true solution for problem graphs and compare it with other methods.This is encapsulated in a quantity referred to as approximation ratio which is discussed in this Section.With regards to the efficiency of obtaining the solution for a given algorithm, it is possible to have a scenario where finding a solution for a particular graph is faster in run-time than compared to another graph of similar size.This indicates that the latter graph is more complex or hard to solve which is an inherent property of the problem graph.However, it is nontrivial to characterize the complexity of arbitrary graphs.Motivated by [27,41], we generalize their hardness parameter to include a broader class of problems using the notion of degenerate sub-spaces.

A. Approximation ratio
In general, the approximation ratio quantifies the worstcase performance of an algorithm for solving a particular problem.In certain cases, this ratio can be evaluated analytically.This is the case for the Goemans-Williamson algorithm [55] solving the Max-Cut problem which is considered to be the best classical approximate algorithm with an approximation ratio of 0.878.Any quantum algorithm outperforming this ratio is suggestive of having an advantage over the classical case.In order to benchmark quantum algorithms, it is more convenient to evaluate this ratio numerically as defined in [56] and give as where C opt is the optimal value of the cost function and C obt is the obtained value of the cost function evaluated through  The table shows a one-to-one correspondence between the degenerate solutions (DS) of the original Max-Cut and MIS problem to the MBGS of the corresponding Rydberg Hamiltonians given by Eq. ( 5) and ( 7) respectively.
our method.In order to evaluate C opt , the exact solution to the given problem should be available a priori.In the Results section, we evaluated R for different graphs.

B. Hardness parameter
For Max-Cut and MIS, there are certain graphs such as weighted ones that are generally more challenging to solve computationally [57].One way to identify the complexity in a graph is by looking at the symmetries in the adjacency matrix [58].However, a physically more intuitive way is to study the degeneracy in the many-body ground state of the interacting spin systems.As a result of mapping the cost functions using local detuning described earlier, there is a one-toone correspondence between the degenerate solutions of the original Max-Cut and MIS problem to the many-body ground state (MBGS) of the corresponding Rydberg Hamiltonians as shown in Table I.The larger the degeneracy for the manybody ground state (solution space), the higher the probability to get to the optimal solution during the dynamics.But if the orthogonal sub-spaces (not belonging to the solution space) also have large degeneracy, then there is a possibility for the solution to get stuck in one of these unwanted sub-spaces.
Using this notion of comparing degeneracy between relevant sub-spaces, a hardness parameter was defined only in the context of the MIS problem [27], HP = Here |MIS | and D |MIS | correspond to the size of the MIS and the degeneracy of the solution space.D |MIS |−1 defines the degeneracy of the sub-space of suboptimal independent sets which has one less arbitrary element (vertex) from the MIS solution set.This hardness parameter is limited as it does not take into account all the other sub-spaces (with sizes < |MIS | − 1) with relatively higher degeneracy.Furthermore, this hardness parameter is not well-defined for the solution to the Max-Cut problem.
Thus in this work, a more general hardness parameter is introduced for both Max-Cut and MIS.This parameter includes information about the degenerate sub-spaces and is defined as where D opt represents the degeneracy of the solution space, C opt is the optimal value of the cost function and D represents the degeneracy of a sub-space.The sum is taken over all orthogonal sub-spaces whose degeneracy is greater than a specific cutoff value.The cutoff value is arbitrarily chosen until the value of the hardness parameter HP shows convergence.
It should be noted that if D opt > D cuto f f then this degenerate configuration is excluded from the numerator sum.The hardness defined for a problem instance is an inherent property of the problem, independent of the algorithm.This is because the hardness here characterizes if one instance of the problem has a more non-convex optimization landscape as compared to the other instance.If a problem has more degenerate subspaces, the optimization landscape will be more non-convex and hence harder to navigate for any classical algorithm.We numerically find that our hardness parameter faithfully captures the complexity of the problem graphs for various cases as shown in the Results section.
V. RESULTS

As mentioned in
Step 2 of Sec.III, although we consider prototypical graphs with a degree not more than five and with restricted connectivity that is realizable in a two-dimensional array of atoms in tweezers, we have solved for a large class of graphs out of which only a few are shown here.These other graphs can e.g.be generated by combining smaller graphs that are known to be solvable through an additional node between them.Refer to Supplemental material S-II for numerical simulation details.
To understand the working principle of our optimal quantum annealing protocol, we consider a relatively simple case of the weighted graph of size 5 and degree 3 for which the Max-Cut and MIS solutions are shown in Fig. 2(a-e) and Fig. 2(f-j) respectively.The optimized shapes of the laser parameters (dotted black symbols for detunings, solid dark red for Rabi frequency) are shown in panels (b,g), plots of fidelity (dashed green) and energy (solid blue) in panels (c,h), plot of ordered energies E i (t) along with basis state contribution in the instantaneous eigenstates (indicated by green color bar) in panels (d,i) and the population of basis states at three different time steps in (e,j).The solutions for both problem graphs are obtained efficiently (∼ 3.5µs) and with high fidelity (∼ 0.99).By turning on the Rabi pulse, one obtains multiple avoided crossings in the many-body energy spectrum.The rate at which the individual detunings are swept control the Landau-Zener transitions across multi-level crossings.This results in a non-trivial path of dynamics for the many-body system to its target state.Let C = [∆ 1 (T ), ∆ 2 (T ), ..., ∆ N (T )] be the set of final local detunings for N atoms determined apriori by the weights on the graph given by Eqs. ( 6) and (8).The ratio between the set of detunings ∆ 1 (T ) : ∆ 2 (T ) : ... : ∆ N (T ) is kept fixed during the entire protocol.Thus a single timedependent parameter ∆ G (t) is defined such that the temporal variation of each individual detuning at each atom j is given as ∆ j (t) = ∆ G (t)∆ j (T ) as shown in Fig. 2 (b,g).This factor ∆ G (t) is initialized such that it starts at a negative value at t = 0 and increases to 1 at t = T , thereby providing the desired set of detuning values at the end of the protocol.The physical implementation of the local-detuning protocol is elaborated in Sec.VI.The quantity E (defined by Eq. ( 10)) is minimized and this is reflected in F (defined by Eq. ( 11)), shown in (c,h).The approximation ratios R (defined by Eq. ( 12)) of the final state are also indicated.For a general Max-Cut problem graph, every cut has at least 2-fold degeneracy because of the symmetry between up and down spins.This is also the case for the chosen graph and thus the Max-Cut problem has more degenerate states compared to the MIS problem.This key feature is reflected in the initial part of the dynamics.The change in E during the initial part of dynamics is synonymous with the change in the instantaneous eigenstates (panel (h) and (i)) for the MIS problem.However, this is not the case for the Max-Cut problem.For Max-Cut protocols, we find that E stays constant for a more extended initial period signifying that the energy of the instantaneous state is close to the high energy manifold of the ĤMax−Cut and remains there longer as a result of the degeneracy.Both for Max-Cut and MIS, the initial state |ggggg⟩ is populated as shown in the left-most panel of (e) and (j).At t = T for Max-Cut (right-most panel of (e)), the degenerate ground states (|egegg⟩ and |gegee⟩) are populated while for the MIS case, a single state |eggge⟩ is populated (right-most panel of (j)).At intermediate times (middle panels of (e) and (j)) due to the non-zero values of the transverse field (Ω 0), multiple states get populated.Specifically, 12 out of 32 in (e) and 9 out of 32 in (j) highlighting the fact that the optimal protocol dynamics across the energy landscape is non-intuitive.This aspect is reinforced in Fig. 3 where we have 15 spins.topology and the exact solution for Max-Cut and MIS.For this graph, there will be in general 15 final detuning values ∆ j (T ) at time T .The set of local detunings will all vary with the same fixed ratios between them.For illustration purposes, we show the variation of the factor ∆ G (t) in (b,f).Since the complexity of the problem is significantly increased, we have a different choice for the initial parameters.In particular, a linear guess is used for ∆ G (t) in (f) while a linear ramp with a flat top is used in (b).Panels (c,g) show the corresponding variation of the expectation value E and fidelity F during the protocol.Another signature of the complexity of this problem graph can be seen in panels (d,h) where we find the significant population of a large number of basis states out of the 2 15  states, especially at inter-mediate times.In panel (d), there are 1380 populated states at t = 3.50µs, and 310 populated states at t = 7.10µs.Similarly, in panel (h), there are 678 populated states at t = 3.50µs, and 1029 populated states at t = 7.10µs.It is remarkable that at the final time, we end up with a small set of degenerate MBGS.This is a direct consequence of the optimizer minimizing E which forces the system to transfer the population to the lower energy levels of the target Hamil-tonian.Despite the population of so many basis states, the optimal protocol ensures efficient transfer to the optimal solution (R ∼ 0.99) making the relevance of energy gaps redundant in the context of multi-level dynamics.One may speculate that this intuition will hold for larger systems.
Fig. 4 compares the performance of our protocol (LOQAL) with fast simulated annealing (SA) for different graphs.More details about SA are provided in the Supplemental material S-IV.In particular, we illustrate the variation of the approximation ratio error 1−R (R is defined by Eq. ( 12)) with varying system size N = 2 − 15 and hardness parameter HP (defined by Eq. ( 13)).In the case of Max-Cut, the approximation ratio error (1 − R) shows similar trends with respect to the system size N and the hardness parameter HP as shown in (a,b).As system size increase, so does the number of degenerate orthogonal sub-spaces (which affects the numerator of HP).But since the degeneracy of the solution space for the Max-Cut is always 2-fold (this affects the denominator of HP), the overall hardness parameter increases with system size.For the MIS problem, the SA gives a general upward trend with oscillations in the approximation ratio error with increasing system size N as seen in panel (c).These oscillations are attributed to the degeneracy of the solution space being linked to the system size which by construction is 2-fold for even N and 1-fold for odd N for the graphs considered in this study.The larger the degeneracy of the solution space, the easier it is for the system to reach one of the solutions.This oscillatory behavior is not reflected in the hardness parameter because our definition of HP takes into account the degeneracy of all the relevant subspaces and thereby successfully characterizes the problem graphs.For the graphs considered in Fig. 4, it clearly shows the quantum algorithm is more robust and performs better than SA.
All results so far did not include any noise in the dynamics.Using Max-Cut as an example, we apply the optimal quantum annealing protocol to a noisy system as is expected in real experiments.Fig. 5 displays the Max-Cut problem of size 10 with noise up to 8% in the laser parameters chosen randomly which is shown as shaded regions across the bold lines (mean value) in (b,d).Two different approaches are adopted to simulate the noisy model.In (b) the noise is added after the optimization of the pulses while in (d) it is added during the optimization procedure and the optimizer adapts to the noisy protocol.The shaded region of the fidelity F is broader than that of E (which is barely visible) and is shown in (c), indicating that slight variations in the parameters bring a small change in E, but drastic changes in F. This is due to the fact that the gap in the lower energy levels of the target Hamiltonian is vanish-ing.That results in the system occupying low-energy excited states close to the MBGS, at the same time dropping the fidelity as it is a quantity highly sensitive to the occupied states.In order to better simulate real experiments, a random error of up to 8% is added to the laser parameters during each run which is shown in (d).This causes the optimization landscape to adapt to the random variations in the parameters.Thus in (d), the shaded region of the pulses is not as broad as in (b).In (e) both the expectation value E and the fidelity F show no significant effect of the noise thereby being more resilient to the laser noise.

VI. PHYSICAL IMPEMENTATION OF THE PROTOCOL
For the numerical simulations shown in this work, we considered Rydberg states 60S 1/2 of the Cs atoms which have a van der Waals coefficient C 6 ∼ 139 GHz • µm 6 and a radiative lifetime of ∼ 234µs [25,59,60].An essential aspect of the protocol is to realize an array of trapped atoms with adjustable inter-atomic distances ranging from 1 − 7µm which should be achievable with the current state-of-the-art optical tweezer technology [43].An ingredient of our protocol is the optimal control of the laser parameters.Apart from one global near-resonant laser that couples the ground state atom to its Rydberg state with Rabi frequency Ω, we have another laser whose intensity will be distributed over the atoms selectively The optimal protocol for the Rabi frequency (solid dark red) and the factor controlling the detunings ∆ G (dashed black) are shown in (b) where noise is added to the laser parameters of the pre-optimized protocol and in (d) where noise is added during the optimization of the parameters.The shaded regions represent the fluctuations in the laser parameters for each run that were chosen from a random distribution.(c) and (e) show the corresponding expectation value (solid blue) of the problem Hamiltonian with respect to the instantaneous state (E) and the fidelity F (dashed green) of the instantaneous state with respect to the ground state for both cases, where R indicates the approximation ratio.
using a spatial light modulator [61].This will provide individual local light shifts for the ground-Rydberg transition thereby implementing specific local detunings.At the end of the experimental sequence, measurements of the distribution of the Rydberg excitations can be performed by fluorescence imaging.The whole process will be repeated multiple times to efficiently calculate and classically minimize the expectation value of the problem Hamiltonian (Eqs.10).
We considered the case of noise in laser parameters in our simulations but the experiments can suffer from a variety of errors, introducing different types of noise to the system which includes motional dynamics, interaction induced dephasing [62], blackbody radiation [63] and state preparation and measurement (SPAM) errors [64]; all of which could lead to a significant drop in readout fidelity.There are ways to address and control errors on the Rydberg platform, including converting them into erasures [63] and transforming complicated error models into Pauli-Z errors through the use of ancillary atoms [65].Future works will involve optimizing for noisy models using our protocol.

VII. CONCLUSIONS AND OUTLOOK
Rydberg atoms can achieve the required scalability with arbitrary connectivity between qubits making them highly desirable platforms to solve optimization problems.The protocol presented in this work is fairly universal in the sense, it can handle both weighted and unweighted graphs within the same framework which is not always obvious for other schemes and possibly could be generalized to other QUBO problems beyond Max-Cut and MIS.Although the focus of this work is on a Rydberg annealer, any quantum device that permits local qubit control along with the global driving of the many-body system can implement our protocol.A promising aspect of the optimized dynamics is that it goes beyond the energy gap dependence that would usually limit the adiabatic protocols for large system sizes, which need to be verified for larger systems.
The motivation for presenting an alternative encoding scheme that is implementable on an optimized Rydberg annealer is to have a more favorable scaling of the required number of qubits for problem graphs with N vertices.For the sizes and complexity of graphs considered here, we do find accurate solutions for both, the Max-Cut and MIS problem graphs with N = 5, 15 vertices that do have a O(N) scaling.However, the topology of the graphs is limited due to the issue of unwanted interactions.Such effects can be mitigated using a three-dimensional arrangement of atoms but further investigations are needed to generalize our encoding scheme for graphs with arbitrary connectivity and degree.Despite having shown the advantage of our protocol with respect to fast simulated annealing, benchmarking our method with other classical algorithms can shed more insights into the performance of our method.A more thorough analysis of optimal dynamics that can adapt to different types of noise and in particular where Bayesian methods [53,54,66] can prove to be useful will be the scope of future work.
where X j ∈ {0, 1}.The sum is over all the vertices V with weights W in the first term and in the second term, the sum is over all the edges E. Similar to the case of Max-Cut, X j is replaced by Z j and the cost function C MIS becomes Now ( j,k)∈E w j w k Z j + Z k contains the sum over all connected to the j th vertex.So for each vertex, the contribution is coming from the neighbors alone, where S j is the set consisting of the neighbours of the j th vertex.The cost function then becomes The first two terms are just constants.If Z i is replaced by σz i , the problem of finding the maximum of C (or minimizing −C) is equivalent to finding the many-body ground state of ĤWMIS , The sum ( j,k)∈E in the second term in Eq. (S15) can be converted to a double sum as follows, e jk w j w k σz where e i j = 1, if and only if there is an edge between vertices i and j, otherwise e i j = 0.The parameter e i j captures the information about the neighbors of the vertex i.The sum j∈S i is over all the neighbors of spin i and is replaced by N j=1, j i e i j .

S-II. NUMERICAL DETAILS OF RYDBERG ANNEALER
The time-dependent Rydberg Hamiltonians provided by Eq. ( 9) is implemented using the sesolve function from QuTip library [68].The objective is to find the many-body ground state of the target Hamiltonians given by Eq. ( 5) for Max-Cut and Eq.(7) for MIS (see main text).The expectation value E of the target Hamiltonian with respect to the instantaneous many-body state is minimized during the optimal quantum annealing protocol.The profile of a global Rabi frequency Ω(t) is optimized in time while the initial and final values remain zero.An initial guess for Ω(t) is provided such that Ω(t) = A(1 − cos(πt/T )) 2 , where T is the total time of the protocol and A is of the order of few MHz.A localized detuning ∆ j (t) on each atom is also optimized in time where each ∆ j (t) at time T is given by Eqs. ( 6) and (8).For a particular graph, the ratios between individual ∆ j (T ) are kept fixed during the entire protocol, hence, optimizing a single pre-factor ∆ G (t) will result in the optimization of individual detunings given by ∆ j (t) = ∆ G (t)∆ j (T ).∆ G takes a negative value at t = 0 such that the many-body ground state of the Hamiltonian at t = 0 is |gg...g⟩, which means all atoms are in the ground state.∆ G is then varied in time and reaches one at t = T , the many-body ground state of the Hamiltonian at t = T corresponds to the solution of Max-Cut/MIS problems.The total time T is divided into sub-parts by selecting 8 points, Ω(t) and ∆ G at these 8 time points are optimized during the optimal quantum annealing to reduce E and are connected by b-splines.After the system reaches its final state, fidelity F is also calculated to measure the accuracy of the obtained many-body state.F represents the probability of finding the system in one of the ground states, which differs from the approximation ratio.The approximation ratio R is also calculated at the end of the protocol to measure the quality of the solution.

S-III. OPTIMAL CONTROL THEORY METHODS
Optimal control theory is a mathematical framework that helps determine the best way to control a dynamical system by finding parameters that extremize a specific objective function [69].It involves solving an optimization problem by adjusting control inputs over time while considering system dynamics and constraints.In physics, it has been used for shaping laser pulses, gate operations, and controlling chemical reactions [70][71][72].For such optimization problems, gradient [73,74] and non-gradient [75,76] methods can be used.Gradient-based methods rely on calculating the gradient of the objective function with respect to the parameters and on updating them iteratively in the direction of the negative gradient to reach an optimal solution.Whereas, non-gradient methods are typically heuristic or evolutionary algorithms that iteratively explore the search space to reach an optimal solution.In this work, we use a combination of BFGS (gradient-based) [77][78][79][80] and Nelder mead (non-gradient based) [81].A description of both methods is given below.
A. − Broyden-Fletcher-Goldfarb-Shanno (BFGS) [77][78][79][80] is one such gradient-based method that approximates the inverse of the Hessian matrix of the objective function using information from the gradients of the function.At each iteration, the BFGS calculates the change in the gradient and uses it to update the current estimate of the inverse Hessian.The new estimate of the inverse Hessian is then used to determine the search direction for the next iteration.The update formula is given by, where H k is the inverse Hessian approximation at iteration k, s k = x k+1 − x k is the difference between the current and previous estimates of the parameters, is the difference between the current and previous estimates of the gradient, ρ k = 1/(y T k s k ), and I is the identity matrix.The BFGS method typically starts with an initial inverse Hessian approximation, H 0 , and iteratively updates the approximation until convergence is reached.The search direction at each iteration is given by d k = −H k ∇ f (x k ), and the step size is determined using a line search method.B. − As for non-gradient based, Nelder-Mead (NM) [81] is one such method, it iteratively searches for the minimum of an objective function.NM relies on exploring the simplex, a geometric figure with n + 1 vertices in n dimensions.At each iteration, the algorithm evaluates the objective function at the vertices of the simplex and then performs a set of operations to update the simplex, such as reflection, expansion, contraction, or shrinkage, based on the values of the function.The method continues until a stopping criterion is met, such as reaching a maximum number of iterations or when the function value change is small.Initialization is a step in which a simplex is defined based on an initial set of n + 1 points in n-dimensional space.The function values at the vertices of the simplex are evaluated.The worst vertex of the simplex (i.e., the vertex with the highest function value) is then reflected through the centroid of the remaining vertices.If the reflected vertex has a lower function value than the second-worst vertex, the simplex is expanded in that direction.Otherwise, the reflected vertex is retained.If the reflected vertex has a higher function value than the worst vertex, the simplex is contracted towards the best vertex.If the contracted vertex has a lower function value than the worst vertex, it is retained.Otherwise, the simplex is shrunk towards the best vertex.If none of the previous steps improve the function value, the simplex is reduced toward the best vertex.Finally, the algorithm terminates when a stopping criterion is met, such as a maximum number of iterations or a slight change in the function value.
BFGS when compared to NM, can converge faster for smooth, convex functions.But at the same time, BFGS may get stuck in local minima or saddle points if the function has many local optima.Consequently, a combination of both methods is used in this work.The methods were applied sequentially in the order of BFGS-NM-BFGS (B-Nel-B), balancing exploration and exploitation.Starting with BFGS, the algorithm can quickly converge to a local minimum, and then NM is used to explore other regions of the objective function and get closer to the global optimal value.Finally, using BFGS can refine the solution and potentially converge to a better optimal value.Parameters for NM and BFGS.− Implementation in the code is done by the optimize.minimizefunction from the python library SciPy [82].For both methods a convergence factor of 10 −4 was set.The first layer of BFGS has a maximum iteration variable which was set to 3 − 6 depending on the size of the problem and was set to 2 for the last layer.In NM, both maximum iterations and maximum number of objective function evolution need to be fixed and were given a value of 300 to get convergence.These values were set based on the competition between the time it takes to evaluate one layer and the quality of the solution.In the case when the system was noisy, three times more the number of runs were required as compared to the noiseless case.

S-IV. SIMULATED ANNEALING
In this section, details of simulated annealing (SA) are discussed, which is used for benchmarking the quantum control method used in this work for the specific graphs.All the simulations for this method are performed by using an optimized SA from the SciPy library [82].
Method: Simulated annealing [83] is an algorithm used to solve optimization problems through a global search approach.It works by simulating the process of heating and cooling a material to reduce defects and minimize energy.With each function call, the algorithm searches for a new solution point in the search space.If the new point has lower energy than the previous point, it's accepted as the new optimal value with a probability of 1.If not, the probability of acceptance depends on the temperature.As the temperature decreases, the algorithm becomes more selective and only accepts better solutions.The algorithm consists of multiple cycles, each defined by a one-time cooling process to lower the temperature from an initial to a final value.While the algorithm is stochastic, it can be optimized through a local search approach to reduce the search space and find an optimal solution.Time for the procedure can further be optimized by including the annealing schedule from Fast Simulated Annealing (FSA) [42], which consists of semi-local searches with occasional long jumps.
Setting up of the numerics for SA: For both MIS and Max-Cut problems, the cost function was defined to provide the optimum solution at the function's minima.In this work for SA, the temperature is lowered from 0.4 to 0.01 under a distorted Cauchy-Lorentz distribution schedule(FSA) [42,84].During the simulations, a limit of 2000 function calls per iteration and 50 cycles was fixed.However, the maximum number of function calls was never reached and all cycles were completed.Each full run consists of multiple iterations, an optimum is reached at the end of each full run, and 50 such full runs were performed for statistics.The probability of success P success and run-time t run to reach the optimum is measured as a function of the number of iterations for a weighted problem graph to find an optimum number of iterations for a single run as shown in Fig S1 (a and b).5000 iterations were chosen for subsequent simulations, as it is enough to decrease the error in the statistics of P success and have a reasonable t run .After fixing the number of iterations, 15 unweighted graphs were chosen of a varying number of vertices from 2 to 16, to study the effect of system size in finding the optimum and scaling of run-time.As shown in Fig S1 (c and d), a general trend of increase in run-time and a decrease in the probability of success is observed for both problems.An oscillatory behavior in the probability of success is also observed with varying system sizes for even and odd numbers of vertices in the graph.This behavior is attributed to the degeneracy of the solution space.If the optimum is highly degenerate, there is an increase in the probability of reaching the optimum stochastically as the solution space is larger.

FIG. 1 :
FIG. 1: Setup of weighted Max-Cut and MIS problems.The figure in panel (a) is the problem graph for Max-Cut with weights w i j on the edges (such that ω 34 > ω 14 > ω 24 > ω 12 ), and the figure in panel (d) is the problem graph for MIS with weights w i on the vertices (such that ω 3 > ω 1 > ω 4 > ω 2 ).Panels (b) and (e) show the solution to the corresponding problems.The dashed curve in (b) indicates the cut, dividing the graph into two sets, red and blue vertices.In (e), red vertices constitute the MIS.Panels (c) and (f) correspond to an atomic configuration where the shaded blue region around each atom indicates the local detuning and the global Rabi frequency used in the setup.Each atom is subjected to specific values of detuning depending on the weights which are encoded in the distance-dependent interactions between the atoms as indicated with the solid-black arrow.Distance between atoms in the Rydberg configuration in the case of the Max-Cut problem follows r 34 < r 14 < r 24 < r 12 and in the case of the MIS problem follows r 34 < r 14 < r 12 < r 24 .

FIG. 2 :
FIG. 2: Left panels (a -e) are for the Max-Cut problem and right panels (f -j) are for the MIS problem.(a) and (f) show weighted prototype graphs of size 5 with corresponding solution graphs.Panels (b) and (g) show optimal protocols for the Rabi frequency (solid dark red) and the local detunings (dotted black symbols) with time.The local detuning of each atom is controlled by varying a single time-dependent parameter ∆ G (t) which is explained in the main text.Maximum and minimum speed in detuning change for the Max-Cut problem is 47.4MHz/µs and 15.4MHz/µs respectively, and for the MIS problem, it is 24.5MHz/µs and 4.7MHz/µs respectively.The expectation value E of the problem Hamiltonian with respect to the instantaneous state (solid blue) and the fidelity F (dashed green) of the instantaneous state with respect to the ground state are shown in panels (c) and (h), where R indicates the approximation ratio.(d) and (i) show the ordered energies of the instantaneous eigenstates, with a color bar indicating the population of the basis states during the protocol.The population of the basis states of the Hamiltonian at final time t = T is shown at three different times during the protocol in (e) and (j).As shown in the middle panels of (e) and (j), 12 states and 9 states out of 32 are populated in the middle of the protocol.The output at the end of the protocol captures all the degenerate states corresponding to all the degenerate Max-Cut/MIS solutions.

Fig. 3
Fig.3demonstrates the flexibility and scaling of our encoding for a more complex case involving a weighted graph of size 15 and degree 5. Panels (a) and (e) show the graph

FIG. 3 :
FIG.3: Solutions to the Max-Cut (a-d) and MIS (e-h) problems for a graph of size 15 and degree 5 using the optimal quantum annealing.Weighted prototype problem graphs with solution graphs for the Max-Cut problem are shown in (a) and are shown for the MIS problem in (e).In (a), vertex 11 has a degree of 5 while in (e), vertex 4 has a degree of 5. (b) and (f) show the optimal protocol for the Rabi frequency depicted by the solid dark red curve and ∆ G (t) (defined similar to Fig.2) depicted by the dashed black curve, with time.The maximum and minimum speeds for detuning change, for Max-Cut protocols, are 28.8MHz/µs and 5.1MHz/µs respectively, and, for MIS protocols are 1.8MHz/µs and 0.8MHz/µs respectively.The expectation value E (solid blue) of the problem Hamiltonian with respect to the instantaneous state and fidelity F (dashed green) of the instantaneous state with respect to the ground state is shown in (c) and (g), where R indicates the approximation ratio.The population of states at three different times (t = 3.50µs, t = 7.10µs, and t = 10.65µs) are shown in (d,h).

FIG. 4 :
FIG. 4: Comparison of optimized simulated annealing (SA) and Localised Optimal-control for Quantum Annealing in a Loop (LOQAL) for Max-Cut (a,b) and MIS (c,d) problem.Approximation ratio error 1 − R with respect to system size N is shown in (a,c), and for the hardness parameter HP in (b,d).

FIG. 5 :
FIG. 5:The figure shows the results using the optimal quantum annealing with noisy protocols to solve the Max-Cut problem for a graph of size 10.The weighted prototype graph for the Max-Cut problem is shown in (a) along with the solution graph.The optimal protocol for the Rabi frequency (solid dark red) and the factor controlling the detunings ∆ G (dashed black) are shown in (b) where noise is added to the laser parameters of the pre-optimized protocol and in (d) where noise is added during the optimization of the parameters.The shaded regions represent the fluctuations in the laser parameters for each run that were chosen from a random distribution.(c) and (e) show the corresponding expectation value (solid blue) of the problem Hamiltonian with respect to the instantaneous state (E) and the fidelity F (dashed green) of the instantaneous state with respect to the ground state for both cases, where R indicates the approximation ratio.
FIG. S1: The figure presents optimized simulated annealing for Max-Cut and MIS problems.Panel (a) shows the plot of the probability of failure (1 − P success ) and panel (b) shows run-time (t run ) in minutes as a function of the number of iterations N iterations (varying from 10 to 10000).In (a,b), a weighted problem graph with 10 vertices and 13 edges is chosen.(c) and (d) show the plot of the probability of failure 1 − P success and run-time t run in minutes as a function of the number of vertices N (varying from 2 to 16) for the unweighted problem graphs.