Hardness of the Maximum Independent Set Problem on Unit-Disk Graphs and Prospects for Quantum Speedups

Rydberg atom arrays are among the leading contenders for the demonstration of quantum speedups. Motivated by recent experiments with up to 289 qubits [Ebadi et al., Science 376, 1209 (2022)] we study the maximum independent set problem on unit-disk graphs with a broader range of classical solvers beyond the scope of the original paper. We carry out extensive numerical studies and assess problem hardness, using both exact and heuristic algorithms. We find that quasi-planar instances with Union-Jack-like connectivity can be solved to optimality for up to thousands of nodes within minutes, with both custom and generic commercial solvers on commodity hardware, without any instance-specific fine-tuning. We also perform a scaling analysis, showing that by relaxing the constraints on the classical simulated annealing algorithms considered in Ebadi et al., our implementation is competitive with the quantum algorithms. Conversely, instances with larger connectivity or less structure are shown to display a time-to-solution potentially orders of magnitudes larger. Based on these results we propose protocols to systematically tune problem hardness, motivating experiments with Rydberg atom arrays on instances orders of magnitude harder (for established classical solvers) than previously studied.

Rydberg atom arrays are among the leading contenders for the demonstration of quantum speedups.Motivated by recent experiments with up to 289 qubits [Ebadi et al., Science 376, 1209(2022)] we study the maximum independent set problem on unit-disk graphs with a broader range of classical solvers beyond the scope of the original paper.We carry out extensive numerical studies and assess problem hardness, using both exact and heuristic algorithms.We find that quasi-planar instances with Union-Jack-like connectivity can be solved to optimality for up to thousands of nodes within minutes, with both custom and generic commercial solvers on commodity hardware, without any instance-specific fine-tuning.We also perform a scaling analysis, showing that by relaxing the constraints on the classical simulated annealing algorithms considered in Ebadi et al., our implementation is competitive with the quantum algorithms.Conversely, instances with larger connectivity or less structure are shown to display a time-to-solution potentially orders of magnitudes larger.Based on these results we propose protocols to systematically tune problem hardness, motivating experiments with Rydberg atom arrays on instances orders of magnitude harder (for established classical solvers) than previously studied.

I. INTRODUCTION
Combinatorial optimization problems are pervasive across science and industry, with prominent applications in areas such as transportation and logistics, telecommunications, manufacturing, and finance.Given its potentially far-reaching impact, the demonstration of quantum speedups for practically relevant, computationally hard problems (such as combinatorial optimization problems) has emerged as one of the greatest milestones in quantum information science.
The physics of the Rydberg blockade mechanism has been shown to be intimately related to the canonical (NPhard) maximum independent set (MIS) problem [8], in particular for unit-disk graphs.The MIS problem involves finding the largest independent set of vertices in a graph, i.e., the largest subset of vertices such that no edges connect any pair in the set; compare Fig. 1 for a schematic illustration.As shown in Ref. [8], MIS problems can be encoded with (effectively two-level) Rydberg atoms placed at the vertices of the target (problem) graph.Strong Rydberg interactions between atoms then prevent two neighboring atoms from being simultaneously in the excited Rydberg state, provided they are within the Rydberg blockade radius, thereby effectively implementing the independence constraint underlying the MIS problem.By virtue of this Rydberg blockade mechanism, Rydberg atom arrays allow for a hardware-efficient encoding of the MIS problem on unitdisk graphs, with the (tunable) disk radius R b ∼ 1 -10µm setting the relevant length-scale [6].
Overview of main results.Recently, a potential (superlinear) quantum speedup over classical simulated annealing has been reported for the MIS problem [12], based on variational quantum algorithms run on Rydberg atom arrays with up to 289 qubits arranged in two spatial dimensions.This work focused on benchmarking quantum variational algorithms against simulated annealing by viewing it as a classical analogue of the adiabatic algorithm, yet left open the question of benchmarking against other state-of-the-art classical solvers.Motivated by this experiment, we perform a detailed analysis of the MIS problem on unit-disk graphs and assess problem hardness using both exact and heuristic methods.We provide a comprehensive algorithmic and numerical analysis, and we demonstrate the following: (i) Typical Schematic illustration of the problem.(a) We consider unit-disk graphs with nodes arranged on a twodimensional square lattice with lattice spacing a and filling fraction ϱ ∼ 80%, and edges connecting all pairs of nodes within a unit distance (illustrated by the circle).For √ 2a ≤ R b < 2a (as considered here), nodes are connected to nearest and next-nearest neighbors resulting in a (quasiplanar) Union-Jack pattern with maximum degree dmax = 8.(b) Our goal is to solve the MIS problem on this family of instances (as depicted here with nodes colored in red in the right panel) and assess the hardness thereof using both exact and heuristic algorithms.
quasi-planar instances with Union-Jack-like connectivity (as studied in Ref. [12]) can be solved to optimality for up to thousands of nodes within minutes, with both custom and generic commercial solvers on commodity hardware, without any instance-specific fine-tuning.(ii) Systematic scaling results are provided for all solvers, displaying qualitatively better runtime scaling for solvers exploiting the quasi-planar problem structure than generic ones.In particular, we find that by relaxing the detailed balance constraint, and considering the low depth regime (both of which are required for analytic runtime lower bounds on SA described in Ref. [13]) our implementation of classical simulated annealing is competitive with the quantum algorithm's performance in Ref. [12].(iii) Conversely, while the definition of problem hardness may be specific to the method used, instances with larger connectivity or less structure display a time-to-solution typically orders of magnitudes larger.(iv) Based on these results, we propose protocols to systematically tune problem hardness (as measured by classical time-to-solution), motivating experiments with Rydberg atom arrays on instances orders of magnitude harder (for established classical solvers) than previously studied.This paper is organized as follows.In Sec.II we first formalize the problem we consider.Next, in Sec.III we describe the algorithmic tool suite with which we address this problem.In Sec.IV we then describe our numerical experiments in detail.Finally, in Sec.V, we draw conclusions and give an outlook on future directions of research.

II. PROBLEM SPECIFICATION
The MIS problem is a prominent combinatorial optimization problem with practical applications in network design [24], vehicle routing [25], and finance [26,27], among others, and is closely related to the maximum clique, minimum vertex cover, and set packing problems [28].

A. Definition
Formally, the MIS problem reads as follows.Given an undirected graph G = (V, E), an independent set S ⊆ V is a subset of vertices of G such that no two vertices in S share an edge in E. The maximum independent set problem is then the task of finding the largest independent set in V.The cardinality of this largest independent set is referred to as the independence number |MIS|.One way to formulate the MIS problem mathematically is to first associate a binary variable x i ∈ {0, 1} to every vertex i ∈ V, such that x i = 1 if vertex i = 1, . . ., N belongs to the independent set, and x i = 0 otherwise.The MIS problem can then be expressed as a compact integer linear program of the form with the objective to maximize the marked vertices while adhering to the independence constraint.Generalizations to the maximum-weight independent set problem are straightforward [25].
A formulation of the MIS problem that is commonly used in the physics literature expresses the integer linear program in Eq. (1) in terms of a Hamiltonian that includes a (soft) penalty to non-independent configurations (i.e., when two vertices in the set are connected by an edge) [8].This Hamiltonian is given by with a negative sign in front of the first term because the largest independent set is searched for within a minimization problem, and where the penalty parameter V enforces the constraints.Energetically, this Hamiltonian favors having each variable in the state x i = 1 unless a pair of vertices are connected by an edge.For V > 1, the ground state is guaranteed to be a MIS, because it is strictly more favorable to have at most one vertex per edge in the set as opposed to both vertices being marked [12].Still, within this framework, the independence constraint typically needs to be enforced via post-processing routines (as is done, for example, in Ref. [12]).Mapping the binary variables x i to two-level Rydberg atoms subject to a coherent drive with Rabi frequency Ω and detuning ∆, one can then search for the ground state of the Hamiltonian H (encoding the MIS) via, for example, quantum-annealing-type approaches using quantum tunneling between different spin configurations [8,12].

B. Problem hardness
The MIS problem is known to be strongly NP-hard, making the existence of an efficient algorithm for finding the maximum independent set on generic graphs unlikely.As such, the MIS problem is even hard to approximate [29], and in general cannot be approximated to a constant factor in polynomial time (unless P = NP).
Here, however, we primarily focus on the MIS problem on unit-disk graphs (dubbed MIS-UD hereafter), given their intimate relation to Rydberg physics [8,12].Our main goal is to empirically assess the hardness of MIS-UD.As schematically depicted in Fig. 1, unit-disk graphs are defined by vertices on a two-dimensional plane with edges connecting all pairs of vertices within a unit distance.The MIS-UD problem appears in practical situations with geometric constraints such as map labeling [30] and wireless network design [24].While approximate solutions to MIS-UD can be found in polynomial time [31], solving the problem exactly is still known to be NP-hard for worst-case instances [12,32].

C. Problem instances and figures of merit
The unit disk (UD) problem instances of interest can be characterized by the number of nodes N , the side length of the underlying square lattice L, and the filling fraction ϱ, with N ≈ ϱL 2 , as schematically depicted in Fig. 1.Following Ref. [12], we focus on single-component (non-planar, but quasi-planar) UD instances with nearest and next-nearest (diagonal) couplings only, resulting in Union-Jack-type graphs with maximum degree d max = 8.Accordingly, these instances consist of (at maximum) n corner ≤ 4 corner, n boundary ≤ 4(L − 2) boundary, and n bulk ≤ (L − 2) 2 bulk nodes, with (at maximum) three, five, and eight neighbors, respectively, and a total of at most |E| max edges, with |E| max = 4L 2 − 6L + 2. Because |E| max ∼ N , the graph density D graph scales as showing that these instances become sparser as the system size N grows.If not otherwise specified, we take ϱ = 80%, as was done in Ref. [12].For comparison, we also run experiments on (unstructured) random Erdős-Rényi (ER) graphs denoted as G(n, m), chosen uniformly at random from the collection of all graphs with n nodes and m edges, or similarly G(n, p) for graphs constructed by connecting nodes randomly with probability p.
To assess and compare the performance of various algorithms (as specified below) we consider the following figures of merit.We use |MIS| to denote the indepen-dence number, while P MIS refers to the probability of observing an (exact) MIS within a fixed number of steps [12].For a given instance, many MIS solutions may be available, with the corresponding number of MIS solutions (i.e., the MIS degeneracy) denoted as D MIS .Similarly, the quantity D MIS−1 refers to the number of first excited states (i.e., independent sets of size |MIS| − 1).As shown in Ref. [12], in the context of simulated annealing, problem hardness may further be specified in terms of the conductance-like hardness parameter with the factor |MIS|D MIS denoting the number of possible transitions from a first excited state into a MIS ground state.Finally, we are interested in the time-tosolution (TTS).For exact methods, TTS refers to the time needed to find the optimal solution (i.e., ground state).While the optimum may be found after time TTS, additional time may be required to provide an optimality certificate, resulting in the time-to-optimality (TTO) time scale, with TTO ≥ TTS.By definition, provable optimality is not available with heuristic methods.Here, we define TTS 99 as the time required to find the exact solution (ground state) with 99% success probability.We can then write TTS 99 as where τ refers to the time of a single run (shot) and is the number of shots (repetitions) needed to reach the desired success probability [33].For small values of P MIS , we have R 99 ≈ 4.6/P MIS , showing that the success probability of a single run P MIS determines the inverse of the time-to-solution for heuristic algorithms.

III. ALGORITHMIC TOOL SUITE
In this section we detail the algorithms used to solve the MIS problem on UD graphs.We distinguish between exact methods (which by design can deterministically find the ground state, typically at the expense of an exponential runtime) and heuristics (which cannot provide an optimality certificate, but may require shorter runtimes).

A. Exact Methods
Sweeping line algorithm.We first consider an exact sweeping line algorithm (SLA) [34] that efficiently exploits the quasi-planar structure of the UD instances considered here.The anatomy of the SLA is schematically illustrated in Fig. 2. The SLA is based on full enumeration and works by sweeping a fictitious line across the two-dimensional plane (e.g., from left to right).Specifically, the algorithm proceeds as follows.We define the boundary as the set of all processed nodes which still share an edge with an unprocessed node.At each step, we track the size of the largest independent set (so far) for each valid boundary variant (set of assigned nodes on the boundary).As the line is swept across the graph, it stops at every node i = 0, . . ., N − 1.We then generate the new variants at step i from those at step i − 1 as follows: • If the variant has a neighboring node of i assigned, we generate only a new variant with i not selected (that is, x i = 0).
• Otherwise we also create a new variant with i assigned as x i = 1 (which increases its independent set).
Note that forward-looking information is not required, and only boundary nodes are relevant for this decision.
Once the new variants have been generated, node i becomes part of the boundary, and we proceed with the next step.When moved across the graph, this recipe generates all valid sets with runtime O(2 N ).However, we can efficiently summarize information that is not relevant to finding the MIS: • Adding new nodes typically removes older nodes from the boundary (because they no longer have a connection to an unprocessed node); c.f. Fig. 2.
• We only need to track the size of the largest independent set for each boundary variant (i.e., for any boundary configuration we can discard any option with equal or smaller number of assignments).
For the UD graphs considered here, the number of variants tracked on any given boundary is limited by the number of valid assignments on the boundary, #MIS.
When processing the nodes in order, the boundaries form continuous (mostly) one-dimensional strips and #MIS(L) ≤ Fib(L + 1) (which can be shown by induction).As a result, the memory requirement for SLA scales as O mem (Fib(L)) = O mem (Fib( √ N )) (to hold the variants at each step).SLA finds the optimal solution after N steps (TTS = TTO), processing all variants at each step with a runtime of , where φ ≈ 1.62 is the golden ratio.This procedure can also be modified to count the degeneracy of the ground and first excited state (i.e., D MIS and D MIS−1 , respectively) by adjusting summarization accordingly.
Branch & bound solvers.We complement our custom SLA solver with commercial solvers based on the branch and bound (B&B) search method.In particular, we use the solver offered by CPLEX [35]; similar results were observed with Gurobi [36].In practice, these solvers are among the de facto go-to tools for many hard, mixedinteger optimization problems.By design, B&B solvers provide upper and lower bounds on the solution, with the difference between these yielding an optimality gap, thereby giving information about the quality of the solution (at any step throughout the algorithmic evolution).Assuming a maximization problem, the lower bound corresponds to the best known feasible solution whereas the upper bound refers to the optimal value for the corresponding relaxed problem in the B&B procedure.In this work we focus on the TTS and TTO figures of merit, which are readily provided by our chosen solvers.We evaluate TTO by enforcing a zero gap between the upper and the lower bound, and TTS by setting the upper bound to be the optimal solution found previously.Thus, the solver terminates successfully as soon as it reaches the optimal solution.Typically, we find TTO is strictly greater than TTS because of additional time required to prove optimality.In order to draw a clear line between the B&B solvers and the heuristic solvers described below, we deactivate the B&B solvers' capability to find feasible solutions heuristically.As such, in practice we expect smaller values for TTS when utilizing these additional features of modern B&B solvers, effectively making the TTS values we report here upper bounds for B&B based performance.Finally, to account for the multithreading capabilities of B&B solvers, we report the process time, i.e., the sum of system and user CPU seconds of each core used during the calculation, rather than the wall-clock time.We checked that the multi-threading overhead does not impact the O complexities inferred from the numerical experiments; see Appendix A 2 for more details.

B. Heuristics
Apart from the exact methods outlined above, we utilize two established (physics-inspired) heuristic algorithms [37], namely simulated annealing (SA) and par-FIG.3: Schematic illustration of the heuristic simulated annealing (SA) solver.The original configuration (a) is overlaid with connectivity statistics in (b): nodes in the set (red), nodes without marked neighbors (white), and nodes with a single neighbor (blue).Potential moves are (i) removal of (red) nodes currently in the set, (ii) addition of currently white nodes to the set, and (iii) swapping a blue node with its red neighbor.Grey nodes have more than one adjacent node in the set and no valid moves.From (b) to (c), one node is added to the set and the statistics are updated accordingly.From (c) to (d) one blue node is swapped with its adjacent red node.
allel tempering (PT).For a schematic illustration see Fig. 3.In these Markov chain Monte Carlo (MCMC) samplers, a random modification to the current solution is proposed at each step of the algorithm and accepted depending on its effect on a specified figure of merit.For the MIS-UD problem at hand, we use the size of the independent set (IS) as the figure of merit to optimize.Any proposed move (update) is then accepted with a probability governed by a temperature parameter T > 0 according to the Metropolis acceptance criterion [38]: That is, moves that increase the size of the IS (i.e., ∆ IS > 0) are always accepted, while those reducing its size are suppressed -at first only marginally at high temperatures (during initial exploration), but then heavily at low temperatures T (during final exploitation).Our Markov dynamics consist of individual additions/removals and swaps of neighboring sites.We ensure the independence criterion is never violated by continuously tracking the full list of valid moves.The random selection from this set is then biased towards adding nodes and performing swaps in order to increase the acceptance rate (while accounting for the shift in energy scales by adapting the cooling schedule).
Simulated annealing (SA) aims to find a highquality solution by starting from a random initial solution and initially high temperature to then gradually lower T → 0 (according to some annealing schedule) until no further improvement is seen [39].This allows the system to first explore the solution space while the temperature is high, but eventually drives the state into a nearby (local, potentially global) optimum.The cooling schedule is optimized to quickly identify this local optimum, and frequent restarts (as specified by the parameter num_restarts) from different initial positions are used to increase the chance of finding the global optimum.We note that our implementation of SA differs from the one presented in Ref. [12] in several ways: (i) Our implementation is not based on a (soft) penalty model as described by Eq. ( 2), but rather involves only moves compatible with the (hard) independence criterion (such that only the feasible space is searched).(ii) Proposal probabilities are biased towards additions and exchanges to increase acceptance.(iii) We use a geometric cooling schedule in combination with frequent restarts, as opposed to a constant low temperature as used in Ref. [12].In particular, we note that our implementation of SA breaks detailed balance, for the sake of improved performance.
Parallel tempering (PT) attempts to efficiently explore the solutions space by simulating several Markov chains concurrently at different temperatures [40][41][42].Exchange moves allow for swapping of states between neighboring chains (in temperature space) such that the best solutions are shuffled to lower temperatures for further local optimization.At the same time, less promising candidates are moved to higher temperatures, where large-scale restructuring is possible.
In the following we will focus on SA, since our implementation of PT did not provide any substantial performance benefits over SA for the problem instances studied here.This is likely because of the relatively fast identification of local minima by SA (within tens to hundreds of sweeps) compared to the mixing time needed to exploit the benefits of PT; as such it is more efficient to restart at a random position than to invest in overcoming local energy barriers.

IV. NUMERICAL EXPERIMENTS
We now turn to our numerical results.We report on the TTS for all exact and heuristic algorithms described above, as a function of system size N and hardness parameter H, with the goal to provide a comprehensive assessment of the hardness of random MIS-UD problem instances.For reference, we also study the MIS problem on similar yet less structured instances (with the same number of nodes and edges), and we provide protocols to systematically tune problem hardness (as measured by TTS) over several orders of magnitude.The classical hardware on which our numerical experiments were run is specified in Appendix A. FIG.4: Time-to-solution (TTS) for the exact solvers.(a) TTS for the exact SLA solver as a function of system size N .For every system size N , 1000 random UD instances with ϱ = 0.8 have been considered.The data fit reasonably to TTS(N ) ≈ cN ϕ √ N , where the basis of ϕ ≈ 1.62 is the theoretical expectation for the Fibonacci sequence.At larger system sizes (N > 500), high memory usage causes slower access times (cache misses), resulting in a substantially larger pre-factor c ′ .(b) TTS for the B&B solver as a function of system size N .For every system size N , 1000 random UD instances have been considered; see Section A 3 for box plot description.Problems with hundreds (thousands) of nodes can be solved to optimality in sub-second (minute) timescales.The solid line is the linear regression over instances whose TTS are in the top highest 2%.The linear regression minimizes the residual sum of squares of log(TTS).

FIG. 5:
Time required to reach 99% success probability (TTS99) for the heuristic SA solver as a function of system size N (i.e., how long the solver should run for a 99% chance of finding the optimal solution).For every system size N , 1000 random UD instances at ϱ = 0.8 filling have been considered; see Section A 3 for box plot description.The solid line is the linear regression over instances whose TTS are in the top highest 2%.The linear regression minimizes the residual sum of squares of log(TTS).

A. Scaling with the Problem Size
We first report on TTS for the MIS-UD problem as a function of system size, given by the number of nodes N ≈ ϱL 2 at fixed density ϱ = 0.8.Note that we have run a few additional experiments for different values of ϱ, with ϱ ≈ 0.8 providing one of the hardest problems for Union-Jack connectivity (as evidenced by the largest TTS).Thus, we primarily focus on ϱ = 0.8, following Ref.[12].As shown in Fig. 4, we find that the MIS-UD problem can be solved to optimality (with both the custom SLA and the generic B&B solvers) for hundreds of nodes in sub-second timescales.Larger instances with up to thousands of nodes can still be solved to optimality within minutes on commodity hardware, without any instance-specific fine-tuning.For the exact SLA solver we infer a runtime scaling of TTS SLA = O(N ϕ √ N ), where ϕ ≈ 1.62 is the golden ratio (expected for the scaling of the Fibonacci sequence) and the notable √ N dependence in the exponent is attained at the expense of exponential memory requirements (as discussed above).For the B&B solver we obtain TTS B&B = O(2 0.0045N ).Similarly, as shown in Fig. 5, we find that UD instances with hundreds of nodes (i.e., L ≤ 25) can typically be solved efficiently in sub-second timescales with the SA heuristic.For the 2% most difficult instances we observe a scaling of TTS 99 = O(2 0.0128N ) for SA.However we also observe a relatively large spread spanning several orders of magnitude in TTS 99 (in particular when compared to the results obtained with SLA), thus motivating a more detailed analysis of problem hardness, as discussed next.

B. Scaling with the Hardness Parameter
Some of the results presented in the previous section display significant instance-to-instance variations in TTS 99 , potentially spanning several orders of magnitude, even for fixed system size N .As argued in Ref. [12], these variations may largely be due to (large) differences in the total number of MIS solutions, given by the ground state degeneracy D MIS , a quantity that can be calculated either with tensor-network methods [43] or within the exact SLA method outlined above (at least for small to intermediate system sizes up to N ∼ 1000).Intuitively, the less degenerate the ground state, the harder it is to hit the global optimum.In particular, based on experiments with up to N = 80 qubits in Ref. [12], a quantum speedup over classical SA has been reported in the dependence of the success probability P MIS on the hardness parameter H that accounts for both the degeneracy of the ground as well as first excited states (denoted by D MIS and D MIS−1 , respectively), as given in Eq. (3).
We now follow Ref. [12] and consider algorithmic performance in terms of this hardness parameter H for values as high as H ∼ 10 8 , complementing existing results based on classical SA [12] with results for the SLA, and B&B solvers.Note that hyperparameter optimization has been performed for our heuristic solvers, although without any instance-to-instance fine-tuning.For direct comparison, the results for the exact SLA as well as the heuristic SA solvers are displayed in Fig. 6 showing a remarkably different behavior.Qualitatively, we find that TTS 99 for the SA solver displays a strong dependence on the hardness parameter H, in line with results reported in Ref. [12].Conversely, virtually no dependence between TTS and H is observed for the exact SLA solver, as expected, thereby demonstrating that the conductance-like hardness parameter H successfully captures hardness for algorithms undergoing Markov-chain dynamics.Alternative algorithmic paradigms such as sweeping line or branch and bound, however, may require a different notion of hardness.Similarly, for the B&B solvers we find that TTS is (weakly) correlated with the hardness parameter H, although mostly because of their common correlation with the system size N .Specifically, using linear regression, we have found that the partial correlation of log 10 (H) and log 10 (TTS) (controlling for system for SA with a fixed depth of 32, for UD graphs selected from the top two percentile of hardness parameter H for each system size L = 13, . . ., 33.Hollow points represent our SA implementation with bias moves in valid configuration space.Power-law fits to the form ∼ H −α are used to extract scaling performance with graph hardness H. size) is smaller than 0.05; see Appendix B 2 for further details.This weak correlation suggests that (similarly to our SLA results) the hardness parameter H does not appear to be a reliable measure of hardness for B&B-type solvers.
Finally, we complement the TTS results above with results for the success probability P MIS as a function of the hardness parameter H for the SA solver, as done in Ref. [12] for both SA and quantum algorithms for instances with hardness of up to H ∼ 10 3 .Our results with hardness of up to H ∼ 10 8 are shown in Fig. 7. Following Ref. [12], for fixed depth (i.e., number of SA sweeps), fits are provided assuming the functional form P MIS = 1 − exp(−CH −α ), where C refers to a positive fitted constant that (in general) could have polynomial dependence on the system size N , and smaller values of α yield larger success rate P MIS .While there are significant variations in the data, on average we observe a scaling P MIS ≈ 1 − exp(−CH −0.66 ) (solid black line), i.e., α = 0.66; additional results with size-dependent depth (not shown) suggest an even higher success probability P MIS but preclude simple fits because of additional sizedependent effects in the data.In particular, a fixed depth of 32 is arguably too small for the largest instances considered here, but was chosen nevertheless for better comparison with results reported in Ref. [12].For comparison, a fit with the exponent α = 1.03 was reported for SA in Ref. [12].If one restricts the analysis to graphs with minimum energy gaps sufficiently large to be resolved in the duration of the (noisy) quantum evolution, the optimized quantum algorithm demonstrated in Ref. [12] was shown to fit best to α = 0.63, i.e., comparable to α = 0.66 as found with our implementation of classical SA.

C. Beyond Union-Jack Connectivity
To provide further context for the results reported above, we now study hardness as we gradually change the topology of the problem instances.Specifically, going beyond Union-Jack-type instances (with √ 2 ≤ r < 2 fixed) as studied so far, we analyze TTS following two protocols by either (i) systematically tuning the blockade radius r = R b /a or (ii) randomly rewiring edges of the graph.While protocol (i) prepares UD graphs (with varying connectivity), protocol (ii) explicitly breaks the UD structure via random (potentially long-range) interactions, ultimately preparing random (structure-less) ER graphs.The results of this analysis are shown in Figs. 8  and 9, respectively.We find that problem hardness (as measured here by TTS for the established B&B solver) can be tuned systematically over several orders of magnitude.
As shown in Fig. 8, we find that the MIS-UD problem is relatively easy for both small (r < 2) and large radii (as expected because the MIS problem is trivial for both edge-less and complete graphs), but significantly harder in between, with a pronounced peak at r = 3.For example, for L = 21 and ϱ = 0.8 (giving N ≈ 350) we observe an increase in TTS over 2 − 3 orders of magnitude for r = 3 instances compared to instances with small or large radii.This behavior appears to be generic, as we have observed similar behavior with the SLA solver (cf.Fig. 12), again showing pronounced peaks in TTS for r = 2, 3, 4; see Appendix B for more details.This observation may be attributed to a density-of-states-like effect: While the boundary size grows monotonically with r, we find that the number of boundary variants peak at r = 2, 3, 4 before decreasing with r as the (feasible) state space becomes smaller, thus correlating the observed TTS behav- ior with the number of boundary variants.This spread in TTS is found to increase with system size, as shown in Fig. 8 (b) for instances with Union-Jack topology (with r = √ 2) as well as instances with r = 3.Alternatively, following protocol (ii) the MIS problem appears to become orders of magnitude harder when randomly rewiring edges (thereby breaking the UD structure), as shown in Fig. 9.In particular, we find that structure-less ER graphs can yield a TTS orders of magnitudes larger than more structured UD graphs (with the same average number of nodes and edges), in agreement with similar results for random UD graphs where vertices are placed randomly (and not on a square lattice) in a two-dimensional box with some fixed density [8].

D. Implementation with Rydberg Arrays
The two protocols outlined above may be implemented in future experiments with Rydberg atom arrays, either by (i) tuning the Rydberg blockade radius R b and/or the lattice spacing a, or by (ii) implementing embedding schemes with ancilla Rydberg chains, as proposed in Refs.[44,45], thus suggesting a potential recipe to benchmark quantum algorithms on instances orders of magnitude harder (for established classical solvers) than previously studied.For example, for today's Rydberg atom arrays (such as the QuEra Aquila device available through Amazon Braket), we estimate that R b /a ∼ 3 (amounting to a maximum degree of d max ∼ 28) should be achievable already today, with potentially even larger values enabled by future hardware improvements.Experiments like these may also provide new insights into effects stemming from the long-range interaction tails associated with the Rydberg interactions.

E. Prospects for Quantum Speedups
With the goal to help identify regimes and system sizes where quantum algorithms could be useful, we now briefly revisit our results in the light of on-going efforts towards quantum advantage.Adopting the taxonomy put forward in Ref. [46], within a larger hierarchy of po-tential quantum speedups, the quantum speedup demonstrated in Ref. [12] could be classified as limited sequential quantum speedup, as it was obtained by comparing a quantum annealing type algorithm over a particular implementation of the classical (sequential) simulated annealing algorithm (that was designed to fulfill detailed balance).Here, we have tried to extend the classical SA benchmarking results (by pushing the hardness parameter up to H ∼ 10 8 ), with a different implementation of SA that breaks detailed balance but shows better performance.While our SA-based scaling results hint at performance similar to the quantum algorithm's performance in Ref. [12], we note that the corresponding exponent shows a dependence on the somewhat arbitrary cut-off in hardness and additional dependence on system size N if the depth is not fixed; the details thereof will be analyzed in future research.Still, within the aforementioned hierarchy of quantum speedups, our analysis points at a potential next milestone, in the form of the experimental demonstration of a (more general) limited non-tailored quantum speedup, by (for example) comparing the performance of the quantum algorithm to the best-known generic classical optimization algorithm.In particular, B&B solvers (as studied here) could be good candidates to fill the role for the latter, and the protocols outlined above motivate potential quantum experiments for such studies.Again, for reference for the MIS-UD problem on Union-Jack instances, empirically for most instances (taken as 98%-th percentile) we can upper-bound the TTS needed by classical B&B solvers through a runtime scaling of TTS = O(1.0031N ) (c.f.Fig. 4), setting an interesting, putative bar for quantum algorithms to beat.

V. CONCLUSION AND OUTLOOK
In summary, we have studied the maximum independent set problem on unit-disk graphs (as it can be encoded efficiently with Rydberg atom arrays [8,12]), using a plethora of exact and heuristic classical algorithms.We have assessed problem hardness, showing that instances with thousands of nodes can be solved to optimality within minutes using both custom and generic commercial solvers on commodity hardware, without any instance-specific fine-tuning.We have also performed a detailed scaling analysis, showing that our implementation of classical simulated annealing is competitive with the quantum algorithm's performance in Ref. [12].Finally, we have devised protocols to systematically tune problem hardness over several orders of magnitude.In future studies, the problem hardness may be tuned even further by generalizing our work from unweighted to weighted graphs, thereby potentially lifting the ground state degeneracy, as could be studied with local detuning control in Rydberg atom arrays.We hope that these protocols may trigger interesting future experiments further exploring the hardness of the MIS problem with Rydberg atom arrays as pioneered in Ref. [12].Thus, for large hardness we expect TTS to scale with the same exponent α as found from the fit shown in Fig. 7.We have further verified this result numerically by excluding small systems with small hardness.for the B&B and SLA solvers in Figs. 8 and 12, respectively.Again, this observation motivates future experiments with Rydberg atom arrays, as the likelihood for a potential quantum speedup may be larger on these instances, thus potentially helping to identify new regimes where quantum algorithms can be useful.

FIG. 1 :
FIG. 1:Schematic illustration of the problem.(a) We consider unit-disk graphs with nodes arranged on a twodimensional square lattice with lattice spacing a and filling fraction ϱ ∼ 80%, and edges connecting all pairs of nodes within a unit distance (illustrated by the circle).For √ 2a ≤ R b < 2a (as considered here), nodes are connected to nearest and next-nearest neighbors resulting in a (quasiplanar) Union-Jack pattern with maximum degree dmax = 8.(b) Our goal is to solve the MIS problem on this family of instances (as depicted here with nodes colored in red in the right panel) and assess the hardness thereof using both exact and heuristic algorithms.

FIG. 2 :
FIG. 2: Schematic illustration of the (exact) sweeping line algorithm (SLA) as applied to the MIS-UD problem.(a) SLA proceeds by sweeping a fictitious line across the graph and tracking all potentially optimal MIS configurations on this boundary, efficiently exploiting the quasi-planar structure of the UD graph.Processed nodes are shown in gray, boundary nodes in blue and unprocessed nodes in yellow.(b) The light blue node from (a) is added to the boundary, while the bottom left blue node is dropped (as a result of not having any more connections to unprocessed nodes).
solved [seconds] Number of Nodes in the Graph N Simulated Annealing Top 2% TTS = O(2 0.0128N )

33 FIG. 6 :
FIG.6:(a) Time-to-solution (TTS) for the exact SLA solver as a function of the hardness parameter H. Virtually no dependence on H is observed, showing that TTS is fully determined by the system size ∼ L 2 .(b) Conversely, for the Markov-chain based SA solver, TTS99 shows a strong correlation with the hardness parameter H, as expected.

33 FIG. 7 :
FIG.7:Estimated success probability PMIS for the heuristic SA solver as a function of the hardness parameter H.Here we plot − log(1 − PMIS) for SA with a fixed depth of 32, for UD graphs selected from the top two percentile of hardness parameter H for each system size L = 13, . . ., 33.Hollow points represent our SA implementation with bias moves in valid configuration space.Power-law fits to the form ∼ H −α are used to extract scaling performance with graph hardness H.

FIG. 8 :
FIG.8:(a) Hardness transition as a function of the disk radius (in units of the lattice spacing) r = R b /a, as given by the time-to-solution (TTS) for the B&B solver, shown here for system size L = 21 and density = 0.8 (i.e., N ≈ 350), with 100 random seeds per radius.(b): TTS as a function of system size N = ϱL 2 for r = √ 2 (blue) and r = 3 (green), the latter referring to the pronounced peak observed in (a).The solid lines are a linear regression fit over 100 instances with TTS in the highest 10%, with corresponding R 2 values of 0.87 and 0.98 for r = √ 2 and r = 3, respectively.Instances with r = 3 appear to be much harder than those with r = √ 2. See Section A 3 for box plot description.

FIG. 10 :
FIG. 10: Coefficient α extracted from time-to-solution (TTS) scaling with system size as TTS = O(2 αN ) for the B&B solver as a function of the lattice filling ϱ.In the main text we have focused on ϱ = 0.8.All results have been taken from the top 2% TTS instances over 1000 random instances.Error bars correspond to the 95% confidence interval.
FIG.12: Time-to-solution (TTS) for the SLA solver as a function of the disk radius r for random UD instances, for system size L = 21 and density ϱ = 0.8 (i.e., N ≈ 350), with 100 random seeds per radius.Similar to our results for our B&B solver, we observe distinct peaks at r = 2, 3, 4.
While the original UD graphs can be solved to optimality in ∼ 10 −2 s, comparable ER graphs (with the same number of nodes and edges) display a TTS orders of magnitudes larger.The red line and the two shaded areas refer to the median TTS over 500 instances for MIS on random ER graphs, the TTS among the 25% and 75% and the minimum and maximum whiskers, respectively.Numerical parameters: L = 21, and ϱ = 80%.