Quantum integration of elementary particle processes

We apply quantum integration to elementary particle-physics processes. In particular, we look at scattering processes such as ${\rm e}^+{\rm e}^- \to q \bar q$ and ${\rm e}^+{\rm e}^- \to q \bar q' {\rm W}$. The corresponding probability distributions can be first appropriately loaded on a quantum computer using either quantum Generative Adversarial Networks or an exact method. The distributions are then integrated sing the method of Quantum Amplitude Estimation which shows a quadratic speed-up with respect to classical techniques. In simulations of noiseless quantum computers, we obtain per-cent accurate results for one- and two-dimensional integration with up to six qubits. This work paves the way towards taking advantage of quantum algorithms for the integration of high-energy processes.


Introduction
In this work, the quantum versions of Monte Carlo algorithms are applied to the problem of integrating elementary-particle cross sections.
In particle physics, integration methods and Monte Carlo programs play a very special role as they are the central link between theory and experiment. On the one hand, they allow the encoding of theoretical predictions including higher-order or beyond-the-Standard-Model effects. On the other hand, by generating theoretical events according to the underlying distribution, they allow a one-to-one correspondence with experimental events. Hence theoretical predictions can be directly compared to experimental measurements in order to get insight into elementary interactions.
For collider experiments such as the Large Hadron Collider (LHC), Monte Carlo simulations are crucial as they simulate all the scattering processes generated in the experiment. It means that considerable computing resources are needed and they are expected to increase further [1,2]. For some analysis, the limited Monte Carlo statistics is even becoming a significant source of uncertainty [3,4]. This calls for a continuous improvement of the performance of such Monte Carlo generators.
It appears therefore particularly timely to apply quantum versions of Monte Carlo algorithms to this problem, given the promising advancements in the industry of quantum devices. The core algorithm of interest for us is Quantum Amplitude Estimation (QAE) [5][6][7][8], that was proven to provide a speedup for the integration of probability distributions, by scaling as O (1/M ), where M is the number of (quantum) samples, as opposed to classical integrators scaling as O(1/ √ M ) with M (classical) samples. In the context of high-energy physics, this would translate into a gain in simulation performance, similarly to what was already assessed in other application fields, and specifically in finance [9][10][11][12].
To apply these techniques, classical data must be loaded into a quantum computer, which is a nontrivial task in terms of computational cost. More precisely, the quantum states that correspond to the data have to be prepared. In order to encode the data, several algorithms and techniques are used in the literature [13][14][15][16][17][18]10].
While the present work is the first application of QAE algorithms to integration, there have been numerous applications of quantum computing to other aspects of collider physics. Such applications have been mainly experiment-oriented: pixel images [19], event topologies [20], event classification [21], Higgs analysis [22], background suppression [23], measurement unfolding [24] or jet clustering [25]. Applications to parton-distribution functions (PDF) have also been carried out by several groups [26,27] in a quantum context. In addition, several investigations of quantum parton shower as well as matrix elements evaluation [28,29] have been carried out [30,28,31]. Finally, in Ref. [32] quantum Generative Adversarial Networks (qGAN) [10] techniques have been used for the purpose of data augmentation.
Here, we focus on two representative cases (e + e − → qq and e + e − → qq W processes) to illustrate the application of QAE. A particularly important point is that the functions to be integrated are significantly more complicated to than usual Gaussian or log-normal distributions. These are typically made of trigonometric functions, polynomials, and logarithms (at least for what concerns the lowest order in perturbative theory). We therefore explore two methods to prepare the quantum states according to the underlying distribution, namely the qGAN [10] and an exact loading [33], respectively. Of particular interest for this application, is the ability to provide correct results when restricting the domain of integration. Finally, we carry out oneand two-dimensional integration of cross sections. For the latter case, we devise a method that is extendable to n dimensions while still allowing the arbitrary reduction of the integration domain. In general, the integrations are accurate at the per-cent level with up to six qubits.
The article is organised as follows: in the first part, the method and tools used for this work are presented. The second part deals with two methods to load the probability distributions. The third one focuses on integrating such probability distributions with quantum amplitude estimate methods. Finally, the last sections contains a brief summary as well as some concluding remarks.

Method and tools
In this section, we first recall some general considerations about Monte Carlo integration and explain how we translate it to our problem. Second, we briefly describe the processes under investigation. Finally, the tools used in the next Sections are presented as well as the numerical input of the cross sections.

General considerations
To start, let us recall some basics of particle physics. A Monte Carlo integration aims at estimating the cross section of scattering processes which can be written schematically as where F is the flux factor, dΦ the phase-space factor [possibly including parton-distribution function (PDF)], and |M| 2 the matrix element squared which encodes the quantum mechanical process. In addition, the phase space (also called integration domain below) can be restricted by the use of so called phase-space cuts which is represented in Eq. (1) by Θ(Φ − Φ c ) which we refer to as the domain function in the following. In particular, the integration is performed over variables that allow to describe the full phase space. While these are not physically observable, they allow the full reconstruction of the event kinematic. In the following, the results are only expressed in terms of these variables that serve as proxies for physical ones. In particular, the domain restriction (or event selection) is only applied to the variables of integration. To obtain a physical restriction of the domain of integration, a simple mapping can be performed. In more general terms, any integral I can be cast into the following form The function f describes the probability distribution, while the function g is the integrand function. In the QAE, f is computed classically, while g is represented by means of a quantum operator. For example, in Ref. [10], the g function is a linear function which represents the payoff. In our case, referring to Eq. (1), we take f = |M| 2 . We then take g = Θ(Φ − Φ c ), so that g is a generalised Heaviside function which only takes the value 1 or 0; such a function is sometimes called the indicator function over the integration domain D, and denoted by χ D or 1 D .  Implementing this procedure on a quantum computer involves in general two main steps: the definition of the quantum states and the integration of the probability distribution. The two approaches that we follow in this work are graphically represented in Fig. 1. The first one is based on an exact loading while the second relies on the qGAN to prepare the quantum states. The details of the implementations are explained in the relevant Sections below.

Particle processes investigated
In order to test our numerical approach with the quantum simulations, we have focused on two simple though non-trivial scattering processes e + e − → qq and e + e − → qq W. In particular, we have not considered hadronic processes as these would require the use of parton distribution functions.
The first process is e + e − → qq. In quantum electrodynamics (QED), this process is rather simple. By parametrising the phase space with two angles, the cross section reads: This means that computing such a process (up to an overall normalisation factor) simply amounts to integrate the function 1 + x 2 on the integration domain [−1; 1] while there is no dependence on φ.
The second one is e + e − → qq W. In this case, we have considered the full electroweak Standard Model and not only QED. Due to the three particles in the final state, this process has 5 variables of integration. These can be chosen as two invariants and three angles and the cross section becomes [34] with s Max 1 = (s 2 − M W ) (s − s 2 ) /s 2 and dΦ 3 = ds 2 ds 1 d cos θ 1 dφ 1 dφ 2 . As in the previous case, one of the φ angle is a trivial integration. The main characteristics of the process are summarised in Table 1.

Software
To check our results, we resorted to an in-house Monte Carlo program, that was used for the computation of various high-energy physics processes before [35][36][37][38]. It is based on the MONACO integration routine which is a modified version of VEGAS [39] which is part of the VBFNLO program [40][41][42]. For the matrix elements, we use either analytical expressions or the matrixelement generator Recola [43,44].
Instead, the results in this article are obtained from the open-source distribution Qiskit [45] which is written in Python. Starting from its libraries, we developed our code to load events, build probability distributions, and calculating integrals. The specific functions used are described below in the relevant Sections. With Qiskit, the IBM Quantum Services offer the possibility to run algorithms on simulated quantum computer as well as test some specific configurations on real quantum chips.

Input parameters
In order to ease reproduction of our results, we provide below the numerical inputs of our simulations. For the centre-of-mass energy, we have used √ s = 500 GeV. The electromagnetic coupling is defined with the help of the G µ scheme [46] which leads to The masses and widths of the massive particles are chosen as [47] M OS Z = 91.1876 GeV, All other fermions are considered massless. The pole masses and widths of the heavy gauge bosons are determined from the measured on-shell (OS) values [48] via 3 Definition of probability distributions A necessary step that enables the usage and exploitation of a quantum algorithm, is the encoding of data into quantum states, by means of a quantum circuit. Today, it is not possible to rely on any quantum native techniques like QRAM [49]. Hence, to solve this potential bottleneck, several approaches were proposed in the literature that allow to encode classical data into quantum states [50]. This is particularly important as the approximation introduced in data loading could affect the quality of the integration. This procedure corresponds to the first steps depicted in Fig. 1.
To investigate it, we have classically generated samples (here 10, 000 events) to be loaded into the quantum state. In particular, we have used two methods: qGAN [10] and an exact loading [33], respectively. Both approaches will be outlined and compared in the following. In this Section, we focus on the simple case of e + e − → qq in QED which amounts to integrate +1 −1 dx 1 + x 2 , meaning that the distribution to be loaded is 1 + x 2 .
We discuss first the qGAN method for our application. A Generative Adversarial Network (GAN) [51], in its classical form, is characterised by the interplay of a generative network and a discriminative network to learn the probability distribution underlying the training data [52]. A qGAN has a similar structure, but the generator is a parametrised quantum circuit (PQC) instead of a classical neural network. This way the generator is trained to load a quantum state, approximating the discretised version of the target distribution. As a consequence, this algorithm belongs to the general class of quantum variational algorithms, namely hybrid algorithms that rely on a continuous interaction between a quantum computer and a classical computer. An initial PQC is defined (called ansatz ) and then using a classical optimiser this circuit is trained iteratively. The update of the parameters is driven by the evolution of the related loss function. There are no general prescriptions about the structure of the variational circuit, so that challenges remain, including the trainability, accuracy and efficiency of any variational quantum circuits. For a general overview of variational algorithms, we refer the reader to Ref. [53].
To apply such a method to the integration of the e + e − → qq cross section, we have loaded the normalised distribution 1 + x 2 , which we define as p(x) = (1 + x 2 ) 3 8 such that +1 −1 dxp(x) = 1, using the implementation of the qGAN in Qiskit [10]. 2 The results obtained are presented in Fig. 2 and Table 2 where several loaded distributions are compared against the true value of the distribution for two cases. The first one is the default loading obtained from default qGAN parameters defined in Qiskit, while the second one is an optimised version of the neural network for this particular functional form obtained after several tests of different variational forms and optimiser parameters.
In both cases, five random seeds have been used to estimate the spread of the loading procedure. From this example, it should be rather clear that an optimisation of the neural network in terms of architecture (rotation gates and entanglement gates) as well as parameters tuning is needed and that a default configuration cannot be used for arbitrary distributions. Specifically, from our study, the best entanglement is the circular one and the best results are obtained with a learning rate of 5.10 −4 and 1.10 −3 for the generator and discriminator, respectively. 3 The improvement in the accuracy of the loading can be observed in Fig. 2 Table 2: Comparison of qGAN loading of the normalised 1 + x 2 distribution for the default learning rates and an optimised one. The results are given for 5 different seeds. The minimum, maximum, and average difference per bin with respect to the true value is provided (in per cent). The root mean squared error from the true value is also given.
In our example, the default qGAN can lead to loading errors up to 40% with an average deviation above 10% per bin. In the case of an optimised neural network, the average accuracy of the loading per bin is significantly better and lies around 5%. To measure the quality of the whole loaded distribution, one can revert to the root mean squared error defined as where i denotes the bins, x i the value of the distributions loaded, and µ i the true value of the distribution. Finally, the better behaviour of the tuned qGAN can be also observed in the relative entropy as a function of the time steps in Fig. 2. The relative entropy S is defined as where P and Q are the output distribution of the quantum generator and the discretised version of the target distribution, respectively. The tuned network shows a much smoother converge than the default one. We now turn to the exact loading which is represented in Fig. 3. Such a technique is an analytical way to initialise complex amplitude on qubit register. Being an exact method, the accuracy of the loaded distribution is obviously better than with the qGAN. This is shown quantitatively in Table 3 where the differences per bin are shown for the best bin, the worst, and the average. For the two qGAN cases above, the best seed has been selected. From Table 3, it is rather clear that the exact loading is performing significantly better in this case. In particular, it shows discrepancies from the truth by no more than 2% and is on average around 1%. This order of magnitude should be kept in mind as the best possible loading accuracy with a sample of 10, 000 events.    Table 3: Comparison of qGAN loading of the normalised 1 + x 2 distribution for the default learning rates and an optimised one as well as the direct loading. The qGAN results are the ones of the best seed in table 2. The minimum, maximum, and average difference per bin with respect to the true value is provided (in per cent). The root mean squared error from the true value is also given.
If instead, one loads the exact distribution which is known analytically (as opposed to a sample generated according to it), there is simply no deviation from the true distribution. This aspect is particularly important as in principle, a closed form of the distribution is not necessarily available. It also means that the quality of the exact loading is directly dependent on the statistics of the sample given as input.
In order to represent the target distributions with N bins, one needs to encode a statevector of size N , and this translates into a number of qubits n such that N = 2 n . In other words, the data resolution is a direct consequence of the number of qubits used. Obviously, the computational complexity is related to such n, as well. As far as the qGAN approach is concerned, and discarding the training process that will be discussed in the next paragraph, the pure loading phase requires O(poly(kn)) gates, where k is the number of layers, which is intrisic in the definition of the ansatz. Assuming that k can be kept under control, qGANs become an efficient data loading technique, and preserve the speedup of the Quantum Amplitude Estimation algorithm for integration [10]. Conversely, for the exact loading algorithm, the number of 2-qubit gates scales as O(2 n ) [33].
In the argument above, the training cost of a qGAN is neglected. This is motivated by the fact that the same distribution is typically used for multiple simulations, and in this case the training process is performed once, so that the training time can be seen as a constant. Nonetheless, it is worth saying that the scaling of the training cost when the distribution size grows, is an open question, whose complexity lies in the unpredictable number of epochs needed to achieve training convergence, in the desired level of approximation, and in the different behavior of various optimisers. The interaction of such hyperparameters on small-scale problems is discussed in Ref. [54].
Given the limited amount of qubits in our study, we could not appreciate the benefits of qGANs in terms of scaling, while on the contrary we had to face the learning, possibly the tuning of the network, and the verification of the result. Moreover, the current absence of analytical estimates for approximations induced by qGANs, limits their applicability for computing arbitrary processes in a quantum Monte Carlo program, especially when probability distributions are known from first principle or vast amounts of classical data representing them are available.

Integration of probability distributions
In this section, we exclusively discuss the integration of probability distributions. While we could also use the qGAN loading, we use exclusively the exact loading method in this Section in order to isolate the integration step from the loading one. In particular, we look at the integration of one-and two-dimensional distributions. This procedure corresponds to the core of the work flow depicted in Fig. 1.
Once the target distribution has been loaded into a quantum channel, the integration is performed through QAE. Assuming efficient data loading, the algorithm achieves a quadratic improvement, compared with classical Monte Carlo simulation. QAE is a very interesting and studied quantum algorithm due to its potential application in different fields such as quantum chemistry, machine learning, finance and high energy physics. QAE is a fundamental routine in quantum computing which generalises the idea behind the Grover's search algorithm, and gives rise to a family of quantum algorithms. The basic idea is that given an operator A that acts as where a ∈ [0, 1] and |Ψ 0 and |exΨ 1 are two normalised states. Quantum Amplitude Estimation (QAE) is the task of finding an estimate for the amplitude a of the state |Ψ 1 : This can be achieved by the definition of a Grover's like operator of the form [5]: where S 0 and S Ψ 1 are reflections about the |0 and |Ψ 1 states, respectively, and phase estimation. This formulation represents the canonical version of QAE which is a combination of Quantum Phase Estimation (QPE) and Grover's Algorithm [55]. On one hand, QPE is theoretically able to achieve exponential speedup, like in the famous Shor's Algorithm for factoring [56], on the other hand its practical implementation in terms of qubits and circuit depth represents an interesting challenge in current technological scenario. Removing the dependence on QPE for a QAE-like routine in a simplified version such that it uses only Grover iterations has been largely studied in the literature. Indeed, there exist different implementations, with respect to the original QAE implementation by Brassard et al. [5], such as the Iterative Amplitude Estimation (IAE) version which does not rely on Quantum Phase Estimation (QPE) as defined in Eq. (11). This is the adopted version for this work which can achieve a provable quadratic speedup over classical Monte Carlo simulation, with a desired asymptotic behaviour in its iterative queries to the quantum computer, reducing the required number of qubits and gates [6]. Additional implementations are the Maximum Likelihood Amplitude Estimation [7,8] which limit resorting to expensive controlled operations.

One-dimensional distribution
As mentioned in the previous section, the direct loading adds no approximation to the probabilities given as input. If such probabilities are obtained through sampling, though, they are in turn approximated. This means eventually that the result of the integration will strongly depend on the quality of the input. To illustrate this, we have use the QAE with samples of different sizes: 1000 events (low statistics), 100,000 events (high statistics), 1M events (very high statistics). In particular, we have made used of the Qiskit functions LinearAmplitudeFunction [15,16], EstimationProblem, and IterativeAmplitudeEstimation. The latter implement and improved version of the original QAE method [6]. The results for the different samples and the loading of the exact distribution are compared to the analytical result in Tables Table 5: Asymmetric integration of the normalised 1 + x 2 probability distribution based on samples with different statistics (low, high, and very high) or the exact probability distribution. The results are compared to the analytical result in per cent. The results are obtained for three qubits. The low, high, and very high statistics refer to 10, 000, 100, 000, and 1 million events, respectively.
these Tables and the following ones, δ[%] = σ−σ truth σ truth in per cent, where σ truth denotes the true analytical integration.
It is particularly visible that the quality of the integration is dependent on the statistics used. For 1 million events, the result of the integration is accurate at around the per-mille level. The loading of the exact distribution, on the other hand, is systematically below half a per mille accuracy. In addition, it is worth emphasising that the relevant statistics for the integration precision is not the one of the full sample but of the sample in the integrated region. This is particularly clear in Tables 4 and 5 where the smaller integration domain have a lower accuracy. This holds true also for the the loading of the exact distribution. It is worth noticing that the relative differences with respect to the true values in Tables 4 and 5 do not necessarily display a scaling behaviour according to the statistics. This is due to the fact that the samples are subject to statistical fluctuation and their central value (as opposed to the error) follow a scaling behaviour only on average and not for every single point. These numbers are particularly useful as they provide an estimate of the error which originates from not knowing the original distribution analytically (as in the 2D case below). In the present case, this error is about few per cent.
While in Tables 4 and 5, the limits of integration corresponds to the eight bins (n = 3 qubits give 2 n bins) on the domain [−1; 1], in Tables 6 and 7 the same exercise is performed with this time integration domains that do not fit the binning of the piecewise definition of the function.
In the present case, only the results of the integration of the high-statistics sample as well as the exact result are provided as a function of the number of qubits. In general, one observes that the results are significantly worse than in the previous case. This is simply due to the ill-defined     Table 7: Asymmetric integration of the normalised 1 + x 2 probability distribution based on a 1 million-events samples as well as the exact probability distribution. The results are compared to the analytical result in per cent as a function of the number of qubits. The high statistics refer to 100, 000 events.
value of the distribution between two bins. By increasing the number of qubits, one observes an improvement of the results until the bin edges correspond to the integration boundaries. Once the bin edges fit the integration boundaries, increasing the number of qubits does not lead to any improvement as the distribution is already best defined within the integration boundaries. This implies that when taking the limit of large numbers of qubits, these artifacts disappear.

Two-dimensional distribution
We now turn to the integration of a two-dimensional function for the case of e + e − → qq W. As it can be seen from Eq. (4), the 3-particles phase space requires the integration over 5 variables.
To simplify the problem while keeping it non-trivial, we integrate over the two invariants s 1 and s 2 . To that end, we take a slice in Eq. (4) by setting cos θ 1 = 0, φ 1 = π/2, and φ 2 = π/2. The cross section then becomes with dΦ 3 = ds 2 ds 1 and M = M e + e − →qq W (cos θ 1 = 0, φ 1 = π/2, φ 2 = π/2). The integration of the cross section therefore amounts to integrate over the variables s 2 and s 1 . The integrand is graphically represented on the left-hand side of Fig. 4  integrated are rather different and more complicated than those that have been tested so far such as Gaussian or log-normal distributions.
To encode the multidimensional distribution of Eq. (12) on qubits register, we revert to the same method as in the one-dimensional case by defining it piecewise. To that end, we introduce a mapping from the two dimensional function to one-dimensional. The way the mapping is performed is represented in the right-hand side of Fig. 4. We opt for this solution instead of simply scanning from left to right and top to bottom, in order to ensure that any physically motivated integration domain restrictions can be mapped to the one-dimensional function in a continuous way, hence minimising the error due to interpolations that are not evident nor considered in previous one-dimensional case. We note that this method is fully general and can be extended to n-dimensional integral. It has also the advantage to be fully flexible and allow for the arbitrary phase-space cut in the integrand. For example in Fig. 4, the blue shaded area represents the restriction of the integration domain. In this case, assuming that the first bin of the 1D function [f →f (X)] is in the top left corner and that it is mapped to S X = [0; 16], break points will be introduced at X = 6, X = 10, X = 14, and X = 15.
For the present experiment, we have produced 100, 000 events according to the two-dimensional distribution over the full integration domain. We reduced the latter by setting the maximum value of s 2 to 20, 000 in order to avoid populating bins with very few events. This leads to a total of 97, 581 events.
The results of our experiments are given in Table 8. In this case, we consider two cases of integrand reduction or cuts: S 1 which implement the cut [62, 500; 187, 500] × [5, 000; 10, 000] and S 2 which corresponds to [80, 000; 150, 000] × [5, 000; 10, 000]. The main difference between these two sub-domains is that the boundaries of S 1 fits the edges of the bins of a 4 × 4 grid while the ones of S 2 do not for the first variable. This explains why in Table 8   reproducing the truth (σ = 0.545833717629457) which is here the classical sample. While in the one-dimensional case, each increase in the number of qubits translates into the halving of the bins, it is not the case here. Indeed, going from 4 qubits to 5, only allows to extend the grid from 4 × 4 to 5 × 5. It explains why for S 1 , while increasing the grid and making the binning finner, the accuracy of the integration does not improve. It only becomes perfect again when the binning is again perfectly fitting the boundaries of the integration domain. This is further exemplified with the case of S 2 where the improvement is not uniform when increasing the grid dimension. For this case, the true value of the cross section is σ = 0.7244852993923. This is due, on the one hand, to the fact that the edges of the second dimension are only matched for the cases of the 4 × 4 and 8 × 8 grids. On the other hand, as seen in the one dimensional case, when the domain of integration does not match the piecewise definition of the function, the result of the integration is uncontrolled. In the present case, the interpolation is such that doubling the number of bins in each dimension does not necessarily increase the precision of the integration (grid 4 × 4 vs. grid 8 × 8).
This implies that only a large number of qubits (implying finer bins) can allow a reliable estimate of the integral. In particular, for the present application which the computation of cross sections in collider experiments, the usual standard for Monte Carlo error is to reach a per mille accuracy. With current technology, this goal could be challenging. Nonetheless, we believe that with the advent of machines with 1000 qubits or more 4 , this is perfectly conceivable. Not only a greater number of qubits is needed but also a greater quantum volume [57] that could allow to run QAE on a quantum computer, where further improvements are required, e.g., longer coherence times and higher gate fidelity.
We note in passing that with our method, we could in principle also sample events according to the underlying distribution as done in Ref. [32]. We defer the study of this aspect to future work. In particular, while the integration of probability distributions with QAE methods has shown to provide a quantum advantage, it is not clear yet if such an advantage can also be observed for the sampling of events.

Conclusion and outlook
This work constitutes the first application of Quantum Amplitude Estimate (QAE) algorithm to high-energy physics. To test its feasibility we have checked two non-trivial elementary processes, namely e + e − → qq and e + e − → qq W.
Complex function appearing in elementary scattering processes can successfully be loaded onto qubits consistently with the results of Ref. [32]. To load the functions we have used two methods, namely: the quantum Generative Adversarial Networks (qGAN) [10] and an exact loading [33]. For our purposes, we have found that the latter one is more appropriate due to its versatility and reliability for what concerns application with a small number of qubits. In particular, it does not require any training nor tuning which makes it very easy to use.
In addition, we have successfully used the QAE algorithm for the integration of the two elementary processes in one and two dimensions, respectively. In particular, we have tested the reliability of the integration when restricting its domain of integration, which would correspond to imposing physical event selection in an experiment. To integrate multi-dimensional functions, we have devised a general method which can be extended to n dimensions.
Following this purely numerical strategy requires large number of qubits in order to be accurate. For our application, we have found that QAE provides per-cent accurate results for one-and two-dimensional integration with up to six qubits. The results support the framework where future physical devices will make quantum computing a viable solution for integrating elementary processes in high-energy physics. An increase in the number of available qubits is critical for the practical application to our domain of study. It should be noted here, though, that other issues emerge in the current era of Noisy Intermediate-Scale Quantum Computers [58]. Indeed, additional challenges originate from the imperfection of present hardware construction, from the limited topological connectivity of qubits, and from the inability to put in place full error correction protocols that would require additional qubits and resources. Practical usage of the algorithms shall therefore be validated and perfected also through the execution on real quantum devices.
This work opens new perspectives for the computation of particle processes with quantum Monte Carlo integration techniques. Following the same method, more complicated processes (with higher multiplicities and hadronic processes) can be investigated.