Quantum Algorithms for Jet Clustering

Identifying jets formed in high-energy particle collisions requires solving optimization problems over potentially large numbers of final-state particles. In this work, we consider the possibility of using quantum computers to speed up jet clustering algorithms. Focusing on the case of electron-positron collisions, we consider a well-known event shape called thrust whose optimum corresponds to the most jet-like separating plane among a set of particles, thereby defining two hemisphere jets. We show how to formulate thrust both as a quantum annealing problem and as a Grover search problem. A key component of our analysis is the consideration of realistic models for interfacing classical data with a quantum algorithm. With a sequential computing model, we show how to speed up the well-known O(N^3) classical algorithm to an O(N^2) quantum algorithm, including the O(N) overhead of loading classical data from N final-state particles. Along the way, we also identify a way to speed up the classical algorithm to O(N^2 log N) using a sorting strategy inspired by the SISCone jet algorithm, which has no natural quantum counterpart. With a parallel computing model, we achieve O(N log N) scaling in both the classical and quantum cases. Finally, we consider the generalization of these quantum methods to other jet algorithms more closely related to those used for proton-proton collisions at the Large Hadron Collider.

Identifying jets formed in high-energy particle collisions requires solving optimization problems over potentially large numbers of final-state particles. In this work, we consider the possibility of using quantum computers to speed up jet clustering algorithms. Focusing on the case of electronpositron collisions, we consider a well-known event shape called thrust whose optimum corresponds to the most jet-like separating plane among a set of particles, thereby defining two hemisphere jets. We show how to formulate thrust both as a quantum annealing problem and as a Grover search problem. A key component of our analysis is the consideration of realistic models for interfacing classical data with a quantum algorithm. With a sequential computing model, we show how to speed up the well-known O(N 3 ) classical algorithm to an O(N 2 ) quantum algorithm, including the O(N ) overhead of loading classical data from N final-state particles. Along the way, we also identify a way to speed up the classical algorithm to O(N 2 log N ) using a sorting strategy inspired by the SISCone jet algorithm, which has no natural quantum counterpart. With a parallel computing model, we achieve O(N log N ) scaling in both the classical and quantum cases. Finally, we consider the generalization of these quantum methods to other jet algorithms more closely related to those used for proton-proton collisions at the Large Hadron Collider.

I. INTRODUCTION
Jets are collections of collimated, energetic hadrons formed in high-energy particle collisions. With an appropriate choice of jet clustering algorithm [1], jets are a robust probe of quantum chromodynamics (QCD) and a useful proxy for determining the kinematics of the underlying hard scattering process. The problem of identifying jets from collision data is a nontrivial task, however, since the jet clustering algorithm must be matched to the physics question of interest. Moreover, it is a computationally intensive task, as it often involves performing optimizations over potentially large numbers of final-state particles.
In this paper, we consider the possibility of using quantum computers to speed up jet identification. We focus on the well-known problem of partitioning an electronpositron collision event into two hemisphere jets, though our results are relevant for other optimization problems beyond high-energy physics. Our main results are summarized in Table I, where the computational scaling is given for N particles in the final state. We show how to improve the well-known O(N 3 ) classical algorithm [2] to an O(N 2 ) quantum algorithm, which includes the cost of loading the classical data into a sequential quantum computing architecture. On the other hand, we also show how to speed up the classical algorithm to O(N 2 log N ), using a clever sorting strategy from Ref. [3], which matches the quantum performance up to log N factors. Finally, using parallel computing architectures, we achieve O(N log N ) scaling in both the classical and quantum cases, albeit for very different computational reasons.
Quantum algorithms have been shown to achieve speedups over classical algorithms [4], resulting, in theory, in time savings which are even more pronounced over large data sets. That said, many proposed quantum al-

Implementation
Time Classical [2] O(N 3 ) -Sec. III A Classical with Sort (using [3]) gorithms for machine learning tasks often omit considerations that would be needed to actually implement them in practice, such as a strategy to interface classical data with a quantum computing architecture. One solution is to assume the availability of qRAM [5], which would let our quantum computer access a classical data set in superposition; however this additional hardware requirement may not be easy to implement in practice. Here, we consider realistic applications of both quantum annealing [6][7][8] and Grover search [9][10][11] to jet finding, including the O(N ) overhead of loading classical collision data into the quantum computer.
The specific jet finding algorithm we use is based on thrust [12][13][14]. Thrust is an event shape widely measured in electron-positron collisions [15][16][17][18][19][20][21][22][23][24][25][26][27]. The optimum value of thrust defines the most jet-like separating plane among a set of final-state particles, thereby partitioning the event into two hemisphere jets. Algorithmically, it poses an interesting problem because it can be viewed in various equivalent ways-such as a partitioning problem or as an axis-finding problem-which in turn lead to different algorithmic strategies.
We note that practical thrust computations typically involve only 10-1000 particles per event, so the current O(N 3 ) classical algorithm [2] is certainly adequate to the task. That said, more efficient jet algorithms are of general interest, for example in the context of active area calculations [28], which can involve up to millions of ghost particles. We also note that the current default jet algorithm at the Large Hadron Collider (LHC) is antik t [29], which already runs in O(N log N ) time [30,31], and it is unlikely that any quantum algorithm can yield a sublinear improvement. On the other hand, anti-k t is a hierarchical clustering algorithm (i.e. a heuristic), whereas thrust is a global optimization problem, and there are phenomenological contexts where global jet optimization could potentially yield superior physics performance [32,33]; see also Refs. [34][35][36][37][38][39][40][41][42][43][44][45][46][47]. Jet finding via global optimization has not seen widespread adoption, in part because of the computational overhead, and we hope the quantum and improved classical algorithms developed here spur more research on alternative jet finding strategies.
Beyond the specific applications to jet finding, we be-lieve that the broader question of identifying realistic quantum algorithms for optimization problems should be of interest to both the particle physics and quantum computing communities. Indeed, we regard thrust as a warmup problem for the more general development of quantum algorithms for collider data analysis. (For other quantum algorithms for collider physics, see Refs. [48,49] for Higgs boson identification, Refs. [50,51] for parton shower generation, and Refs. [52,53] for track reconstruction.) Because collider data is classical (and will likely remain so for the foreseeable future), understanding the limitations imposed by data loading is essential to evaluate the potential of quantum algorithms to speed up or improve data analysis pipelines. At the same time, it is important to assess potential classical improvements to existing collider algorithms, and the sorting strategy of Ref. [3] is an important example of how new classical strategies can sometimes match the gains from quantum computation.
Turning now to an extended outline of this paper, our quantum algorithms build on existing classical strategies to compute thrust. In Sec. II, we define thrust in its various equivalent manifestations, as both a partitioning problem and an axis-finding problem. Then in Sec. III, we review classical algorithms for computing thrust based on a search over reference axes. As already mentioned, the best known result in the literature requires O(N 3 ) time [2]. We show how to improve it to O(N 2 log N ) using a sorting strategy inspired by SISCone [3], which appears to have no quantum analog (see Sec. V D).
The first quantum method we consider in Sec. IV involves formulating thrust as an quadratic unconstrained binary optimization (QUBO) problem, which can then be solved via quantum annealing [6,7]. This comes from viewing thrust as a partitioning problem and then considering the brute force enumeration of all candidate partitions. See Refs. [54,55] for other studies of quantum annealing for clustering with unique assignment.
The core results of this paper are in Sec. V, where we describe quantum algorithms for computing thrust based on Grover search [9]. Although naively Grover search offers a square root speedup over any classical search algorithm, in practice Grover search cannot yield sublinear algorithms. The reason is that data loading over a classical database of size N requires O (N )  The precise speed up achievable in our Grover search strategy depends on the assumed quantum computing paradigm. We implement two models for retrieving and processing the classical data, based on the abstract operations LOOKUP and SUM. The sequential computing model requires O(1) qubits and results in an O(N 2 ) thrust algorithm. Here we use O(·) to mean that we neglect factors that are polylog in N . The parallel computing model requires O(N ) qubits and results in an O(N log N ) thrust algorithm. Both computing models are applicable to any general problem where the size of the search space scales like O(N α ) with α ≥ 2, which are precisely the problems that can typically be sped up with a realistic application of Grover search.
In Sec. VI, we assess whether or not there is any quantum advantage for hemisphere jet finding. Formally, if one has read access to O(N ) classical bits but only write access to O(log N ) bits, then one cannot implement the classical sorting strategy in Sec. III C. In that case, there is a quantum advantage for both sequential and parallel computing models. With write access to O(N log N ) classical bits, though, classical sorting is possible, and the asymptotic performance of our classical and quantum algorithms are identical (up to log N factors) in both the sequential and parallel case. This equivalence appears to be special to algorithms like thrust where the search space scales like O(N 2 ), and we speculate that larger search spaces might benefit from Grover speedups even if classical sorting is possible.
Finally in Sec. VII, we briefly consider generalizations of our results to jet algorithms more closely related to those used at the LHC. We consider jet function maximization [43][44][45], showing that, with suitable modifications, it can be written in QUBO form for quantum annealing. We consider stable cone finding in the spirit of SISCone [3], showing how a single-jet variant we dub SingleCone is amenable to quantum search. We also comment on quantum multi-jet finding motivated by the XCone algorithm [32,33]. We conclude in Sec. VIII with some broader lessons about quantum algorithms for collider physics.

II. DEFINITION OF THRUST
We start by defining thrust [12][13][14], noting that it has multiple equivalent definitions that suggest different algorithmic strategies, as shown in Fig. 1. Thrust can be viewed as a partitioning problem, which lends itself naturally to quantum annealing. Thrust can alternatively be viewed as an axis-finding problem, which we can frame as a quantum search problem. Both definitions of thrust can be stated in terms of operator norms, and through this lens, they are in fact dual to each other.

A. Thrust as a Partitioning Problem
An intuitive formulation of thrust (though not exactly the original one [12,13]) is to separate the particles into a partition H L ∪ H R such that momenta on each side are as "pencil-like" as possible. That is, we seek to maximize the quantity where the second equality follows from momentum conservation. The quantity known as "thrust" corresponds to the maximum obtainable value: The factor of 2 in Eq. (2) is conventional such that 1/2 ≤ T ≤ 1, where T = 1 corresponds to a perfectly pencil-like back-to-back configuration and T = 1/2 is an isotropic event.
There is an equivalent geometric formulation of Eq. (2) due to Ref. [56]. Consider sequentially summing the three-momenta { p i } to form a closed polygon. Each sequence yields a different polygon, and computing thrust is equivalent to maximizing twice the diagonal of the polygon over all possible polygons, normalized by the circumference of the polygon. The diagonal splits the polygon into two halves, which yield the partition H L ∪ H R . The particles in H L are said to be in the "left hemisphere jet" and the particles in H R are said to be in the "right hemisphere jet".
This definition immediately suggests a naive, bruteforce classical strategy for computing thrust. We can enumerate all O(2 N ) possible partitions (which can be reduced to O(2 N −1 ) using momentum conservation), and then we sum the momenta in each to determine the maximum, resulting in an O(N 2 N ) algorithm. This is the version of thrust we will use for the quantum annealing formulation in Sec. IV, which corresponds to attacking the problem using quantum brute force.

B. Thrust as an Axis-Finding Problem
An alternative definition of thrust is as an axis-finding problem, which is a bit closer to the historical definition [12,13]. Letn be a unit norm vector and define Thrust can then be determined by the maximum value of T (n) overn: The optimaln is known as the thrust axis: To gain some intuition for why Eqs. (3) and (5) are equivalent, note that once we find the thrust axisn opt , we can partition the particles into those withn opt · p i > 0 and those withn opt · p i < 0. (It is an interesting bit of computational geometry to show thatn opt · p i can never be exactly zero for a finite number of particles.) Said another way, the plane normal ton opt partitions the event into left and right hemispheres. Any partitioning other than hemisphere partitioning will decrease the value of thrust in Eq. (2), so the optimal partition will be defined by a plane. Because of this equivalence between axis finding and plane partitioning, the thrust objective is sometimes written as where the Heaviside theta function picks out particles in just one hemisphere. Note that the optimal partitioning plane is not unique, since there can be multiple planes that yield the same partition. We can exploit this fact to find a computationally convenient partitioning plane, defined by a normal reference axisr. This reference axis will in general be different from the thrust axisn opt but nevertheless yield the same value of thrust via Eq. (2). Specifically, once the optimal partition is known via a reference axis, the thrust axis can be determined from the total three-momentum in the left hemisphere: We will use this reference axis approach for the classical thrust algorithms in Sec. III and for the quantum search strategies in Sec. V.

C. Duality of Thrust Definitions
Using the formalism of operator norms, we can show that these two definitions of thrust are in fact dual to each other.
Let M : V → W be a map from V = R m with norm · α to W = R n with norm · β . The operator norm of M , known as the induced α-to-β norm, is defined as That is, we search over all vectors v in V with norm 1 and find the maximum norm for the vector M v in W.
The case when α and β are both the usual L 2 norm corresponds to the largest singular value of M , but in general M α→β can be NP-hard to estimate [57]. By duality, we can rewrite this as where y is in W * , the vector space dual to W, defined as W * = R n with dual norm · β * . Thus, the α-to-β norm of M is the same as the β * -to-α norm of M T .
In the context of thrust, we are interested in the following norms for a vector v ∈ R n : These are known, respectively, as the 1-norm, 2-norm, and sup-norm. By Hölder's inequality, the space of vectors endowed with the p-norm is dual to the space of vectors endowed with the q-norm, where 1 p + 1 q = 1. In particular, the 1-norm is dual to the sup-norm, and the 2-norm is dual to itself. Now consider the matrix M ij = ( p i ) j , whose rows are the N three-momenta and whose columns are the p x , p y , and p z components. This is a map from R 3 to R N . Letting α = 2 and β = 1, the induced 2-to-1 norm of M is We recognize the last term as the numerator of T (n) in Eq. (4). Since the denominator of T (n) is independent of n, this is equivalent to the definition of thrust via axisfinding in Sec. II B. Thus, thrust takes the form of an induced 2-to-1 norm problem. By duality, with β * = ∞, thrust can alternatively viewed as an sup-to-2 norm problem: (15) This corresponds to the definition of thrust via partitioning in Sec. II A, since setting s i = −1 denotes flipping the orientation of vector p i relative to the partitioning plane, while setting s i = 1 retains the orientation of p i . Therefore, we see that the problem of computing thrust in particle physics is in fact a special instance of the more general problem of computing induced matrix norms. While there exist choices of α and β for which efficient algorithms for computing M α→β exist for arbitrary M , it is believed that the general problem of computing the induced 2-to-1 norm and that of computing the induced ∞-to-2 norm are both NP-hard [58][59][60]. This suggests that thrust is an excellent testbed to explore possible gains from quantum computation.

D. Alternative Duality Derivation
There is alternative language to understand this thrust duality that will be useful for the generalizations in Sec. VII. This approach is based on Ref. [61], which showed that different jet finding strategies can sometimes be derived from a common meta-optimization problem.
Consider a partition H (not necessarily defined by a plane) with total three-momentum Our analysis is based on the following objective function that depends on both a choice of partition and a choice of axis: where λ is a Lagrange multiplier to enforce that the axiŝ n has unit norm. At this point, P andn are completely independent entities, andn does not play any role in determining the partition H.
For fixed P , we can optimize O( P ,n) overn: Plugging this into Eq. (17) yields which is (half) of the thrust numerator in Eq. (2). For fixedn, we can optimize O( P ,n) over P (or equivalently, over the partition H): Plugging this into Eq. (17) yields which is (half) of the thrust numerator in Eq. (7). Since the order of optimization is irrelevant to the final optimum, this again shows that the two thrust definitions are dual. Either way, the maximum value of the objective function will be: which, following Ref. [56], is just the maximum achievable polygon diagonal.

III. CLASSICAL ALGORITHMS
We now describe the best known classical algorithm for thrust in the literature, which requires O(N 3 ) time, and then show how it can be improved to O(N 2 log N ) using a sorting technique from Ref. [3]. We start by assuming a sequential classical computing model in this section, and end with a brief discussion of parallel classical computing.

A. Plane Partitioning via a Reference Axis
The best known classical thrust algorithm [2] uses the reference axis approach discussed at the end of Sec. II B. 1 This is the thrust algorithm implemented in Pythia as of version 8 [62]. The key realization is that, because of Eq. (8), one only needs to search over inequivalent plane partitions. Two particles are sufficient to determine a separating plane, so there are O(N 2 ) inequivalent plane partitions to consider. For each partition, determining More specifically, for each pair of particles p i and p j , one determines a reference axisr ij normal to the plane spanned by them:r Then, each particle p k is either assigned to the hemisphere H ij ifr ij · p k > 0 or ignored ifr ij · p k < 0. Cases wherer ij · p k = 0 are ambiguous, and we provide a general strategy to deal with this in Sec. III B below. At minimum, we have to treat the cases where k = i or j, which requires testing 2 × 2 = 4 possibilities for whether or not p i and/or p j should be included in H ij , for a total of 4N (N − 1) partitions. (This can be reduced by a factor of 2 using momentum conservation, sincer ij andr ji define the same hemispheres.) The final hemisphere jets are determined by the partition that maximizes Note that, in general, none of the O(N 2 ) reference axes considered will align with the actual thrust axis. Nevertheless, the partitions defined byr opt andn opt will be identical. (In the idealized case of infinitesimal radiation everywhere in the event, all possible separating planes would be considered, sor opt would then equal n opt .) Once the optimal partition is known, the thrust axis itself is determined by Eq. (8).
In terms of computational complexity, for a fixed hemisphere H ij , it takes O(N ) time to compute the hemisphere three-momentum in Eq. (16). The thrust denom- can be precomputed since it is independent of the partition. Once P ij and T denom are known, though, it only takes O(1) time to determine the value of T ij : where we used Eq. for each partition, leading to the O(N 3 ) scaling. In Sec. III C, we can improve on this runtime by iteratively updating P ij in a special order.

B. Doubling Trick
To simplify the thrust algorithm, it is convenient to artificially double the number of particles. Starting from N three-momenta, we create a list of length 2N by including both p k and its negative − p k . Because p k and − p k can never be in the same hemisphere, and because of the momentum conservation relation in Eq. (2), this doubling trick has no effect on the value of thrust. It does, however, provide us with a convenient way to deal with the four-fold ambiguity above, since we can now always define the hemisphere H ij to include particle i.
To deal with cases wherer ij · p k = 0 (i.e. any time three or more particles are coplanar), we offset the reference axis byr and then take the formal → 0 limit. Specifically, if where the factor of 1 2 compensates for the artificial doubling.
We will use this doubling trick repeatedly in this paper, though not for quantum annealing in Sec. IV where it is counter-productive. To simplify the description of the algorithms, we will leave implicit the treatment of all r ij · p k = 0 cases via Eq. (26). It is worth mentioning that an alternative way to deal with coplanar configurations is to offset the momenta by a small random amount, but we find the doubling trick to be more convenient in practice since it avoids the four-fold ambiguity automatically.

C. Improvements via Sort
The O(N 3 ) algorithm can be further improved to run in time O(N 2 log N ). 2 This can be achieved using a strategy from SISCone proton-proton collisions, whereas our interest here is in electron-positron collisions, but the same basic strategy still applies.
The goal of SISCone is to find conical jet configurations J where the enclosed particles are within a distance R from the cone axisn J . Moreover, these cone jets must be stable, meaning that the jet three-momentum P J = i∈J p i is aligned with the cone axisn J . Like thrust, SISCone involves solving a partitioning problem where the naive brute-force approach requires O(N 2 N ) time. Like for thrust, one can reduce the naive runtime for SISCone to O(N 3 ) using the fact that two points lying on the circumference of a circle are sufficient to determine the cone constituents. There is an eight-fold ambiguity in the cone assignments, which we discuss further in Sec. VII C.
The key insight of Ref. [3] is that one need not recompute P J for all O(N 2 ) candidate cones. Ignoring the eight-fold ambiguity, let the candidate cones be labeled by i and j. For fixed i, one can define a special traversal order for j such that only one particle enters or leaves the cone at a time. There are N particles labeled by i, and for fixed i, sorting over j takes O(N log N ) time. After the initial O(N ) determination of P J for the first j values in the sorted list, updating the value of P J for each j iteration only requires O(1) time, since you need only add the momentum of a point entering the cone or subtract the momentum of a point leaving the cone. Thus, the final algorithm is O(N 2 log N ).
We can apply exactly the same sorting strategy to the computation of thrust, as shown in Fig. 2. The reason is that the reference axisr ij depends only on the cross prod-uct p i × p j . This means that for fixed p i , we can choose an ordering of the p j such that the partitions induced by {r i1 , . . . ,r iN } are specified by a single sortable parameter. To see this, it is convenient to transform to a coordinate system where p i points in the z direction, i.e. p i = |p| i (0, 0, 1). For any j = i, we can write p j in spherical coordinates as p j = | p j | (sin θ j cos φ j , sin θ j sin φ j , cos θ j ), where θ j is the polar angle and φ j is the azimuthal angle. Thenr ij = (− sin φ j , cos φ j , 0), so the partition is indeed determined by the single parameter φ j , independent of θ j . Specifically, particle k is in hemisphere H ij if is positive. This implies that 0 < φ k − φ j < π, where azimuthal angle differences are calculated modulo 2π. Furthermore, because of the doubling trick, there is a simple way to determine which particles are in the partition. With the doubling, there are 2N possible choices for i, and by Eq. (26) we know that the doubler − p i cannot be in the same partition as p i . Using the above coordinate system, we can sort the remaining 2N − 2 vectors according to their φ coordinates, i.e. so that 0 ≤ φ j1 ≤ φ j2 ≤ · · · ≤ φ j 2N −2 < 2π. (In cases where two particles happen to have identical values of φ j , their relative ordering does not matter for the argument below, as long as the doublers are also put in the same order.) Crucially, for a particle at position a in this sorted list, its doubler (which is π away in azimuth) must be at position a + N − 1. To see why, note that a hemisphere either contains a particle or contains its doubler, so there must be exactly N particles in each hemisphere. Particle i is already accounted for, meaning that any candidate partition must contain N − 1 entries from the sorted list. Since the sorted list is ordered by azimuth, and since the partitioning is determined by azimuth alone via Eq. (28), the N − 1 elements from position a to position a + N − 2 inclusive must be in a common partition, and the doubler must be the next one on the list. Therefore, candidate thrust partitions always take the form: (Because we are searching over all partitions, we can ignore the subtlety from Eq. (26) and just assume that j a is always contained in H i,ja .) These observations allow us to construct an O(N 2 log N ) algorithm for thrust. The outer loop involves iterating over all 2N choices for i. The inner loop involves the following O(N log N ) algorithm. We perform the sorting procedure above for fixed i, which takes O (N log N ) time. For the first element in the sorted list, we determine the partition H i,j1 using Eq. (29) with a = 1. We can readily compute P i,j1 via Eq. (27) in time O(N ), and then compute the associated thrust value via Eq. (25) in O(1). For the subsequent 2N − 3 elements of the sorted list, we step through them one by one, updating the partition from H i,ja = {i, j a , j a+1 , . . . , j a+N −2 } to H i,ja+1 = {i, j a+1 , j a+2 , . . . , j a+N −1 }. In doing so, we need to subtract p ja and add p j a+N −1 (which is the same as − p ja by the doubling trick), leading to the update step: where one has to remember the factor of 1 2 in Eq. (27). From the updated momentum, we recompute the associated thrust value via Eq. (25) in O(1) time. The total time from stepping through the 2N − 3 partition momenta is O(N ), so the inner loop is dominated just by the initial O(N log N ) sorting step. The maximum T ij over all i and sorted j determines the final hemisphere jets.

D. Parallel Classical Algorithm
The sorting algorithm above requires O(N 2 log N ) operations. In a model with a single CPU and randomaccess memory, this corresponds to time O(N 2 log N ) as well. We can also consider parallel computing models in which the N words of memory are accompanied by N parallel processors; see Ref. [63] for more discussion of these models. In this case, we will see that a runtime of O(N log N ) can be achieved. For simplicity, we do not consider the general case in which the number of parallel CPUs and the amount of memory can be varied independently, nor will we discuss the varying models of parallel computing in Ref. [63].
We briefly sketch here how the sorting strategy in Sec. III C can be sped up with parallel processors. There are three main computational bottlenecks: iterating over all particles i (C iter ), sorting over particles j for fixed i (C sort ), and determining the hemisphere constituents over each j for fixed i (C hemi ), leading to a runtime of O C iter (C sort + C hemi ) . For sequential classical computing, we found C iter = O(N ), C sort = O(N log N ), and C hemi = O(N ). A parallel computer cannot improve on C iter , but there are parallel computing algorithms for sorting [64] and partial sums [65] that would allow us to achieve C sort = C hemi = O(log N ), leading to a O(N log N ) runtime. We will compare the quantum and classical parallel architectures in Sec. VI.

IV. THRUST VIA QUANTUM ANNEALING
The first quantum algorithm we describe is based on quantum annealing [6,7]. In a quantum annealer such as the D-Wave system [8], the solution to an optimization problem is encoded in the ground state of a target Hamiltonian. Such a Hamiltonian takes the form of an Ising model: where each of the N Ising spins s i ∈ {−1, +1} corresponds to a qubit, and the {h i } and {J ij } correspond to programmable weights and couplings between qubits, respectively.
Equivalently, under the transformation s i = 2x i −1, we can frame the optimization problem as a QUBO problem, where the objective function takes the form for x i ∈ {0, 1}. Note that the fact that i, j are now summed with repeated indices and the fact that x 2 i = x i allow us to absorb the linear terms into the quadratic terms.
For the thrust problem, it is convenient to first define the three-momentum of a candidate partition as where x i = 1 if particle p i is in the partition and x i = 0 otherwise. Following Eq. (25), the thrust of this partition is given by Because of the square root factor, this is not a QUBO problem, but since the optimal partition is the same for any monotonic rescaling of T ({x i }), we can optimize the squared relation: which now takes the form of the QUBO problem in Eq. (32), as desired. Finding the ground state of −T ({x i }) 2 (note the minus sign) is the same as determining thrust.
The space usage of a quantum annealing algorithm is O(N ), corresponding to one qubit for each x i . The annealing time required depends on the spectral gap of the particular Hamiltonian, and we leave the question of determining the spectral gap of the thrust objective function to future work.

V. THRUST VIA QUANTUM SEARCH
We now describe a quantum algorithm for thrust based on Grover search. We first describe the algorithm in terms of two abstract operations, LOOKUP and SUM, both of which perform data loading in superposition. Then, we describe two computing models for loading the classical data into quantum memory: the sequential

A. Algorithm Overview
Our quantum thrust algorithm is based on the quantum maximum finding algorithm of Dürr and Høyer [11], which returns the maximum element of an unsorted array with K elements in O( √ K) time, assuming quantum query access to the array. This algorithm is itself a generalization of Grover search [9].
In this context, quantum query access means that for an array A[1], . . . , A[K], we can efficiently perform a unitary operation U such that along with its inverse U † . The first register, containing |i , should have dimension at least K, so that |1 , . . . , |K are each orthogonal states of the register, and the second register should be large enough to store the values A[i]. Note that Eq. (36) does not fully specify the unitary U A since it does not specify its action when the second register is not initially in the state |0 . One possible way to define U A fully is to have U |i |x = |i |x + A[i] with addition defined over an appropriately sized finite ring such as Z n 2 , but this is not necessary for applications such as in Refs. [9,11]. Quantum query access to an array A is more demanding than simply having A stored on disk, as we will discuss below.
Recall that Grover search finds one marked item out of an array of K items, assuming the ability to reflect about the marked item. Ref. [10] further extends Ref. [9] to find one marked item when there are t > 1 marked items, assuming the ability to reflect about the multiple marked items. Generic Grover search then consists of the following steps: When the number t of marked items is unknown, Ref. [10] employs an exponential searching algorithm that guesses the number of marked items, increasing the guess by a constant factor each time. This is a probabilistic algorithm that performs a measurement for each guess, finding a solution in overall expected time O( K/t).
The maximum finding algorithm of Ref. [11], summarized in Fig. 3, is based on this probabilistic exponential searching algorithm. It keeps track of the current best maximum seen so far and considers marked states to be those that have a larger array entry value than the current maximum. It employs the Grover-based exponential searching algorithm of Ref. [10] for an unknown number of marked states, performing measurements to obtain the maximum with probability at least 1/2. If desired, we can improve the success probability to 1 − η with η > 0, at the cost of an extra O(log 1/η) factor, by performing O(log 1/η) rounds of the algorithm.
Our quantum thrust algorithms are then a direct application of quantum maximum finding, but now to an array with K = O(N 2 ) entries corresponding to the choice of separating plane. To deal with the four-fold ambiguity, we use the doubling trick of Sec. III C, including each original vector p k and its negative − p k in the list of three-momenta to obtain a search space of size K = 4N 2 . Our problem, now, is to find the maximum value of T ij with i and j each ranging over 2N possible indices. This requires us to be able to load the momentum vectors corresponding to each array index, which means that the quantum algorithm must have some means of accessing the classical data.

B. Data Loading Considerations
We can describe our quantum thrust algorithms in terms of two abstract operations, LOOKUP and SUM. Their implementation will be described in Sec. V C for the sequential model and Sec. V E for the parallel model. Beyond thrust, these operations are quite general in their application to loading classical data into quantum algorithms.
Note that our search space is of size O(N 2 ), while data loading over N items in a classical database takes time O(N ). Therefore, we can conceptualize our quantum speedup as resulting from being able to perform data loading over the superposition of search space items. It is important here that the set of search space items is not the same as the set of data points. In general, any application of Grover search over a search space of size O(N α ) with α ≥ 2 will result in a square root speedup, whereas for α < 2, the cost of the algorithm will be dominated by the O(N ) data loading cost.
The LOOKUP operation is queried with one index corresponding to a given particle, returning the momentum corresponding to that index: Note that the second register, initialized as | 0 , has to be large enough to store the three-momenta to the desired (qu)bit accuracy. To make U LOOKUP unitary, we define U LOOKUP |i | q = |i | q + p i for general vectors q, where the addition is done modulo some value larger than the 1. Randomly pick an index j and set curr max = j.  with K elements [11]. The number of Grover steps is chosen at random, since this is a search over an unknown number of marked times [10]. maximum momentum encountered in the problem. To deal with pairs of particles, we can call U LOOKUP twice on different registers to map |i |j | 0 | 0 → |i |j | p i | p j . This LOOKUP operation will be used to determine all O(N 2 ) reference axesr ij , taking O(N ) time in the sequential model and O(log N ) time in the parallel model.
The SUM operation returns the sum over all momenta, possibly with a transformation f ( p; c) applied to each momentum vector: where c represent possible control qubits. From a given reference axisr ij , SUM will be used to calculate the value of T ij . It is crucial that calculating T ij for fixed i and j takes the same runtime as LOOKUP, i.e. O(N ) for sequential and O(log N ) for parallel. Notably, a wide class of collider observables can be computed in linear runtime [66], even those that would naively scale like a high polynomial power. Using LOOKUP and SUM, our quantum thrust algorithm is described in Fig. 4. As with standard Grover search, we need to be able to reflect about the initial state and the marked states, namely those whose corresponding values of thrust are larger than the best maximum seen so far. To identify the marked states, we compute thrust for each choice of separating plane, using LOOKUP and SUM to interface the quantum algorithm with the classical data. We uncompute intermediate steps of our calculations using standard methods (e.g. Section 3.2 of Ref. [4]) to make sure that, after computing T ij , the system can be reflected about the initial state.
Let C LOOKUP be the asymptotic cost of LOOKUP and C SUM be the asymptotic cost of SUM. The runtime of this algorithm is O N (C LOOKUP + C SUM ) since there is a O(N ) outer Grover search loop, while the inner loop is dominated by one application of LOOKUP and one application of SUM. Note that the computation of the initial guess for the maximum, T mn , can be performed in O(N ) time classically, while preparation of the initial state and reflection about the initial state can each be performed in O(log N ) time, the time required to perform a Hadamard gate.

C. Sequential Computing Model
The first computing model we consider is one in which one gate, classical or quantum, can be executed per time step. We should think of the classical computer as controlling the overall computation. In a single time step, it can either (a) perform a classical logic gate, (b) choose a quantum gate or measurement, or (c) read a word from the input (e.g. a single momentum). Another way to think about this model is that we measure cost by the circuit size, i.e. the total number of gates.
While fault-tolerant quantum computers are expected to require parallel control to perform error correction, there are still plausible models in which the cost of the computation will be proportional to the number of logical gates. One possibility is that the cost is dominated by generating magic states or by long-range interactions. Another possibility is that we are using a small quantum computer without fault tolerance, but in an architecture such as a one-dimensional ion trap, where the available gates are long-range and cannot be parallelized. Under this sequential model, the operations LOOKUP and SUM each take O(N ) time and require O(log N ) qubits. Specifically, LOOKUP requires a register of size O(log N ) to store the query index i, along with a register to store the requested three-momentum p i . It operates by performing a sequential scan through all N items in the classical database to fetch and return p i . More con-1. Randomly pick indices m, n and set curr max = (m, n). iii. Let iter count = iter count + grov steps. iv. Repeat grov steps times:

Compute
A. Call subroutine COMP T to compute Tij: B. Reflect about states with Tij > Tcurr max with a phase factor: Subroutine COMP T: 1. Load pi, pj using LOOKUP: 2. Calculate the reference axis viarij = ( pi × pj)/| pi × pj|:

FIG. 4. Our
Grover-based quantum thrust algorithm, written in terms of the abstract LOOKUP and SUM operations. The symbols | 0 , |0 , and |0 refer to initial states for a three-momentum, normalized axis, and real number, respectively. Note that we have applied the doubling trick from Sec. III B, such that each p k has its negative − p k in the set of three-momenta. Cases whererij · p k = 0 are treated implicitly via Eq. (26). A key difference compared to Fig. 3 is that the quantity to maximize, Tij, is calculated quantumly via the COMP T subroutine. cretely, in O(1) time, we can perform U LOOKUP,i , defined by One might wonder whether the runtime of the quantum thrust algorithm could be reduced from O(N 2 ) to O(N 3/2 log N ), using the same strategy that we used in Sec. III C to reduce the classical thrust algorithm time from O(N 3 ) to O (N 2 log N ). The answer is yes, in principle, but it would require a computing model beyond the sequential one.
Recall that two points define the partitioning plane, and after selecting the first point, we could sort the second point according to a special traversal order. This allowed us to avoid the O(N ) cost of re-summing the momenta for each candidate plane. Quantum algorithms require Ω(N log N ) time for sort [67], which means that they cannot be used to speed up this part of the classical algorithm. In principle, though, we could still obtain a Grover square root speedup when searching over the The challenge here is that to perform quantum sort, all of the data needs to be stored somehow in quantum memory, which goes beyond the sequential computing model above where only one data point is ever accessed in a given time step. We leave to future work the design of a quantum computing architecture suitable for loading and sorting data from a classical database.
Assuming that such a sort-friendly architecture exists, one might ask about the origin of the O(N 2 log N ) to O(N 3/2 log N ) speed up. Such an improvement is only possible since the strategy in Sec. III C converts thrust into a structured search problem [68,69], which evades the naive bounds on quantum search performance. Of course, no matter the degree of structure, we can never do better than the O(N ) cost to examine each data point once.

E. Parallel Computing Model
The parallel computing model reduces the time usage of the sequential model at the expense of additional space usage. 3 Under this model, the operations LOOKUP and SUM each take O(log N ) time but require O (N log N ) qubits.
An abstract version of this model is the standard quantum circuit model, in which on N qubits we can perform up to N/2 two-qubit gates on as many disjoint pairs of qubits as we like. A controlling classical computer with the same parallelism can also be used to process the measurement outcomes and feed the results back in to the quantum computer. To implement this in an actual quantum computer, we would need to assume long-range connectivity but not all-to-all connectivity. For example, Brierley [70] describes how connecting each qubit to four other qubits is enough to simulate full connectivity with O(log N ) time overhead. In what follows, we neglect any O(log N ) or other factors from converting the abstract circuit model to a concrete architecture.
Parallel data retrieval requires first pre-loading all N database items into the O(N ) qubits. This can be done in O(1) time, since it requires only parallel copy (or CNOT) operations from the classical bits onto the qubits. (Even a cost of O(N ) at this stage would not change the asymptotic runtime, so one could also consider input models in which the data could only be accessed sequentially, such as tape storage.) This results in the state (40) Note that this is not the same as qRAM [5], since we are loading the classical data into a product state once, and not assuming any kind of query access to the data. Now, given our pre-loaded data, we can perform

VI. IS THERE A QUANTUM ADVANTAGE?
Starting from the previously best known O(N 3 ) classical algorithm on a sequential computer, we found an improved O(N 2 log N ) classical algorithm and an O(N 2 ) quantum algorithm. Because these scalings are identical up to a log N factor, one might wonder if there is any real quantum advantage for the task of hemisphere jet finding.
Formally, there is a quantum advantage if we make a rather restricted assumption about the computing model. The sequential quantum computing model in Sec. V C only requires read access to the O(N ) classical dataset, whereas the sorting strategy in Sec. III C requires write access to O(N log N ) classical bits. Thus, if one restricts the computing model to have write access to only O(log N ) classical bits, then the classical sorting strategy cannot be implemented. In that case, the best classical algorithm would be the O(N 3 ) one from Ref. [2], which would be bested by our O(N 2 ) quantum algorithm.
For any realistic application of thrust, this computing model is overly limited, since data from a single collider event can easily be read into random-access classical memory. On the other hand, it is not possible to read in the entire LHC dataset into memory, and indeed some collider datasets are only stored on tape drives. For this reason, there may be interesting quantum advantages for clustering algorithms that act on ensembles of events (instead of on ensembles of particles in a single event). See Ref. [71] for recent developments along these lines.
For the parallel computing models, there is no formal limit with a quantum advantage, since we need O(N ) (qu)bits with read-write access in both the quantum and classical cases. Note that the speed up in the classical and quantum cases come from rather different sources.

VII. GENERALIZATIONS
In this section, we discuss how to apply the quantum algorithms from Secs. IV and V to jet identification methods that generalize thrust. These algorithms are more closely related to the ones used at the LHC, since they involve a jet radius parameter R.
We start with algorithms that divide the event into a single jet and an unclustered region, as in Fig. 5. (For the thrust problem, R = π/2 and the unclustered region is the opposite hemisphere.) We then mention strategies to identify multiple jets. To simplify the discussion, we continue to use the (p x , p y , p z ) coordinate system for electron-positron collisions, noting that the methods below can be adapted to the standard proton-proton coordinate system of transverse momentum (p T ), rapidityn R FIG. 5. Partitioning an event into a stable cone jet of radius R and an unclustered region. This is the same as Fig. 1 when R = π/2.

A. SingleCone
The generalizations we consider are all based on or inspired by the analysis of Ref. [61], which showed that the thrust duality in Sec. II D holds for a one-parameter family of jet finding algorithms. No matter which dual formulation is used, we refer to this jet finding strategy as SingleCone, since it finds a single stable cone jet of radius R.
To match the literature, we use four-vector notation in this section. The four-momentum of a particle is where the energy E i ≡ p 0 i = p 2 i + m 2 i depends on the mass m i of particle i. The four-momentum of a candidate partition H is where E ≡ P 0 is the total energy of the partition. A light-like axis is given by withn 2 = 1. We contract indices with the mostly minus metric: The SingleCone jet finder is based on maximizing the following objective function [61]: where λ is again a Lagrange multiplier, and we maximize over both the choice of partition and the choice of axis. The second line makes it clear that R = π/2 returns the thrust objective function in Eq. (17).
Performing the same manipulations as in Sec. II D, the optimum axis (for fixed partition) is Since the optimum axis is aligned with the jet threemomentum, this is an example of a stable cone algorithm; see Sec. VII C below. The reduced SingleCone objective function is which is an example of a jet function maximization algorithm [43][44][45]. The optimum solution partitions the event into a clustered region H and an unclustered region (the complement of H). This definition of the problem naturally lends itself to quantum annealing in Sec. VII B. Doing the dual manipulation, the optimum partition (for fixed axis) is Writing the Heaviside theta function requirement in three-momentum language, we see that for massless particles (E i = | p i |), the jet constituents are those within an angular distance R of the jet axis. For R = π/2, this yields the thrust hemisphere regions. The reduced SingleCone objective function is now where we dropped the Lagrange multiplier term for compactness. The second term in Eq. (50) is an example of an N -jettiness measure [72][73][74] with N = 1, whose minimum yields the XCone jet algorithm [32,33]. This definition of the problem naturally lends itself to quantum search in Sec. VII C.

B. Jet Function Maximization
In the jet function maximization approach of Refs. [43][44][45], the goal is to optimize P µ for a global jet function.
The original jet function from Ref. [43] can be written as where the jet mass is In the limit M E, this matches the reduced Single-Cone objective function of Eq. (47), though they yield different optimal jet regions for finite-mass jets.
Since jet function maximization is a kind of partitioning problem, it is natural to try to write these objective functions in QUBO form. However, the original jet function from Eq. (51) is not quadratic since it involves a 1/E factor, and the SingleCone function in Eq. (47) is not quadratic since | P | involves a square root. Thus, these cannot be rewritten as QUBO problems without some kind of modification.
In the analysis of Sec. IV for thrust, we got around this issue by squaring the thrust objective function, which nevertheless yielded the same partitioning solution. This approach does not work in this more general case because of non-quadratic cross terms.
What we can do, however, is square the SingleCone objective in Eq. (47) but only keep the lowest non-trivial term in the M E limit. 4 (Squaring and expanding Eq. (51) yields the same result.) This gives the following QUBO objective function: where again x i ∈ {0, 1}. Taking R = π/2 in Eq. (53) then recovers the thrust (squared) problem. It is interesting that Eq. (53) has the same form as the generalized jet functions in Ref. [44] (Ref. [45]) with n = 2 (α = 2). This objective function corresponds to a QUBO problem and can thus be solved on a quantum annealer. It will, however, generally yield a different solution compared to SingleCone. Unlike SingleCone, which yields perfectly conical jets for massless particles via Eq. (49), this QUBO jet finder has an effective jet radius that depends on the mass of the jet [44,45]. Quadratic objective functions were also explored in Ref. [47] for jet clustering at the LHC. In future work, we plan to characterize the general phenomenological properties of jets identified using QUBO objectives.

C. Stable Cone Finding
Stable cone algorithms search over candidate jet regions of radius R and select ones that are "stable" [34,75], meaning that the center of the jet region aligns with the jet momentum. As shown in Eqs. (46) and (49), Sin-gleCone is an example of a stable cone algorithm, which is closely related to SISCone [3].
It is worth emphasizing two key differences between SingleCone and SISCone. First, SingleCone finds a single jet, where as SISCone finds all stable cones, and a separate split/merge step is needed to determine the final jet regions. That said, it is possible to run SISCone in progressive removal (PR) mode, where one finds the most energetic stable cone, removes the found jet constituents, and repeats the SISCone procedure on the unclustered particles. In this way, SISCone-PR acts like an iterated application of SingleCone. Second, SingleCone finds the jet region with the largest value of Eq. (47) (= E − O(M 2 /E)), whereas SISCone-PR would typically take the stable cone with the largest plain energy E. As we will see below, though, it is still possible to develop quantum algorithms for stable cones with alternative jet hardness sorting schemes.
It is straightforward to implement the SingleCone algorithm (a.k.a. SISCone-PR with Eq. (47) ordering) via quantum search. Just as two points define a partitioning plane, two points are enough to determine a cone region of radius R [3]. (This is true up to an eightfold ambiguity, which is twice that of the thrust case because the two candidate cones are not complements of each other as they are for hemispheres.) We can use the LOOKUP operation to determine all O(N 2 ) candidate reference axes (which are not the same as the jet axes, but yield the same partitions). We can then use SUM to calculate Eq. (47) for a fixed reference axis, since finding P µ for the particles in the candidate jet region is a linear operation. We finally use Grover search to find the partition that maximizes Eq. (47), and we are guaranteed that the found cone jet will be stable via Eq. (46). This algorithm now has the identical structure to thrust, with the same asymptotic scaling as in Table I, taking us from a classical O(N 3 ) algorithm (without sort) to a quantum O(N log N ) algorithm (with parallel data loading).
Note that the quantum maximum finding algorithm only returns one maximum element of an array, so we cannot use it to speed up an algorithm for identifying all stable cones. We can, however, use it to find one stable cone with a different objective function from Eq. (47). For example, to implement SISCone-PR with standard energy ordering, we can use a subroutine consisting of two SUM operations in series. The first SUM determines P µ for the candidate jet region, while the second SUM findsP µ for all particles within a radius R of P µ . This subroutine would return P 0 if P µ =P µ , while it would return 0 if P µ =P µ . One would then use Grover search to find the maximum subroutine output, with the same asymptotic quantum gains as in the SingleCone case.

D. Multi-Region Optimization
Typical collider studies involve more than one jet per event, so it is interesting to ask whether these quantum methods can be adapted to the multi-jet case. As already mentioned, one can use a PR strategy to identify multiple jet regions, so finding M jets just requires M iterations of the algorithms above. Except in specialized circumstances, the number of desired jet regions does not grow with N and is at most O(1/R 2 ), so the runtime of SingleCone-PR would scale linearly with M . That said, we are interested in simultaneously optimizing the jet regions as in XCone algorithm [32,33], in order to treat the overlapping jet regions in a more sophisticated way than just PR.
The QUBO objective in Eq. (53) can be easily generalized to the M -jet case using O N (M + 1) qubits, suitable for quantum annealing. Instead of a binary assignment of each particle to the clustered or unclustered region, we can do a one-hot encoding with M + 1 qubits per particle to indicate their assignment to one of the M jet regions or to the unclustered region. Specifically, let x ir ∈ {0, 1} for i ∈ {1, ..., N } and r ∈ {0, 1, ..., M }. We assign x i0 = 1 if particle i is in the unclustered region, x ir = 1 if particle i is in jet region r for r ∈ {1, ..., M }, and x ir = 0 otherwise. We then add a penalty term to the objective function such that, for fixed i, x ir = 1 for only one value of r.
The multi-jet QUBO objective function is Here, there is a copy of Eq. (53) for each of the M jet regions, taking the schematic form of O = − i,j Q ij x i x j . The coefficient of the penalty term must be taken to be Λ 2 > N max ij Q ij to ensure that it is never favorable for a particle to be assigned to more than one jet region. Because Eq. (54) is quadratic in the momentum, it will not have the same behavior as XCone (which has a linear objective function), though we expect the results to be similar for well-separated jets of comparable energies. This objective function does not penalize empty jet regions, so it might be interesting to run this algorithm with a large value of M to let the number of non-empty jet regions be determined dynamically.
Compared to the single-jet case, the multi-jet case will likely be more difficult to implement on currently available quantum annealing hardware. Previous numerical studies [54] have shown that clustering problems that use multiple qubits to implement one-hot encoding are prone to errors. The reason is that on annealing hardware, qubit couplings have a maximum dynamic range, which in turn limits the effectiveness of the Λ penalty term. In practice, this means that annealers often output a fuzzy assignment rather than a hard assignment to one cluster. We would also like to argue that this problem is conceptual in origin. The search space of the single-jet QUBO problem is 2 N , whereas the search space of the multi-jet QUBO problem is 2 M N . However, the QUBO quantum search space contains many extra unphysical states, since the actual (non-QUBO) search space is size M N = 2 N log M . While the most natural way to address this would be to use qudits with d = M instead of qubits, such hardware is not currently available.
Turning to the quantum search case, finding M conical jet regions naively requires searching a space of O(N 2M ), with the added complication of needing to treat overlapping jet regions. We are unaware of any classical approach to this problem apart from brute force, though one expects an O(N 2M +1 ) algorithm for the XCone objective should be feasible, though it likely requires a more sophisticated treatment of reference axes. (The current implementation of XCone in FastJet contrib 1.041 [31,76] only finds a local minimum starting from suitable seed axes.) Using quantum search with sequential (parallel) data loading, one might hope that this could be improved to O(N M +1 ) (O(N M log N )), though one would have to generalize the LOOKUP and SUM operations to deal with the multi-jet case. At minimum, LOOKUP would have to load the momenta into 2M registers (to label the candidate partitions), and SUM would have to have M distinct outputs (for each of the M jet regions). Even with quantum gains, this is computationally daunting, motivating future studies of multi-jet algorithm whose computational complexity grows only polynomially with M .

VIII. CONCLUSIONS
In this work, we demonstrated how quantum computers could be applied to a realistic collider physics problem, which requires interfacing a classical data set with a quantum algorithm. We focused on maximizing thrust to identify hemisphere jets, but the quantum methods developed here are relevant to a broader range of optimization and cluster-finding problems. The asymptotic performance of our quantum annealing and quantum search algorithms are summarized in Table I. We found a way to improve the previously best known O(N 3 ) classical thrust algorithm to an O(N 2 ) sequential quantum algorithm. Along the way, we found an improved O(N 2 log N ) classical algorithm, based on the sorting strategy of Ref. [3]. Both the quantum and improved classical algorithms can be implemented on parallel computing architectures with asymptotic O(N log N ) runtime. Formally, we found a quantum advantage, but only when assuming a computing model with read access to O(N ) (qu)bits but write access to only O(log N ) (qu)bits.
Going beyond thrust, we briefly generalized our quantum methods to handle structurally similar jet clustering algorithms. These involve maximizing an objective function with a radius parameter R, which partitions the event into a conical jet region and an unclustered region. While we focused on electron-positron collisions, it is known how to adapt these methods to proton-proton collisions [45,47,61]. In future work, we plan to investigate the phenomenological performance of these "quantum friendly" jet algorithms at the LHC, to assess whether they offer improved physics performance relative to hierarchical clustering schemes like anti-k t .
The main take home message from this work is that the overhead of data loading must be carefully accounted for when evaluating the potential for quantum speedups on classical datasets. In many ways, optimization-based jet algorithms are an ideal platform to think about quantum algorithms for collider physics, since these problems tend to involve searching over a large space of possibilities, O(N α ) with α ≥ 2, and therefore benefit from Grover search methods. By contrast, even though the number of events in a collider data sample (N events ) is usually much larger than the number of final-state particles in a jet, typical collider tasks like filling a histogram involve O(N events ) operations, such that data loading is already the limiting factor. On the flip side, this motivates further quantum investigations into classically O(N 2 events ) data manipulation strategies, such as the metric space approach recently proposed in Ref. [71], since they might be reducible to O(N events ) quantum algorithms under suitable circumstances. We also note that Grover search is limited to a square-root speedup on unstructured search, whereas collider data has additional structures like symmetries and heuristics which might lead to further quantum gains.