A quantum algorithm for track reconstruction in the LHCb vertex detector

: High-energy physics is facing increasingly computational challenges in real-time event reconstruction for the near-future high-luminosity era. Using the LHCb vertex detector as a use-case, we explore a new algorithm for particle track reconstruction based on the minimisation of an Ising-like Hamiltonian with a linear algebra approach. The use of a classical matrix inversion technique results in tracking performance similar to the current state-of-the-art but with worse scaling complexity in time. To solve this problem, we also present an implementation as quantum algorithm, using the Harrow-Hassadim-Lloyd (HHL) algorithm: this approach can potentially provide an exponential speedup as a function of the number of input hits over its classical counterpart, in spite of limitations due to the well-known HHL Hamiltonian simulation and readout problems. The findings presented in this paper shed light on the potential of leveraging quantum computing for real-time particle track reconstruction in high-energy physics.


Introduction
Progress in the field of experimental Particle Physics is often limited by statistical precision.For this reason, near-future upgrades to collider experiments, such as the Large Hadron Collider (LHC), are investing in an order of magnitude increase of data collection rates or instantaneous luminosity.While this brings exciting physics targets within experimental reach [1][2][3], it comes with a drastic increase in the number of particles in the final state to be detected.These particles traverse the various sensor layers of the detector and generate hits, which have to be reconstructed into the original particle trajectories in real-time for data storage selection and further processing.
As the complexity of particle track reconstruction scales with the number of sensor-hits to the power 2-3 [4,5], depending on the algorithmic approach and geometry of the detector, time constraints and reducing fake rates provide challenges that are still to be overcome for High Luminosity-LHC (HL-LHC).Various approaches are under study, such as parallel reconstruction on GPUs [6] as well as the addition of picosecond timing information to sensors [7].In this paper, we will explore the use of quantum algorithms for particle track reconstruction in the Upgrade LHCb vertex detector (VELO) [8].
Quantum Computing (QC) represents one of the current frontiers of Computer Science and has seen rapid developments on both the hardware and the software sides over the past decade.QC has gained the attention of the High Energy Physics community [9], offering the possibility of overcoming the current computational limitations of algorithms on simulation [10,11], data analysis [12][13][14] and reconstruction [15].
Several studies exist on the potential advantage of QC algorithms for track reconstruction.A complexity study [16], shows some degree of potential speedup in quantum search algorithms using a seeding/following approach.Alternatively, a Quadratic Unconstrained Binary Optimisation (QUBO) formulation of the tracking problem and proof-of-principles have been built using quantum annealing [17] or a Variational Quantum Eigensolver (VQE) [18] to find solutions, obtaining reasonable tracking performance numbers but making no promise on HL-LHC conditions or timing improvements.Considering the reconstruction of tracks from hits as a clustering problem, one deals with a global algorithm, processing all the hits at once.An example of such an implementation was done at the LUXE experiment [19], where positrons trajectories have been reconstructed using a VQE to solve a QUBO problem, and a quantum Graph Neural Network (GNN) [20,21].A similar approach, based on a hybrid quantum-classical GNN, has been tested on the TrackML dataset, which emulates HL-LHC conditions [22].In both cases, a practical implementation on quantum hardware is limited by the number of input hits that the algorithm can process, due to the circuit becoming too deep to be compatible with current devices.
In this paper we adopt the QUBO formulation of the tracking problem, based on an custom Ising-like Hamiltonian, and we solve it in an alternative way, using a matrix inversion method.The well-known Harrow-Hassadim-Lloyd (HHL) algorithm [23] is employed for matrix inversion of a Hamiltonian embedded in a quantum circuit.We show that the preparation of the input quantum state can be executed efficiently on a quantum device and that the matrix to be inverted is highly sparse and well-conditioned, matching the requirements of HHL.We benchmark our algorithm on LHCb simulated data and test a quantum implementation on toy model events.When the hypotheses for an efficient implementation of HHL will be met, this approach promises an exponential speedup over classical counterparts.

Track reconstruction in LHCb vertex detector 2.1 The LHCb vertex detector
The LHCb detector [24,25] is a single-arm forward spectrometer that covers the pseudo-rapidity1 range 2 <  < 5.The reconstruction use case investigated is based on the LHCb VELO tracker [8].The VELO consists of a series of vertically oriented pixel sensors arranged along the beam-line (on the -axis), positioned close to the collision point of the LHC, and measuring the  coordinates of passing particles (see Figure 1).The collisions occur in a 5 cm long region along the -axis, with detection layers surrounding the beam-line axis  =  = 0.The VELO is placed upstream to the dipole magnet, in a region practically free from magnetic field, leading to straight tracks and simplifying the track model when considering only hits in the VELO.

Track reconstruction algorithms
The reconstruction of particle trajectories in high-energy physics experiments at a hadron collider is a computing-intensive task due to the high multiplicity of particles in the final state.
Approaches to track reconstruction can be classified into two categories: local and global algorithms.Global algorithms process all the hits in an event simultaneously, aiming to provide a globally optimal solution to the track-finding problem; this category includes approaches based on the Hough transform methods [26], Hopfield neural networks [27,28], etc.On the other hand, local algorithms process the hits in a sequential way, starting from a set of track seeds: as a consequence, 1The pseudo-rapidity is defined as  = − log tan   2 , where  is the angle between the three-momentum of the particle and the beam axis.
the result of the algorithm may depend on the choice of the starting seeds and processing order of the hits.This class includes Kalman filter techniques [29,30].
The choice of which algorithm fits best with a specific use-case depends on the specifics of the detector, the operating conditions and the available computing resources.Generally, global methods tend to scale better in time complexity, while for space complexity local algorithms may be more suitable.Both types of methods have been widely adopted in particle physics experiments [28,31].
The pattern recognition problem used to find particle trajectories is increasingly approached using a heterogeneous computing model.Where before the task was fully CPU-based, recent developments add GPGPU-based algorithms in LHCb.
The current state-of-the-art VELO track reconstruction algorithm is such a GPGPU-based algorithm called Search-by-Triplet [4] exploiting a high degree of parallelism.Track candidates are built from hits in three consecutive layers close to the primary proton-proton interaction point, which are then projected in  to find hits in the following detector layers (track following).The quantum algorithm investigated in this paper, instead, is a global method considering all hits at once, and the performance of Search-by-Triplet will act as a baseline comparison.
3 Tracking with Ising-like Hamiltonian

Ising-based methods
In this paper, we present a quantum algorithm for track reconstruction based on the minimisation of an Ising-like Hamiltonian.The algorithm takes, as an input, a collection of hits left by a set of particles over the several layers of the detector.The algorithm will reconstruct each particle trajectory as a collection of track segments (also called doublets) that belong to that track.The first step consists in constructing all the possible doublets of hits in consecutive modules of the detector: each doublet is associated with a binary variable   ∈ {0, 1} where 1 if the doublet is part of a track 0 otherwise .Therefore, the track reconstruction task consists of determining the vector S = { 1 ,  2 , ...,   } containing the correct activation state of every doublet, i.e. 1 if it is part of any track or 0 otherwise.Ising-based approaches to this problem consist in finding a suitable quadratic Hamiltonian in the form such that the solution to the problem coincides with its ground state.This represents a generalised Ising Hamiltonian with couplings described by the real symmetric matrix  (the Ising interaction matrix) and an external field described by the real vector b (bias vector).In general, finding the optimal value of S is a QUBO problem and requires a number of computations of H that scales exponentially with the size of S [17].The explicit representation of (3.1), i.e. the coupling matrix  and bias vector b, must be chosen in such a way that the minimum of H corresponds to the solution of the tracking problem: in this paper, we adopt a modified version of the Denby-Peterson Hamiltonian.

Denby-Peterson Hamiltonian
The Denby-Peterson (DP) Hamiltonian [32,33] is an Ising-like Hamiltonian in the form (3.1), which is often used in track reconstruction algorithms.It favours the reconstruction of smooth and long tracks while preventing non-physical bifurcations.Let   be the oriented segment stemming from hit  to hit .The DP Hamiltonian can be expressed in the form highlighting the following three main components: Angular term where   is the angle between the segments   and   , as depicted in Figure 2c;   and   are the length of the segments   and   respectively.This term favours aligned triplets (i.e.where   is small) with short doublets. is a constant integer: large values of  penalise less straight triplets in favour of aligned ones.

Bifurcation term
is weighted by the Lagrangian multiplier  ≥ 0. This term penalises, with a constant penalty +, bifurcations in the trajectories like the ones shown in Figure 2b.
Occupancy term where  hits is the number of hits in the detector.This term constrains the number of active doublets to be roughly equal to the number of hits in the event and it is weighted by the Lagrangian multiplier  ≥ 0.
The constants ,  and  are free parameters of the model and must be fixed in advance according to an appropriate hyper-parameter search.
The DP Hamiltonian has been designed to reconstruct curved trajectories, minimising the amount of bifurcations in the reconstructed tracks.However, as the LHCb VELO operates in a region practically free of magnetic field, tracks are approximately straight lines, therefore we expect the angular term to have the dominant contribution, relatively to the other two terms, which can be eventually neglected.Moreover, as the quantum implementation for this algorithm is expected to scale better when number of non-zero couplings in the Hamiltonian is small, it can benefit from a simplification of the terms described above.For these reasons, in this work, we use a modified version of the DP Hamiltonian, removing the bifurcation and occupancy terms (i.e. and  are fixed to 0), and using a simplified angular term: with This term, compared to (3.3), removes the dependence on the segment length and replaces the smooth cosine behaviour with a step function of width arccos(1 − ).The constant  represents an angular tolerance for the model: only consecutive segments that form an angle smaller than arccos(1 − ) are coupled by the angular term, while the others are left uncoupled.The parameter  allows to deal with trajectories that are not perfectly straight.Therefore, its numerical value can be related to the physical characteristics of the detector, such as its geometry and resolution.It also allows taking into account multiple-scattering effects, i.e. when a particle is slightly scattered by the interactions with the detector material.Figure 3 shows the decay behaviour of cos  () for different values of the exponent , compared to the step function  with  = 10 −5 .

The quantum algorithm
The DP Hamiltonian has been employed in track reconstruction algorithms based on the optimisation of Hopfield networks [27,28,32,33]: nevertheless, these approaches become unfeasible in high hit multiplicity conditions due to the scaling of the size of the network.Several quantum versions of the Hopfield network are available in the literature, trying to exploit the unique features of quantum computation to achieve an advantage [34][35][36][37][38].Our quantum algorithm follows a similar approach to Rebentrost et al. [38], in which the optimisation of a Hopfield network is turned into a matrix inversion problem and solved with the Harrow-Hassadim-Lloyd (HHL) algorithm (also known as quantum algorithm for linear systems of equations), which has the potential of achieving an exponential advantage with respect to its classical counterpart [23].

Hamiltonian Relaxation and local minimum
The first step consists of relaxing the original problem, where we allow the binary-valued spin   in the Ising-like Hamiltonian (3.1) to have any possible real value   ∈ R, instead of only 0 or 1.This way, the Hamiltonian is a differentiable function of the variables   , and this allows us to calculate the gradient.The drawback of this approach is that we allow each   to have a non-physical state, such as a negative or fractional value: to solve this problem, a threshold operation is performed in the final stage of the algorithm, which will ensure valid final states.Using the fact that  is symmetric, the gradient of H with respect to S has the following expression Now we look for a S that satisfies ∇H (S) = 0, the necessary condition for being a local minimum.This implies solving the square system of linear equations In general a solution for (4.2) is not guaranteed to exist (i.e. if det  = 0).To overcome this, we look for a least-squares solution for this system: where || • || 2 denotes the Euclidean 2-norm.In general, (4.3) may have multiple solutions: in that case, we select the one with the smallest 2-norm.This problem has a unique solution that is provided by the Moore-Penrose pseudo-inverse matrix [39]  + applied to the vector b: If det  ≠ 0 then the pseudo-inverse  + and the inverse  −1 coincide.

Regularization terms
In order to deal with the relaxed version of the problem, we have to introduce two regularisation terms in the Hamiltonian (3.6), which becomes Spectral term This term is purely quadratic and contributes to the diagonal of the matrix .The purpose of this term is to shift the eigenvalues spectrum of  by a constant  to make it positive-definite and avoid numerical instabilities.However,  cannot be chosen arbitrarily large as this would result in a polynomial increase of the run-time [38].

Gap term
Due to the presence of linear and quadratic contributions, this term contributes to both the matrix  and the vector b, by a constant .The purpose of this term is to enforce a gap in the spectrum of the solution vector S, suppressing the occurrence of ambiguous states that are not close to 0 or 1.This allows to identify a stable threshold value  and discretise the output values as The numerical values of the hyper-parameters have been chosen according to the testing conditions of the algorithm.The parameter , which appears in the angular term (3.7), has been fixed to match the average multiple scattering angle in a LHCb VELO module [40].The spectral parameter  has been fixed such that the eigenvalues spectrum of  stays positive.It has been observed that the value of gap parameter  is correlated with the threshold , but does not directly affect the performance of the algorithm: therefore  has been fixed to 1.0 and the corresponding threshold  has been determined as the average gap midpoint observed on events generated with a toy model, described in Section 5.1.The numerical values of all hyper-parameters are reported in the table below.To solve the system of linear equations (4.2) we use the HHL quantum algorithm [23,38].Under proper assumptions, the HHL algorithm is able to solve a  ×  sparse linear system of equations in O ( 2 log ) run-time, where  is the condition number of the matrix [41], achieving an exponential advantage in the problem size over classical algorithms [23].This advantage is achieved thanks to a powerful subroutine called Quantum Phase Estimation (QPE), which is also used in other quantum algorithms that show an exponential advantage, such as Shor's algorithm [42].It requires the matrix  to be square and its dimension to be a power of 2. If that is not the case, in practice, it is always possible to pad  to round up its dimension to the closest power of 2. Similarly, the b vector can always be padded with a suitable number of ones and normalised.

Angular Spectral Gap Threshold
Therefore we will henceforth assume  to be a power of 2 and ||b|| = 1.QPE also requires the matrix  to be sparse, in order to be implemented efficiently.The Hamiltonian presented in this work has been designed to match the sparsity condition in practice: a more detailed study about the sparsity dependence on the event size can be found in Appendix A. The HHL run-time complexity has a quadratic dependence on the condition number : a study of the behaviour of the condition number with respect to the number of particles and detection layers has revealed that  can be upperbounded irrespective of the problem size and therefore does not degrade the exponential advantage of HHL.More details can be found in Appendix B. The algorithm requires a register of   = log 2  qubits |Φ⟩  for storing the vector b and the solution vector S.An additional register |Φ⟩  of   qubits is required for QPE, according to the required precision on the estimated eigenvalues:   is chosen as follows However, since  is upper-bounded, 1 +   is dominant for  ≥ 4. Finally, the matrix inversion step requires an additional single-qubit register |Φ⟩  .The global quantum state is represented as the tensor product of these 3 registers The algorithm is implemented by the quantum circuit shown in Figure 4 and consists of the following steps: 1. Initialisation: Each qubit is initialised to the |0⟩ computational state 2. State preparation: the entries   of the vector b have to be loaded into the amplitudes of a quantum register |Φ⟩  .
where |⟩ represents the -th computational basis vector.This requires a suitable unitary   to construct the state |⟩ from the initial state |0⟩, namely |⟩ =   |0⟩.This is usually a strong limitation of the HHL algorithm as no efficient arbitrary state preparation is currently known.In this application, the vector b gets its only contribution from the gap term(4.7),which is a constant.Therefore, including the normalisation constant , b has the following structure b = 1
The corresponding state |⟩ can be constructed from the initial state |0⟩ by applying a single Hadamard gate to each of the   qubits.

Quantum Phase Estimation:
The QPE subroutine is applied to the operator  =    and will determine its phases with   -bit precision, where   is the number of ancilla qubits.QPE employs multiple controlled executions of integer powers of ,  2 ,  3 , ....In order for this to be done efficiently, one needs a quantum circuit that implements the Hamiltonian simulation operator  () =    .Such an efficient implementation for this approach is not available yet and therefore represents the current limitation of this algorithm.

Matrix inversion:
The matrix  is inverted in its own eigenbasis, by taking the reciprocal of each eigenvalue.This is accomplished by a controlled rotation of the ancilla qubit |Φ⟩  .The algorithm ensures the successful inversion of the matrix if the ancilla is measured in the state |1⟩.If |Φ⟩  is found to be in state |0⟩, then the state is discarded and steps 1,2,3 and 4 are repeated until |1⟩ is measured.

Uncomputing:
The ancilla registers |Φ⟩  are cleaned up by uncomputing them: in practice, the QPE subroutine is inverted and applied to the |Φ⟩  register to bring them back to their initial state.This step is necessary to avoid that, after the computation, the floating state of the ancilla registers interferes with the measurement of the qubits where the result is stored [43].
6. Post-processing and measurement: At this stage, the register |Φ⟩  contains the entries of the solution vector The measurement of all the amplitudes of the final state, at this point, would require at least  repetitions of the algorithm, voiding the exponential advantage.One can employ quantum post-processing techniques on the final state and measure the expectation value of an observable ⟨|O|⟩.The identification of such an observable is currently out of the scope of this work, and state-vector analyses have been performed using a quantum simulator.The quantum algorithm has been implemented [44] using the Qiskit open-source Python library [45].

Figure 4:
Quantum circuit that implements the algorithm.In red, the |Φ⟩  , which encodes the b input vector and the S output vector.In green, the   -qubit ancilla register used for QPE.In blue, the additional ancilla qubit used for the matrix inversion.

Results
The algorithm is validated by assessing the track reconstruction performance using simulated events.Due to the lack of an efficient Hamiltonian simulation subroutine, processing LHCb simulated events using HHL is currently unfeasible: the HHL implementation available inside the Qiskit library produces circuits that are too deep to be run on current quantum hardware and to be simulated efficiently.For this reason, we have opted for a hybrid validation procedure: track reconstruction performances have been measured on LHCb simulated events, replacing the HHL algorithm with a classical solver, while the quantum approach has been tested using much smaller and cleaner events generated with a toy model.In the full simulations, the LHCb VELO geometry is simulated including multiple scattering of the particles on the detector material, as well as realistic hit detection efficiency and resolutions.

Detector toy model
Toy simulations are used for R&D purposes of the algorithm: they implement a flexible detector geometry and allow us to study the algorithm in ideal conditions, i.e. perfectly straight tracks and no spurious hits, varying the number of detection layers and particle multiplicity.The detector is made of a tunable number of layers placed parallel to the -plane at regular distances from the origin.Particle trajectories are lines stemming from the origin, generated uniformly within the geometrical acceptance of the detector modules.Hits are collected from the intersections between the track and the detector layers.For the toy simulations ideal detector performance, in particular, no scattering, 100% detection efficiency and perfect hit resolutions are assumed.The code implementing the detector toy model is available in Ref. [44].

Track reconstruction performance on LHCb simulated data
To benchmark the algorithm, we have used a selection of LHCb simulated events of   → , being a decay channel of major interest for the LHCb physics programme.The simulations used in this study have been obtained generating   collisions at centre-of-mass energy √  = 13.6 TeV from proton bunches crossing in the 2022-2025 operating conditions of LHC, using Pythia [46] with a specific LHCb configuration [47].Decays of hadronic particles are simulated by EvtGen [48], in which the final-state radiation is generated using Photos [49].The interaction of the particles with the detector and its response are implemented using Geant4 [50] as described in Ref. [51].We selected events that leave up to 1000 hits in right-half of the VELO and these have been processed using the described algorithm where a classical solver has been used to solve the system (4.2).To validate the result, the list of active segments is processed into a set of tracks, each one being the set of hits it contains.The display of an example event and the related reconstructed tracks are shown in Figure 5.The reconstructed tracks are then validated using the generated, true information and the tracking performances of the algorithm are evaluated in terms of the following metrics: • Track-finding efficiency  track : the fraction of correctly found tracks: where  corr track is the number of correctly reconstructed tracks and  acc gen is the number of generated particles that are in the acceptance of the detector volume.The acceptance is taken to be the fraction of generated tracks that leave hits in at least 3 detection layers of the right-side of the VELO.For a track to be considered correctly reconstructed, at least 70% of the hits of the reconstructed track should be correctly assigned.
• Track-finding fake rate  fake : the fraction of found tracks that are not associated to a generated particle: where  fake =  all track −  corr track is the number of reconstructed tracks that are not associated with a real track.
• Hit-purity  pure : the fraction of hits of a track that are correctly assigned to the same generated particle associated with the track: where  all hit are the total number of hits assigned on a reconstructed track and  correct hit are the correctly assigned hits.
• Hit-efficiency  eff : fraction of the hits left by a generated particle that are correctly associated: where  gen hit is the total number of hits left by the generated particle.
Figure 6 shows the performance of the algorithm, in terms of hit-efficiency, hit-purity and trackfinding efficiency, as a function of the pseudo-rapidity , the momentum  and the transverse momentum  T of the generated particles.The algorithm shows an average track-finding efficiency of around 97%, with degraded performance at the boundary of the detector acceptance ( below 2.5 and above 4.5), and for particles with low momentum (below 4 GeV) as their trajectories are affected more by multiple-scattering effects.The hit-efficiency is high (around 98%) and shows a constant behaviour, except for the high pseudo-rapidity region where a few hits of high hit-multiplicity tracks are missed with minor impact on the track-finding efficiency.The hit-purity is also high (above 98%) with values above 99% in the low pseudo-rapidity region, where tracks have a low hit-multiplicity.The integrated fake tracks rate has been found to be 4.3%.These results are comparable to the performance of the current state-of-the-art algorithm, Search-by-Triplet [4].Nevertheless, the latter has still an edge, indicating that our algorithm can still be further improved with a more refined hyper-parameter search, which goes beyond the purpose of this work.
The algorithm using a classical solver has been also tested on toy events, generated with the model described in Section 5.1, containing up to 1000 hits over 26 detection modules and resulted in a perfect (i.e.100%) segment finding efficiency, with a segment purity of 100% (93%) when setting the angular parameter to  = 10 −9 (10 −5 ), indicating excellent performance of our new algorithm under ideal conditions.

Qiskit implementation
To validate the usage of the HHL subroutine, the algorithm has been tested using events generated with a detector toy model, described in Section 5.1.In order to test the implementation using the HHL algorithm, a number of particles ranging from 2 to 5 are generated.For each generated event, the matrix  and the vector  of (4.2) are calculated for the list of hits left in the detector modules.The HHL implementation provided by the Qiskit Aqua Python package [45] has been used to generate a circuit that solves the system of linear equations.The circuit is then executed on the Qiskit Aer noise-less simulator and results are retrieved by performing a state-vector analysis of the final quantum state.As an example, the largest (in terms of the number of doublets) reconstructed event is shown in Figure 7. Additionally, circuits are transpiled (optimisation level 3) to fit the ibm_hanoi 27-qubit r5.11 Falcon quantum processor.This allows us to get realistic circuit sizes for the current implementation of the algorithm on present devices.The number of qubits required, the total depth and number of 2-qubit gates are shown in Table 1, together with the size of the input events.Under these conditions, the algorithm was able correctly identify all the tracks in each toy event, with hit-purity  pure = 100% and hit-efficiency  eff = 100%.On the downside, the obtained circuit depths range from a few thousands to several millions suggesting that this implementation is not yet in a state to be executed on current circuit-based quantum hardware: the Qiskit implementation   The Qiskit implementation has been tested for several event sizes, in terms of the number of layers and particles.For each configuration, the number of qubits and depth of the resulting circuit is reported, as well as the number of 2-qubit gates required on ibm_hanoi.

Layers Particles Doublets
of HHL computes of the operator    classically and performs a Pauli decomposition of the result.Only then the resulting unitary matrix is implemented as a quantum circuit.This approach does not take into account the sparsity pattern of the matrix A for the exponentiation and, as a result, the computation has a high classical computational cost and results in very deep circuits, nullifying the advantage provided by HHL.Nevertheless, several theoretical results suggest that, given the high sparsity of the matrix , its condition number behaviour and structure, such an efficient implementation is possible, and it will be investigated in future works [52][53][54][55].

Conclusions
In this paper, we have presented a new algorithm for the reconstruction of charged particle tracks based on the minimisation of an Ising-like Hamiltonian H .The ground state of H is obtained using a relaxation procedure that allows turning the problem into a sparse system of linear equations.This approach to the problem is suitable to be implemented on a circuit-based quantum computer using the HHL algorithm to solve the system of linear equations.If all the hypotheses for its efficient implementation are met, this novel approach has the potential to lead to an exponential advantage over classical algorithms.The studies of the sparsity and condition number scaling show that the algorithm presented in this paper produces sparse systems of linear equations that are very well-conditioned when scaling up the size of the problem, therefore it matches the requirements to achieve an efficient implementation of the HHL algorithm.
Our new Hamiltonian approach has been validated using events coming from simulated   →  events in the LHCb vertex detector: the algorithm has shown a performance close to the state of the art, with an average reconstruction efficiency above 96% and an average fake rate below 5%, using a classical sparse matrix inversion method.In addition, we have presented a quantum implementation of the same algorithm using the Qiskit library and we have tested it on smaller events generated with a toy model.In all the cases, the correct solution has been retrieved by the algorithm.
An analysis of the size of the resulting quantum circuits clearly shows that the deployment on current quantum hardware is unfeasible at present.This result is strongly limited by the current implementation of Hamiltonian simulation in the QPE subroutine of HHL, which produced extremely deep circuits: this behaviour is expected as the Qiskit implementation of HHL does not make any assumption on the structure of the matrix  and therefore the QPE subroutine becomes extremely expensive.In general, it is hard to achieve efficient simulation if no assumptions are made on the structure of the Hamiltonian.Nevertheless, due to the structure of our matrix , there is strong theoretical support that this should be possible in our use-case [52][53][54][55].Future studies on improvements to this section of the algorithm will bring drastic enhancements to the size of the resulting quantum circuits, possibly enabling the execution of the full pipeline on a quantum device.
The HHL algorithm encodes the solution vector as amplitudes of the final state.Therefore, the read-out of the activation state of each segment requires a state tomography, which can be accomplished by using O (2  ) - being the number of qubits -repeated execution of the circuit, voiding the exponential advantage.As a consequence, quantum post-processing of the final state is needed to fully exploit the capabilities of HHL.Quantum clustering techniques could be employed to collect all the active segments and group them into tracks which can then be extracted without full-state tomography.A different approach could exploit the amplitude amplification techniques to increase the amplitude associated with active segments up to a point when the solution can be retrieved with a limited number of shots.
This work opens new possibilities in the application of Quantum Computing to the problem of track reconstruction in HEP.Our new approach, based on the HHL algorithm, has the potential to provide an exponential speed-up over classical counterparts if the conditions discussed above are met.The design of efficient read-out and Hamiltonian simulation subroutines is paramount and will be addressed in future works.

Code availability
The code associated with this work is available at Ref. [44].Any update will be made available as a new version at the same Zenodo DOI.

B Condition number studies
The condition number  of the real symmetric matrix  is defined as the ratio between the largest and the smallest (in modulo) eigenvalues of  ( ) = | max | | min | , and it plays a major role in the complexity and convergence of several numerical linear algebra methods [41].The HHL algorithm run-time complexity has a quadratic dependence on the condition number that can potentially nullify the exponential advantage when ill-conditioned matrices are involved.In this section, we study the behaviour of the condition number ( ) as a function of the number of particles in the event and layers in the detector.
Figure 10a clearly shows that the largest and smallest eigenvalues have no dependence on the number of particles in the event, thus the condition number is constant.On the other hand, studying the behaviour as a function of the number of layers, Figure 10b shows that the largest(smallest) eigenvalue has asymptotically increases(decreases) towards 5 (1).Therefore, the condition number increases with the number of layers but it can be safely bounded from above:  < 5.The bound holds irrespective of the event/detector size, allowing us to ignore the dependence of the run-time complexity on .This property makes our approach to track reconstruction an excellent candidate for an HHL application.

Figure 1 :
Figure 1: Schematic representation of the VELO.The 52 detection modules are arranged on the so-called A-side (left side) and C-side (right side).The interaction region is located in the area with the highest density of modules and depicted in red.

Figure 2 :
Figure 2: (a) Event display with all the possible doublets represented: the ones that are part of a real trajectory are highlighted while the others are left in grey.(b) Diagram showing an example of a track bifurcation.(c) Diagram explicitly representing the angle   between two adjacent doublets   and   .

Figure 3 :
Figure 3: Plot of cos  () for several values of , compared with the step function  (, )

4. 3
Solving the system using quantum computingThe next step consists of solving the linear system (4.2) where the coupling matrix  obtains contributions from each term in the Hamiltonian, including the hyper-parameters =  ang +  spec +  gap ,while the only non-zero contribution to the bias vector b comes from the gap term b = b gap = (1, 1, 1, ..., 1) .

Figure 5 :
Figure 5: Event display of the hits (in red) coming from a collision event in half of the VELO together with the tracks reconstructed using the classical linear Ising solver presented in this work (in black).

Figure 6 :Figure 7 :
Figure 6: Track reconstruction performances in bins of ,  and  T

Figure 8 :
Figure 8: Sparsity of the  matrix (a) as function of the number of particles in the event (generated with a 10-layer detector) and (b) as function of the number of layers in the detector (fixing the number of particles in the event to 25), for a regularised and not regularised Hamiltonian.

Figure 9 :
Figure 9: Sparsity of the  matrix as function of the number of particles in the event, for several values of the angular tolerance term .This has been obtained by processing toy model events with a 10-layer detector.Particles have been generated in a polar angle range of [0, 1.8 • ] to increase the particle density in the event.

Figure 10 :
Figure 10: Condition number , largest eigenvalue  max and smallest eigenvalue  min of  (a) as a function of the number of particles in toy model event (generated with a 10-layer detector) and (b) as a function of the number of layers in the toy model detector (the number of particles in the event is fixed to 25).