Digital-analog quantum computation with arbitrary two-body Hamiltonians

PHYSICAL REVIEW RESEARCH 6, 013280 (2024) Digital-analog quantum computation with arbitrary two-body Hamiltonians Mikel Garcia-de-Andoin ,1,2,3,*Álvaro Saiz ,4,5,†Pedro Pérez-Fernández ,5,6Lucas Lamata ,4,6 Izaskun Oregi ,2and Mikel Sanz 1,3,7,8 1Department of Physical Chemistry, University of the Basque Country UPV/EHU, Apartado 644, 48940 Leioa, Spain 2TECNALIA, Basque Research and Technology Alliance (BRTA), Astondo Bidea Edificio 700, 48160 Derio, Spain 3EHU Quantum Center, University of the Basque Country UPV/EHU, Barrio Sarriena s/n, 48940 Leioa, Spain 4Departamento de Física Atómica, Molecular y Nuclear, Universidad de Sevilla, Avenida de la Reina Mercedes s/n, 41012 Sevilla, Spain 5Departamento de Física Aplicada III, Universidad de Sevilla, Camino Descubrimientos s/n, 41092 Sevilla, Spain 6Instituto Carlos I de Física Teórica y Computacional, Universidad de Granada, Avenida de Fuente Nueva s/n, 18071 Granada, Spain 7IKERBASQUE, Basque Foundation for Science, Plaza Euskadi 5, 48009 Bilbao, Spain 8Basque Center for Applied Mathematics (BCAM), Alameda Mazarredo 14, 48009 Bilbao, Spain (Received 15 December 2023; accepted 16 February 2024; published 14 March 2024) Digital-analog quantum computing is a computational paradigm which employs an analog Hamiltonian resource together with single-qubit gates to reach universality. Here, we design a new scheme which employs an arbitrary two-body source Hamiltonian, extending the experimental applicability of this computational paradigm to most quantum platforms. We show that the simulation of an arbitrary two-body target Hamiltonian of nqubits requires O(n2) analog blocks with guaranteed positive times, providing a polynomial advantage compared to the previous scheme. Additionally, we propose a classical strategy which combines a Bayesian optimization with a gradient descent method, improving the performance by ∼55% for small systems measured in the Frobenius norm. DOI: 10.1103/PhysRevResearch.6.013280 I. INTRODUCTION When quantum computing was originally proposed [1,2], it was envisioned as a way of simulating the dynamics of a quantum system employing another controllable system. This set the foundations of what we now call analog quantum computing (AQC) [3]. A different approach was introduced when Deutsch et al. proposed the concept of a quantum gate [4], which finally led to the digital quantum computing (DQC) paradigm. One of the main advantages of AQC is the robustness of the simulation. Quantum control techniques have been developed in the last decades, providing a high fidelity and further protecting the dynamics against different sources of errors [5,6]. Despite their robustness, AQC is strongly limited by the dynamics of the system, making it difficult to implement most dynamics of interest. In contrast, DQC is performed through the sequential application of quantum gates in a discrete manner, mimicking classical computations. It is proven that any unitary can be decomposed with arbitrary precision *[email protected] †[email protected] Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. in terms of single-qubit gates (SQGs) and at least one entangling two-qubit gate (TQG) [7]. One of the main features that DQC provides is the possibility of applying quantum error correction (QEC) techniques [8]. In the current noisy intermediate-scale quantum (NISQ) era [9], the qubits and the gates available are noisy and prone to errors. Thus, the only hope of reaching fault-tolerant quantum computing is through the application of sophisticated QEC techniques [10] once we fulfill the requirements for the quantum threshold theorem [8,11,12]. The digital-analog quantum computing (DAQC) paradigm was proposed as a way of combining the robustness of AQC with the versatility of DQC [13,14]. The main idea behind DAQC is employing the natural interaction Hamiltonian of a system as an entanglement resource. By alternating the evolution under this Hamiltonian (analog blocks) and the application of SQGs (digital blocks), one can simulate an arbitrary target Hamiltonian. Here, we can distinguish two kinds of approaches. If the interaction Hamiltonian is turned off during the application of the digital blocks, we call this approach a stepwise-DAQC (sDAQC) circuit. Otherwise, if for practical purposes the system Hamiltonian is always on, and the SQGs are performed on top of this dynamics, we call this approach a banged-DAQC (bDAQC) circuit. Interestingly, although this introduces a systematic error, this scales better than main error sources found in quantum computers [15]. It has already been experimentally proven that DAQC is a suitable paradigm for the NISQ era, for instance, in the implementation of a variational quantum algorithm in a system with up to 61 qubits [16]. 2643-1564/2024/6(1)/013280(14) 013280-1 Published by the American Physical Society MIKEL GARCIA-DE-ANDOIN et al. PHYSICAL REVIEW RESEARCH 6, 013280 (2024) In order to enhance the range of quantum platforms suitable for DAQC, we must extend the techniques to arbitrary resource Hamiltonians. Previously, the resource Hamiltonian was the aforementioned Ising Hamiltonian. Additionally, the construction of an arbitrary target Hamiltonian was performed through a two-step procedure: first by transforming the source Hamiltonian into an adequate ZZ Hamiltonian [17], then by employing sequences of SQGs to build an arbitrary Hamiltonian [14]. A different step-by-step construction for arbitrary qudit Hamiltonian simulation was proposed, in this case for the analogue quantum computing paradigm [18]. This technique also requires a step-by-step construction, which involves additional ancillary qubits for the encoding. However, the question of systematically performing the transformation in a single step for DAQC and two-body qubit Hamiltonians was still open. In this article, we provide an explicit construction for a DAQC protocol to approximate the evolution under an arbitrary target two-body Hamiltonian by evolving under another arbitrary two-body resource Hamiltonian up to a certain Trotter error. By means of the Trotter-Suzuki formula, we argue that by repeating this sequence nTtimes one can reduce the error exponentially with nT. The tools we develop in this article allow for a practical realization of DAQC schedules in faulty hardware, in which spurious couplings prevent us from approximating the system as an Ising Hamiltonian. We also solve the problem of the negative analog blocks times that limited the implementability of previous protocols. Additionally, we introduce a classical Bayesian optimization technique to find optimal angles for the SQGs. Taking into account the depth limitations of quantum circuits in the NISQ era, our objective is to maximize the fidelity of the circuit employing a fixed amount of digital-analog blocks. We show how it is possible to employ DAQC schedules with a low number of digital-analog blocks that achieves fidelities compared to a systematic approach with a higher count of blocks. The rest of the article is organized in the following manner. In Sec. II, we review the previous protocol employing ZZ Hamiltonians. In Sec. III, we present the new protocol, which extends it to an arbitrary two-body source Hamiltonian, and discuss the error scaling. Then, in Sec. IV, we introduce an optimization technique for approximating arbitrary dynamics employing a fixed number of blocks, and illustrate the technique for a particular problem. Finally, in Sec. Vwe conclude with some final remarks. II. WARMUP: DAQC PROTOCOL FOR ZZ HAMILTONIANS As a warmup, let us review the previous protocol for simulating the dynamics during a time Tof a target ZZ all-to-all (ATA) Hamiltonian, HT,ZZ = n  i<j gi,jσz iσz j,(1) by employing a source ZZ ATA Hamiltonian, HS,ZZ = n  i<j hi,jσz iσz j,(2) for {gi,j,hi,j}∈R. ... X t1,2   X  X t1,3 X t4,5 X X X X X X X X FIG. 1. sDAQC circuit for the Ising ZZ ATA Hamiltonian for five qubits. For simulating an arbitrary evolution, we sandwich several analog blocks with a couple of Xgates applied to all combinations of two qubits. The blocks labeled with ti,jrepresent the unitary evolution under the source Hamiltonian for the time ti,j;thisise−iti,jHS. For achieving this, in Refs. [14,19] the authors proposed a universal protocol, pictorially shown in Fig. 1. It consists in sandwiching each analog block with two Xgates, applied to a different pair of qubits each time. Effectively, this changes the sign of all couplings in which only one of the qubits is selected. Noticing that all terms of the Hamiltonian commute with each other, we have that the Trotter formula is exact, so we can write UT=e−iT HT = i<j exp −iti,j <m (−1)δi+δim+δj+δjm h,mσz σz m,(3) where δij is the Kronecker delta and ti,jis the analog time for the corresponding analog block. Note that, throughout this work, we are considering ¯h=1. We can rewrite the equation more conveniently as a linear system of equations M t=T−→ g/h,(4) where the elements of the matrix M(i,j),(,m)= (−1)δi+δim+δj+δjm represent the effective signs of the couplings between qubits (i,j) in every analog block (, m),  tis the vector of times of each analog block, and (g/h)(i,j)is a vector of the proportion between the target and the source coupling strengths between the qubits (i,j). If there is a missing coupling in both source and target Hamiltonians, gi,j=hi,j=0, then we remove the corresponding element from the vector −→ g/hand the corresponding row in M.Ifthe coupling is missing but the target coupling is nonzero, then the Hamiltonian cannot be simulated directly, as one would need implement a SWAP strategy to simulate the Hamiltonian [19]. It can be proven that the matrix Mis nonsingular for all numbers of qubits except 4, so we can obtain an exact simulation of the desired dynamics employing this schedule. III. EXTENSION OF DAQC TO ARBITRARY TWO-BODY HAMILTONIANS The proof that almost any entangling two-body Hamiltonian, together with SQGs, can be employed to simulate the dynamics of another two-body Hamiltonian was shown in Ref. [20], but no constructive method was provided. From 013280-2 DIGITAL-ANALOG QUANTUM COMPUTATION WITH … PHYSICAL REVIEW RESEARCH 6, 013280 (2024) FIG. 2. sDAQC circuit for an arbitrary ATA Hamiltonian for five qubits. The protocol for simulating any target Hamiltonian can be divided into two steps. In the first step (top), we select each of the n(n−1)/2 pairs of qubits in our system. In the second step (bottom), for each pair of qubits, we apply all the nine possible combinations of the Pauli gates. The number of analog blocks for this protocol scales quadratically with the number of qubits, O(n2). now on, we will refer to the source Hamiltonian as HS= n  i<j μ,ν∈{x,y,z} hμ,ν i,jσμ iσν j,(5) and the target Hamiltonian as HT= n  i<j μ,ν∈{x,y,z} gμ,σ i,jσμ iσν j,(6) where σμ iis the Pauli operator μacting on qubit iand {gμ,ν i,j,hμ,ν i,j}∈R. Then, the objective is to obtain a circuit to simulate the evolution of HTfor a time T,U=e−iT HT. The first step of the proof is to note that it is possible to decouple a pair of qubits from the rest in an n-qubit system. Then, by employing a 36-step digital-analog protocol, it can be proven that any two-qubit interaction can be simulated. The proof for universality can be obtained by extending this to every coupling in the target Hamiltonian, and repeating the circuit nTtimes for simulating a time =T/nTin each Trotter step. However, this protocol gets convoluted as the number of qubits in the system, n, increases. In general, the circuit requires O(n3nT) analog blocks, with an error of ε∼ O(n2T2/nT). As it stands, the question of obtaining a more efficient protocol is still open. For the general case, we can extend the ideas reviewed in Sec. II to arbitrary two-body Hamiltonians. The first step for our protocol is to select each pair of qubits {i,j}of our system. Now, instead of just applying an Xgate to both of them, we will apply all nine possible choices of pairs of gates, {XX,XY,XZ,YX,YY,YZ,ZX,ZY,ZZ}. This is illustrated for a simple example in Fig. 2. Applying this to every pair of qubits, we effectively change the sign of some of the couplings, generating a nonsingular system of 9n(n−1)/2 equations. As the number of equations coincides with the number of variables and the number of parameters to define an arbitrary two-body Hamiltonian, this protocol is optimal in the number of digital-analog blocks. Unfortunately, in the general case the Hamiltonian does not commute with itself. This means that, unlike the warmup case, the Trotter formula is not exact, and thus, if we want to achieve an arbitrarily small error, we need to employ more Trotter steps. As a quick sketch of the proof, we write the problem as a system of equations similar to the one in Eq. (4), M(n) t=T−→ g/h.(7) The matrix for an n-qubit system can be constructed recursively as a block matrix, M(n)=A(n)P(n) Q(n)M(n−1),(8) where the blocks A(n), Q(n), and P(n) can be constructed systematically by taking into account the change of signs of the effective Hamiltonian terms after sandwiching them by Pauli gates. Then, by using the properties of this matrix and the definition of the formal determinant, we can prove that it is nonsingular. Further details for the proof of the universality of the protocol are given in Appendix A. With this result, we can then employ the same results as in the original work by Suzuki [21] to argue that an arbitrarily small error can be attained. A. Analysis of the errors In order to obtain a bound for the maximum error of this protocol, we can resort to the original error analysis of the Suzuki-Trotter formula [21]. Since we have proven that the sum of the effective Hamiltonians in each block is exactly the target Hamiltonian, we can employ the formula for the (nT,1) approximant, ε=UT−US=     e−iT HT− k e−itk nTH(k) SnT     ⩽2 nT ktkH(k) S2 e nT+2 nTktkH(k) S,(9) where H(k) Sis the effective Hamiltonian in the kth analog block and ·is the Frobenius norm defined as A=√AA†.We will employ this norm for matrices throughout the text. Since in each block only the sign of some Pauli string terms in the Hamiltonian changes, the norm is the same for all blocks H(k) S=HS. With this, we have that the error is bounded by the sum of the times of the analog blocks. If we assume that we have a correct protocol in which all the times of the analog blocks are positive, we can rewrite ε⩽2 nT t2 AHS2e nT+2 nTtAHS,(10) where tA=ktk. 013280-3 MIKEL GARCIA-DE-ANDOIN et al. PHYSICAL REVIEW RESEARCH 6, 013280 (2024) The total analog time will be lower bounded by the norms of HSand HTand the time Tof the simulation, tA⩾ THT/HS. This corresponds to a situation in which both Hamiltonians are proportional to each other, HS∼HT.Let us now study a general case. Assume an optimal protocol in which the total time is minimized over all possible protocols. In this case, tAis upper bounded by the weakest source coupling and the strongest target coupling as tA⩽C(n)Tmax(i,j,μ,ν)gμ,ν i,j min(i,j,μ,ν)hμ,ν i,j ,(11) where C(n)⩾1 is a constant depending both on the protocol and the system size. In the case a coupling hμ,ν i,jis zero we take it out of the corresponding element from both  gand  h and the corresponding row from the matrix M. In general, the constant C(n) heavily depends on the protocol, but it does not monotonically grow with the number of qubits. For example, we find that the original protocol in Ref. [14] is upper bounded by C(n)=n(n−1)/(n−4)(n−5) ⩽15 when n>5fora balanced solution ( t∼(1,...,1)). Similarly, the protocol in Ref. [19] for a nearest-neighbor Hamiltonian has C(n)⩽3/2 for any system size. Employing this result, we see that we can make the error arbitrarily small by increasing the number of Trotter steps of the protocol. However, this would only work for sDAQC circuits. When working on the bDAQC paradigm, increasing the number of Trotter steps would increase the bang error in the protocol, as the time to apply the SQGs remains the same. This means that there is a point at which increasing the number of Trotter steps actually reduces the fidelity of the circuit. However, this analysis should be performed case by case, as it heavily depends on the problem and the system characteristics. As a rule of thumb, the time for the shortest analog block in the circuit should be at least two orders of magnitude above the time for applying a SQG, min(tk)102tSQG. B. A note about negative times When computing the times of the analog blocks in Eq. (4), one can obtain a solution comprising some negative times. Simulating the evolution over a negative time would require to completely flip the sign of all the terms in the corresponding H(k) S. However, this cannot be done in general. For instance, one can straightforwardly prove that one cannot do this for the three-qubit all-to-all connected system by exhausting all possible combinations. As a suggestion for solving this problem, it was originally proposed to add one extra analog block, without it been sandwiched by any SQG. However, this only solves the problem in some particular cases. Approximate solutions to similar problems have been proposed [22], but the question of finding a systematic solution was still open. Here we propose a method for constructing DAQC protocols in which the times are all positive. This method is highly inefficient in single-qubit gate counts, but useful for proving the existence of such a solution. The problem is exactly the same as in Eq. (7) but, instead of a square matrix M(n), we will employ all combinations of Pauli gates plus the identity {1,X,Y,Z}to construct a matrix with 4ndifferent columns Mi. Then, we map the problem to a non-negative least-squares (NNLS) problem, for which we then employ the Algorithm NNLS to obtain a positive solution [23]. However, we first need to prove that a positive solution exists. Here we roughly sketch the proof for this claim. First, we note that the columns Micorrespond to the vertices of a polytope. Second, we prove that there is a strictly positive solution for the homogeneous system −→ g/h= 0, with  t=4−n 1. Last, we build a hypersphere centered in  0 with small radius with a strictly positive vector   s=4−n 1+ ε, from which we can reach to any possible problem −→ g/hby scaling with a positive number, −→ g/h=Mλ(4−n 1+ ε)forλ>0. Even though we are using an exponential number of blocks for the proof, numerically we observe that the solution contains only ∼9n(n−1)/2 nonzero elements. The full proof is provided in Appendix B. IV. CLASSICAL OPTIMIZATION OF THE DAQC SCHEDULE The protocol discussed in the previous section is a systematic method to obtain an arbitrary target Hamiltonian using another arbitrary Hamiltonian as a source. Its implementation involves a number of digital-analog blocks that grow quadratically with the number of qubits. With the current limitations of NISQ devices, long circuits can accumulate large experimental errors, so it becomes necessary to find a trade-off between the accuracy of the theoretical approximation and the required experimental resources. With the goal of reducing the number of digital-analog blocks, we propose a classical optimization strategy to find a set of SQGs sandwiching Kanalog blocks such that the digital-analog schedule is as close as possible to the ideal evolution. This way, we propose an optimization problem where the parameters to be optimized are the times of the analog blocks, tk, and the parameters of an arbitrary SQG, R(θ,φ,λ)=⎛ ⎝ cos θ 2−eiλsin θ 2 eiφsin θ 2ei(λ+φ)cos θ 2⎞ ⎠,(12) where {θ,φ,λ}∈[0,2π) are the rotation angles of the SQG. The cost function we minimize is the Frobenius distance between the target evolution UTand the circuit with the optimized parameters UC, similar to the calculation in Eq. (9). By employing the Frobenius distance between the unitaries as a proxy for the fidelity, we can test the approach in general, without imposing any assumptions about the initial state of the system or without expensive Haar integral calculations. For simulating the circuits, we have employed two techniques. In one, we simulated the exact evolution of an ideal quantum computer. Since the cost of computing the exact evolution under a Hamiltonian scales exponentially with the number of qubits, n,O(23n), we have tested a less resourcedemanding method as well. By means of the first-order Trotter expansion, we can approximate the evolution as UT=e−it n i<jμ,ν∈{x,y,z}hμ,ν i,jσμ iσν j ≈Uappx = n  i<j e−it μ,ν∈{x,y,z}hμ,ν i,jσμ iσν j⊗1r(i,j),(13) 013280-4 DIGITAL-ANALOG QUANTUM COMPUTATION WITH … PHYSICAL REVIEW RESEARCH 6, 013280 (2024) where 1r(i,j)is the identity matrix for the subspace of all qubits except iand j. Using this approximation, we can reduce the cost of calculating the matrix exponential. Additionally, since every term is a sparse matrix with sparsity ∼1−2−n, we can employ efficient functions for the matrix products. Techniques involving matrix product states (MPSs) or matrix product operators (MPOs) could be useful for extending this classical optimization strategy to larger systems [24,25]. A. The parameter space The complete parameter space is given by the parameters of the SQG applied to qubit iof the kth analog block, {θi k,φi k, λi k}∈[0,2π), and the evolution time of the kth block, tk.We have restricted the evolution time of the analog blocks to be 0⩽tk⩽THP/KHS, where Kis the number of analog blocks, to avoid both negative and times much larger than the total time Tof the target evolution. As the number of qubits increases, this quickly leads to a wide and hard-toexplore parameter space, which makes convergence slow in the optimization process. We found a good compromise in reducing the kind of SQG applied on each block to just two, one Rk(θk,φ k,λk) applied to all even qubits and an additional R k(θ k,φ k,λ k) applied to all odd qubits. These kind of oddeven SQG layers are the same as those required to compute the Trotterized evolution in many DAQC systems [26], which gives a good starting point and sufficient flexibility for the optimization. Calculating the evolution under the source Hamiltonian for all analog blocks is a resource-intensive task. As a mean of simplifying the computational cost and reducing the number of variables for the optimization, we have performed tests in which we fix the analog-block times to a fraction of the total evolution time of the target evolution, T, divided by the total number of analog blocks, K. Even though in certain cases one can achieve better results by having the evolution time tkas an optimization parameter, fixing tk=T/Kresults not only in a faster computation, but in a much faster and consistent convergence as well, as we show in Sec. IV C. B. Optimization protocol Despite the compromises previously described, the optimization landscape is complex and shows multiple local minima. Moreover, evaluating the cost function for a new set of parameters is computationally very expensive. There are many strategies to tackle an optimization of this kind of “black-box” function, such as genetic [27]orswarm[28] algorithms. In this work, we propose a combined strategy of Bayesian optimization and a gradient descent algorithm. The popular gradient descent consists of evaluating a point, computing the gradient of the function at that point, and tuning the parameters following the slope of the function. It is a fast and efficient tool to exploit local minima, but falls short when searching for a global minimum and is highly dependent on the initial point of the optimization. Meanwhile, Bayesian optimization is especially advantageous to efficiently explore unknown and computationally expensive functions. Bayesian optimization treats the black-box function as a random function and considers a prior upon it. Then, it evaluates only the function, and not its derivative, to compute an acquisition function, which usually is based in the expected improvement or the probability of improvement. Finally, it uses this acquisition function to find the next point to evaluate and updates the prior for the next iteration. For an in-depth read on Bayesian optimization, we refer the readers to Refs. [29,30]. As mentioned, the gradient descent is highly dependent on its initial point, as it usually converges towards the nearest local minimum, not the global one. A common way of dealing with this problem is performing several gradient descent optimizations by changing the initial point, such as with random or grid searches on the parameter space [31]. The aforementioned Bayesian method allows us to minimize the number of initial points by means of a guided search through its acquisition function. The combined strategy leads to a much faster convergence to the global minimum than random or grid searches, achieving faster and better results. C. Example: Nearest-neighbor Hamiltonian A simple case of simulating the dynamics of the onedimensional XY model is presented as a test for this strategy. Let us assume a target nearest-neighbor Hamiltonian of the form HT=g n−1  i=1 σx iσx i+1+σy iσy i+1.(14) Let us also assume that the source Hamiltonian to which we have access is HS= n−1  i=1 hxz iσx iσz i+1+hzx iσz iσx i+1+hzz iσz iσz i+1.(15) This problem was also studied within the DAQC paradigm in Ref. [17]. In the following, we will discuss two cases: the homogeneous case, in which all coupling strengths are equal to the target Hamiltonian coupling strength, hk i=h∀i,k, and the inhomogeneous case. For the inhomogeneous case, we propose a more realistic scenario in which the coupling strengths of the system follow a Gaussian distribution around the value of the target coupling strength ˜ hk i=N(g,σ), where σis the variance of the distribution. To simulate Eq. (14)usingEq.(15) in the homogeneous case, the general approach can be simplified into applying a single-qubit πrotation around the xaxis to all the qubits, R(n) x(π)=⊗ n i=1σx, so that the crossed terms (σxσzand σzσx) change signs. Then, since the rotated and system Hamiltonian do not commute, we employ the Trotter approximation R(n) x(π)†e−iHST 2nTR(n) x(π)e−iHST 2nTnT≈e−iHZZT,(16) where HZZ =gn−1 i=1σz iσz i+1and nTis the number of Trotter steps. The second step consists in sandwiching the resulting HZZ evolution in either a Rx(π/2) or Ry(π/2) to rotate all qubits from σz ito σy ior σx i, respectively. This, of course, requires an additional Trotterization. This approximation can be made arbitrarily precise by increasing the number of Trotter 013280-5 MIKEL GARCIA-DE-ANDOIN et al. PHYSICAL REVIEW RESEARCH 6, 013280 (2024) steps at the cost of increasing the number of required digitalanalog blocks (four blocks per Trotter step). As previously described, we can use two or more analog blocks sandwiched between arbitrary rotations on even and odd qubits, and find the optimal rotation angles to achieve the best approximation. Since the set of all arbitrary rotations includes the specific rotations used for the Trotter approximation, the optimization must result in values necessarily equal to or better than the Trotter approximation for the same amount of blocks. In Fig. 3, we show a comparison between the performance of a Trotter approximation and the proposed Bayesian plus gradient descent optimization with N=6 qubits. The chosen coupling strengths for the studied case were hk i=g=1∀i,k in the homogeneous case and ˜ hk i=N(1,0.175) for the inhomogeneous one. Additionally, we performed the optimization process both with fixed evolution times and with tkas an optimization parameter. For each point, we computed 20 runs, using only a ten-step optimization for the Gaussian process. Fixed-time optimizations achieve convergence in these ten steps, as evidenced by their small error bars. Larger error bars show a higher variance in the achieved value for the cases with time as a parameter, meaning a lack of convergence and the need for a longer optimization with more steps to reach the global minimum. Comparing the results obtained with the Trotter formula we see an improvement on the Frobenius distance. In particular, we see an improvement of ∼68% and ∼78% for the four-block cases in which the time is not a parameter for the inhomogeneous case and homogeneous case, respectively. When we include the analog block times as a parameter, we only obtain mean improvements of ∼45% and ∼65%, respectively. Although leaving the time as a parameter can lead to better results, the improved consistency and convergence rates of fixed-time optimization, along with a much lower computational cost, makes the fixed-time optimization the most desirable approach. Homogeneous and random inhomogeneous couplings perform comparably. As expected, the case with homogeneous couplings has a slight advantage, as the target Hamiltonian is also homogeneous. This is good evidence that, as long as the coupling strengths associated with the system and the target Hamiltonians are similar, the approach is suitable for arbitrary two-body Hamiltonians. Most importantly, all optimization approaches achieve a significantly smaller Frobenius distance and, therefore, a better fidelity than the traditional first-order Trotter approximation for the same number of analog blocks. This way, by employing this classical optimization strategy on the digital layers (SQGs) of the simulation, one can achieve the same or better theoretical precision in the approximation, saving up experimental resources. With current NISQ devices, this could directly lead to reduced experimental errors. Additional nonextensive testing has been done in combining this optimization strategy with classical approximate simulations of the system. As the cost of classically calculating the unitary evolution under a Hamiltonian scales exponentially with the number of qubits, we need to find an alternative with a better scaling and which still can provide some advantage. With this goal, we employed the firstorder Trotter formula for simulating the analog blocks, which FIG. 3. Frobenius distance between the optimized circuit and the exact evolution vs the number of digital-analog blocks. (a) We solve the problem of the nearest-neighbor Hamiltonian for six qubits. We show two cases, one in which the couplings in the system Hamiltonian is homogeneous, and the other where the coupling strengths are in a normal distribution around the values of the former. We also distinguish two versions for the optimization process, which fixes the length of the analog blocks or leaves them as a variables. As a baseline, we represent the distances obtained with first-order Trotterization up to nT=19 Trotter steps. The shown values indicate the mean Frobenius distance and the error bars show the interquartile range. The inset on top provides an extended view in the coordinate axis for a higher number of blocks in logarithmic scale. (b) We repeat the procedure for a single run with a different classical simulation method, in which we employ a first-order Trotter approximation for computing both the circuit and the target. reduces the cost of the calculation when combined with sparse matrix multiplication algorithms. Here, instead of exactly calculating the evolution of our circuit, we employed the first-order Trotter formula for approximating the evolution under both the target and the system Hamiltonians. These 013280-6 DIGITAL-ANALOG QUANTUM COMPUTATION WITH … PHYSICAL REVIEW RESEARCH 6, 013280 (2024) approximations, UT,appx and US,appx, respectively, are employed during the Bayesian optimization process. For checking the performance one would get when using the result of the optimization for running the quantum circuit, we calculate the Frobenius distance with the exact expressions. This way, we can directly compare the results obtained in both tests. We conclude that the Trotter approximation limits our capacity for optimizing the circuit. For the lowest amount of blocks for the Trotterization (four blocks), we can obtain ∼54% (∼51%) of improvement for the inhomogeneous (homogeneous) problems when we discard the time as an optimizable parameter, and ∼56% (∼55%) improvement when we include it. This confirms the higher expressibility of the optimization model when we include the time as a parameter, at the cost of increasing the time to solution of the classical algorithm. However, for shallow circuits, the improvement of the fidelity is worth the optimization process. n particular, the fixed-time optimization gives us a competitive Frobenius distance for circuits with five or less blocks compared to the Trotter decomposition. As it is also shown, we can simulate a system with a single analog block while maintaining the fidelity of the evolution obtained with four blocks employing the Trotter formula, Uexact −UTrotter=11.1, Uexact −Uappx=10.3for the inhomogeneous problem, and Uexact −Uappx=10.5for the homogeneous. V. CONCLUSIONS In this article, we provide a set of versatile tools to implement DAQC circuits in realistic scenarios, where the terms in the Hamiltonians of both the system and the target do not commute by pairs. Following on the universality of DAQC, we have proposed a protocol for simulating an arbitrary twobody Hamiltonian by employing another arbitrary two-body Hamiltonian. This was achieved by sandwiching the analog blocks with pairs of SQGs picked from the Pauli gate set. This systematic method is universal, and provides an optimal number of digital-analog blocks Qthat scales quadratically with the number of qubits for all-to-all Hamiltonians, Q⩽ 9n2. This strategy employs less blocks compared to previous approaches, which required Q⩽36n3steps. Additionally, we have proposed a new method for obtaining DAQC protocols in which all analog block times are positive and, thus, are implementable in practice. With the same objective, we have explored a less resourceintensive approximation. We have proposed a classical strategy in which all possible SQG rotations are considered as digital blocks. Solving this optimization problem, we have shown that one can achieve an accuracy comparable to or even better than large sequences of Trotter steps employing significantly less analog blocks. The benefits are twofold, as not only is the accuracy of the systematic error decreased, but a reduced amount of analog blocks also reduces the experimental error. In this regard, we have proposed a classical optimization strategy consisting in a Bayesian optimization method that explores the parameter space to find the optimal initial point of a gradient descent. This process is done without any assumption about the initial state, as the information given to the routine is just the target unitary evolution. For a low number of Trotter steps, we have shown that a precision comparable to that of regular Trotter approximations can be achieved by using only a fraction of the experimental resources. The two approaches we have proposed in this work provide new tools for the implementation of DAQC circuits in NISQ devices. These new protocols for working with arbitrary two-body Hamiltonians pave the way for generating and scaling experimental implementations of this paradigm. The next steps for widening DAQC should search for an optimal solution for the negative analog block times problem, and further extend the protocols for k-body Hamiltonians. The exact method involves O(3kn2) digital-analog blocks in the worst case, but ∼9n2in numerical simulations. This suggests that approximate solutions, similar to the one proposed in this article, should be explored for obtaining implementable DAQC circuits. The code employed for the numerical experimentation is available under reasonable request. ACKNOWLEDGMENTS The authors thank M. Reichert, P. Rodríguez, and M. Capela for useful discussions and valuable feedback. A.S.C., P.P.-F., and L.L. acknowledge financial support from the Spanish Grants No. PID2019-104002GB-C21, No. PID2019-104002GB-C22, No. PID2022-136228NBC21, and No. PID2022-136228NB-C22 funded by Ministerio de Ciencia e Innovación/Agencia Estatal de Investigación MCIN/AEI/10.13039/501100011033, FEDER “A Way of Making Europe,” Consejería de Conocimiento, Investigación y Universidad, Junta de Andalucía, European Regional Development Fund (ERDF) under Project No. US-1380840, Grant Groups FQM-160 and FQM-177, and the project PAIDI 2020 with Reference No. P20_01247 and No. P20_00617 funded by the Consejería de Economía, Conocimiento, Empresas y Universidad, Junta de Andalucía (Spain). The authors acknowledge resources supporting this work provided by the CEAFMC and Universidad de Huelva High Performance Computer (HPC@UHU) funded by ERDF/MINECO Project No. UNHU-15CE-2848. M.G.d.A. and M.S. acknowledge support from the EU FET Open project EPIQUS (899368) and HORIZON-CL4-2022-QUANTUM01-SGA Project No. 101113946 OpenSuperQPlus100 of the EU Flagship on Quantum Technologies, the Spanish Ramón y Cajal Grant No. RYC-2020-030503-I, Grant No. PID2021-125823NAI00 funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe” and “ERDF Invest in your Future,” Basque Government, through Grant No. IT1470-22 and from the IKUR Strategy under the collaboration agreement between Ikerbasque Foundation and BCAM on behalf of the Department of Education of the Basque Government. M.G.d.A. acknowledges support from the UPV/EHU and TECNALIA 2021 PIF contract call. This work has been financially supported by the Ministry for Digital Transformation and of Civil Service of the Spanish Government through the QUANTUM ENIA project call–Quantum Spain project, and by the European Union through the Recovery, Transformation and Resilience Plan–NextGenerationEU within the framework of the “Digital Spain 2026 Agenda.” 013280-7 MIKEL GARCIA-DE-ANDOIN et al. PHYSICAL REVIEW RESEARCH 6, 013280 (2024) APPENDIX A: UNIVERSALITY OF THE PROTOCOL We want to solve the problem of simulating an arbitrary n-qubit two-body ATA Hamiltonian HT= n  i<j μ,ν∈{x,y,z} gμ,ν i,jσμ iσν j,(A1) with a DAQC schedule employing another arbitrary two-body Hamiltonian, HS= n  i<j μ,ν∈{x,y,z} hμ,ν i,jσμ iσν j,(A2) where {gμ,ν i,j,hμ,ν i,j}∈R. We are assuming that if we have a nonzero gμ,ν i,jterm, we will have a nonzero hμ,ν i,jcoupling in the source Hamiltonian. The protocol we propose is the following: select each of the possible pairs of qubits i,j, and apply all the combinations of {x,y,z}gates, which are the corresponding Pauli gates. This changes the signs of the effective couplings, according to a (±1) matrix, which we will call M(n). Now, we want that our DAQC schedule simulates our original problem; this is e−iT HT≈ n  i<j μ,ν∈{x,y,z} exp −itμ,ν i,j n  i<j μ,ν∈{x,y,z} ×M(n)μ,ν,μ,ν i,j,i,jgμ,ν i,jσμ jσν j.(A3) To prove that this protocol is universal we have to prove that we can always find a set of times for the duration of the analog blocks, (tμ,ν i,j). In our proof, we focus on the definition of universality provided in Ref. [20] for its simplicity. However, universality can also be proven using a complementary definition, for example, the one proposed in Refs. [18,32]. In this case, conditions i, ii, iii, and iv from Ref. [18] are fulfilled by construction, and condition iv can be obtained by controlling the Trotter error following the discussion in Sec. III A. 1. Notation and definitions We define the notation employed for the proof. We employ a single index for labeling pairs of qubits, b(i,j,n)=n(i−1) −i(i+1)/2+j,1⩽i<j⩽n. (A4) Whenever we refer to a label of a coupling or a pair of qubits, we will be using this labeling convention. We also give a formula for the inverse of this indexing [33]: i(b,n)=n−n(n−1) −2b+2+1 2, j(b,n)=b−n(i(b)−1) +1 2i(b)(i(b)+1).(A5) Additionally, we define an indexing method for the different selection of SQGs for a pair of qubits, f(μ, ν). These indices are given from 1 to 9 according to the order in the following list: {xx,xy,xz,yx,yy,yz,zx,zy,zz}. For example, the pair of SQGs {y,z}has the index f(y,z)=6. The elements M(n)μ,ν,μ,ν i,j,i,jcan be expressed as a matrix by employing the previous labeling of the pair of qubits, which can be then mapped to a matrix M(n)μ,ν,μ,ν i,j,i,j= M(n)g(i,j,μ,ν),g(i,j,μ,ν). The rows and columns of this matrix are defined as g(i,j,μ,ν)=9(b(i,j)−1) +f(μ, ν).(A6) The times for each analog block, tμ,ν i,j, can be expressed as a column vector by employing the same labeling as in Eq. (A6), tμ,ν i,j=tg(i,j,μ,ν). We also define a new column vector −→ g/h, which contains the information about the couplings in the hardware and the couplings we want to simulate; again, we employ the same labeling as in Eq. (A6) such that (−→ g/h)μ,ν i,j= (−→ g/h)g(i,j,μ,ν). Now, the problem of finding the times for the analog blocks can be expressed as M(n)−→ t=T−→ h/g.(A7) Let us write the full M(n) matrix as a block matrix, M(n)=A(n)P(n) Q(n)M(n−1),(A8) where A(n)isan(n−1)×(n−1) block matrix, and M(n) an n(n−1)/2×n(n−1)/2 block matrix. Here, each row labels a coupling in HSand each column labels the pair of qubits where the SQGs are applied. For example, the coupling between qubits 1 and 4 would be addressed in row I= b(1,4) =4, and the blocks in which the SQGs are applied on qubits 2 and 3 correspond to the column J=b(2,3) =n. Each of these terms is defined also as block matrices (which we call from now on submatrices). These blocks are 9×9 matrices, that can be one of the four {M2,M1.1,M1.2,M0}, according to the following formula: MI,J=⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ M2if I=J M1.1if i(I)=i(J)ori(I)=j(J) M1.2if j(I)=i(J)orj(I)=j(J) M0else. (A9) For a detailed construction of an M(n) matrix see Appendix A4. Let us define each of these submatrices: (i) M2: This is the matrix that represents the signs of the effective couplings μ, ν when we apply a SQG to each of the two qubits i,j. The columns represent each pair of gates, in the order given by f(μ, ν). Equivalently, each column represents the sign of each effective coupling between the pair of qubits i,j, with the same order as the columns. Then, 013280-8 DIGITAL-ANALOG QUANTUM COMPUTATION WITH … PHYSICAL REVIEW RESEARCH 6, 013280 (2024) we have the following matrix: M2= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1−1−1−11 1−11 1 −11−11−11 1−11 −1−11 11−111−1 −11 1 1−1−1−11 1 1−11−11−1 1 −11 11−1−1−11 11−1 −11 1−11 1 1−1−1 1−11 1−11−11−1 11−111−1−1−11 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .(A10) (ii) M1.1: In this case we only apply a SQG in the first qubit of the pair, M1.1= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1−1−11−1−1 1 −1−1 −11−1−11−1−11−1 −1−11−1−11−1−11 1−1−1 1 −1−1 1 −1−1 −11−1−11−1−11−1 −1−11−1−11−1−11 1−1−1 1 −1−1 1 −1−1 −11−1−11−1−11−1 −1−11−1−11−1−11 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .(A11) (iii) M1.2: In this case we apply one SQG on the second qubit, M1.2= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 111 −1−1−1−1−1−1 111 −1−1−1−1−1−1 111 −1−1−1−1−1−1 −1−1−1111−1−1−1 −1−1−1111−1−1−1 −1−1−1111−1−1−1 −1−1−1−1−1−1111 −1−1−1−1−1−1111 −1−1−1−1−1−1111 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .(A12) (iv) M0: In this case neither SQG is applied on any of the two qubits. This is the trivial case in which the sign of the effective couplings does not change: M0= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 111111111 111111111 111111111 111111111 111111111 111111111 111111111 111111111 111111111 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .(A13) For example, we apply an xgate to qubit 1 and a ygate to qubit 3. We want to know how the sign of the coupling yz changes between qubits 2 and 3. In this case, we are applying a SQG only to the second qubit, so we have to look at the M1.2submatrix. Now, the SQGs we have applied are xy,so this corresponds to the second column, and the yz coupling to the sixth row. Looking at the matrix, we see that in this case we have a change of sign of the effective coupling, and thus we have a −1 in the corresponding matrix element. We have some properties with these subblocks: (1) All the {M2,M1.1,M1.2,M0}commute in pairs. (2) M1.1M1.2=M2M0=M0. (3) M1.1M0=M1.2M0=(−3)M0. (4) (M0)2=9M0. (5) (M1.1)2=−3M2M1.1. (6) (M1.2)2=−3M2M1.2. Also, we have some properties related to how these subblocks appear in the matrix: (1) M(n)I,I=M2. (2) M(n)I,J={M1.1,M1.2,M0},∀I= J. (3) M(n)I,J=M0⇔M(n)J,I=M0. (4) M(n)1,J={M2,M1.1,M0}. (5) M(n)n(n−1)/2,J={M2,M1.2,M0}. (6) The number of subblocks that fulfill M(n)I,J= M(n)J,I=M1.1is the same as the number of subblocks fulfilling M(n)I,J=M(n)J,I=M1.2. 013280-9