QFitter -- A Quantum Fitting Framework Applied to Effective Field Theories

The use of experimental data to constrain the values of the Wilson coefficients of an Effective Field Theory (EFT) involves minimising a $\chi^2$ function that may contain local minima. Classical optimisation algorithms can become trapped in these minima, preventing the determination of the global minimum. The quantum annealing framework has the potential to overcome this limitation and reliably find the global minimum of non-convex functions. We present QFitter, a quantum annealing method to perform EFT fits. Using a state-of-the-art quantum annealer, we show with concrete examples that QFitter can be used to fit sets of at least eight coefficients, including their quadratic contributions. An arbitrary number of observables can be included without changing the required number of qubits. We provide an example in which $\chi^2$ is non-convex and show that QFitter can find the global minimum more accurately than its classical alternatives.


I. INTRODUCTION
Determining model parameters by adjusting a theoretical hypothesis until it best fits experimental data is commonly referred to as fitting. Fitting algorithms use optimisation methods which aim to find the model's parameters that result in the smallest value of a characteristic function, like the negative log-likelihood function, χ 2 , which incorporates all experimental observables of interest simultaneously. As such procedure allows an interpretation of experimental data in terms of a given theoretical model, fitting algorithms are at the core of phenomenological studies of new physics.
Effective Field Theories (EFTs) provide a modelindependent framework to tension experimental data with the imprints of new physics scenarios, assuming that any new particles are significantly heavier than the Standard Model particles. The free parameters of an EFT are the coefficients of the local operators that appear in its Lagrangian. They are known as the Wilson coefficients. The EFT formulation ensures that the Wilson coefficients can capture any effects of the heavy particles at low energies compared to their masses.
Because of this, EFT fits have become an essential tool in interpreting experimental data in the current context, in which no resonant new physics has been observed at the LHC. Subsets of the coefficients of the EFT for the non-linear realisation of the electroweak symmetry have been fitted to the relevant experimental data in Refs. [1,2]. Similarly, fits for the linear realisation have been performed in Refs. [3][4][5][6][7][8][9][10], for example.
If non-linear contributions of the Wilson coefficients are considered in the fit, the χ 2 function to be minimised may develop local minima which differ from the global one. Furthermore, local minima can also occur when some nuisance parameters are included. This may prevent classical optimisation algorithms from reaching the values of the Wilson coefficients that fit the observed experimental values optimally. Indeed, such algorithms usually obtain the solution by incrementally improving a given point in parameter space, which can lead to them becoming trapped in a local minimum.
In this work, we construct a formulation of EFT fitting calculations that can be implemented in these devices. An arbitrary number of observables can be included in the fit without changing the required number of qubits. Furthermore, general polynomial dependence of the theoretical predictions on the Wilson coefficients can be implemented in this formulation, with non-linear terms encoded by auxiliary qubits. In practice, we find that current quantum annealing devices can be used to perform fits with at least 8 Wilson coefficients and quadratic dependence of observables on them. Using practical examples, we show that the quantum approach generates the best-fit estimator for the coefficients accurately and more consistently than classical methods in a non-convex problem.
The rest of this paper is organised as follows. In Section II, we present the QFitter method for fitting EFT coefficients using quantum annealing. In Section III, we provide the results of applying QFitter to three different problems, including a fit in which the χ 2 function to be minimised is non-convex. We compare QFitter to several classical approaches. We give our conclusions in Section IV.
where s is a free parameter that can be controlled externally; A(s) and B(s) are continuous functions such that A(0) > 0 = B(0) and A(1) = 0 < B(1); and H 0 is a Hamiltonian whose ground state is known in advance. The solution to the problem is obtained by preparing the system in the ground state of H 0 and changing s continuously from s = 0 at an initial time t i to s = 1 at a final time t f . The function s(t) described by the time evolution of s is known as the schedule. The adiabatic theorem ensures that, if the change in s is sufficiently slow, the system is likely to end in the ground state of the target Hamiltonian H 1 .
A concrete realization of this method is transversefield quantum annealing, which has been implemented in real-world devices. In this realization, the system can be viewed effectively as a collection of qubits (that is, quantum systems with two independent states), with the Hamiltonians H 0 and H 1 given by where h i and J ij are adjustable parameters, andσ i x,z are the x, z Pauli matrices acting on the ith qubit. To perform a computation using transverse-field quantum annealing, one needs to encode its result as the ground state of the Ising Hamiltonian H 1 . The D-Wave devices provide a physical implementation of this setup, in which not all J ij couplings can be set to a non-vanishing value. In the state-of-theart Advantage system architecture, there are more than 5000 available qubits, but each one is coupled only to 15 others. To find the ground state of Ising models with a higher degree of connectivity, several qubits are chained together with large coupling to act as a single qubit with more connections. The mapping between the abstract Ising model Hamiltonian to be minimized and the one implemented in the physical device is known as an embedding.
Once an embedding has been found, and the schedule s(t) and h i , J ij parameters are set, the annealer is typically run several times to reduce the effects of external noise. Then, the final state with the least energy obtained from the different runs is selected. The number of runs is referred to in the context as the number of reads.

B. QUBO formulation
The eigenstates of the quantum Ising Hamiltonian H 1 correspond to the states of its classical analogue, whose Hamiltonian is with the σ i being classical variables taking the values σ i = ±1. Thus, the problem solved by the transversefield quantum annealers can be viewed equivalently as finding the set of values for the σ i such that H classical is minimized: This problem can be solved both using quantum annealing and classical algorithms, such as simulated annealing. Quantum annealing has been shown to be more consistent in finding the ground state of some non-convex functions [14]. In Section III, we will compare the performance of quantum annealing to several classical methods for EFT fits.
A useful reformulation of the classical Ising Hamiltonian minimization problem is obtained by making use of the binary variables τ i = (σ i + 1)/2, whose possible values are 0 or 1. In terms of them, the problem can be expressed as the minimization of a homogeneous quadratic polynomial: This is known as a Quadratic Unconstrained Binary Optimization (QUBO) problem. We will refer to the function L = τ i Q ij τ j to be minimized in such a problem as the loss function.
where C −1 is the inverse covariance matrix. The maximum-likelihood estimator for the coefficients can be obtained by minimizing χ 2 .
In any EFT, the theoretical predictions O will thus be polynomials in c. For simplicity, we restrict ourselves here to quadratic polynomials although the method that we present can be straightforwardly extended to higher degrees. Nevertheless, quadratic polynomials are sufficient for most current EFT applications. In particular, dimension-8 and dimension-6 squared contributions in the Standard Model EFT can be dealt with in this setting.
We thus have a χ 2 function to be minimized that is a quartic polynomial in real variables c i . To turn this into a QUBO, we need to re-formulate it as a quadratic polynomial in binary variables. Binarization can be achieved by means of a binary encoding [15,29]: with L i , U i and N being fixed parameters, and the c a (c) into a quadratic function of the binary variables, which means that χ 2 is a quartic polynomial in them.
The reduction from a fourth-degree polynomial to a second-degree one can be made by means of auxiliary binary variables c In this way, the c i c j factor that appears in the last term of Eq. (6), can be written as a linear function in the binary variables Replacing Eq. (7) in the linear term of Eq. (6), and Eq. (8) in the quadratic one, O (th) a (c) and χ 2 become linear and quadratic in the binary variables, respectively.
In order for this procedure to work, the constraints c need to be enforced somehow. This can be done by means of a constraint Hamiltonian: a quadratic function P (x, y, z) of binary variables x, y, z that achieves its minimum P = 0 if and only if xy = z. In particular, we use the function Now, to construct the loss function, we add together χ 2 , viewed as a quadratic function of the binary variables, as described above, and the constraint Hamiltonians for all the relations between them: When the coefficient λ is large enough, the minimum of L will correspond to the minimum of χ 2 over the set of values of the binary variables that satisfy the constraints. This concludes the re-formulation of the problem as QUBO: the maximum likelihood estimator for the Wilson coefficients can be obtained by minimizing a function quadratic function of binary variables L.
The QUBO problem we have obtained can be solved by the usual methods, including quantum and simulated annealing. We will show in Section III that, in practice, quantum annealing is the most effective of the two for this purpose, especially when χ 2 is a non-convex function of the coefficients c i . The result of any of these methods will be a set of values of the binary variables c  7) and (8).
The number n bin of binary variables employed in this formulation is highly relevant for the practical applications using quantum annealing devices since, together with the connectivity of the QUBO, it will control the number of qubits in the final embedding. n bin is controlled by the number M of Wilson coefficients involved in the fit and the number N of binary variables used in the encoding for each of them. In general, one has The first term is the number of binary variables required to encode the coefficients. The second one is for their products: there are N (N + 1)/2 different combinations of the form c i c j , and at most M 2 different choices for m and n given (i, j). In many situations, only certain combinations of coefficients appear in the quadratic terms in Eq.
Concerning the connections, an entry of the Q matrix of the QUBO is non-vanishing if and only if there is one observable to which the corresponding coefficients (or product of coefficients) contribute simultaneously. Since there are observables to which many coefficients may contribute, the degree of connectivity is typically high. Finally, notice here that the number of observables does not play any role in determining n bin . That means that the number of observables that can be included in the fit is not limited by the available number of qubits when the QUBO is embedded in physical annealing devices.

D. Zooming
As discussed in Section II A, the total number of available qubits in state-of-the-art devices is in the order of several thousand. However, for problems with a high degree of connectivity that require embedding several physical qubits per abstract qubit in the original problem, the effective number of available qubits is much lower, reaching around 100 for fully connected QUBOs. In an EFT fit, the number M of coefficients is fixed, so the only way to reduce the number of qubits, according to Eq. (11) is to reduce the number N of binary variables per coefficient. Thus, the bound in the number of qubits translates into a bound in N , which limits the precision that can be achieved for the best-fit values of the coefficients.
We overcome this limitation using a zooming process, similar to the one presented in Ref. [30]. The idea is to perform the fit in several steps, which we call epochs. Each epoch consists of a quantum annealing run for a refined version of the QUBO problem of the previous run, in which the range [L i , U i ] for each coefficient c i is updated. Given a parameter 0 < f < 1, the zoom factor, the update rule for the lower and upper limits of the range is where c i denotes the value of the corresponding coefficient obtained from the previous epoch. The effect of this transformation is to reduce the length of the [L i , U i ] range by a factor f , while centering it around c i . The zooming process thus allows achieving any desired precision at the price of introducing a classical update step between quantum annealing runs. We stress that its use is necessary for embedding in current quantum annealing devices because of the limited amount of available qubits. In future devices with a larger quantity of qubits, it might be possible to use a larger number N of binary variables per coefficient, thus making it possible to reach the relevant precision in a single quantum annealing run.

A. EWPO fit
As a first and simple test of the method, we perform a fit of the S and T oblique parameters to the electroweak precision observables. We use the set of such observables defined Ref. [31], with the values, uncertainties and correlations given in Refs. [32][33][34][35][36]. The S and T parameters can be defined in this context as where v is the Higgs vacuum expectation value, g and g are the gauge coupling constants for the SU (2) and U (1) factors of the electroweak symmetry, and c φW B and c φD are the Wilson coefficients for the following dimension-6 operators in the Standard Model EFT: Here, φ is the Higgs doublet, W a µν and B µν are the SU (2) and U (1) field-strength tensors, D µ is the corresponding covariant derivative, and σ a are the Pauli matrices.
We choose the following hyperparameters for the fit: N = 5 binary variables per coefficient, ranges for them given by c φW B ∈ [−0.005, 0.005] TeV −2 and c φD ∈ [−0.03, 0.03] TeV −2 , no zooming, 100 reads, and a linear anneal schedule s(t), with annealing time t f − t i = 50 µs. A visual representation of the Q matrix of the corresponding QUBO problem is shown in Fig. 3. We implement this in the D-Wave Advantage system4.1 quantum annealer. The total computing time spent in the quantum annealer for the fit is thus 5 ms. Since there is no zooming procedure involved, the fit is performed purely quantum mechanically in this case, with no classical updating steps. As the result, we obtain the central values S = 0.05 and T = 0.09, close to the values obtained in similar fits in the literature [37].
As the next step, we generate the profile of the minimal ∆χ 2 = χ 2 − χ 2 min as a function of the parameter S. To do so, we perform the fit while keeping the value of S fixed. The ranges of for the coefficients are   in a range around its central value. We use the same procedure to produce the profile for ∆χ 2 as a function of T . The results are shown in Fig. 1. From these profiles, we can derive the following 1σ-intervals for the parameters: which are again close to those obtained in Ref. [37]. The total quantum annealing time employed in generating the profiles is 0.2 s.

B. Higgs fit
We now use our method to perform a larger fit, including 8 Wilson coefficients and quadratic terms in the theoretical predictions for observables. The Lagrangian we consider is a subset of the dimension-6 Standard Model EFT Lagrangian, given by: where y t and y b are the top and bottom Yukawa couplings, m W is the mass of the W boson, and g S is the strong coupling constant.
We use the set of observables described in Ref. [6]. They correspond to projections of Higgs production and decay processes at the High-Luminosity LHC. We consider only inclusive measurements in gluon fusion, vector boson fusion, and the production associated with an electroweak gauge boson, two top quarks or a jet. The decays are into the following final states: γγ, W + W − , ZZ, µ + µ − , τ + τ − and bb.
We compute the theoretical predictions for the observables, up to quadratic order in the coefficients. We find that c H , c W , c HB and c HW have quadratic contributions to some observables. We only consider contributions of the form c 2 i for observables where a linear parametrization results in a poor description of the parameter dependence 1 . The only coefficients that have non-negligible contributions at the quadratic level are c W , c HB and c HW . We assume that the measured values of the coefficients are the same as their (dimension-4) Standard Model values. We use the QUBO formulation of the fit to generate the individual ∆χ 2 profiles for the coefficients.
There are several choices of hyperparameters that give rise to similar results. However, we find a greater consistency among different trial runs with the following set of values: N = 2 binary variables per coefficient, 20 zooming epochs, a zoom factor of f = 80%, 200 reads in each epoch, and t f − t i = 100 µs of annealing time for each read, with a linear anneal schedule s(t) ∝ t − t i . The ranges we use for the coefficients are centred around the origin (L i = −U i ), with the following upper limits: The matrix Q for the corresponding QUBO problem is shown in Fig. 3. Again, we use the D-Wave Advantage system4.1 quantum annealer to perform the fits. The total quantum annealing time per fit is 0.4 s. We show our results in Fig. 2. For the coefficients c g , c γ , c H , c u3 and c d3 , we obtain a ∆χ 2 profile with quadratic shape, which is symmetric under c i → −c i , as only have linear contributions to the observables. For c W , c HB and c HW , we observe an asymmetry under c i → −c i , generated by the corresponding quadratic terms.

C. Non-convex Higgs fit
Finally, we consider a modified scenario in which the χ 2 function has a local minimum close to c i = 0 and a (much lower) global one away from it. We thus have to minimize non-convex loss functions, a problem in which the quantum approach might provide an advantage over its classical counterparts [14].
To generate this setup, we consider a scenario in which several observables have been measured to be away from their Standard Model values. The observable that plays the most important role in generating the non-convexity is µ pp→W H→W Zγ , which we set to −0.3, with total uncertainty to be ∆µ pp→W H→W Zγ = 0.01. For the rest of the introduced deviations from the Standard Model values, we keep the projected experimental uncertainties while setting the measured values as follows: • For gluon fusion, vector-boson function, ttH and jH processes with W W or ZZ in the final state, we take µ i = 0.8.
• For any process with Zγ is in the final state, we set µ i = −6.
• For all the processes with Higgs production in association with a vector boson, we choose µ i = 0.2.
We keep the remaining observables in the fit at their Standard Model values. Concerning the coefficients, we only allow c HW , c W , c g and c γ to vary, while fixing the rest of the coefficients to zero.  Fig. 3.
We embed this QUBO in the D-Wave Advantage system4.1 quantum annealer, fixing c HW at 20 equally-spaced values in its range to generate the ∆χ 2 profile for it. The results are given in Fig. 4. This suggests that χ 2 indeed has 2 local minima, one at c HW −0.01 and a lower one at c HW −0.05. We check that the gradients are small at these two points and that the Hessian is definite positive. As a quantitative measure of the smallness of the gradients, we find that the following equation is satisfied: We conclude that we have found an excellent approximation to two of the local minima of χ 2 . As shown in Appendix A, the form of the χ 2 function that we consider here implies that we can have at most two local minima, so we have identified all of them. With χ 2 ≥ 0, this implies that the minimum at c HW −0.05 is the global one.
We now compare different methods for finding the maximum-likelihood estimations of the four coefficients. This corresponds to a search for the global minimum of χ 2 . The classical methods that we consider are: • The MIGRAD algorithm provided by the Minuit code [38]. This type of variable-metric algorithm works by incrementally improving an approximation of the error matrix and the best-fit point.
• Simulated annealing. This method also proceeds by iterativaly improving a point in parameter space. In each step, the algorithm moves to a random nearby point if the new point has a lower χ 2 . If the new point has a higher χ 2 , it moves to it with probability e −∆χ 2 /T , where ∆χ 2 is the difference in χ 2 between the new and the old points, and T is a parameter that decreases along the run, known as the temperature. The algorithm explores the vicinity of local minima while changing T depending on the number of failures to find a better minimum in the previous steps. This may allow the algorithm to escape local minima in some cases. We chose the initial value of T to be 10 6 with a minimum value of 10 −6 , an adaptive speed of 1 and a step size of 0.01. The maximum number of steps is chosen to be 10 5 .
• Finally, we also consider the same algorithm as for the quantum annealing approach, with the QUBO formulation and a zooming process, replacing the quantum annealing runs with simulated annealing. In this case, the parameter points to be updated by the algorithm are the sets of 0 or 1 values of the binary variables. The random nearby point is obtained by flipping a random variable. The number of steps is usually measured here in terms of sweeps, which correspond to as many updating steps as binary variables are present in the problem.
We refer to the methods not using the QUBO formulation as standard-formulation methods. We show a comparison of the best-fit values of the coefficients and χ 2 obtained from all the methods in Table I. The classical methods are initialized to a point with c W = c g = c γ = 0. Inspecting the χ 2 values given in the Table I, one can see that classical algorithms that start at the point c i = 0 typically get trapped in the local minimum nearby. This problem is not present in the quantum annealing approach. The minimum with χ 2 = 135 obtained by Minuit for initial point c HW = −0.05 corresponds to the one from the quantum annealer. The small differences in the minimum value of χ 2 and the best-fit parameters are because we have neglected some quadratic contributions in the formulation of the QUBO problem, which has no consequences on the general shape of the χ 2 function, with a global minimum around c HW = −0.05 and a local one close to c HW = 0.
The QUBO-formulation simulated annealing results depend strongly on the schedule for the temperature. For large starting temperatures, there is no dependence on the starting point. The results also vary considerably for fixed annealing parameters from run to run. We find the most consistent results with a schedule that exponentially increases the inverse-temperature parameter β from 1 × 10 −5 to 10 in 1 × 10 6 steps, performing one sweep per step. The results fall into three classes: (A) those on the "wrong" side of the barrier, with χ 2 4000; (B) those on the correct side, with χ 2 < 1000; and (C) those for which the zooming gets stuck at χ 2 > 10 7 . We perform 40 runs and find that 13% of the runs belong to class A, 82% to class B, and 5% to class C. The results in class (B) also present a considerable variation. We show an arbitrary example in Table I.
In Table I, we also provide the total time spent performing the optimization for each method. Again, we run Minuit, the standard-formulation simulated annealing and the QUBO simulated annealing in an Apple M1 processor. Since quantum annealing is performed on a dedicated device, the numbers cannot be compared directly. However, we note that quantum annealing requires orders of magnitude less time to perform this task.

IV. CONCLUSIONS
We have presented QFitter, a quantum annealingbased method for fitting EFT coefficients to experimental measurements. The χ 2 is encoded as a QUBO problem which can be directly embedded in the currently-available quantum annealers. The required number of qubits depends on the number of coefficients to be fitted, the nonlinear terms in their contributions to observables (which require auxiliary qubits), and the precision to which they are to be determined. The number of observables included does not affect this.
Since physical annealers only provide a limited amount of qubits, the practical implementation of QFitter can only be done for a limited number of coefficients, precision and non-linearities. We have used a zooming algorithm, in which the precision is increased iteratively through several annealing runs, to overcome this limitation partially. With this setup, we have found that fitting problems involving at least eight coefficients and their quadratic dependencies can be embedded in current quantum annealing devices.
Finally, we have tested the performance of QFitter with three examples. The first two, the EWPO and the Higgs fit, involve a convex χ 2 function. The quantum approach gives comparable results to the classical ones here. We have then modified the χ 2 for the Higgs fit to make it non-convex. By comparing with several classical algorithms, we have found that the quantum one is the one that ends in the global minimum most consistently with a considerable gain in processing time.