Variational quantum eigensolver for causal loop Feynman diagrams and directed acyclic graphs

We present a variational quantum eigensolver (VQE) algorithm for the efficient bootstrapping of the causal representation of multiloop Feynman diagrams in the Loop-Tree Duality (LTD) or, equivalently, the selection of acyclic configurations in directed graphs. A loop Hamiltonian based on the adjacency matrix describing a multiloop topology, and whose different energy levels correspond to the number of cycles, is minimized by VQE to identify the causal or acyclic configurations. The algorithm has been adapted to select multiple degenerated minima and thus achieves higher detection rates. A performance comparison with a Grover's based algorithm is discussed in detail. The VQE approach requires, in general, fewer qubits and shorter circuits for its implementation, albeit with lesser success rates.

We present a variational quantum eigensolver (VQE) algorithm for the efficient bootstrapping of the causal representation of multiloop Feynman diagrams in the Loop-Tree Duality (LTD) or, equivalently, the selection of acyclic configurations in directed graphs.A loop Hamiltonian based on the adjacency matrix describing a multiloop topology, and whose different energy levels correspond to the number of cycles, is minimized by VQE to identify the causal or acyclic configurations.The algorithm has been adapted to select multiple degenerated minima and thus achieves higher detection rates.A performance comparison with a Grover's based algorithm is discussed in detail.The VQE approach requires, in general, fewer qubits and shorter circuits for its implementation, albeit with lesser success rates.

I. INTRODUCTION
In recent years, there has been a tremendous progress in achieving highly-precise theoretical predictions for particle colliders, which has been possible because of the development of new techniques in Quantum Field Theories (QFT) as well as improved (classical) hardware.For instance, from the theory point of view, the calculation of multiloop scattering amplitudes and Feynman integrals was optimized by applying different sophisticated techniques.Some of these theoretical advancements include Mellin-Barnes transformations [1][2][3][4][5], algebraic reduction of integrands [6][7][8][9][10][11][12][13], integration-by-parts identities [14,15], contour deformation assisted by neural networks [16] and sector decomposition [17][18][19][20], among several other highly efficient methods 1 .Special emphasis was recently put in the development of four dimensional methods [21][22][23], with the purpose of achieving a seamless combination of algebraic/analytic strategies and numerical integration directly in the four physical dimensions of the space-time.
Overcoming the current precision frontier will certainly require to push these methods to their limits.
Future collider experiments require challenging theoretical predictions from full calculations at order n in the perturbative expansion (N n LO), with n ≥ 3.One foreseeable complication is directly related to the appearance of several scales in multiloop multileg Feynman integrals.These objects are analytically known for specific processes and kinematic configurations up to two loops with a limited number of scales and very few examples of three-loop amplitudes are starting to pop-up [24,25].At this point, it is very unlikely that all the required ingredients will become analytically available and the development of novel numerical approaches seems unavoidable.
Given present and future needs, a major progress in solving (at least) two challenges is necessary.On one side, new computational techniques must be developed, exploiting the fundamental properties of QFT and the mathematical concepts behind scattering amplitudes.On the other side, these new methods must be implemented in efficient event generators, capable of overcoming the bottlenecks of the currently available hardware.
Regarding the first challenge, we focus on the Loop-Tree Duality (LTD) formalism [26][27][28][29][30][31][32][33][34][35][36][37], which exhibits very attractive mathematical properties and a manifestly causal physical interpretation.Multiloop scattering amplitudes are transformed in LTD into the so-called dual amplitudes by integrating out one component of each of the loop momenta through the Cauchy's residue theorem.In this way, physical observables are expressed in terms of Euclidean integrals that combine loop and tree-level contributions (as well as renormalization counter-terms) to achieve a fully local cancellation of singularities [38][39][40][41][42][43][44].Besides the possibility of a local regularization, including the simultaneous cancellation of infrared and ultraviolet singularities, the LTD formalism fully exploits causality in QFT.A recent reformulation [45][46][47][48][49][50][51][52][53][54][55][56][57] showed that multiloop scattering amplitudes and Feynman integrals can be represented in terms of a subset of cut diagrams that generalize the well-known Cutkosky's rules [58].In particular, a manifestly causal representation can be directly obtained by decomposing the original Feynman diagrams into binary connected partitions in the equivalence class of topologies defined by collapsing propagators into edges (or multi-edges) [51,59].As a result, a well defined geometrical algorithm is available and a strong connection between causality and directed acyclic graphs can be established [60].
The other challenge is related to an efficient calculation of Feynman integrals, overcoming current hardware limitations.In this direction, the development of novel strategies for classically hard problems based on quantum algorithms (QAs) is gaining momentum across different areas.For instance, there are several ideas that exploit the potential speed-up of quantum computers, such as database querying through Grover's algorithm [61], the famous Shor's algorithm for factorization of large integers [62] or Hamiltonian minimization through quantum annealing [63].In the context of particle physics, QAs are often applied to solve problems related to lattice gauge theories [64][65][66][67][68][69].Recent applications for high-energy colliders include jet identification and clustering [70][71][72][73][74][75][76], determination of parton densities (PDFs) [77], simulation of parton showers [78,79], anomaly detection [80], and integration of elementary particle processes [81].This list is rapidly growing, since QAs are suitable for several uses, especially those involving minimization problems.
With this panorama in mind, the purpose of this article is to explore the application of QAs for unveiling the causal structure of multiloop scattering amplitudes and Feynman integrals, or equivalently acyclic configurations of directed graphs.Our strategy consists in exploiting the properties of the adjacency matrix in graph theory to build a Hamiltonian which weights the cost of different momentum flow configurations.Explicitly, causal configurations are those with minimum energy, so we select them by identifying the minima of the Hamiltonian.This identification is implemented through a Variational Quantum Eigensolver (VQE) [82][83][84], namely a hybrid quantum-classical algorithm that seeks the ground state of a given operator.
The classical problem of identifying or counting directed acyclic graphs is #P-hard, see Ref. [85].Explicitly, given a graph G = (V, E) with V vertices and E edges, there are 2 |E| possible directed graphs that must be checked for cycles.Efficient classical algorithms have been proposed in the literature: e.g., Ref. [86] gives an example of a Fixed Parameter Tractable algorithm, whose scaling depends on a single structural parameter of the graph, see also references therein.This algorithm consists of a Binary Decision Diagram algorithm [87], which performs asymptotically in O(2 p 2 w /4 ) running time per solution, with the path-width p w always smaller than the input size of the graph, becoming closer to it for dense graphs (i.e., maximally connected graphs).Instead, the motivation of our paper is to encode the complexity of the graph by means of a Hamiltonian, i.e. a cost function to be minimized in order to find the directed acyclic configurations.Since VQE is expected to offer a speed-up in minimization problems compared to purely classical algorithms, we explore up to what extent this speed-up remains in the detection of directed acyclic graphs.
The outline of this paper is the following.In Sec.II, we present a brief introduction to LTD, and we offer a description of its manifestly causal representation.Then, in Sec.III, we discuss the geometrical aspects of the causal configurations, presenting explicit examples in App. A. After that, we introduce the Hamiltonian approach in Sec.IV, putting special emphasis in the connection with the adjacency matrix of directed acyclic graphs.We offer a definition of the loop Hamiltonian and its classical reconstruction algorithm in Secs.IV A and IV B, respectively.We discuss subtleties of the encoding of vertex registers in App.B. Then, we discuss the implementation of the VQE in Sec.V, presenting an explicit example for a representative two-loop topology in Sec.V A. The full list of Hamiltonians for the topologies studied in this article is given in App. C. Right after in Sec.V B, we discuss an improved strategy based on multiple runs of the VQE to collect, step by step, all the possible causal solutions.In Sec.VI, we carefully compare the VQE approach with the Grover's based algorithm described in Ref. [60], paying attention to the resources required for a successful implementation in (real) quantum devices.Finally, conclusions and further research directions are presented in Sec.VII.

II. LOOP-TREE DUALITY AND CAUSALITY
To reach highly precise theoretical predictions, it is necessary to deal with multiloop scattering amplitudes and the corresponding Feynman integrals.In the Feynman representation, the most general L-loop scattering amplitude with P external particles is given by where q i with i ∈ {1, . . ., n} are the momenta flowing through each Feynman propagator, G F (q i ) = (q 2 i − m 2 i + ı0) −1 , and N represents a numerator, whose specific form depends on the topologies of the different diagrams, the interaction vertices and the nature of the particles that propagate inside the loops.Regarding the momenta of the internal particles, they are linear combinations of the primitive loop momenta associated to the integration variables in the loop, ℓ s with s ∈ {1, . . ., L}), and the external momenta, p j with j ∈ {1, . . ., P }.Given a Feynman diagram, we can group the different internal momenta into sets: two lines are said to belong to the same set if their propagators involve the same linear combination of primitive loop momenta [46].This is useful for achieving an efficient classification of diagrams according to their causal structure [47][48][49].
The denominator of Eq. ( 1) characterizes the singular structure of the scattering amplitude, and singularities provide important information to simplify loop calculations.In this direction, the Loop-Tree Duality (LTD) [26,27] makes use of Cauchy's residue theorem (CRT) to reduce the dimensionality of the integration domain, enabling the possibility of transforming it into an Euclidean space.By means of an iterated application of CRT, we can get rid of one integration variable per loop: this is equivalent to say that the LTD representation of a L-loop amplitude is obtained by cutting (or setting on-shell) L internal lines.The effect of cutting a line is the modification of the infinitesimal complex prescription, leading to the so-called dual propagators [26,27].It is important to highlight that this modified prescription plays a crucial role to preserve causality, as we will discuss later.
The traditional formulation of LTD leads to the socalled dual representation, in which there is one contribution for each possible connected tree obtained by setting on-shell L internal propagators.If we study the structure of the denominators of each dual term, we immediately realize the presence of divergences that do not correspond to any physical configuration.These unphysical singularities are spurious, and they vanish when we add all the dual contributions together.This is the so called manifestly causal representation within LTD [46,88].
Before moving on, let us recall one crucial fact about the Feynman propagators.If we define the positive on-shell energy then the propagator can be written as This decomposition suggests that each Feynman propagator encodes a quantum superposition of two on-shell modes, with positive and negative energy respectively.The so-called causal representations are obtained by consistently aligning these modes.Furthermore, this superposition also motivates the identification of each internal propagator with a qubit, as we will explain later.
As carefully discussed in Refs.[46,89], the calculation of the nested residues within the LTD formalism leads to a manifestly causal representation of multiloop multileg scattering amplitudes.Thus, it can be shown that Eq. ( 1) is equivalent to with i,0 , h σ(i) = ±1, and the integration measure in the loop three-momentum space.It turns out that Eq. ( 4) only involves denominators with on-shell energies, added together in samesign combinations within the so-called causal propagators, where σ(i) includes information about the partition p of the set of on-shell energies and the corresponding orientation of the energy components of the external momenta, k p,0 .Each causal propagator is in a one-to-one correspondence with any possible threshold singularity of the amplitude F , which contains overlapped thresholds.These are known as entangled causal thresholds: σ(i) indicates the set of causal propagators that can be simultaneously entangled.The set Σ in Eq. ( 4) contains all the combinations of allowed causal entangled thresholds, which involves overlapped cuts with aligned momenta flow.
An important advantage of the representation shown in Eq. ( 4) is the absence of spurious nonphysical singularities.In fact, since causal propagators involve same-sign combinations of positive on-shell energies, the only possible singularities are related to IR/UV and physical thresholds.More details about the manifestly causal LTD representation and its benefits can be found in Refs.[47,48,50,51].

III. GEOMETRIC INTERPRETATION OF CAUSAL FLOWS
In order to identify the causal configurations of a multiloop Feynman diagram in a quantum device, it is necessary to translate the problem into a suitable language.The geometrical formulation of the manifestly causal LTD representations [51] provides a clear and intuitive interpretation.In the following, we briefly describe the most relevant features of this formalism, recalling some basic concepts and definitions as presented in Refs.[51,59,60].
Any scattering amplitude can be represented starting from Feynman diagrams built from interaction vertices and internal lines (or propagators) connecting those vertices.Regarding causality, those lines that connect the same vertices can be merged into a single edge, leading to reduced Feynman graphs because the only allowed con-figurations which are causal are those in which all the momentum flows of these propagators are aligned in the same direction.Then, reduced graphs are built from vertices and edges, as described in Refs.[50][51][52], and their causal structure turns out to be equivalent to that of the dual representation of the original Feynman diagrams [49].The number of loop integration variables in the LTD representation given by Eq. ( 4) is directly related to the topological independent loops present at the level of the Feynman diagram.However, the reduced Feynman graph has a fewer number of graphical loops, also called eloops, since bunches of propagators connecting the same vertices were collapsed to a single edge.Furthermore, it turns out that vertices, edges and eloops are the only required ingredients to obtain the causal representation of any multiloop scattering amplitude [50,51].
Specifically, given a reduced Feynman graph, the associated causal propagators 1/λ ± p correspond to the set of connected binary partitions of vertices.Also, they can be graphically interpreted as lines cutting the diagrams into two disconnected pieces, in a full analogy with Cutkosky's formulation [58].Then, the representation in Eq. ( 4) can be constructed from all the possible causal compatible combinations of k causal propagators, the so-called causal entangled thresholds.The number k is known as the order of the diagram and defined as k = V − 1 [46,50,51].It naturally induces a topological classification of families of Feynman diagrams, or equivalently, reduced Feynman graphs.The family with k = 1 is known as Maximal Loop Topology (MLT) [46] and only involves two vertices; a Next-to-Maximal Loop Topology (NMLT) has three vertices, and so on.
It was shown that causal entangled thresholds can be identified by imposing geometrical selection rules.These rules are deeply connected to the algebraic formulation presented in Refs.[50,90], and they establish that: 1.When considering all the k causal thresholds that can be simultaneously entangled, all the edges can be set on shell simultaneously.This is equivalent to impose that the causal entangled thresholds depends on the on-shell energies q (+) i,0 of all the edges.2. Two different causal propagator λ i and λ j can be simultaneously entangled if they do not cross each other.In other words, this means that the associated partition of vertices do not intersect or are totally included one in the other.
3. The momentum of the edges that crosses a given binary partition of vertices must be consistently aligned, i.e., momentum must flow from one partition to a different one.
A careful explanation of these causal rules and their geometric interpretation was presented in Ref. [51], including pedagogical examples and practical cases of use.By imposing these three rules on the set of all the possible combinations of k thresholds, we construct Σ, namely the set of all causal entangled thresholds leading to Eq. ( 4).
Furthermore, as previously shown in Ref. [60], the third condition is equivalent to ordering the internal edges in such a way that the reduced Feynman diagram is a directed acyclic graph.This means that there cannot be closed cycles, allowing the information to propagate consistently from one partition of the diagram to the other; this condition is necessary for having a causal-compatible partition.For this reason, given a reduced Feynman graph, it is crucial to efficiently detect all the associated directed acyclic graphs, since they are a vital ingredient to reconstruct the LTD causal representation.To clarify this discussion, we present an explicit application example in App. A.

IV. HAMILTONIAN FORMALISM FOR CAUSAL-FLOW IDENTIFICATION
The first proof-of-concept of a quantum algorithm for Feynman loop integrals was presented in Ref. [60], managing to successfully unfold the causal configurations of multiloop Feynman diagrams.The implementation follows a modified Grover's quantum algorithm to identify the presence of directed acyclic configurations, representing the causal solutions, over different subloops of the multiloop topologies.
The aim of this Section is to define a Hamiltonian whose ground state corresponds to a superposition of all the acyclic configurations of a given multiloop Feynman diagram2 .We first discuss the construction of such Hamiltonian based on direct inspection of the multiloop topology, and then we discuss a more general approach based on the adjacency matrix of a graph, an object that concentrates all the information regarding the orientation of edges and the vertices that they connect within a reduced Feynman graph (as defined in Sec.III).We also give a more formal presentation, relying on concepts and notations from graph theory.

A. Cycle detection via Hamiltonian optimization
Let us consider a generic undirected graph G = (V, E), with V a set of distinct vertices and E ⊂ V × V a set of edges (also known as links).Given the graph G, one can build a set D G of 2 |E| possible directed graphs where we associate a specific direction to each link.We will denote the subset of directed acyclic graphs as In order to formulate this problem as a Hamiltonian optimization, we first fix the encoding for all the possible (classical) solutions as elements of the computational basis as follows.We start by establishing a conventional orientation to each edge, selecting a specific element from the set of directed graphs G 0 ∈ D G .This set is associated to the state |00 . . .0⟩ ∈ H E , mapping each qubit to an edge as done in Ref. [60]: the single qubit state |0⟩ e for the edge e ∈ E corresponds to a link with the same orientation of e in the conventional graph G 0 , while the state |1⟩ e would correspond to the opposite orientation.The problem can then be cast in the search of all the possible states of the computational basis of edge orientations associated to directed acyclic graphs D G , whose span generates a Hilbert space After these definitions, we need to identify the set Γ G0 of all the possible oriented cycles.In this notation, the subscript G 0 means that we will use the reference graph to specify the orientations of those cycles.Here, cycle refers to simple cycles (i.e. a closed loop in which a vertex appears only once).If the set Γ G0 of directed loops is already known, one can build the so-called loop Hamiltonian, which is a Hamiltonian on the space H E , diagonal in the computational basis described above.We proceed by penalizing oriented cycles using diagonal projectors, i.e.
where s(e; G 0 ) evaluates to 0 or 1 according to the orientation of e ∈ γ relative to the conventional graph orientation G 0 , and the π symbol identifies the projector operators on single edge registers, with Z the Pauli matrix in the z-axis.
To provide an explicit example, let us consider Fig. 1.It depicts a graph G 0 (i.e., G equipped with a conventional orientation), where one can identify by direct inspection the set of all the possible directed loops: where the presence or absence of the bar identifies the orientation (i.e., s(e, G 0 ) = 0, s(ē, G 0 ) = 1).Then, the loop Hamiltonian has the following form: Notice that a term in H G acts trivially on a non-indicated edge, i.e.
with I i corresponding to the identity operator on H i .This loop Hamiltonian will evaluate to 0 for any solution in D G of graphs without oriented cycles (i.e., Ker(H G ) = H E ), while the other directed graphs are associated to a strictly greater than 0 integer eigenvalue.Moreover, the same oriented graph could have more than one closed cycle, and this Hamiltonian counts the number of these oriented cycles.

B. Classical algorithm for reconstructing the loop Hamiltonian
The adjacency matrix of a graph provides a powerful tool to further characterise the associated multiloop topology, and it can be used when the loop set is not straightforward to compute.Classically, the adjacency matrix A = ||a|| of a directed graph of N = |V | vertices (labelled in an arbitrary manner) is defined such that a i1i2 = 1 (i 1 , i 2 = 1, . . ., N ) if there is a directed edge from i 1 to i 2 , and 0 otherwise.Note that a cycle of length ℓ along the vertices i 0 , . . .
, the following properties are equivalent: 1. the graph of which A is the adjacency matrix is acyclic; A is equivalent to a triangular matrix, where equivalence has a precise algebraic meaning, corresponding to the existence of a relabelling; 4. the eigenvalues of A are all zero, which directly implies tr(A k ) = 0 for any positive integer k.
We can exploit these equivalent properties to build suitable candidates for Hamiltonians.For this purpose, let us consider a graph G = (V, E) and a corresponding conventionally oriented graph G 0 = (V, Ē).Then, we define an extended Hilbert space of vertices and edges H V E ≡ H V × H E , and we build a linear operator acting on it and encoding the information about the oriented adjacency matrix: where σ ± = (X ± i Y )/2 are the ladder operators.This construction can be interpreted as a sum of hopping terms from the tail to the head vertices of each oriented edge e, with the edge projector operators π s e keeping track of the orientation relative to the conventional graph.Thus, joining two adjacent edges would give nontrivial contributions only if the orientation is the same for both edges (i.e., if the head of one is joined with the tail of the other), and that is ensured by the properties of the ladder operators on the corresponding vertex register.Therefore, the n-th power of the A operator will have non-zero diagonal terms in the vertex side only if all the ladder operators are completely contracted and only projectors π s appear.On the edge side, these terms identify oriented cycles (or disjoint products of oriented cycles), so that a possible loop Hamiltonian, whose kernel is the span of directed acyclic graph states H E , can be written as follows: where the trace runs only over the vertex space, and M G is the maximal length of loops in G, i.e., M G = |V |.An oriented cycle going across j vertices corresponds to a term in this sum when n = j.The loop Hamiltonian built in this way has the same terms as the loop Hamiltonian in Eq. ( 7), plus some terms representing disjoint products of oriented cycles (e.g., two oriented cycles going respectively across j 1 and j 2 vertices and corresponding to a term in the sum above when n = j 1 + j 2 ).Also, it includes coefficients in front of the different products of projectors π, which are related to topological properties of the graph under consideration (see, e.g., Eq. ( 16) below).
In any case, the ground state will be associated with the subspace of directed acyclic graphs.For the practical implementation, we will rely on the normalized version of H G where all the coefficients in front of the π operators are exactly 1.
It is worth mentioning that, to derive a loop Hamiltonian, in principle any positive-definite function of A can be used.For example, would have a different spectrum, but exactly the same kernel as the loop Hamiltonian in Eq. ( 13).Given the relation tr(A j ) = N i=1 λ j i , where λ i (i = 1, . . ., N ) are the eigenvalues of the adjacency matrix, one sees that cycles are deeply connected to spectral properties of graphs.So, the problem we have at hand consists in knowing whether any of these eigenvalues is non-zero, in which case cycles would be present.Since the Hamiltonians defined above are related to the eigenvalues, they are independent of the specific labelling of the graph in consideration.
To conclude this section, note that a discussion about efficient encoding methods for the vertex space is presented in App.B, having in mind the possibility of a Hamiltonian reconstruction by means of a purely quantum algorithm.Also, we would like to mention that the number of terms in the loop Hamiltonian does not depend on the number of external legs.As shown in App.C, more complex topologies require more terms because each term of the Hamiltonian penalizes the presence of cycles.Thus, more loops leads to more possible cycles to be checked, and this translates into more terms in the Hamiltonian, but more external legs do not generate additional cycles.Computing the Hamiltonian in a fully quantum way, as we discuss in App.B, could allow to speed-up the codification of the characteristics of the graph for the posterior minimization with a quantum device.

V. VQE IMPLEMENTATION AND NUMERICAL RESULTS
We now discuss the application of the Variational Quantum Eigensolver (VQE) algorithm [82,83] to find directed acyclic graphs given a fixed topology.Fixing a positive orientation for an arbitrary link (the first one for example), the |E| − 1 remaining qubits encode the direction of the other links, i.e., the states considered have the form |ψ⟩ ⊗ |0⟩ with the restricted loop Hamiltonian H E/{e0} ≡ H E | s(e0,G0)=0 acting only on |ψ⟩.This is possible without loss of generality, since the solutions to the unrestricted problem can be found by adding a copy of the solutions to the restricted problem with all links flipped: where X ⊗(|E|−1) |ψ j ⟩ corresponds to applying a tensor product of NOT gates to the first |E| − 1 bits to invert the directions of the links/edges described in |ψ j ⟩.
To explain how the proposed algorithm works, let us consider again the topology shown in Fig. 1, where |E| = 5.The unrestricted Hamiltonian obtained from Eq. ( 13) is where we point out the presence of nontrivial coefficients in front of each term, that results from having identical terms in Eq. ( 13).As explained before, the ground state will have always 0 energy and it will be associated to the set of directed acyclic graphs, independently of the spectrum of H E .For the sake of simplicity, from now on, we will work with the normalized Hamiltonian, eliminating the positive coefficients seen in Eq. ( 16).Then, we can fix the orientation of the first qubit in agreement with the initial choice which implies with |ψ⟩ ∈ {|ψ ′ ⟩ ⊗ |0⟩}.This is equivalent to set π 0 0 ≡ I 0 and remove all the terms proportional to π 1 0 .Explicitly which acts over the reduced edge space.The computational basis of the restricted problem represents the 2 |E|−1 = 16 possible directed graphs, of which only 9 are acyclic (i.e., the solutions to our minimization problem).
In practice, all the 9 solutions of the reduced problem can be found by applying the VQE where we minimize ⟨H E/{e0} ⟩ with any of the Hamiltonians built as discussed in the previous sections.As argued above, the space of solutions corresponds to the intersection between the canonical basis and the eigenspace associated with the eigenvalue 0. A single run of the VQE converging to ⟨H E/{e0} ⟩ = 0 (within precision), will, in general, result in a superposition of states of the computational basis in Ker H E/{e0} .In analogy to the Grover's approach, we collect solutions whose probability appears above a certain threshold λ.However, in this case, there is no guarantee that all possible solutions are found in the superposition of a single run.In the following, we will explain a possible strategy to overcome this difficulty, based on running the VQE multiple times.

A. Example: two-eloop case
Now, we concentrate again our attention to the explicit example of Fig. 1.As a first step, we proceed to translate Eq. ( 18) to the basis of Pauli operators using the definition of the projectors given in Eq. ( 8).The result is Then, we setup the VQE using the package qiskit.algorithms 3together with the function RealAmplitudes to generate the Ansätze.In this 3 http://qiskit.org/experiment, we use the optimizer COBYLA [91], with maxiter={1024,10240}, and the simulation is implemented in the aer simulator with a fixed number of shots (by default shots=1024).Notice that iterations refers to the number of times that the classical optimization algorithm is called.This has to be distinguished from shots, which is the number of times the quantum circuit is executed and measured.In any case, we observe that more iterations allow to achieve a better convergence of the optimization whilst more shots are associated with a more precise determination of the mean value of the Hamiltonian.
In Fig. 2, we show the probability of the different configurations contributing to the approximate ground state.The yellow bars are associated to acyclic (or causal) configurations, that the algorithm is aimed to identify.For the present topology, there are 9 acyclic configurations in the reduced space, namely: In this particular example, the simulation leads to the identification of 3 solutions, for both maxiter= 1024 and maxiter= 10240.The extremely good approximation to the exact energy of the ground state implies in practice that there was no contamination with noncausal configurations.However, this is not guaranteed, as we will see in subsequent examples along the article.Furthermore, we notice that a single VQE run (i.e., the set of quantum and classical steps discussed so far, which includes a certain number of shots and iterations) is generally not enough to collect all the solutions of the problem.Since there are several local minima, the minimization algorithm could be stuck in one of them and achieve convergence without exploring the full set of possibilities.The situation is even more complicated because there are multiple degenerated global minima, whose contribution to the ground state can be re-scaled arbitrarily without changing the energy.This implies that the random selection of the initial points might produce very different probability distributions from run to run: in other words, if we re-run the experiment, it is highly probable that the distributions in Fig. 2 change substantially, selecting a different subset of causal solutions.For these reasons, we define the success rate as which also takes into account the possibility of misidentification of causal states.For the particular example presented in Fig. 2, we have r success = 0.89 and r success = 0.78, respectively.If we increase the number of iterations, we find the results in Fig. 3: we have r success = 0.56 for maxiter= 102400 and r success = 1 for maxiter= 1024000 iterations.In the first case, one of the states, |0110⟩, is noncausal, which leads to the drastic decrease of the success rate 4 .
With the purpose of reducing the probability of collecting incorrect states, we can select only those states contributing to the approximated ground state with a probability above a certain threshold λ.Since the VQE is optimized to find the states with lower energy, it is reasonable to assume that there are much higher probabilities to include acyclic rather than noncausal configurations in the approximated ground state.Thus, if the solution found by VQE is then we propose the threshold λ to be i.e. the average minus a term proportional to the standard deviation of the probabilities |c j | 2 , or the minimum of such probabilities.This value of λ corresponds to the red line in Figs. 2 and 3. We would like to emphasize that this definition is empirical, and motivated by the exploration of the numerical results.Also, if we want to increase the certainty that only causal states are selected, we can take higher values of λ.This could eventually lead to a lower success rate (since more states are rejected).
To balance this situation, we can run multiple times the VQE procedure and collect solutions in each step until we fulfill a certain convergence criterion, discussed in the next section.
To conclude this section, we would like to mention 4 Alternative definitions for the success rate could be examined.
Here we stick to this one in order to penalize the misidentification of causal states.
that different versions of the optimizers and Ansätze were used.In particular, besides COBYLA, we performed the optimization with SPSA [92] and NFT [93], achieving similar results.However, COBYLA turned out to run faster on our computer 5 .Regarding the Ansatz, we relied on RealAmplitudes and EfficientSU2, both circuits being part of qiskit.circuit.library.Again, the success rates for the topology A of Fig. 1 were rather similar, although RealAmplitudes involves half of the parameters of EfficientSU2.However, the use of RealAmplitudes manifests into slightly higher values of the Hamiltonian mean value.Moreover, as we will see in the next section, it was necessary at times to switch to EfficientSU2 to avoid missing solutions, while also avoiding misidentifications.

B. Multiple-run VQE approach
In order to find the complete set of solutions, one can set up a sequence of VQE runs where the Hamiltonian used contains additional terms which suppress the solutions found at each previous step of the sequence.The proposed theoretical algorithm is the following: 1. We start with the reduced Hamiltonian H (0) = H E/{e0} and apply the VQE algorithm using an Ansatz.We can choose as an initial point the uniform superposition of all the configurations, i.e.
with {|ϕ i ⟩} the canonical basis of the reduced edge space H E/{e0} .FIG. 3: Approximation of the ground state found by VQE with COBYLA optimizer, corresponding to the two-eloop diagram in Fig. 1.We present the probability of the different states contributing to the ground state after 102400 (left) and 1024000 (right) iterations, with 1024 shots.The red line is a threshold to limit the contamination by noncausal configurations.

The VQE minimization is executed till it finds a
state with energy E 1 < μ or until the maximum number of iterations, maxiter, is reached 6 .In case the energy condition is not fulfilled, we run the VQE again till it converges to E 1 < μ (which is guaranteed by the presence of, at least, 1 causal configuration).The approximation to the ground state is expressed as 3. We pick up the subset S 1 = {|ϕ j ⟩ | |c (1) j | 2 > λ}, which corresponds to elements with a high probability of belonging to Ker H E/{e0} .Then, we define a penalization term to exclude these potential solutions from the next run.Explicitly, we consider with b (1) l positive coefficients to be set appropriately.Steps 1-3 define the single-run VQE, used in Sec.V A.

We define the modified reduced Hamiltonian
H (1) = H (0) + Π (1) , (27) and re-run the VQE.If it does not find any state such that E 2 < μ within the given allowed number of iterations, then the algorithm stops.The set of solutions remains S = S 1 .

If a state
such that E 2 < μ is found, then we pick up the set j | 2 > λ}, build the associated penalization term Π (2) and add S 2 to the set of solutions.

Repeat
Step 4 and 5 till the algorithm stops when the mean energy is greater than μ.
This approach is rather general and can be modified by adjusting several parameters: we can choose different energy thresholds μ for each step, different Ansätze, probability thresholds, modify the number of shots and/or classical iterations for each VQE run, etc..Although in the next section, we will discuss some of these observations, we defer a more detailed study for future works.Also, we would like to point out that this procedure does not guarantee collecting all the solutions.Due to the presence of so-called barren plateaus [94], the VQE could get stuck around local minima which do not belong to the true ground-state.This is a shared drawback of several minimization strategies (including VQE), in which the presence of large regions in the parameter space with rather small gradients prevents an efficient optimization [95].In the next section, we discuss improvements to overcome barren plateaus and increase the success rate, keeping a very low error rate.Also, in the VQE literature [96], it is normally highlighted the importance of choosing a suitable Ansatz that properly represents at least the low-lying eigenstates of the Hamiltonian.In the context of quantum chemistry, this is usually possible since the Hamiltonians represent physical systems and analytical approximations or mathematical properties of the solutions are known.By choosing an appropriate Ansatz it would be possible to achieve a faster convergence and partially avoid the barren plateaus.However, in our case, the situation is not straightforward and identifying an optimal Ansatz is far from trivial.The choice of the Ansatz in general should have parametrized rotational gates and one or more layers of entanglement to be able to connect the qubits, if the entanglement entropy of the system is not zero.So, for finding convergence in this particular problem, it is enough to have a generic circuit structure.Our approach can be considered as a first exploration, and further investigations in this direction will be performed in the future.Now, let us consider again the example of Fig. 1, keeping the two test scenarios (i.e. the ones with 1024 and 10240 iterations, respectively).The starting point is the set of solutions collected from Fig. 2. We set the detection threshold λ according to Eq. ( 23), and we iterate the VQE excluding the solutions found in the precedent step, as in Eq. ( 27).In Figs. 4 and 5, we present the probabilities of the different states contributing to the approximated ground state, for 1024 and 10240 iterations respectively.In each run, we collect the solutions above the threshold only if the energy is below μ = 0.1.In both cases, we manage to successfully identify the 9 causal states, after 5 and 6 runs, for maxiter= 1024 and maxiter= 10240 respectively.No misidentification took place because the selection threshold turned out to be rather rigorous.By lowering λ, the procedure converges faster, but there is an increased risk to pick up a noncausal state.
Before concluding this section, let us briefly discuss the complexity of this approach.Again, the multi-run VQE strategy involves a hybrid classical-quantum setup.On one side, there is a polynomial computational cost for calculating the Hamiltonian, which is O(|V | 4 ) based on Eq. ( 13), i.e., distinct |V | rows and |V | columns have to be multiplied, at a cost of |V | per multiplication, and this has to be repeated a total of |V |−1 times when n = |V |.Then, this Hamiltonian has to be evaluated in a quantum system, leading to a potential speed-up since all the expectation values (i.e., one for each term; see, e.g., Eq. ( 19)) are simultaneously computed on a superposition of many configurations [84].After that, a classical parameter optimization is performed with standardised algorithms (such as COBYLA or NFT) and the procedure is iterated.

VI. COMPARING HAMILTONIAN OPTIMIZATION WITH GROVER'S ALGORITHM
After a careful study of the Hamiltonian minimization through VQE, we carry out in this Section a performance comparison with Grover's based algorithm [60].We rely on the set of six representative topologies shown in Fig. 6.These topologies allow us to understand how the quantum circuits' complexity scales with the number of vertices and edges.In Appendix C, we show the explicit form of the Hamiltonians in H E for the aforementioned topologies, so that a classical verification of the procedure explained in Sec.IV can be implemented.
For each topology, we study the set of causal/acyclic configurations comparing with the output of a classical algorithm 7 .We used the success rate defined in Eq. ( 21) to quantify the accuracy of the identification of causal states.Then, we compare: 1. the number of qubits required to run the algorithm; 2. the quantum depth of the transpiled circuits; 3. the total number of executions needed to get the full solution.
The first point is important since most of the current quantum devices have a (relatively) small number of qubits.
Of course, we expect that this limitation will be overcome in the (near) future, but still it is a tangible drawback.
In the same direction, physical bottlenecks to ensure the stability of the qubits prevents to implement very deep circuits.So, from the point of view of the hardware requirements, avoiding a fast growth of these two aspects is crucial to ensure the feasibility of the quantum approach.Finally, the third point implies more execution time, since re-running the circuit implies preparing the state and repeating the measurements.However, this is not a hard limiting factor as the previous two properties, although it penalizes the performance.A further comment is needed to explain the definition of quantum depth and executions in this comparison.Regarding the quantum depth, it is related to the number of operations or quantum gates that are involved within the circuit.Each specific device has a set of fundamental gates, and any other gate/operator is defined by a combination of the fundamental set.As a consequence, the real depth of the circuit is hardware-or simulator-dependent, so we need to fix a convention to compare Grover's and VQE approaches.Since the calculations are implemented in Qiskit, we use the function qiskit.compiler.transpileand aer simulator as target device.Then, we estimate the quantum depth through the .depth()method.For the case of Grover's algorithm, the depth depends not only on the complexity of the topology but also on the number of iterations: for the topologies studied in this article, we can achieve a perfect reconstruction with only one iteration [60].Note that in this case, the complexity of the oracle and the diffuser are included.The quantum depth for VQE is related to the complexity of the Ansatz and the Hamiltonian.In fact, for the VQE, we estimate the depth by applying the .depth()method to the circuit composed by one term of the Hamiltonian applied to the Ansatz.So, in general, VQE involves shorter circuits, as discussed in e.g.[84].
About the concept of executions, it is directly related to the number of times that we need to measure the output.In the case of Grover's algorithm, this quantity exactly agrees with the number of shots [60,97].However, VQE requires a deeper analysis.We need to remember that the VQE pipeline mixes a classical and a quantum algorithm.The quantum circuits depend on parameters, that are adjusted after performing the measurement in order to minimise a certain cost function [84].For each choice of parameters, we define a quantum circuit and execute s shots (controlled by the parameter shots within the backend).Then, we allow the algorithm to modify the parameters r times (until the cost function evaluates under a certain threshold), which is what we call an iteration (controlled by maxiterations within the classical optimizer routine).So, each run of the VQE requires s × r measurements.On top of that, the multiple-run VQE algorithm explained in Sec.V must be repeated k times, adding in each repetition a certain number of penalization terms to the Hamiltonian.As a consequence, we will need to make s × r × k executions to complete the identification of the winning states (or causal configurations).For these reasons, in general, VQE involves a larger number of executions compared to Grover's approach.
In Tables I and II, we present the number of qubits, executions and the quantum depth of the circuits in-volved in the implementation of the Grover's algorithm and the single-run VQE approach, respectively, for all the topologies shown in Fig. 6.In the case of Grover, the success rate is 100 % because we are using the optimal number of executions explored in Ref. [60], which guarantees in practice a perfect casual/noncausal discrimination, in spite of requiring more qubits and larger coherence.For the single-run VQE, we estimate the success rate by using Eq. ( 21).We rely on COBYLA optimizer and RealAmplitudes Ansatz which represents less quantum resources with respect to Grover's implementation but involves more shots, classical iterations and the efficiency is noticeably smaller.This is because, in the experiments that we executed to obtain the inferior values in the success rate ranges of Tab.II, we are collecting all the states in the approximate ground-state as solutions of our problem.However, if the energy is not 0 (or compatible with 0 within machine error), this ground-state might be contaminated with noncausal configurations; this will reduce the success rate since we are penalizing the misidentification with Eq. ( 21).This illustrates the need for defining a proper threshold to retain only the true minima from the approximated ground-state.Then, we tested the performance of the multiple-run VQE algorithm.We considered different scenarios: 1. Setup 1 : COBYLA optimizer with 1000 iterations, 1000 shots, RealAmplitudes with a random parameter initialization as Ansatz and the threshold λ given by Eq. ( 23).
2. Setup 2 : NFT optimizer with 1000 iterations, 1000 shots, RealAmplitudes with a random parameter initialization as Ansatz.The selection threshold for the i-th run is given by    6, using a single run.We indicate the number of executions as shots times iterations.The success rate rsuccess range contemplates the worst case scenario (no threshold, all the states are selected) and the best one (optimal threshold for distinguishing causal from noncausal).Thus, these results illustrate the need for defining a threshold to avoid misidentifications.
FIG. 6: Representative topologies used for the performance comparison of Grover's and VQE implementations.We considered the following configurations: (A) two eloops with four vertices; (B) three eloops with four vertices; and (C) four eloops with five vertices.Also, we include the complete set of four-eloop six-vertex topologies: (D) t-channel; (E) s-channel; and (F) u-channel.
compatible with 0, then all the states are solutions).
3. Setup 3 : NFT optimizer with 1000 iterations, 1000 shots, EfficientSU2 using the approximated solution found in the previous run as initial point of its subsequent run.The selection threshold is defined as in Setup 2, and we repeat 3 times the execution of the VQE if E i > μ in the i-th run.
Setup 1 corresponds to the configuration tested in Sec.V B. We do not present a table with explicit results for complex topologies (such as D, E of F) because the success rates are rather low, even below O(10 %).This is due to several withdraws of this naive implementation, namely some limitations of COBYLA optimizer and the initial Ansatz.For this reason, we tried other configurations, tweaking the parameters of the multi-run VQE to achieve higher success rates.In Tabs.III and IV, we present the results for Setup 2 and Setup 3 respectively.We observe that NFT identifies approximations to the real ground-state with lower energy but with less configurations in superposition w.r.t.COBYLA.In fact, with COBYLA the identification rate could reach up to O(50%), but there are several incorrect states, which leads to very low success rates, see Eq. ( 21).In contrast, the combination of the NFT optimization and the improved selection threshold from Eq. ( 29) has an extremely low misidentification rate: in all the topologies tested, the error rate was absent, but each run collects a smaller number of solutions.The difference in the performance is due to the nature of the optimization algorithm.NFT is an optimization method specifically designed for quantum-classical hybrid algorithms based on parameterized quantum circuits, whilst COBYLA is a more general approach envisaged for problems where the derivative of the objective function is unknown.Since it assumes certain properties of the cost function, NFT allows to reach configurations with lower energy by performing fewer measurements.Even if we were unable to test the performance of these two approaches in real quantum devices, we expect NFT to perform better for the multi-run VQE since it is also claimed to be more effective in the presence of noise (as well as SPSA).Still, in general, it is not clear from first principles which is the best choice of optimizer and, most of the times, this is found by empirical exploration.Still, in presence of barren plateaus, the implementation in Setup 2 stops before identifying all the possible solutions, since the algorithm gets stuck around local minima.To overcome this limitation, in Setup 3, we consider a more flexible Ansatz and select an optimized value for the initial point.Also, we implemented a routine to randomly kick the parameters when the energy is higher than μ = 1 (which means a non-global minimum was found).In this way, we avoid a premature interruption of the multi-run VQE, which keeps collecting solutions, with null error rate.

VII. CONCLUSIONS AND OUTLOOK
The computation of accurate theoretical predictions for current and future high-energy particle colliders requires efficient ways to deal with scattering amplitudes at higher orders in the perturbative expansion.In the present article, we put forward a quantum strategy to re-cast the selection of causal configurations of multiloop Feynman integrands into a minimization problem.For this purpose, we exploit the geometrical approach to the causal Loop-Tree Duality (LTD) representation of Feynman integrals [51] and relate causality with the identification of directed acyclic graphs [60].
We rely on the concept of adjacency matrix in graph theory to define the so-called loop Hamiltonian, whose different energy levels correspond to graph configurations with different number of cycles.By construction, the ground-state is directly associated to the subspace of states related to directed acyclic graphs.In this way, given a multiloop Feynman diagram, we fix a reference orientation of the internal propagators and then define the associated loop Hamiltonian.By construction, this loop Hamiltonian has mean value 0 when evaluated on acyclic states.Therefore by minimizing the energy, we automatically detect the subset of causal configuration.
After explaining the construction of the loop Hamiltonian, we proceed to solve with quantum minimization algorithms.We use the Variational Quantum Eigensolver (VQE) within the Qiskit framework to test the different configurations.In particular, we develop a multi-run VQE strategy that allows to achieve higher detection rates.The VQE approach is compared against the Grover's detection procedure implemented in Ref. [60].Our results indicate that even if VQE requires less quantum resources than Grover's algorithm (i.e. it involves fewer qubits and shorter circuits), still a larger number of executions must be performed in order to achieve high success rates.However, since VQE is a hybrid algorithm that employs classical optimization, it allows to use classical hardware and save quantum resources.Given the current limitations in quantum hardware, VQE approach allows to study complex multiloop topologies that would require a prohibitive number of stable qubits within Grover's implementation.
In conclusion, along this article, we have presented a proof-of-principle for the detection of causal configurations using a loop Hamiltonian implemented in VQE.We have demonstrated the feasibility of our approach, paving the road to achieve an efficient causal reconstruction of multiloop Feynman amplitudes through quantum algorithms using less hardware resources.We have presented a multiple-run VQE algorithm with much better success rates than a naive single-run VQE.We stress that we avoid the misidentification of solutions, reaching null error rates in the simulators.Still, the success rates achieved are not as satisfactory as the ones achieved with a Grover-based algorithm and further investigation is required to develop a better strategy through quantum minimization.In any case, we successfully tested, for the first time, a VQEbased strategy to solve a problem with highly-degenerated ground-states, which constitutes a hard stress-test for this kind of minimization algorithms.A promising road to enhance the success rate and reduce the number of VQE runs is Conditional-Value-at-Risk (CVaR) [98], which has already shown improved performance on certain practical problems, such as the allocation of flight gates [99,100].Also, different strategies must be explored regarding optimal state-preparation, since several results point towards a significant impact of the state-preparation in the performance of the minimization [101,102].In particular, it would be highly interesting to test our algorithms in quantum annealers, which could lead to an improved performance and much higher success rates with a lower consumption of quantum resources.We defer these improvements and studies for future investigations.
Examples of partitions that are not causally entangled are: 1. {1, 2, 5} is not allowed because e 4 remains uncut, i.e. condition 1 is not fulfilled.
3. {1, 2, 3} is not allowed because e 2 can not be oriented in a consistent way.Explicitly, if we take λ + 1 , then we have λ − 2 and λ − 3 ; but λ − 2 implies S(e 2 , G 0 ) = 0 and λ − 3 implies S(e 2 , G 0 ) = 1 (following the notation introduced in Sec.IV).These 3 cases are depicted in Fig. 7. Finally, we would like to close this Appendix by reminding that conditions 1-3 are equivalent to the following ones: 1.The edges must be oriented in such a way that cyclic loops are not present, i.e. we need to identify all the possible directed acyclic graphs obtained from the original reduced Feynman graph.
2. Then, we need to dress the acyclic graphs with k = V − 1 simultaneous causal propagators, in such a way that the conditions 1-2 from the list in Sec.III are fulfilled.
This was the strategy first used in Ref. [60] to bootstrap the causal representation, and justifies the importance of developing an efficient algorithm to detect directed acyclic configurations for a given topology.

Appendix B: Compact encoding of vertex registers
Having in mind the extension of the VQE algorithm presented in this paper to complex multiloop topologies, it is relevant to think about suitable scaling techniques.As we already know from Ref. [60], the number of qubits required to implement the computation scales very fast with the complexity of the underlying graphs.For this, it is important to identify strategies that allow us to reduce the qubit consumption.
One source of proliferation of qubits is associated with the description of the vertex space V .The encoding used for the vertex registers in the discussion above is of the one-hot kind, i.e., each classical vertex state for a specific vertex with label v ∈ {0, . . ., |V | − 1} is associated to a unique string with all 0's except for a single 1 on the v-th position.For example, the vertices in Fig. 1 are associated to the following elements of the computational basis: With this encoding, the oriented link e (i,j) : v i → v j is represented by an hopping term σ − vi σ + vj π 0 e (i,j) .However, it is simpler (and more scalable) to consider a compact encoding using q = ⌈log 2 (|V |)⌉ qubits instead of the |V | necessary in the one-hot case8 .Returning to the example in Fig. 1, one can use a possible compact encoding with just q = 2 qubits, as the following: The hopping terms for the adjacency matrix can be computed, as for the one-hot encoding, by looking at the oriented link e (i,j) : v i → v j , which will contain a mix of ladder (at least one) and projector operators.Using Fig. 1, the oriented link e (1,0) : v 1 → v 0 is associated with the non-Hermitian operator After computing all the terms for the operator representing the adjacency matrix A, one can build the loop Hamiltonian by tracing out the vertex register as in Eqs. ( 13)-( 14) (or any equivalent positive semi-definite matrix with H E as kernel).This is guaranteed by the fact that the only terms contributing to H G after computing the trace are those in A n that do not contain ladder operators (i.e., diagonal terms), because these can be obtained only if the products of the n hopping terms connect the same vertex as input (e.g., |v i ⟩⟨v i | V ).In other words, the terms present in A n are those representing a closed chain, starting and ending in the same vertex.

Appendix C: Loop Hamiltonians for representative multiloop topologies
In this appendix, we include the formulae for the loop Hamiltonians for the multiloop topologies described in Fig. 6. Results for diagram A were discussed in detail in Sec.V A, so here we present only expressions for the remaining topologies.Explicitly, we have for diagram B

3 FIG. 1 :
FIG. 1: Example of a diagram with two eloops and four vertices representing external particles (Topology A).Edges are directly numbered, and vertices are explicitly indicated as vi.

1 :FIG. 2 :
FIG.2: Approximation of the ground state found by VQE with COBYLA optimizer, corresponding to the two-eloop diagram in Fig.1.We present the probability of the different states contributing to the ground state after 1024 (left) and 10240 (right) iterations with 1024 shots.The red line is a threshold to limit the contamination by noncausal configurations.

4 FIG. 4 :
FIG.4: Iterated application of the VQE with COBYLA optimizer, corresponding to the two-eloop diagram in Fig.1.We present the probability of the different states contributing to the ground state after 1024 iterations.States in yellow (light-blue) represent causal (noncausal) solutions.The red line is a threshold to limit the contamination by noncausal configurations, described in the text.The algorithm converges after 4 runs.

FIG. 5 :
FIG.5: Iterated application of the VQE with COBYLA optimizer, corresponding to the two-eloop diagram in Fig.1.We present the probability of the different states contributing to the ground state after 10240 iterations.States in yellow (light-blue) represent causal (noncausal) solutions.The red line is a threshold to limit the contamination by noncausal configurations, described in the text.The algorithm converges after 5 runs.
if E i > 10 −8 , and λ i = 0 otherwise (i.e. if the mean energy of the approximated ground-state is

TABLE I :
Performance of the Grover-based algorithm applied to the representative multiloop topologies depicted in Fig.6

TABLE II :
Performance of the Hamiltonian minimization through VQE applied to the representative multiloop topologies depicted in Fig.

TABLE III :
Performance of the Hamiltonian minimization through VQE applied to the representative multiloop topologies depicted in Fig.6, using the iterative approach explained in Sec.V. We used 1000 iterations and 1000 shots in each VQE run, with the Setup 2.Even if the success rate was below 100 %, no misidentification took place.

TABLE IV :
Performance of the Hamiltonian minimization through VQE applied to the representative multiloop topologies depicted in Fig.6, using the iterative approach explained in Sec.V. We used 1000 iterations and 1000 shots in each VQE run, with the Setup 3. A major improvement compared to Setup 2 was reached, keeping a null error rate.