1 Introduction↩︎

Interacting atomic, molecular and solid-state spin systems are the fundamental constituents of various quantum processors. This includes in particular trapped ions [1], [2], Rydberg atoms [3], [4], and superconducting circuits [5], [6]. In these systems, one typically combines the ‘native’ interaction (resource) Hamiltonian \(H_R\) with external control fields, in order to engineer a certain quantum dynamics described by a unitary operator \(U\). This is of key interest for building a quantum circuit in a quantum computer, where \(U\) represents a sequence of quantum gates, but also for quantum simulation, where the aim is to engineer the dynamics \(U=\exp(-i\tilde{H}T)\) of a many-body Hamiltonian \(\tilde{H}\). The goal of this article is to present a complete framework for quantum circuit and Hamiltonian engineering of generic quantum spin systems. Importantly, we show that this framework is at the same time flexible, i.e. it can implement arbitrary pairwise interactions, and robust with respect to various relevant experimental imperfections.

In recent years, we have seen multiple proposals and experimental realizations of protocols where external driving fields engineer quantum dynamics from a spin system. For instance, one can physically move the spins to make them interact or decouple in a programmable fashion. While this was demonstrated in remarkable pioneering experiments in Rydberg atom arrays [7], coherently moving the spins on short timescales, as well as exporting this method to other quantum platforms, are still significant technical challenges. We will focus here on another approach that consists in applying external field ‘pulses’ to the spins, in order to engineer an effective dynamics.

There are two main families of engineering protocols using pulses. First, in digital-analog quantum simulation/computing (DAQC) [8]–[13], one aims to realize a large class of unitary operators \(U\) given a single resource Hamiltonian \(H_R\), see in particular Ref. [11], [13] for numerical pulse sequence design to engineer universal quantum computation in spin systems. Note that, instead of addressing the spins individually, one can also use multiple spins to implement each qubit and drive them with a global field [14]. While the DAQC approach is theoretically supported by theorems stating that we can always design a pulse sequence realizing a unitary operation \(U\) to arbitrary precision using a fixed entangling Hamiltonian [15]–[17], strong experimental requirements are imposed. In particular, one requires either to use pulses that are perfectly calibrated and quasi-instantaneous, as, in general, the interplay between errors arising from discretization of the dynamics and pulse calibration errors hinders a practical implementation [18], requiring either error mitigation techniques [19] or additional pulses [20]. The second family of protocols is related to the idea of Floquet engineering [21]–[26]. Instead of aiming at universal quantum computation, here the effective dynamics is realized with a sequence of few pulses, designed mainly analytically and in a model-dependent fashion. Importantly, when the pulse sequences are global, i.e. when all qubits are subject to the same pulses, one can make the resulting dynamics robust to experimental imperfections, such as external noise fields and errors due to imperfect pulses [27]–[30]. The simplicity and the robustness features have made this second approach appealing for quantum simulation experiments [31]–[37].

In this work we demonstrate a protocol that merges the advantages of both DAQC and Floquet engineering, i.e. we provide a single experimental recipe to implement universal quantum processors based on arbitrary two-body interactions, using analytical pulse sequences that are inherently robust to experimental imperfections. This makes our proposal accessible in current setups, like Rydberg atom arrays and trapped ions. The use of Walsh functions [38]–[40] is at the heart of our proposal. This set of discrete orthogonal signal functions allows us to easily parametrize single qubit rotations both in space and time, and to obtain a simple Hamiltonian dynamics corresponding to the restriction of the resource Hamiltonian \(H_R\) on a certain graph structure. Based on simple graph decomposition arguments, we then show how to program multiple sequences of pulses to obtain an arbitrary two-body target Hamiltonian \(\tilde{H}\). This allows in turn to realize arbitrary quantum computation, using for instance parallel two-qubit quantum gates. Similar pulse sequence design ideas based on Walsh functions have been previously considered to engineer Ising couplings [41]–[44]. Importantly, we also show that the orthogonality property of Walsh functions offers ‘direct’ analytical strategies for decoupling the system from unwanted sources of errors, without adding additional pulses. This is in contrast to previous DAQC protocols where these errors have been shown to be detrimental [18], and makes our proposal promising for realizing quantum computation. We analytically derive these robust pulse sequences in the case of imperfect pulses, as well as in presence of background fields.

In Sec. 2 we outline the general protocol and discuss how it can realize arbitrary quantum computation and simulation, along with a resource analysis in terms of number of pulses needed and errors, and finally present a controlled approximation scheme to reduce the number of pulses. In Sec. 3 we translate robust strategies that were developed in the context of global pulses [27] to the present case of locally-addressable pulses, and we derive analytically the robustness conditions for Walsh pulse sequences. We then provide three examples showing the versatility of the protocol: engineering dynamics of quantum spin models beyond the resource Hamiltonian (Sec. 4), the implementation of a surface code (Sec. 5), and quantum optimization algorithms (Sec. 6). Finally, in Sec. 7 we discuss in detail the possibility of implementing our protocol in state-of-art Rydberg atom arrays and trapped ions experiments.

2 Introducing Walsh sequences↩︎

In this section we present our protocol in a way that is thought to be self-contained for most readers. We start by reviewing the pulse sequences approach to Hamiltonian engineering, in particular average Hamiltonian theory (AHT) and Trotter errors. We then introduce the first key result of our work: pulse sequences parametrized with Walsh functions result in a simple structure of the engineered Hamiltonian. This provides us with an experimental recipe, using multiple such pulse sequences, to engineer arbitrary target Hamiltonians (or quantum circuits).

2.1 Pulse sequences for Hamiltonian engineering↩︎

In this subsection we review the theory of Hamiltonian engineering based on pulses sequences. Let us consider for the moment a system consisting of \(N\) spins \(1/2\) (i.e., qubits) interacting via a static Hamiltonian \(H_R\) (which we refer to as resource Hamiltonian), on which we can also apply pulses represented by instantaneous single-qubit gates \(P = \bigotimes_{i=1}^N p_i\). Note already that in an experimental setup, such gates can be approximately realized with external laser pulses, which are not instantaneous. We study the action of realistic pulses in Sec. 3, and we derive conditions under which the instantaneous pulse approximation is accurate. Commonly used pulses consist of Pauli rotations \(p_i = R_{O}^{\alpha} = e^{-i\alpha O/2}\), \(O \in \{X,\;Y,\;Z\}\) being a Pauli operator acting on the qubit \(i\). For \(\alpha = \pm \pi\), they amount to directly applying the operator \(O_i\). Several works about dynamical Hamiltonian engineering focus on the case in which \(p_i = p_j \;\forall i,j\), which means that the pulse is acting in the same way on every qubit. We will refer to these pulses as global pulses.

In this work, we will focus instead on local \(\pi-\)pulses, where \(p_i \in \{1, X, Y, Z\}\) and we can have \(p_i \neq p_j\). Let us now consider a pulse sequence \(\{P^{(k)}\}=\{\bigotimes_{i=1}^N p_i^{(k)}\}\) of length \(n\), which is repeated after a period of time \(\tau\). We partition the period \(\tau\) in \(n\) intervals of length \(\tau/n\), and we apply the pulse \(P^{(k)}\) at the beginning of the interval \(k\), and the inverse pulse \((P^{(k)})^{-1}\) at the end of the interval \(k\) (see Fig. 2). The unitary evolution to which the system is subject over a period is \[\tilde{U}(\tau) = \prod_{k = 1}^n (P^{(k)})^{-1} e^{-i H_R \frac{\tau}{n}} P^{(k)} \label{eq:Utau}\tag{1}\] where the product is applied right-to-left (\(\prod_{k=1}^n A_k = A_n...A_2 A_1\)). Note that in an experiment, one may consider to merge the two pulses \((P^{(k-1)})^{-1}\) and \(P^{(k)}\), and simply implement \((P^{(k-1)})^{-1}P^{(k)}\) at the time \(\frac{(k-1)}{n}\tau\). However, it will be easier for us to derive robustness conditions in Sec. 3 keeping those two pulses as distinct objects. Using the property of matrix exponentials \(A^{-1} e^B A = e^{A^{-1}BA}\), we obtain \[\label{eq:productformula} \tilde{U}(\tau) = \prod_{k = 1}^n e^{-iH^{(k)}\tau/n},\;H^{(k)} = (P^{(k)})^{-1} H_R P^{(k)}\tag{2}\] where \(H^{(k)}\) are called the toggling-frame Hamiltonians. The evolution of the system over a period is then described by a product formula. Product formulas are connected to the evolution under a static Hamiltonian using the Floquet-Magnus expansion in small \(\tau\to 0\) [45]–[47]. At the leading order, this expansion yields the average Hamiltonian theory (AHT), in which \[\tilde{U}(\tau) = e^{-i\tau\tilde{H}} + O(\tau^2),\;\tilde{H} = \frac{1}{n} \sum_{k=1}^n H^{(k)}\] and \(\tilde{H}\) is the average Hamiltonian. To simulate a Hamiltonian dynamics up to a time \(T\), usually the same product formula is applied for a fixed small time \(\tau\), and then repeated \(T/\tau\) times, \(T/\tau\) being an integer. The unitary error \(\mathcal{E}(\tau, T) := ||(\tilde{U}(\tau))^{T / \tau} - e^{-iT\tilde{H}}||\) can be quantified as [48] \[\label{eq:lloyd95formula} \mathcal{E}(\tau, T) = \frac{\tau T}{2n^2} \Big{|} \Big{|} \sum_{k<l} [H^{(k)}, H^{(l)}] \Big{|} \Big{|} + O(\tau^2 T)\tag{3}\] where we take \(|| \cdot ||\) to be the spectral norm, i.e. the largest singular value of the operator. Eq. 2 is also referred to as the first-order \(p = 1\) Trotter formula for \(\tilde{H}\) [49]. It is possible to engineer our pulse sequence to realize higher orders Trotter formulas, which feature a smaller error in simulating \(\tilde{H}\). For example, the \(p=2\) Trotter formula requires a pulse sequence of length \(n' = 2n\) (and the time duration of an interval is \(\tau/n' = \tau/(2n)\)) with \(H^{(n + k)} \leftarrow H^{(n - k + 1)}\) for \(k<n\), yielding \(\mathcal{E} = O(\tau^2 T)\) [49], [50]. \(p>2\) Trotter formulas require to evolve the system with the negative resource Hamiltonian \(-H_R\), which is not easily accessible in quantum simulation platforms.

2.2 Hamiltonian engineering from Walsh functions↩︎

In this subsection we start explaining the core ideas for our proposal. Using the formalism of Walsh functions, we design pulse sequences \(\{P^{(k)}\}\) that can realize target evolutions \(\tilde{H}\). From now on, we will consider a general \(XY\) model as a resource Hamiltonian, i.e. \[\label{eq:H95resource} H_R = \sum_{i<j} (J^X_{ij} X_i X_j + J^Y_{ij} Y_i Y_j) \\\tag{4}\] This covers in particular two very common types of interactions spin systems: Spin-exchange interactions \(J^X_{ij}=J^Y_{ij}\) as for instance in dipolar Rydberg atom arrays, and Ising interaction \(J^Y_{ij}=0\) as between trapped ions. More details about experimental implementation of such Hamiltonians combined with local pulses are given in Sec. 7. While our analytical formulas will be general, we will consider for all examples power-law decaying spin-exchange interactions \[\label{eq:LR95coupling} J^X_{ij} = J^Y_{ij} = - \frac{J}{r_{ij}^{\alpha}} \;.\tag{5}\] We now proceed by noting that, for our resource Hamiltonian Eq. 4 , and using Pauli pulses \(p_i^{(k)}\in\{1_i,X_i,Y_i,Z_i\}\) the expression of the toggling-frame Hamiltonians \(H^{(k)}\) simplifies to \[\begin{align} H^{(k)} &= \sum_{O = X, Y}\sum_{i<j}J_{ij}^{O}((p_i^{(k)})^{-1} O_i p_i^{(k)}) ((p_j^{(k)})^{-1} O_i p_j^{(k)}) \nonumber \\ &= \sum_{O = X, Y}\sum_{i<j} J_{ij}^{O} \left( s_{i,O}^{(k)} s_{j,O}^{(k)} O_i O_j \right). \end{align}\] Here \(s_{i,O}(k) = 1\) (\(-1\)) if \(p_i\) commutes (anticommutes, respectively) with \(O_i\). Importantly, for each qubit \(i\), it is easy to verify that we can always program the pulses \(p_i^{(k)}\) to choose the two signs \(s_{i,X}^{(k)}\), \(s_{i,Y}^{(k)}\) \[\begin{align} \label{eq:walsh95rules} p_i^{(k)} &= 1_i\text{ if } s_{i,X}^{(k)}=s_{i,Y}^{(k)}=1 \\ &= X_i\text{ if } s_{i,X}^{(k)}=-s_{i,Y}^{(k)}=1 \\ &= Y_i\text{ if } -s_{i,X}^{(k)}=s_{i,Y}^{(k)}=1 \\ &= Z_i\text{ if } s_{i,X}^{(k)}=s_{i,Y}^{(k)}=-1. \end{align}\tag{6}\]

For a sequence of \(k=1,\dots,n\) pulses, we can write the average Hamiltonian as \[\label{eq:Haverage} \tilde{H} = \sum_{i<j} ((s_{i,X}|s_{j,X})J^X_{ij}X_i X_j + (s_{i,Y}|s_{j,Y})J^Y_{ij}Y_i Y_j) \;.\tag{7}\] with the scalar product \((s_1|s_2)= \frac{1}{n}\sum_{k=1}^{n} s_1^{(k)}s_2^{(k)}\). The expression Eq. 7 brings us to the key part of our proposal: As we can choose the signs \(s_{i,X/Y}\) arbitrarily, let us define \(s_{i,X}^{(k)} = w_{x_i}^{(k)} \;,\;s_{i,Y}^{(k)} = w_{y_i}^{(k)}\) where \(w_a^{(k)}\in \{-1,1\}\), with \(a\) integer and \(k\in \{1,\dots,n\}\), are Walsh functions, which feature the orthonormality property \((w_a|w_b) = \delta_{a,b}\) [38]–[40] (see App. 9 for more details). In the case of Ising Hamiltonians, the same results hold by designing sequences of \(X\) pulses using Walsh functions [41]–[44]. We obtain the final form of \(\tilde{H}\) which we will use throughout this work \[\label{eq:single95walsh95h} \tilde{H} = \sum_{i<j} \left(G^X_{ij}J^X_{ij}X_iX_j + G^Y_{ij}J^Y_{ij}Y_iY_j\right)\tag{8}\] where \(G^X_{ij}=\delta_{x_i, x_j}\) (\(G^Y_{ij}=\delta_{y_i, y_j}\)). This means that we can propose a general experimental recipe to program the pulses in such a way that two spins \(i\) and \(j\) interact via \(X_iX_j\) couplings (or not) by assigning them the same Walsh index \(x_i=x_j\) (\(x_i\neq x_j\), respectively). The same result applies independently to \(Y_iY_j\) couplings. This is illustrated in Fig. 1. We show next the universality

and flexibility of this approach (and later its robustness).

It will be useful for what follows to use the language of graphs: \(G^X_{ij}\) and \(G^Y_{ij}\) are the adjacency matrix of the graph \(\mathcal{G}_1^X\) (\(\mathcal{G}_1^Y\)). We can view these two graphs as made with qubits as vertices, which are connected by a link if their Walsh indices match. Formally, the links of these graphs are \[\label{eq:graph95equation} \begin{align} E(\mathcal{G}_1^X) = \{(i,j)\;|\;x_i=x_j\},\\ E(\mathcal{G}_1^Y) = \{(i,j)\;|\;y_i=y_j\} \;. \end{align}\tag{9}\] Note that we have obtained a mapping between a graph structure of interactions and a local pulse sequence which realizes it. We will refer to each pulse sequence designed in this way as a Walsh sequence (or a Walsh pulse sequence, as in the title of the article). An explicit example for \(N=2\) qubits is shown in Fig. 2 a)-c). We choose \(x_1 = x_2\) in order to couple the two qubits with respect to the \(XX\) interaction, and \(y_1 \neq y_2\) to decouple them with respect to the \(YY\) interaction. The choice of Walsh indices \((x_i, y_i)\) for each qubit associates two Walsh functions to it, from which we derive the resulting local pulse sequence \(p_i^{(k)}\). Finally, we note that the length \(n\) of a Walsh sequence, i.e. the number of pulses to be applied in an experiment, is equivalent to the size of the Walsh function with the largest index, which is (see App. 9) \[\label{eq:length95pulse95seq} n = 2^{\lceil \log_2 (\text{max}_i(x_i, y_i)+1)\rceil} \;.\tag{10}\] The longest Walsh sequence is obtained in the case all the qubits are required to be decoupled, for which \(x_i=y_i=i\), which results in \(n=O(N)\) (in particular, \(N\le n < 2N\)).

2.3 Universal quantum computing from Walsh sequences↩︎

In this section we illustrate how to engineer a target spin Hamiltonian using \(Q\ge 1\) Walsh sequences. From the definition in Eq. 9 , we note that the interaction graphs \(\mathcal{G}_1^X\) (\(\mathcal{G}_1^Y\)) cannot be arbitrary. For instance, one can readily see that it is not possible to couple two atoms \(i_1\) and \(i_2\) to a third atom \(i_3\), without coupling \(i_1\) to \(i_2\). In order to engineer more general average Hamiltonians, one thus needs to consider a pulse sequence made by alternating \(q=1,\dots,Q\) Walsh sequences of length \(n_q\), each one with a period \(\tau_q\). The global length of the resulting pulse sequence will be \(n = \sum_{q=1}^Q n_q\), which thus scales at most as \(n_q \le O(NQ)\).

In AHT, \(Q\) Walsh sequences results in the unitary operator [49], [50] \[\label{eq:multiple95walsh} \tilde{U}(\tau) = \prod_{q=1}^Q \tilde{U}^{(q)}(\tau_q) = e^{-i\tau \sum_q \frac{\tau_q}{\tau} \tilde{H}^{(q)}} + O(\tau^{p+1}) \;.\tag{11}\] where the effective Hamiltonian \(\sum_q (\tau_q/\tau) \tilde{H}^{(q)}\) can now be engineered for any possible spin Hamiltonian with two-body interactions, and thus in particular any quantum circuits of two-qubit gates (which is universal for quantum computing). While we report the complete proof of universality in App. 10, let us sketch here the main idea.

First, assuming here for the rest of this subsection \(J_{ij}\neq 0\), and considering a target Hamiltonian of the form \[H_{\mathrm{T}} = \sum_{i<j}(\tilde{J}^X_{ij} X_i X_j + \tilde{J}^Y_{ij} Y_i Y_j ),\] with \(\tilde{J}^X_{ij},\tilde{J}^Y_{ij}\ge 0\), we define \(g^O_{ij} := \tilde{J}^O_{ij} / J^O_{ij}\), which we decompose as \[\label{eq:decomposition} g^O_{ij}=\sum_{q} \frac{\tau_q}{\tau} [G^O_q]_{ij}\tag{12}\] where \((G^O_q)_{ij}\in \{0,1\}\) are adjacency matrices of graphs whose structure is defined by Eq. 9 . Therefore, we can implement each contribution \(\propto [G^O_q]_{ij}\) with a Walsh sequence.

The decomposition Eq. 12 is always possible, and can be performed in various ways. The simplest situation arises when \(g^{O}_{ij}\) is homogenenous (\(g^{O}_{ij} \in \{0, 1\}\) \(\forall i, j\)). This corresponds to various physically relevant scenarios both in quantum simulation and quantum computing, as we will illustrate below. In this case, we can simply decompose \(g^{O}_{ij}\) in graphs of degree 1 with adjacency matrices \([G^O_q]_{ij}\) and \(\tau_q=\tau\), i.e. in each such graph a spin \(i\) is coupled to at most another spin. Based on Eq. 8 , this corresponds to choosing for each \(q\) a Walsh sequence \(x_i=x_j\) (\(y_i=y_j\)) whenever \([G^X_q]_{ij}=1\) (\([G^Y_q]_{ij}=1\)). As shown in Sec. 4, 5, this decomposition can be directly ‘seen’ analytically in various relevant examples, when the number of Walsh sequences is \(Q = O(1)\). Alternatively, as shown in App. 10, such decomposition can also be achieved in full generality by either

1. starting from a decomposition of the complete graph [51]–[53], resulting in \(Q=O(N)\) sequences.

2. employing a numerical algorithm (see Algorithm 7 in the appendix) to obtain a smaller number of Walsh sequences, i.e. \(Q = O(d_\text{max})\), \(d_\text{max}< N\) being the degree of the graph with adjacency matrix \(g_{ij}^O\).

When \(g^O_{ij}\notin \{0,1\}\) is not homogenous instead, the graph \([G^O_q]_{ij}\) needs to be simulated for different times \(\tau_q\), with possible repeated links. For instance, a strongly inhomogenenous and connected \(g^O_{ij}\) can require to engineer as many Walsh sequences as the number of interacting links, resulting in the maximal value \(Q=(N-1)N/2\). This would imply that a random target Hamiltonian would require up to \(\sum_{q=1}^Q n_q \lesssim O(N^3)\) pulses. At the same time, a target Hamiltonian ‘commensurate’ with the resource Hamiltonian, e.g. \(g_{ij}^O\) with \(O(1)\) different weights, will result in the much more favorable scaling of \(\sum_{q=1}^Q n_q \lesssim O(Nd_{\text{max}})\) pulses.

Finally, in order to extend the construction above to more general couplings, it is sufficient to add a set of \(Q\) pulses \(\{P_{\text{set}}^q\}\), where \(P_{\text{set}}^q\) is applied before the \(q-\)th Walsh sequence, and \(\prod_{q=1}^QP_{\text{set}}^q = I\). More details about this construction are given in App. 10, and examples are given in the following sections.

To conclude, we briefly discuss the connection between Walsh pulse sequences and quantum circuits. In fact, any layer of a quantum circuit made of two-qubit gates acting on different pairs of qubits can be realized by the dynamics of a single pair-wise interacting Hamiltonian \(\tilde{H}\). Since the connectivity of such a Hamiltonian will be represented by a graph of degree 1, this set of gates can be implemented in parallel via Walsh sequences (or a single Walsh sequence, in some cases). Each one of this Walsh pulse sequences will require \(n = O(N)<2N\) pulses, and they will be applied in total \(T/\tau\) times, where \(\tau\) fixes the Trotter error, and \(T\) is the longest evolution time needed by a two-qubit gate to be realized. An explicit example will be discussed in Sec. 5. Finally, let us note that allowing for groups of \(N_g>2\) qubits to interact might be used to realize \(N_g-\)qubit gates. This implies that not only Walsh pulse sequences can find application in realizing quantum gates in systems with collective interactions, but also that these interactions could be exploited to realize novel gates, using our framework.

2.4 Trotter errors↩︎

We now assess the Trotter errors from realizing the average Hamiltonian \(\tilde{H}\) with \(Q\) Walsh sequences, each one with a period \(\tau_q\). Starting from Eq. 3 , we derive a bound for unitary Trotter errors for one-dimensional systems with power-law decaying interactions as in Eq. 5 , which reads for \(p=1\) (see App. 11) \[\label{eq:trotter95bound} \begin{cases} \mathcal{E}(\{\tau_q\}, T) \leq a_{\alpha} N(JT)(J\sum_q\tau_q) \;\text{if} \;\alpha > 1 \\ \mathcal{E}(\{\tau_q\}, T) \leq b_{\alpha} N^{3 - 2\alpha}(JT)(J\sum_q \tau_q) \;\text{if} \;\alpha < 1 \end{cases}\tag{13}\] where \(a_{\alpha}, b_{\alpha} = O(1)\) are functions of only \(\alpha\). Therefore, to simulate a dynamics with a Trotter error \(\mathcal{E} \leq \epsilon\), we can fix the period of the total cycle to be \(J\sum_q\tau_q = \kappa_{\epsilon}N^{-1}(JT)^{-1}\) for \(\alpha > 1\), or \(J\sum_q\tau_q = \kappa_{\epsilon}N^{2\alpha - 3}(JT)^{-1}\) for \(\alpha < 1\), \(\kappa_{\epsilon} = O(\epsilon)\) being a constant that can be tuned to obtain the desired error value. We believe this can be generalized easily to higher spatial dimensions and Trotter orders \(p>1\), using the techniques presented in Ref. [49]. This estimate serves as a performance guarantee for a generic Hamiltonian engineering protocol, and shows that the required number of number of pulses is polynomial with system size.

2.5 Reducing the number of pulses with a cut-off distance↩︎

We conclude this section by presenting a controlled approximation scheme which allows to reduce the number of pulses \(n\) in Walsh sequences. This is important in light of an experimental implementation since upon fixing the physical limit for pulse repetition frequency \(\tau/n\), decreasing \(n\) allows us to reduce \(\tau\), and therefore Trotter errors. In particular, let us consider the situation in which we have a target Hamiltonian with finite-range interactions, i.e. it exists \(r = O(1)\) such that \(\tilde{J}^O_{ij} = 0\) for \(r_{ij}>r\) for \(O = X, Y\).

Let us consider the associated Walsh sequence with length \(n = O(N)\). The approximation which we rely on consists in choosing a cut-off distance \(\Lambda_r\) for which we assume \(J^O_{ij} \approx 0\) for \(r_{ij}>\Lambda_r\) in the resource Hamiltonian. For one-dimensional systems, it is possible to prove that, under this assumption, the Walsh indices can be chosen to be bounded as \(x_i,y_i< \Lambda_w\), with \[\Lambda_w := \text{max}_i (x_i, y_i) < 2\Lambda_r,\] which corresponds to the largest index of the ‘compressed’ Walsh sequence, of which the length is given by Eq. 10 . A proof of this statement and a Walsh indices assignment algorithm are given in App. 12. We believe that a similar approach can also be adapted to higher dimensions.

Let us now quantify the errors made by the cutoff approximation of the resource Hamiltonian. We have \[(\tilde{U}(\tau))^{T/\tau} = e^{-iT(\tilde{H} + H_E)} + O(\tau^p T)\] where \[H_E = \sum_{i, j, |i-j|>\Lambda_r}(J^X_{ij}\delta_{x_i,x_j}X_i X_j + J^Y_{ij}\delta_{y_i,y_j} Y_i Y_j)\] is a perturbation to our target Hamiltonian \(\tilde{H}\). For small enough \(\tau||H_E||\), the unitary error induced by \(H_E\) in a single Trotter step is \[\tilde{U}(\tau) - e^{-i\tau\tilde{H}} = \tau H_E e^{-i\tau\tilde{H}} + O(\tau^2)\] and therefore if \(|| H_E || = O(\tau)\) the residual error will be still proportional to \(p=1\) Trotter errors. For \(J^O_{ij} = J/|i-j|^{\alpha}\), \(\alpha>3\), we provide an analytical guarantee that the errors given by the choice of a cut-off can be made arbitrarily small. In particular, we find that \(||H_E|| = \epsilon \tau\) is ensured by choosing \(\Lambda_w = O((N/\epsilon)^{\frac{2}{\alpha - 1}})\). We provide a numerical example for the cut-off scheme in Sec. 4 for \(\alpha = 0.2,\;3\), finding that for \(\alpha = 3\) the cut-off errors can be made much smaller that Trotter errors.

3 Walsh sequences are robust against experimental imperfections↩︎

The formalism of Walsh sequences presented in the previous section provides us with analytical pulse sequences to decouple/couple selectively groups of qubits to engineer specific Hamiltonian terms. Let us show now that the use of Walsh functions, based on similar ‘decoupling strategies’, also allows us to suppress the effect of experimental imperfections. This step is essential to implement reliable quantum computation and simulation based on pulse sequences, as alternating global Hamiltonian evolution for short times with single-qubit imperfect pulses can yield a quick build-up of errors [18]. In the following, we provide analytical evidence of robustness of our protocol with respect to rotation angle errors, finite pulse duration, and presence of background magnetic fields. In particular, the latter situation will be useful to draw connections between our formalism and the research topic of dynamical decoupling.

3.1 Rotation angle errors: From single averaging to double averaging↩︎

Let us first consider errors in the rotation angle applied by the pulses during one Walsh sequence. This kind of errors can be related to a miscalibration of the single-qubit pulses [18], or to over- or under-illumination of certain qubits by the applied laser [35], [36]. We model a pulse with rotation angle errors as \[\begin{align} p_i^{(k)} = e^{-is_i(\pi + \delta_i)[O]_i/2}, \end{align}\] where \(s_i = \pm 1\) when the pulses are implemented as a \(\pi\) rotation (\(-\pi\) rotation, respectively), \(\delta_i\) is a single-qubit rotation angle error, and \([O]_i \in \{1_i, X_i, Y_i, Z_i\}\) represents the pulse we wish to apply on the \(i-\)th qubit (i.e., \(p_i^{(k)}=[O]_i\) in absence of errors \(\delta_i=0\)). In the following, we will take advantage of the freedom to choose the type of pulses \(s_i\) to design robust Walsh sequences.

As shown in App. 13, the effect of rotation errors under Walsh sequences leads at order \(O(\delta)\) to a remarkably simple expression of the average Hamiltonian \(\tilde{H}+\tilde{H}_{\text{RA,err}}\), with \[\begin{align} \tilde{H}_{\text{RA,err}} &= \sum_{O = X, Y}&\sum_{i<j}J^O_{ij}\left(\delta_i s_i \mathcal{O}_i^j O_j + \delta_j s_j O_i \mathcal{O}_j^i \right), \label{eq:HRAerr} \end{align}\tag{14}\] and \[\begin{align} \mathcal{X}_i^j =& \frac{1}{n}\sum_{k=1}^n \frac{1 - w_{x_i}^{(k)}}{2}w_{x_j}^{(k)}\left(\frac{1 - w_{y_i}^{(k)}}{2} Y_i - \frac{1 + w_{y_i}^{(k)}}{2}Z_i \right) \\ \mathcal{Y}_i^j =& \frac{1}{n}\sum_{k=1}^n \frac{1 - w_{y_i}^{(k)}}{2}w_{y_j}^{(k)} \left(\frac{1 + w_{x_i}^{(k)}}{2} Z_i - \frac{1 - w_{x_i}^{(k)}}{2}X_i \right) \;. \end{align}\]

The simple form of Eq. 14 suggests to build a robust sequence using a second averaging strategy [27], [54] based on varying with time the sign functions \(s_i\). Let us first recall that in the previous section, we considered repeating the same Walsh sequence identically \(T/\tau\) multiple times to propagate the Hamiltonian until a certain evolution time \(T\). To make this evolution robust, we will realize instead the evolution \[\underbrace{ (\tilde{U}(\tau))^{T/\tau} } _{\text{single averaging}} \rightarrow \underbrace{\prod_{l=1}^{T/\tau}\tilde{U}^{(l)}(\tau)} _{\text{double averaging}}.\] where the sign \(s_i^{(l)}\) becomes a function of the index \(l\), and we define them to be Walsh functions as \(s_i^{(l)} = w_{e_i}^{(l)}\). Here we choose \(e_i>0\). Note that this second averaging does not require the introduction of additional pulses in our protocol. In this case, averaging the expression Eq. 14 of the error \(\tilde{H}_{\text{RA,err}}\) over the \(T/\tau\) evolutions, one obtains that the average error Hamiltonian at order \(O(\delta)\) vanishes exactly, i.e. \(\sum_{l=1}^L\tilde{H}_{\text{RA,err}}^{(l)} = 0\), \(L\) being the periodicity of the sign function with the largest Walsh index \(e_i\). This demonstrates how a simple choice of pulsing angles can make Walsh pulse sequences resilient to calibration errors. We report in App. 13 numerical data illustrating the use of robust Walsh sequences in the presence of rotation angle errors.

3.2 Double averaging finite pulse duration errors↩︎

So far we have considered the qubits to be not evolving with the resource Hamiltonian during the application of the pulses. While this assumption is certainly valid in ion chains, where laser mediated interactions can be switched on and off, the situation may be different in Rydberg atom arrays in the ‘frozen gas’ regime, where interactions are always on, see Sec. 7. At the same time, it might be desirable to keep the evolution under the resource Hamiltonian on even during the application of pulses, to avoid ramp-up and ramp-down errors [11], [18]. Let us consider thus the effect of the finite duration of a pulse and define the single-qubit Hamiltonian term \[H_p = \frac{\pi}{2t_p}\sum_i s_i [O]_i\] generating the pulse \(P^{(k)} = e^{-it_pH_p}\), where \(t_p\) is the time duration of a pulse, and we assume to a have a square-wave pulse. So, during the interval \(k\), the system is subjected to the following unitary evolution \[U^{(k)}(\tau) = e^{-it_p(H_R - H_p)}e^{-i\left(\frac{\tau}{n} - 2t_p \right)H_R}e^{-it_p(H_R + H_p)} \;.\] In the limit of an infinitely short pulse \(t_p=0\), we recover the idealized limit of the previous section, c.f. Eq. 1 . Upon defining the finite pulse duration error parameter \(\epsilon_{\text{FP}}=2nt_p/\tau\), which amounts to the time portion of the interval \(k\) during which pulses are applied, we obtain the following average Hamiltonian (see App. 14 and Refs. [27], [55]) \[\begin{align} \tilde{H}_{\text{FP}} &= \left(1 - \epsilon_{\text{FP}} \right)\tilde{H}\;+ \nonumber \\ &+\epsilon_{\text{FP}}\frac{1}{n}\sum_{k=1}^n \frac{1}{t_p}\int_0^{t_p}dt\;(Q^{(k)}(t))^{-1} H_R Q^{(k)}(t) \end{align}\] where \(Q^{(k)}(t) = e^{-itH_p}\) can be interpreted as a partially applied pulse. Notice that we can repeat the above treatment for the case of pulses with different wave-shapes.

As for rotation errors, the expression of the average Hamiltonian takes a simple form when parametrizing the pulses with Walsh functions (under the assumption here that \(x_i\neq 0\), \(y_i\neq 0\)), see App. 14. When using a second averaging with time changing signs \(s_i^{(l)}=w_{e_i}^{(l)}\), and \(e_i\neq e_j\neq 0\) we first obtain \[\begin{align} \tilde{H}'_{\text{FP}} =\left(1 - \frac{5 \epsilon_{\text{FP}}}{8} \right)\tilde{H} + \frac{3\epsilon_{\text{FP}}}{8}H_R \;. \end{align}\] At this point, it is possible to further deform the Walsh sequence to simulate \(\tilde{H}\) in AHT. This is achieved by 1) realizing the interval \(k=1\), the ones in which \(H^{(k=1)} = H_R\), for a reduced time \(\frac{\tau}{n} - \frac{3nt_p}{4}\), and 2) correcting the total simulation time \(T\) as \(T \rightarrow T\left(1 - \frac{5 \epsilon_{\text{FP}}}{8}\right)^{-1}\). These robustness conditions, and time deformations, completely remove any contribution in AHT of the finite pulse duration effects, and at the same time they remove the dominant contribution in AHT for rotation angle errors. To corroborate our analysis, a numerical example is given in Sec. 4 and App. 13, showing both the correction of finite-pulse duration effect alone, and together with rotation angle effects.

3.3 Suppression of background magnetic fields and connections to dynamical decoupling↩︎

As a last source of errors, we consider single-qubit noise in the resource Hamiltonian. To decouple a qubit system from single qubit noise (or other sources of unwanted interaction) is referred to in the literature as dynamical decoupling [56], [57]. Recently, dynamical decoupling has been connected to Hamiltonian engineering [28], [35], [36], in order to increase accuracy in quantum simulation protocols beyond the simulation of the resource Hamiltonian. Here, we show that this idea further extends to the case of universal quantum evolutions treated here.

We consider a resource Hamiltonian which also has space-dependent unwanted single-qubit fields \[H_R' = H_R + H_{\text{ext}}\] with \[H_{\text{ext}} = \sum_i (h^x_i X_i + h^y_i Y_i + h^z_i Z_i)\] which can model disorder fields, or slowly-varying noise. After a Walsh sequence, we obtain again a simple expression (see App. 15) \[\tilde{H}_{\text{ext}} = \sum_i ( \delta_{0, x_i} h^x_i X_i + \delta_{0, y_i} h^y_i Y_i + \delta_{x_i, y_i} h^z_i Z_i) \;.\] In this case, there is no need for double averaging: Full decoupling can be then achieved after a single Walsh sequence by 1) not using the Walsh function \(w_0^{(k)}\) and 2) using different Walsh indices for the \(X\) and \(Y\) operators on the same qubit. This can be always be achieved by enumerating the \(y_i\) starting from the number after the largest \(x_i\), which doubles the number of Walsh indices (hence the length of the resulting Walsh sequences). Note that this robustness condition can be applied simultaneously with the ones described in the two previous subsections.

4 Application to Hamiltonian engineering↩︎

We are now ready to illustrate examples of our approach. In what follows, we consider a long-range resource Hamiltonian, with couplings defined from Eq. 5 . We start with the implementation of the Ising model with nearest-neighbor couplings, used as a paradigmatic example to benchmark the performances of our approach. This is for instance related to the problem of reducing interactions from long-range to short-range in quantum simulators [58]. This is also useful in quantum computing, as the model allows us to prepare graph states, in particular cluster states [59], [60].

Our target model is defined as \[\label{eq:NN95Ising} H_{\text{Ising}} = -J\sum_{i=1}^{N-1} X_i X_{i+1}.\tag{15}\] This corresponds to the two following rescaled interactions defined in Sec. 2.3 \[g_{ij}^X = \delta_{i, j\pm 1},\;g_{i j}^Y = 0 \;.\] which represent a homogenenous graph of degree 2, and can be straightforwardly decomposed in two graphs of degree \(1\) \[[G_{q = 1}^X]_{2i-1, 2i} = 1,\;[G_{2}^X]_{2i, 2i+1} = 1\] while all the other terms are zero. Finally, \([G_{1}^X]_{i, j}\) can be realized by choosing the Walsh indices \(x_i = \lfloor (i-1)/2 \rfloor\), and \([G_2^X]_{ij}\) by choosing \(x_i = \lfloor i/2 \rfloor\). \([G_{q}^Y]_{i, j}\) instead are always realized with \(y_i = i-1\). In Fig. 3a-c) we show an example of the decomposition and Walsh indices assignment for \(N = 4\).

We benchmark numerically the protocol in Fig. 3d)-g). In particular, we assess Trotter errors numerically studying the state fidelity \(\mathcal{F}=|\langle\psi_{\text{cluster}}|\psi_{\text{sim}}\rangle|\) for \(T = \pi / (4J)\) and initial state \(\ket{\psi_0} = \bigotimes_i \ket{0}_i\). This choice outputs cluster states \(\ket{\psi_{\text{cluster}}}\) in the case of perfect Ising dynamics (\(\tau = 0\)) [59], [60]. The output state fidelity is related to the unitary Trotter errors bounded by Eq. 13 , since the latter is a bound for the trace distance between the simulated and the target states [8], which for pure states coincides with the following fidelity error [61] \[||(\tilde{U}(\tau))^{T/\tau} - e^{-iT\tilde{H}}|| \geq \sqrt{1 - \mathcal{F}} \;. \label{eq:fide95trotter}\tag{16}\] To check the output state fidelity, we numerically simulate the pulse sequence protocol using exact numerical methods, with size varying from \(N=8\) to \(N=16\), and various periods \(\tau\). We consider the two cases \(\alpha = 3\) and \(\alpha = 0.2\). The results are shown in Fig. 3d)-g). In panels d,e), we plot the fidelity error per qubit \((1 - \mathcal{F})/N\), for pulse sequences that correspond to \(p = 1\) (round markers) and \(p = 2\) (triangular markers) Trotter formulas. For the behavior with the timestep \(\tau\), we observe a power-law decay of the fidelity error per qubit \(\sim \tau^{2p}\), which is for \(p=1\) consistent with the upper bounds of Eqs. 13 16 . For the behavior with the size \(N\) instead, in the case \(\alpha = 3\) and \(p = 1\) the fidelity error per qubit appears to be independent of the size \(N\), whereas for \(\alpha=0.2\) it increases with \(N\). To better understand the scaling with \(N\), we consider the case \(p = 1\), where can expect scalings of at most \(\sim N^2\) for \(\alpha=3\) and \(\sim N^{2(3 - 2\alpha)}=N^{5.2}\) for \(\alpha =0.2\), based on the results of Eq. 16 and Eq. 13 . By performing a data collapse, we observe a scaling of the fidelity error of \(\sim N\) for \(\alpha =3\) and \(\sim N^{2.6}\) for \(\alpha=0.2\).

This shows that the Trotter bound predicted is loose in this case.

In the case presented, the length of the pulse sequence considered was \(n=N\), which was required to cancel completely the \(YY\) component of the resource Hamiltonian. As seen in Sec. 2.5, we can reduce the pulse sequence length drastically using a cut-off distance \(\Lambda_r\) above which we can neglect interactions. Using the Walsh indices \(x'_i = x_i\;\text{mod}\;\Lambda_r\) (\(y'_i = y_i\;\text{mod}\;\Lambda_r\)), we obtain the fidelity shown in Fig. 4b). As expected, for the short-range case \(\alpha = 3\), the introduction of a cut-off \(\Lambda_r<N\) induces a very modest error \(1 - \mathcal{F} \sim 10^{-5}\) for the worst case considered \(N = 14\), \(\Lambda_w = 8\). For the long-range case \(\alpha = 0.2\), the notion of the cut-off distance is meaningless, since, as soon as \(\Lambda_r<N-1\), the fidelity error becomes independent of \(\tau\) and of order one.

Finally, we check numerically the effect of pulse induced errors, and the performances of the robustness conditions presented in Sec. 3. We plot the results in Fig. 4c). While the black lines represent the fidelity error for perfect pulses, the red and blue lines represent the fidelity error when pulse errors are present (see Apps. 13-14 for the effects studied separately). We model the rotation angle errors by sampling a qubit-dependent rotation angle error of \(\delta_i\), which we sample from a uniform distribution in the interval \([-2\epsilon_{RA} J, 2\epsilon_{RA} J]\). At the same time, we consider a finite duration for the pulses of \(t_p = \frac{\tau}{2n}\epsilon_{\text{FP}}\). The fidelity errors are plot for various values of the errors \(\epsilon_{\text{RA}} = \epsilon_{\text{FP}} = \epsilon\), \(N = 6\) qubits and \(\alpha = 1.2\). We observe that the introduction of the robustness conditions using the double averaging strategy discussed in Sec. 3 (blue lines) allows to reduce the fidelity error of several orders of magnitude compared to a single averaging (red lines).

Note that this example for the Ising model can be straightforwardly extended to the situation of translation-invariant \(XXZ\) type Hamiltonians on arbitrary spatial geometry with finite connectivity. This can for example be done by implementing separately the \(XX\), \(YY\) and \(ZZ\) interactions on degree 1 graphs. This idea can be applied in particular to build Floquet Kitaev models, which are relevant for topological quantum computing with non-abelian anyons [62], [63].

5 Application to parallel quantum circuits and quantum error correction↩︎

As discussed in Sec. 2.3, Walsh pulse sequences can find application in realizing quantum gates in analog quantum processors. This can be done by effectively deleting interactions until qubits are interacting in pairs, realizing two-qubit gates. Let us illustrate this powerful connection between quantum circuits and Walsh sequences in the context of quantum error correcting circuits.

We consider the surface code, which is routinely realized with a set of physical qubits placed on a square lattice encoding the logical qubit, and ancilla qubits used to perform error detection [64]. The building block of the code is the measurement of stabilizer operators \(\bigotimes_{i \in \text{plaq}(X^{(j)})} X_i\) (\(\bigotimes_{i \in \text{plaq}(Z^{(j)})} Z_i\)) involving the \(4\) (or less) physical qubits connected to the ancilla qubit in the center, which is dubbed as \(X^{(j)}\) (\(Z^{(j)}\)) (by abuse of notations). The standard circuit realizing the \(X\) stabilizer measurement with CNOT gates is shown in Fig. 5a). Such a circuit naturally maps to a sequence of two-qubit rotations parametrized by \([R_{XX}^{\pi/2}]_{ij}=\exp(iX_i X_j\pi/4)\) (see App. 16 for a quick derivation). Note that these rotations can be realized using pair-wise Ising interactions. In QEC, such measurements in the \(X\) basis are complemented by measurements in another basis, which we choose here to be the \(Y\) basis (instead of the traditional \(Z\) basis). The circuit for \(Y\) measurements is shown in Fig. 5b).

Let us assemble these different circuits blocks to show the implementation of a minimal surface code made of \(4\) physical qubits, and \(3\) ancilla bits, namely the \(7-\)surface code. The circuit realizing the \(7-\)surface code in terms of Ising interactions is shown in Fig. 5c). One realizes sequentially the three (commuting) measurements circuits involving two \(X^{(i)}\) and one \(Y\) ancilla qubits. In order to obtain an efficient mapping of Ising evolution that can be implemented withing our framework in parallel, we take advantage of reordering strategies that have been proposed for instance in the context of superconducting surface codes [65]. This consists in permuting the order of commuting gates. For the \(7-\)surface code, we then obtain the circuit in Fig. 5c). Note that for larger codes, we can obtain in a similar way circuits that are also of depth \(4\). In each of these \(4\) layers shown in Fig. 5c), every qubit is subject to at most one gate, thus we can implement all these gates in parallel. In particular, within our framework, each layer corresponds to a different set of independent \(XX\) and \(YY\) Ising interactions. We first consider an implementation of the \(7-\)surface code where the qubits have the same geometry as the ones in the surface code, and are interacting via power-law decaying interactions. In this setting, each layer can be realized with a single Walsh sequence. For instance, we show the Walsh sequence required to implement the first layer in Fig. 5d). We benchmark numerically this scenario for \(\alpha = 3\) by plotting the expectation values of the stabilizer operators, which should be \(\langle O \rangle = \pm 1\), following the ancilla measurement performed after the evolution with Walsh sequences. The results, averaged over \(64\) initial Haar random states for the data qubits, are reported in Fig. 5e), where we observe errors decaying with \(\sim \tau^4\), as expected for \(p=2\) Trotter formulas, and an accuracy \(|\langle O \rangle| \sim 1 - 10^{-6}\) already for \(\tau / n \sim 10^{-2}\).

We also address the implementation of a \(7-\)surface code starting from a one-dimentional interacting qubit chain. In this geometry, certain layers might require \(Q>1\) Walsh sequences. This is because the interaction rescaling graph for each block has inhomogeneous weights. For instance, as shown in Fig. 5f), the first block has \(g^X_{2,X^{(2)}} = 1\) and \(g^Y_{1,Y} = 3^{\alpha}\). Therefore, it needs two different Walsh functions, the first one implementing \(g^X_{2,X^{(2)}} = 1\) and \(g^Y_{1,Y} = 1\), for a time \(\tau_{q = 1} = \tau\), and a second one implementing \(g^X_{2,X^{(2)}} = 0\) and \(g^Y_{1,Y} = 1\), for a time \(\tau_{q=2} = (3^{\alpha} - 1)\tau\). The total number of Walsh sequences needed to implement the surface code in this geometry is \(Q = 5\). In Fig. 5g) we report numerical results for this implementation in a long-range interacting chain with \(\alpha = 0.2\), mimicking a trapped-ion chain. We observe a very similar behaviour to the 2d case initially considered. Note finally that our protocol can be naturally applied to more general quantum error correcting codes, such as the LDPC codes [66], which are raising recent interest.

6 Application to quantum optimization algorithms↩︎

While our first two examples could be implemented with typically few \(O(1)\) Walsh sequences due to the low connectivity of \(\tilde{H}\), let us address a final situation where one requires dense coupling matrices between qubits. This is the case for example in quantum optimization algorithms, where the solution to a complex computational problem is encoded in the ground of the Hamiltonian realized by a quantum simulator, see e.g. [67]–[70] with Rydberg atoms. In particular, we consider quantum optimization schemes with digital quantum annealing (dQA) [71]–[73]. dQA is a quantum algorithm that allows to approximate the ground state of a target Hamiltonian \(\tilde{H}\) by starting from a product state which is the ground state \(\ket{\psi_0}\) of a single qubit Hamiltonian \(H_0\) such that \([\tilde{H}, H_0]\neq 0\), and evolving the system as \[\label{eq:dQA95walsh} U_{dQA} = \prod_{k=1}^K e^{-i (K - k) \tau H_0} e^{-i\tau \tilde{H} k} \;.\tag{17}\] When \(K\tau \ll 1\) this approach reduces to continuous-time quantum annealing, for which we know that, for large total annealing times \(K^2\tau \gg 1\), (with requirements related to the spectral gap of the instantaneous Hamiltonian [74]) the dynamics is adiabatic and we can prepare in good approximation the ground state of \(\tilde{H}\).

We demonstrate dQA with Walsh functions considering the MaxCut problem that consists in finding the ground-state of an Ising Hamiltonian \[\label{eq:maxcut95ham} \tilde{H} = J\sum_{i<j} G_{ij}X_i X_j\tag{18}\] where \(G_{ij}\) is the adjacency matrix of the target graph \(\mathcal{G}\). For simplicity, we consider, as a resource Hamiltonian, Eq. 5 with \(\alpha = 0\), meaning \(J_{ij} = J\) \(\forall i\neq j\). In this case, the interaction rescaling graph for the \(XX\) interaction is exactly \(g^X_{ij} = JG_{ij}/(-J) = -G_{ij}\), while for the \(YY\) interaction \(g^Y_{ij} = 0\). Therefore, to implement \(\tilde{U}(\tau)\), we first decompose \(\mathcal{G}\) in degree 1 graphs using Algorithm 7. Then, for each degree 1 graph \(\mathcal{G}_q'\) we design a Walsh sequence implementing it for the \(XX\) interaction, while we decouple all the \(YY\) interactions by choosing \(y_i\neq y_j\) \(\forall i\neq j\). To flip the sign of the resulting Hamiltonian, it is sufficient to apply the pulse \(P_{\text{set}}^q\) before and after each Walsh sequence \(P_{\text{set}}^q = \bigotimes_{i, (i, j)\in E(\mathcal{G}_q')}Y_i\), which flips the sign of one of the two \(X\) operators on each link of the graph. We choose the initial state to be \(\ket{\psi_0} = \bigotimes_i \ket{0}_i\), and therefore \(H_0 = -\sum_i Z_i\). The evolution under \(H_0\) can be implemented with global pulses.

We numerically benchmark this protocol by computing the energy difference \(\langle\tilde{H}\rangle - E_{\text{gs}}\) between the prepared state and the ground state of \(\tilde{H}\) as a function of the discretization time \(K\tau\) and the total annealing time \(K^2\tau\). The numerical results are reported for the two graphs \(\mathcal{G}_1\) and \(\mathcal{G}_3\) in Fig. 6. It is possible to see that the discretization error (due to \(K\tau > 0\)) quickly saturates after \((K\tau)^{-1} \ge 2\), and so in that regime only increasing the total annealing time \(K^2 \tau\) reduces the energy difference. In particular, the energy difference reduces exponentially with \(K^2 \tau\). To further validate our method, we plot the amplitudes of the prepared state \(\ket{\psi_{\text{sim}}}\) in the \(x\) basis, which give us the probability of a certain bitstring \(s\) as measurement outcome. For parameters giving a energy difference \(\langle \tilde{H} \rangle - E_{\text{gs}} \simeq 10^{-2}\), we observe that the ‘wrong’ solutions for the MaxCut problem have at most a probability \(|\langle s | \psi_{\text{sim}}\rangle|^2 \sim 10^{-3}\), which is two order of magnitude less likely than the correct outcomes. Note that we can straightforwardly adapt the present example to different quantum optimization circuits, such as the quantum approximate optimization algorithm (QAOA) [74], [75].

7 Implementation↩︎

Finally, we briefly comment on the experimental realizations of our protocol. While we believe our protocol can also be easily implemented with solid-state qubits, we focus here for concreteness on atomic qubits.

7.1 Rydberg atoms↩︎

Rydberg atoms trapped in optical tweezer arrays [3], [76], [77] provide a natural experimental platform for the implementation of the proposed framework to universal quantum computing and simulation, featuring long-range interactions as well as fast, local pulse operations. The first proposed implementation builds on the natural emergence of the resource Hamiltonian \(H_R\) in the interaction between two opposite-parity Rydberg states, which interact via resonant dipole-dipole interactions. These interactions have been utilized for quantum simulations of the \(XY\)-model using tweezer-trapped atoms [33], [34], [78], [79]. The required local single-qubit gates can be constructed in the Rydberg-qubit basis by a combination of two global microwave rotations \(R_X^{\pi/2}\) about the X-axis by an angle of \(\pi/2\) combined with three local phase shifts, which correspond to rotations about the Z-axis, \(R_Z^{\theta_i}\), as [80] \[\begin{align} R_{O}^{\theta} = R_Z^{\theta_1} R_X^{\pi/2} R_Z^{\theta_2} R_X^{-\pi/2} R_Z^{\theta_3}. \end{align}\] This Euler decomposition [81], [82] of an arbitrary single-qubit rotation \(R_{O}^{\theta}\) has been realized experimentally in a chain of trapped ions [83]. The required local \(\pi\)-pulses about the \(X\)-axis (\(Y\)-axis) are given by the choice \((\theta_1, \theta_2, \theta_3)=(\pi/2, \pi, -\pi/2)\) (\((\theta_1, \theta_2, \theta_3)=(-\pi, \pi, \pi)\)) respectively. Local phase rotations by controllable angles can be experimentally implemented through spatial addressing of individual atoms with acousto-optic deflectors, using off-resonant coupling to lower-lying states [84]. This induces light-shifts selectively on one of the two dipole-coupled Rydberg states, leading to the required phase rotations in the qubit basis. In experiments, typical interaction strengths can reach up to \(1\,\)MHz at distances of approximately \(10\,\mu\)m [34]. We note that our pulse schemes are robust against small variations in the interaction strengths of approximately \(2\%\) arising from the finite spread of typically \(30\,\)nm of atoms trapped in the motional ground state of typical optical tweezer or lattice potentials. To suppress finite-pulse-duration errors, see Sec. 3.2, the single-qubit manipulations must be significantly faster than the timescale set by the interaction strength in \(H_R\). Experimentally, this criterion can be fulfilled, with both light-shifts and also global microwave couplings exceeding \(10\,\)MHz in state-of-the-art experiments [34], with a straightforward scaling perspective of at least one order of magnitude faster pulses. The associated timescales with interaction and single qubit gates are significantly shorter than typical Rydberg state lifetimes on the order of \(100\,\mu\)s, which can be significantly increased e.g., using circular Rydberg states [85]–[87].

In an alternative approach, the qubits can be encoded in atomic ground or low-lying metastable states. Here, single-qubit gates can be performed using single- or two-photon coupling with coupling strengths in the MHz-regime [7]. The resource Hamiltonian can be generated using off-resonant Rydberg dressing, where the Rydberg character is optically admixed to the ground state. Rydberg-dressed interactions have been realized and benchmarked in optical lattices [88]–[91], optical tweezers [92], [93] and bulk gases [94]. The optical switchability of the interactions in this case allows for a direct realization of the required pulse-schemes with alternating global application of the resource Hamiltonian and single-qubit rotations, thereby eliminating the need to compensate for finite pulse durations. Furthermore, controlling the dressed admixture, the timescale set by the resource Hamiltonian can be tuned, such that also more complex single-qubit pulse sequences can be engineered. During the application of the single-qubit gates, the resource Hamiltonian can also be turned off, such that the effective decoherence due to finite Rydberg lifetimes is minimized. To generate the required \(XY\)-type interactions in \(H_R\), two-color dressing schemes have been proposed and experimentally implemented [93], [95]–[97].

7.2 Trapped ions↩︎

With trapped ions, the resource Hamiltonian amounts to a long-range Ising model (\(J^Y_{ij}=0\)) [2], where each ion encodes a spin using two long-lived electronic states. These interactions decay as a power-law with tunable exponent \(\alpha\). Since \(H_R\) is here laser-mediated, the time-evolution can be interspersed efficiently with single-qubit pulses [25], [28], [35]. Finally, these pulses can be programmed in a local way using experimentally-demonstrated optical addressing techniques, see e.g. Ref. [83].

7.3 Polar molecules↩︎

Arrays of ultracold polar molecules have been considered for a long time as candidates to perform quantum processing tasks [98], [99]. Molecules feature rich internal structure, long coherence times [100], [101], and long-range dipolar interactions. Recently, single molecules have been directly laser cooled [102] or associated from single atoms [103] before trapping in tweezers. Recently, such tweezer-trapped single molecules also have been entangled successfully via dipolar interactions [104], [105]. We argue that ultracold polar molecules trapped in optical tweezers or lattices can be used to realize universal quantum processing with Walsh pulse sequences. In particular, it is possible to encode a spin in the rotational states. Together with the dipolar interactions, we obtain a resource Hamiltonian similar to Eq. 4 . The spin coupling depends strongly on the exact implementation and molecular species, with values as of \(J^X_{ij} = J^Y_{ij} \sim 2\pi \times 210 \text{ Hz}\) realized recently in an optical lattice of spacing \(a\sim 540 \text{ nm}\) [106] or \(J^X_{ij} = J^Y_{ij} \sim 2\pi \times 43 \text{ Hz}\) in optical tweezers at \(a\sim1.93 \,\mu\)m [104]. In lattices, Hamiltonian engineering with global pulse sequences has been already demonstrated [106]. The use of rotational states implies the need of global microwave fields in the GHz-range to perform \(X\) pulses. The experimental toolbox also allows for applying local AC Stark shifts for local molecular control away from magic trapping conditions for the qubits. This indicates that the building blocks to implement Walsh pulse sequences in array of polar molecules are currently within reach.

8 Conclusions and outlook↩︎

This work proposes a complete and robust framework for building universal quantum processors in quantum spin systems with existing quantum technology. This builds on the use of Walsh functions to parametrize in a flexible and robust way pair-wise interactions in average Hamiltonian theory. Our work points to future interesting directions.

First, while we can effectively switch on/off interactions of the resource Hamiltonian, we cannot with the present proposal increase its connectivity. Thus, it will be interesting to combine our approach with strategies based on moving atoms [7] to achieve programmable connectivity in a robust and time-optimal fashion. An alternative strategy would be to use ancilla qubits realizing, withing our framework, multiple SWAP gates, to effectively enhance the connectivity.

Furthermore, the possibility to ‘isolate’ interacting groups of spins can find further applications beyond Hamiltonian engineering. In particular, this can be used for robust Hamiltonian learning in which interacting pairs are isolated [107], and realizations and couplings of identical quantum systems for entanglement measurements [7], [108]–[110]. Our protocol should also find extensions beyond qubits, in particular for fermionic/bosonic atoms in optical lattices [111]. In this case, a Hubbard model evolution can be used as resource Hamiltonian, and the dynamics can be enriched with local pulses (Stark-shifts or beam splitter operations) that can be described in AHT. It would be interesting in particular to see if Walsh functions can be used to build fermionic quantum processors [112] from analog systems with dipolar interactions, such as fermionic magnetic atoms [113].

9 Walsh functions↩︎

A central tool in the design of our protocol are Walsh functions. Walsh functions are an orthonormal family of binary and discrete functions of time, introduced in harmonic analysis [38], [39] and became popular in signal theory [40]. In our work, we use a specific representation of Walsh functions, which relates them to Hadamard matrices. Concretely, we define the unnormalized Hadamard matrix as \[H_1 = \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\] and, recursively, a family of generalized Hadamard matrices as \[H_q = H_1 \otimes H_{q-1} = \begin{pmatrix} H_{q-1} & H_{q-1} \\ H_{q-1} & -H_{q-1} \end{pmatrix}\] where we intend the tensor product as Kronecker product. We then define Walsh functions from the rows of a generalized Hadamard matrix \(H_q\), i.e. \[w^{(k)}_a = [H_q]_{a,k}\] with \(a = 0, 1, ..., 2^q - 1\). Since \(a < 2^{q}\), in order to choose high Walsh indices \(a\) we need to increase \(q\) accordingly. From this, we obtain that the periodicity of a Walsh function is \[n = 2^q,\;q = \lceil \log_2 (a) \rceil\] which translates to the length of pulse sequences designed with Walsh functions in Eq. 10 . The orthonormality of Walsh functions is then guaranteed by the unitarity and the symmetry of the Hadamard matrices: since \([H_q^2]_{ij} = n\delta_{ij}\), and \(H_q^T = H_q\), then \[(w_a|w_b) := \frac{1}{n}\sum_k w^{(k)}_a w^{(k)}_b = \delta_{ab} \;.\] A straightforward property used multiple times in the text can be now derived: since \(w^{(k)}_0 = 1\), we obtain \[\frac{1}{n}\sum_k w^{(k)}_a = \frac{1}{n}\sum_k w^{(k)}_0 w^{(k)}_a = (w_0|w_a) = \delta_{0a} \;.\]

10 Proof of universality for two-body Hamiltonians↩︎

In this section, we prove that any two-body spin Hamiltonian can be realized in AHT using multiple Walsh sequences. Let us assume at first that we want to prove the possibility of engineering any possible Hamiltonian of the form \[\label{eq:app95univ95Hamiltonian} \tilde{H} = \sum_{i<j} (\tilde{J}^X_{ij}X_iX_j + \tilde{J}^Y_{ij}Y_iY_j)\tag{19}\] using Eq. 4 as resource Hamiltonian and multiple Walsh sequences. We recall that a single Walsh sequence provides the average Hamiltonian \[\tilde{H}^{(q)} = \sum_{i<j} \left([G^X_q]_{ij}J^X_{ij}X_iX_j + [G^Y_q]_{ij}J^Y_{ij}Y_iY_j\right)\] where \([G^X_q]_{ij}\) (\([G^Y_q]_{ij}\)) is the adjacency matrix of the graph \(\mathcal{G}^X_q\) (\(\mathcal{G}^Y_q\)), the graphs being defined by \[\label{eq:app95graph} \mathcal{G}^X_q = \{(i,j)\;|\;x^q_i=x^q_j\},\;\mathcal{G}^Y_q = \{(i,j)\;|\;y^q_i=y^q_j\}\tag{20}\] with \(x_i^q\) and \(y_i^q\) the Walsh indices chosen in the \(q-\)th Walsh sequence. Since the two graphs \(\mathcal{G}^X_q\) and \(\mathcal{G}^Y_q\) can be parametrized independently, we can focus on the \(XX\) part of the Hamiltonian. Alternating \(Q\) Walsh sequences, each one having a different duration \(\tau_q\), realizes a unitary dynamics \(\tilde{U} (\tau) = e^{-i\tau\tilde{H}} + O(\tau^{p+1})\), where \[\tilde{J}^X_{ij} = \frac{1}{\tau}J^X_{ij}\sum_q \tau_q [G^X_q]_{ij} \;.\] We then define the interaction rescaling as \[\label{eq:app95int95rescaling} g_{ij} := \tilde{J}^X_{ij} / J^X_{ij} = \sum_q c_q [G^X_q]_{ij}\tag{21}\] where \(c_q := \tau_q/\tau\), and we assume here \(J_{ij}\neq 0\). Seeing \(g_{ij}\) as the adjacency matrix of a weighted graph \(\mathcal{G}\), we note that this definition amounts to a weighted decomposition of \(\mathcal{G}\) in subgraphs \(\mathcal{G}^X_q\) respecting Eq. 20 . Therefore, to prove that we are able to engineer any Hamiltonian of the form of Eq. 19 amounts to prove that \(\mathcal{G}\) can be any graph. The proof is done in three steps of increasing complexity:

1. if \(g_{ij}\in \{0, 1\}\), then we can take \(c_q = 1\;\forall q\) in Eq. 21 , and prove that any unweighted graph \(\mathcal{G}\) can be decomposed in subgraphs respecting Eq. 20 ;

2. if \(g_{ij}>0\) with more general values, it is possible to extend the above proof by decomposing each weighted subgraph in a set of unweighted subgraphs with coefficients \(c_q>0\);

3. if \(g_{ij}\) can also be negative, we know from Eq. 21 that multiple Walsh sequences are not enough to prove universality (\(c_q>0\) by definition), but we will prove that it is sufficient the introduction of few single-qubit pulses in between Walsh sequences.

Finally, we will show that the introduction of single-qubit pulses in between Walsh sequences allows us to engineer two-body interactions between arbitrary spin operators.

10.1 Case I: \(g_{ij}\in\{0, 1\}\)↩︎

We first consider the case in which \(g_{ij}\in \{0, 1\}\). For power-law decaying couplings described by Eq. 5 , this could represent the case of engineering a model with only nearest-neighbor couplings with coupling strength strength \(J\), such as in Sec. 4. Further couplings are admitted only if their coupling strength is \(J/r_{ij}^{\alpha}\). Another relevant case is the one in which the resource Hamiltonian has \(\alpha = 0\), and the target interaction is \(J_{ij} = G_{ij}J\), i.e. the coupling strength is homogeneous and solely defined by a graph structure, such as in Sec. 6. In all of these cases, we take \(\tau_q = \tau \;\forall q\), obtaining \[g_{ij} = \sum_q (G^X_q)_{ij} \;.\] Since \(g_{ij}\in\{0, 1\}\) can be seen the adjacency matrix of a target unweighted graph \(\mathcal{G}\), this equation amounts to a graph decomposition in terms of graph respecting Eq. 20 . In particular, graphs of degree \(1\) fulfill this condition, using \(x_i=x_j\) only for the connected pairs. As any graph can be decomposed in graphs of degree \(1\), the decomposition of \(g_{ij}\) is always possible. One can proceed in different ways:

1. We can always decompose a graph \(\mathcal{G}(V, E)\) in terms of all its links \(E\) \[g_{ij} = \sum_{(k, l)\in E} \delta_{k, l} g_{ij}\] where each degree 1 graph corresponds to a single link. This approach is however inefficient, since a densely connected graph can decomposed in up to \(O(N^2)\) subgraphs;

2. Any complete graph \(K_N\) (with even \(N\)) can be decomposed in terms of \(N/2\) open (Hamilton) paths \(P_N^q\) [51]–[53] with \(N^2/2\) assignments. In turn, any open path \(\mathcal{P}_N^q\) can be decomposed in two degree 1 subgraphs \(\mathcal{P}_N^q = \mathcal{G}_q' + \mathcal{G}_q''\) by alternating links along itself (\(\mathcal{G}_q'\) consisting of the odd links, \(\mathcal{G}_q''\) consisting of the even links). Having decomposed any complete graph in terms of degree 1 graphs, we have achieved the decomposition of an arbitrary graph as well \[g_{ij} = \sum_{q} [G_q']_{ij}g_{ij} + \sum_{q} [G_q'']_{ij}g_{ij} \;.\] Note that this decomposition results in at most \(N\) subgraphs, which is optimal for complete graphs, but not for all graphs. Additionally, its runtime \(O(N^2)\), is comparable with listing all links (see [53] for the explicit algorithm). For \(N\) odd, a similar construction holds, resulting in at most \(N + O(1)\) subgraphs;

3. Algorithm 7 achieves a decomposition of any graph \(\mathcal{G}\) (of maximum degree \(k_{\text{max}}\)) in a set of \(Q = O(k_{\text{max}})\) graphs of degree 1, i.e. \[g_{ij} = \sum_{q} [g_q^1]_{ij} \;.\] This algorithm has a runtime of order \(O(N^3 k_{\text{max}})\). While classically it takes more time than the previous approach, the number of resulting subgraphs is smaller than \(N\) for a large number of graphs (according to their sparsity). This is less expensive for the resulting quantum simulation, requiring a smaller number of Walsh sequences to engineer the same Hamiltonian. For instance, we apply this method in Sec. 6, where we consider the realization of Hamiltonians with complex interaction graphs.

10.2 Case II: \(g_{ij}>0\)↩︎

We now consider the case in which the target Hamiltonian has a different geometry than the resource Hamiltonian (i.e. it requires \(g_{ij}\in \mathbb{R}\)). For the moment, we restrict to the case in which \(g_{ij}>0\). We want to prove that any positive-weighted adjacency matrix can be written as \[\label{eq:app95weighted95deco} g_{ij} = \sum_q \frac{\tau_q }{\tau}[G^X_q]_{ij}\tag{22}\] where \((G^X_q)_{ij}\) are unweighted adjacency matrices of graphs \(\mathcal{G}^X_q\) of degree 1. As a first step, following one of the procedures given in the previous subsection, we decompose the graph \(\mathcal{G}\) in terms of inhomogeneously weighted degree 1 graphs \[\label{eq:app95unweighted95deco} g_{ij} = \sum_p [g_p']_{ij}\tag{23}\] where the weights of the links \(E_p\) in the subgraphs are the same as the ones of the links in \(E\). Then, for each degree 1 graph \(\mathcal{G}_p'\), we call the smallest weight \(c_{(p, 1)}\), and we define the decomposition \[[g_p']_{ij}^{(1)} = [g_p']_{ij} - c_{(p, 1)}[g_{(p, 1)}'']_{ij}\] where \(\mathcal{G}_{p, 1}''\) is an unweighted graph with the same topology as \(\mathcal{G}_p'\). Iterating this decomposition on the resulting graph \([g_p']_{ij}^{(1)}\) for all the values of the weights, we obtain \[\label{eq:app95deg1indeg195deco} [g_p']_{ij} = \sum_k c_{(p, k)} [g_{(p, k)}'']_{ij}\tag{24}\] with \(c_{(p, k)} > 0\). Plugging this equation in Eq. 22 , we obtain \[g_{ij} = \sum_{(p, k)} c_{(p, k)} [g_{(p, k)}'']_{ij} \;.\] Taking \((p, k) = q\), and \(c_q = \tau_q / \tau\), we have shown that it is possible to reach the decomposition of Eq. 23 , where \(\{\mathcal{G}_q^X\}\) are graphs of degree 1. We can follow this same procedure when designing pulse sequences: for each weighted graph of degree 1, we will require as many Walsh sequences as there are different weights. While this case requires more Walsh sequences than the one of the previous subsection, following the decomposition of Eq. 24 will require a total time \(\sum_k\tau_{(p, k)} = \tau \cdot \text{max}_{(i, j)\in E_p}(g_{ij})\), like in the latter case.

10.3 Case III: arbitrary \(g_{ij}\) and two-body Hamiltonians↩︎

Having demonstrated how to realize any possible \(g_{ij}>0\) using Walsh sequences, we now address the more general case in which we can also have \(g_{ij}<0\). From the previous proof, if each degree 1 subgraph we obtain in Eq. 22 can take values \([g^X_q]_{ij}\in \{-1, 0, 1\}\), then \(g_{ij}\) can take any possible real value. We now show that, to obtain \([g^X_q]_{ij}\in \{-1, 0, 1\}\) it is not necessary to add more Walsh sequences, but just a single pulse at the beginning and at the end of the \(q-\)th Walsh sequence. To begin with, we recall that the \(q-\)th Walsh sequence realizes the product formula \(\tilde{U}^{(q)}(\tau) \simeq e^{-i\tau{H}^{(q)}}\), with average Hamiltonian \[{H}^{(q)} = \sum_{i<j}( |[g^X_q]_{ij}| J_{ij} X_i X_j + |[g^Y_q]_{ij}| J_{ij} Y_i Y_j ) = \sum_{i<j}( [\hat{h}^{XX}_q]_{ij} + [\hat{h}^{YY}_q]_{ij} ) \;.\] Assuming that the graph \(\mathcal{G}^X_q + \mathcal{G}^Y_q\) is also of degree 1 (if it is not the case, it can be decomposed in two degree 1 graphs), the product formula can be written as \[\tilde{U}^{(q)}(\tau) \simeq \bigotimes_{(i, j)\in E^X,\;(k, l)\in E^Y} e^{-i\tau[\hat{h}^{XX}_q]_{ij}} e^{-i\tau[\hat{h}^{YY}_q]_{kl}} \;.\] We now divide the links in two subsets \(E^X = E^X_+ + E^X_-\), where \((i, j) \in E^X_+\) (\((i, j) \in E^X_-\)) if \([g^X_q]_{ij}>0\) (\([g^X_q]_{ij}<0\)). We define the setting pulse \(P^{\text{set}} = \bigotimes_i p^{\text{set}}_i\), where:

1. \(p^{\text{set}}_i = X\) if \((i, j)\notin E^X_-\) and \((i, j)\in E^Y_-\);

2. \(p^{\text{set}}_i = Y\) if \((i, j)\in E^X_-\) and \((i, j)\notin E^Y_-\);

3. \(p^{\text{set}}_i = Z\) if \((i, j)\in E^X_-\) and \((i, j)\in E^Y_-\);

4. \(p^{\text{set}}_i = 1\) otherwise.

Let us stress that the assumption that the underlying subgraph \(\mathcal{G}^X_q + \mathcal{G}^Y_q\) is of degree 1 is crucial, as it allows us to define the ‘set’ pulse sequence in terms of individual links. Recalling that \[Z_i e^{-iJ\tau X_i X_j} Z_i = e^{-iJ\tau (Z_i X_i Z_i) X_j} = e^{+iJ\tau X_i X_j}\] applying the setting pulse results in the effective unitary evolution \[P^{\text{set}}\tilde{U}^{(q)}(\tau)P^{\text{set}} \simeq \bigotimes_{(i, j)\in E^X_+,\;(k, l)\in E^Y_+} e^{-i\tau[\hat{h}^{XX}_q]_{ij}} e^{-i\tau[\hat{h}^{YY}_q]_{kl}} \bigotimes_{(i', j')\in E^X_-,\;(k', l')\in E^Y_-} e^{i\tau[\hat{h}^{XX}_q]_{i'j'}} e^{i\tau[\hat{h}^{YY}_q]_{k'l'}} = e^{-i\tau\tilde{H}^{(q)}}\;.\] where the new average Hamiltonian is now our target Hamiltonian \[\tilde{H}^{(q)} = \sum_{i<j}( [g^X_q]_{ij} J_{ij} X_i X_j + [g^Y_q]_{ij} J_{ij} Y_i Y_j ) \;.\] At this point, we have proven the possibility to engineer, with \(Q\) Walsh sequences and at most \(Q+1\) additional pulses (which can be absorbed into the Walsh sequences anyway) any Hamiltonian of the form Eq. 19 . This proof can be extended by considering the case in which our target Hamiltonian is a general two-body spin Hamiltonian of the form \[\tilde{H} = \sum_{m=1}^{M}\sum_{i<j} J_{ij}^{(m)} O_i^{(m, 1)} O_j^{(m, 2)}\] where \(O_i^{(m, 1)}\) and \(O_i^{(m, 2)}\) are two sets of spin operators (Pauli operators or rotations of Pauli operators). To achieve this result, we repeat the construction of this subsection, but splitting the links in \(m\) subsets \(E^X = \sum_m E_m^X\). Then, for \((i, j)\in E_m\), we can choose \(p_i^{\text{set}}\) and \(p_j^{\text{set}}\) such that \[p_i^{\text{set}} X_i (p_i^{\text{set}})^{-1} = O_i^{(m, 1)},\;p_j^{\text{set}} X_j (p_j^{\text{set}})^{-1} = O_j^{(m, 2)} \;.\] For instance, for \(O_i^{(m,1)}=Z_i\), we apply the Hadamard gate \(p_i^{\text{set}}=H\). This gives \[(i, j)\in E^X_m \rightarrow p_i^{\text{set}}p_j^{\text{set}} e^{-iJ\tau X_i X_j} (p_i^{\text{set}}p_j^{\text{set}})^{-1} = e^{-iJ\tau O_i^{(m, 1)} O_j^{(m, 2)}} \;.\] Note that for \(M = 2\), \(O^{(m, 1)} = -X\), \(O^{(m, 2)} = X\), and similarly the \(YY\) part, the generalization reduces to the previous case. We have thus proven that, with the same resources as for engineering arbitrary Hamiltonians as in Eq. 19 , it is possible to engineer arbitrary two-body spin Hamiltonians. Note that this result implies the possibility to realize universal quantum computing, for instance using as gate set arbitrary phase gates \(\exp(-i\theta X_iX_j)\) and single qubit rotations.

11 Trotter errors in Walsh sequences↩︎

We discuss Trotter errors obtained by engineering a Hamiltonian \(\tilde{H}\) using a pulse sequence designed with Walsh functions. In a general \(p = 1\) Trotter formula, to realize a Hamiltonian \(\tilde{H}\) for a time \(T\), we run a set of Hamiltonians \(\{H^{(k)}\}\) each one for a time \(\tau_k\), with \(\frac{T}{\sum_k \tau_k}\) cycles and \(\tilde{H} = \frac{\sum_k \tau_k H^{(k)}}{\sum_k \tau_k}\). The Trotter errors for such a general setting can be upper bounded, at order \(O\left(\left(\sum_k\tau_k\right) T \right)\), as [48] \[\label{eq:app95general95lloyd} \left|\left|\left(\tilde{U}\left(\sum_k\tau_k \right)\right)^{T / \sum_k\tau_k} - e^{-iT\tilde{H}}\right|\right| \leq \frac{1}{2}\frac{T}{\sum_k \tau_k} \Big{|} \Big{|} \sum_{k<l} [\tau_k H^{(k)}, \tau_l H^{(l)}] \Big{|} \Big{|} + O\left(\left(\sum_k\tau_k \right)^2 T \right)\tag{25}\] which reduces to Eq. 3 for a single Walsh sequence (\(\tau_k = \tau / n\)). We consider the spectral norm \(||O||_{\infty}:= \text{max}_{v} \frac{||Ov||}{||v||}\) (which for finite Hilbert spaces coincides with the largest singular value of the operator \(O\)), which is 1 for the identity and any Pauli operator.

At first we focus on the \(Q = 1\) case, in which the average Hamiltonian \(\tilde{H}\) is obtained from a single Walsh sequence, which reduces to Eq. 3 . be simplified through the following upper bound \[\mathcal{E}(\tau, T) := ||(\tilde{U}(\tau))^{T / \tau} - e^{-iT\tilde{H}}|| = \frac{\tau T}{2n^2} \Big{|} \Big{|} \sum_{k<l} [H^{(k)}, H^{(l)}] \Big{|} \Big{|} + O(\tau^2 T) \leq \frac{\tau T}{2n^2} \sum_{k<l} ||[H^{(k)}, H^{(l)}]|| + O(\tau^2 T) \;.\] The toggling-frame Hamiltonians obtained in applying Walsh sequences to these systems are of the form \[H^{(k)} = \sum_{i<j}( J_{ij}^X w_{x_i}^{(k)}w_{x_j}^{(k)} X_i X_j + J_{ij}^Y w_{y_i}^{(k)} w_{y_j}^{(k)} Y_i Y_j)\] where \(w_{x_i}^{(k)} = \pm 1\), \(w_{y_i}^{(k)} = \pm 1\). The commutator between two toggling-frame Hamiltonians is \[\begin{align} \nonumber [H^{(k)}, H^{(l)}] & = \sum_{i<j,\;i'<j'} ( J_{ij}^X J_{i'j'}^Y w_{x_i}^{(k)} w_{x_j}^{(k)} w_{y_{i'}}^{(l)} w_{y_{j'}}^{(l)} [X_i X_j, Y_{i'}Y_{j'}] + J_{ij}^Y J_{i'j'}^X w_{y_i}^{(k)} w_{y_j}^{(k)} w_{x_{i'}}^{(l)} w_{x_{j'}}^{(l)} [Y_i Y_j, X_{i'}X_{j'}]) = \\ & = \sum_{i<j,\;i'<j'} J_{ij}^X J_{i'j'}^Y (w_{x_i}^{(k)} w_{x_j}^{(k)} w_{y_{i'}}^{(l)} w_{y_{j'}}^{(l)} - w_{x_i}^{(l)} w_{x_j}^{(l)} w_{y_{i'}}^{(k)} w_{y_{j'}}^{(k)}) [X_i X_j, Y_{i'}Y_{j'}] \;. \end{align}\] Using the triangle inequality, we can bound the norm of these commutators in the following way \[||[H^{(k)}, H^{(l)}]|| \leq \sum_{i<j,\;i'<j'} J_{ij}^X J_{i'j'}^Y |w_{x_i}^{(k)} w_{x_j}^{(k)} w_{y_{i'}}^{(l)} w_{y_{j'}}^{(l)} - w_{x_i}^{(l)} w_{x_j}^{(l)} w_{y_{i'}}^{(k)} w^{y_{j'}}_k| \cdot ||[X_i X_j, Y_{i'}Y_{j'}]||\] and further, noting that \((w_{x_i}^{(k)} w_{x_j}^{(k)} w_{y_{i'}}^{(l)} w_{y_{j'}}^{(l)} - w_{x_i}^{(l)} w_{x_j}^{(l)} w_{y_{i'}}^{(k)} w_{y_{j'}}^{(k)}) \in \{-2, 0, 2\}\), we obtain \[\label{eq:app95walsh95bound} ||[H^{(k)}, H^{(l)}]|| \leq 2 \sum_{i<j,\;i'<j'} J_{ij}^X J_{i'j'}^Y ||[X_i X_j, Y_{i'}Y_{j'}]||\tag{26}\] independently of \(k\) and \(l\). We have now obtained a bound on the norm of a commutator which does not take into account the structure of the specific Walsh sequence we are implementing, and therefore applies to the Trotter errors inherited from any Walsh sequence.

We now quantify the bound obtained from Eq. 26 for the case of a one-dimensional long-range spin-exchange model, where \[\label{eq:app95LR95coupling} J_{ij}^X = J_{ij}^Y = \frac{J}{|i-j|^{\alpha}} := J_{|i-j|}\tag{27}\] First, we re-arrange the sum \[\sum_{i<j,\;i'<j'} J^X_{ij} J^Y_{i'j'} ||[X_i X_j, Y_{i'} Y_{j'}]|| = \frac{1}{4} \sum_{i\neq j,\;i'\neq j'} J^X_{ij} J^Y_{i'j'} ||[X_i X_j, Y_{i'} Y_{j'}]|| \;.\] Only the terms in the sum in which only one between \(i\) or \(j\) is equal to \(i'\) or \(j'\) will be non-zero, and in particular they will have spectral norm 2 (for the property \([X, Y] = 2iZ\), and \([X\otimes X,Y\otimes Y]=0\)). Since \(J^X_{ij} = J^Y_{ji} = J_{|i-j|}\), we obtain \[\label{eq:app95err1} \sum_{i<j,\;i'<j'} J^X_{ij} J^Y_{i'j'} ||[X_i X_j, Y_i Y_j]|| = \sum_{i\neq j, i\neq j'} J^X_{ij} J^Y_{ij'} - \sum_{i \neq j} J^X_{ij} J^Y_{ij} \leq \sum_{i\neq j, i\neq j'} J_{|i-j|} J_{|i-j'|} \;.\tag{28}\] We can rewrite the obtained sum, and upper-bound it again \[\label{eq:app95err2} \sum_{i\neq j, i\neq j'} J_{|i-j|} J_{|i-j'|} = 4 \sum_{i<j, i<j'} J_{|i-j|} J_{|i-j'|} = 4\sum_{r_1=1}^{N-1} \sum_{r_2=1}^{N-1} (N- \text{max}(r_1, r_2))J_{r_1}J_{r_2} \leq 4\sum_{r_1=1}^{N-1} (N- r_1)J_{r_1} \sum_{r_2=1}^{N-1} J_{r_2}\tag{29}\] where we have used that \(r_1\) and \(r_2\) run from \(1\) to \(N - 1\), but at the same time \(i \leq N - r_1\) and \(i \leq N - r_2\). At this point, to simplify the computation of the sums, we use the general property of monotonously decreasing functions (coming from the definition of the Riemann integral) [114] \[\label{eq:app95integral} \sum_{k = a}^b f(k) \leq f(a) + \int_a^b \text{d}x f(x) \;.\tag{30}\] The second factor in Eq. 29 admits the following integral upper bound \[\sum_{r_2=1}^{N-1} J_{r_2} \leq J \left(1 + \int_1^{N-1} \text{d}r\;r^{-\alpha} \right) = J\frac{1}{1 - \alpha}((N-1)^{1-\alpha} - \alpha)\] and similarly does the first factor \[\label{eq:app95integral2} \sum_{r_1=1}^{N-1} (N - r_1)J_{r_1} \leq J \left( \frac{1}{(1 - \alpha)(2 - \alpha)}(N - 1)^{2 - \alpha} + \frac{1}{1 - \alpha}(N - 1)^{1 - \alpha} - \frac{\alpha}{1 - \alpha}N - \frac{(1-\alpha)}{2 - \alpha}\right) \;.\tag{31}\] Multiplying these two results yields then an upper bound to Eq. 28 . If one takes the large \(N\) limit of the product it is possible to note that, depending on the value of \(\alpha\), two possible leading terms arise. For \(\alpha>1\), one obtains \[\label{eq:app95trotter95bound95g1} \sum_{i\neq j, i\neq j'} J_{|i-j|} J_{|i-j'|} = 4J^2\left(\frac{\alpha}{\alpha - 1}\right)^2 N + o(N) =: 2J^2a_{\alpha}N\tag{32}\] while for \(\alpha < 1\) \[\label{eq:app95trotter95bound95l1} \sum_{i\neq j, i\neq j'} J_{|i-j|} J_{|i-j'|} = 4J^2\frac{1}{(1 - \alpha)^2 (2 - \alpha)}N^{3 - 2\alpha} + o(N^{3 - 2\alpha}) =: 2J^2b_{\alpha}N^{3-2\alpha} \;.\tag{33}\] Therefore, any quantum simulation over a total time \(T\) with \(p = 1\) Walsh sequences with period \(\tau\) will have a Trotter error bounded by \[\label{eq:app95final95trotter} \begin{cases} \mathcal{E} \leq \frac{\tau T}{2 n^2} \frac{n(n-1)}{2}\cdot 4J^2 a_{\alpha} N \leq a_{\alpha} N(JT)(J\tau) \;\text{if} \;\alpha > 1 \\ \mathcal{E} \leq \frac{\tau T}{2 n^2} \frac{n(n-1)}{2}\cdot 4J^2 b_{\alpha} N^{3 - 2\alpha} \leq b_{\alpha} N^{3 - 2\alpha}(JT)(J\tau) \;\text{if} \;\alpha < 1 \end{cases}\tag{34}\] which implies the bound of Eq. 13 , in the case \(Q = 1\).

This result can be straightforwardly extended to cycles comprised of \(Q\) Walsh sequences. For each Walsh sequence, we have a different period \(\tau_q\) and number of steps \(n_q\). Therefore, Eq. 25 becomes \[\label{eq:app95multi95lloyd} \mathcal{E}(\{\tau_q\}, T) = \frac{T}{2\sum_q \tau_q} \sum_{(q, k) < (q', l)} \frac{\tau_q \tau_{q'}}{n_q n_{q'}}||[H^{(q, k)},H^{(q', l)}]||\tag{35}\] where \(H^{(q, k)}\) is the \(k-\)th toggling frame Hamiltonian of the \(q-\)th Walsh sequence, and \((q, k)<(q', l)\) means \(k<l\) if \(q=q'\), or \(q<q'\) for arbitrary \(k, l\). The bound Eq. 26 still holds, giving, for the previously considered example Eq. 27 \[||[H^{(q, k)},H^{(q', l)}]|| \leq 2\sum_{i\neq j, i'\neq j'} J^X_{ij} J^Y_{ij} ||[X_i X_j, Y_{i'}Y_{j'}]|| \leq 2\sum_{i\neq j, i\neq j'} J_{|i-j|} J_{|i-j'|} \;.\] Plugging this bound in Eq. 35 we obtain \[\mathcal{E}(\{\tau_q\}, T) = \frac{T}{\sum_q \tau_q} \left( \sum_q \sum_{k<l} \frac{\tau_q^2}{n_q^2} + \sum_{q<q'}\sum_{k, l} \frac{\tau_q \tau_{q'}}{n_q n_{q'}} \right)\sum_{i\neq j, i\neq j'} J_{|i-j|} J_{|i-j'|} \le \frac{\left( \sum_q \tau_q \right) T}{2} \sum_{i\neq j, i\neq j'} J_{|i-j|} J_{|i-j'|} \;.\] Now, plugging in the formula above the bounds of Eq. 32 , 33 , we obtain expressions corresponding to Eq. 34 with the substitution \(\tau \rightarrow \sum_q \tau_q\), i.e. \[\label{eq:app95final95trotter95multi} \begin{cases} \mathcal{E}(\{\tau_q\}, T) \leq a_{\alpha} N(JT)(J\sum_q\tau_q) \;\text{if} \;\alpha > 1 \\ \mathcal{E}(\{\tau_q\}, T) \leq b_{\alpha} N^{3 - 2\alpha}(JT)(J\sum_q \tau_q) \;\text{if} \;\alpha < 1 \end{cases}\tag{36}\] which we report in Eq. 13 in the main text.

12 Approximation scheme for shorter Walsh sequences↩︎

In this section we investigate the unitary error due to imposing a cut-off \(\Lambda_w\) to the Walsh indices. We will focus on one-dimensional systems, even though we believe this approach can be extended to higher dimensional systems. We consider the case in which the target Hamiltonian \(\tilde{H}\) has finite range interactions, i.e. there exists a distance \(r = O(1)\) such that \(\tilde{J}_{ij}=0\) if \(|i-j|>r\). In general, the Walsh sequence realizing the dynamics of \(\tilde{H}\) in AHT will have a length \(n = O(N)\). However, the length can be reduced by assuming that interactions in the resource Hamiltonian \(H_R\) are negligible after a certain distance \(\Lambda_r\geq r\) (for resource Hamiltonians with couplings Eq. 5 it is the case if \(\alpha\) is large). In this case, we do not require using different Walsh indices to decouple two qubits that are separated by distance larger than \(\Lambda_r\). Formally, this means that Eq. 9 becomes \[\label{eq:app95graph95cutoff} E(\mathcal{G}_1^X) = \{(i,j)\;|\;x_i=x_j\;\text{and}\;|i-j|<\Lambda_r \},\;E(\mathcal{G}_1^Y) = \{(i,j)\;|\;y_i=y_j\;\text{and}\;|i-j|<\Lambda_r \} \;.\tag{37}\] Fixing the distance cut-off \(\Lambda_r\) provides us with a strategy to reduce the pulse sequence length (see Fig. 8):

1. order the values of Walsh indices \(x_i\);

2. divide the \(N\) qubits in the subregions \(A_l\), such that \([l\Lambda_r + 1, (l+1)\Lambda_r]\in A_l\), in such a way that links can only be present between \(A_l\) and \(A_{l+1}\);

3. assign sequentially (the first \(x_i\) to the 1st qubit) the Walsh indices to \(A_1\) and \(A_2\);

4. assign sequentially the Walsh indices to \(A_3\) starting from the beginning (i.e. from those in \(A_1\)), but excluding those already assigned in \(A_2\)

5. repeat, where in any \(A_l\) are excluded the Walsh indices in \(A_{l-1}\).

This protocol ensures that the shortest residual interaction acts on a distance \(\Lambda_r + 1\). Note that, if the region \(A_l\) is not coupled with \(A_{l-1}\), we can assign the Walsh indices there without constraints. The highest Walsh index assigned in this way will be \(\text{max}_i(x_i) \leq 2\Lambda_r\). This can be proven from the fact that, given \(X_l := \text{max}_{i \in A_l}(x_i)\), we have \(X_{l} \leq \Lambda_r + X_{l-1} - X_{l-2}\) (together with \(X_1\leq \Lambda_r\)), as we can use the Walsh indices in \(A_{l-2}\), but not those in \(A_{l-1}\). At the same time we can apply the condition, \(X_{l-1} \leq \Lambda_r + X_{l-2}\), which further bounds the inequality above. The two conditions together give \(X_{l} \leq \Lambda_r + X_{l-1} - X_{l-2} \leq 2\Lambda_r\), as claimed.

Let us now address the effect of the residual interactions \(|i-j|\ge \Lambda_r\), which we neglected in the first paragraph. The Walsh sequence \((x_i, y_i)\) will not realize exactly (in AHT) the dynamics of the target Hamiltonian \(\tilde{H}\), but instead, in a single Trotter step \[\tilde{U}(\tau) = e^{-i\tau(\tilde{H} + H_E)} + O(\tau^2)\] where \[H_E = \sum_{i, j, |i-j|>\Lambda_r}(J^X_{ij}\delta_{x_i,x_j}X_i X_j + J^Y_{ij}\delta_{y_i,y_j} Y_i Y_j) \;.\] which consists in an additional contribution to the error \(||(\tilde{U}(\tau))^{T/\tau} - e^{-iT\tilde{H}}||\). We now provide an estimate for this contribution at the level of a single Trotter step. First, we note that working in AHT allows us to regroup toggling-frame Hamiltonian arbitrarily, so \[\tilde{U}(\tau) = e^{-i\tau\tilde{H} + H_E} + O(\tau^2) = e^{-i\tau\tilde{H}}e^{-i\tau H_E} + O(\tau^2) = e^{-i\tau\tilde{H}}\left(1 + \sum_{l=1}^{\infty} \tau^l H_E^l \right) + O(\tau^2) \;.\] The norm of the error induced by \(H_E\), beside Trotter errors, is then \[||e^{-i\tau\tilde{H}}\left(1 + \sum_{l=1}^{\infty} \tau^l H_E^l \right) - e^{-i\tau\tilde{H}}|| = ||e^{-i\tau\tilde{H}} \sum_{l=1}^{\infty} \tau^l H_E^l || \leq \sum_{l=1}^{\infty} \tau^l ||H_E||^l \;.\] To ensure that the error Hamiltonian \(H_E\) does not hinder the simulation of the target Hamiltonian \(\tilde{H}\), we would like the errors induced by \(H_E\) to be comparable or smaller than the Trotter errors. Therefore, we want them to be of order \(O(\tau^2)\). This can be obtained by fixing the norm of \(H_E\) to be \[\label{eq:app95precision95cutoff} ||H_E|| = \epsilon \tau\tag{38}\] with \(\epsilon = O(1)\) to be a parameter tunable to make the error smaller or comparable to Trotter errors. This translates in conditions on the cut-off distance \(\Lambda_r\) and the exponent \(\alpha\) of the interactions decay. To obtain them, we evaluate a bound to \(||H_E||\) for the one-dimensional case, of which we consider an example in Sec. 2.5. From the triangle inequality \[||H_E|| \leq \sum_{i=1}^{N-1} \sum_{r=\Lambda_r + 1}^{N-i} J_r ||X_i X_{i+r} + Y_i Y_{i+r}|| = \sum_{r=\Lambda_r + 1}^{N-1}(N - r)J_r \;.\] We now estimate the sum on the RHS using the integral bound Eq. 30 , from which we obtain \[\sum_{r=\Lambda_r + 1}^{N-1}(N - r)J_r \leq \frac{J}{(1-\alpha)(2-\alpha)}(N-\alpha+1)(N-1)^{1-\alpha} - \frac{J}{1-\alpha}N(\Lambda_r +\alpha)(\Lambda_r +1)^{-\alpha} - \frac{J}{2-\alpha}(\Lambda_r + \alpha - 1)(\Lambda_r+1)^{1-\alpha} \;.\] By taking the large \(N\) limit, we notice that \(||H_E||\) has a different leading term depending on the value of \(\alpha\). For \(\alpha < 1\) we have \[||H_E|| \leq \frac{J}{(1-\alpha)(2-\alpha)}N^{2-\alpha} + O(N)\] which is a bound that does not depend on the cut-off \(\Lambda_r\) and thus cannot guarantee the criterion Eq. 38 . We deduce that, in this regime, imposing any cut-off to Walsh indices would lead to a dramatic error in the Hamiltonian engineering. This can be also understood intuitively, by noticing that, for \(\alpha<1\), even far away qubits are strongly interacting, therefore coupling them would lead to a much different dynamics than having them decoupled. For \(\alpha>1\) we have instead \[||H_E|| \leq \frac{J}{(\alpha-1)}N(\Lambda_r +\alpha)(\Lambda_r +1)^{-\alpha} + O(1) \;.\] If we consider \(\tau = O(N^{-1})\) (to have performance guarantees with respect to Trotter errors, see Eq. 34 ), then for \(\alpha>3\) this bound guarantees the existence of a value of \(\Lambda_w = o(N)\) satisfying Eq. 38 . In particular, \(||H_E|| = \epsilon / N\) if \(\Lambda_w = O((N/\epsilon)^{\frac{2}{\alpha - 1}})\). Note that Trotter errors, from our numerical simulations, can have a weaker scaling with respect to \(N\) for interesting quantities than what can be estimated by the unitary error, leading this approximation to work also for smaller \(\alpha\).

Therefore, we conclude that, for \(\alpha\) large enough, it is possible to set a cut-off \(\Lambda_w = o(N)\) to reduce the length of Walsh sequences without introducing further errors. In particular, for the one-dimensional case, we have provided analytical guarantee that this is the case for \(\alpha>3\). In Sec. 2.5 we report a numerical example indicating that a finite \(\Lambda_w\) will create errors comparable or smaller than Trotter errors even for \(\alpha = 3\), but not for \(\alpha < 1\).

13 Rotation angle errors↩︎

In this section, we derive the effect of rotation angle errors in the implementation of Walsh sequences, and derive conditions under which Walsh sequences are robust against this kind of errors. In particular, we study their effect in AHT. If the rotation angle errors are of order \(O(\delta)\), we find that there is an error in AHT of order \(O(\delta)\). The robustness condition we find allows to reduce this AHT error to \(O(\delta^2)\).

13.1 Effect on the average Hamiltonian↩︎

We model rotation angle errors by assuming that each \(s_i\pi-\)pulse in a certain direction (where \(s_i = \pm 1\)) becomes a \(s_i(\pi + \delta_i)-\)pulse instead, \(\delta_i\) being a single-qubit error, i.e. \[\begin{align} P^{(k)} = \bigotimes_i e^{-is_i(\pi + \delta_i)[O]_i/2} \end{align}\] where \([O]_i\) is the operator with respect to which we want to pulse on a given qubit. This mimics the case in which the pulsing time is not well calibrated on the qubits. The exact Walsh sequences are recovered for \(\delta_i = 0\), and any choice of \(s_i\). The resource Hamiltonian considered is \[H_R = \sum_{i<j} (J_{ij}^X X_i X_j + J_{ij}^Y Y_i Y_j) \;.\] To study the effect of the errors, we consider the first pulse of the interval \(k\) to be \(P^{(k)}\) as defined above, and the second one to be \((P^{(k)})^{-1}\). The effect on local \(X_i\) and \(Y_i\) operators is \[\begin{align} X_i^{(k)} := (p_i^{(k)})^{-1} X_i p_i^{(k)} &= \frac{1 + w_{x_i}^{(k)}}{2}X_i + \frac{1 - w_{x_i}^{(k)}}{2}(\cos(\pi + \delta_i)X_i + s_i\sin(\pi + \delta_i)\tilde{X}_i^{(k)}) = \\ &= w_{x_i}^{(k)} X_i - \delta_i\frac{1 - w_{x_i}^{(k)}}{2}(s_i \tilde{X}_i^{(k)} - \frac{1}{2}\delta_i X_i) + O(\delta_i^3) \\ Y_i^{(k)} := (p_i^{(k)})^{-1} Y_i p_i^{(k)} &= \frac{1 + w_{y_i}^{(k)}}{2}Y_i + \frac{1 - w_{y_i}^{(k)}}{2}(\cos(\pi + \delta_i)Y_i + s_i\sin(\pi + \delta_i)\tilde{Y}_i^{(k)}) = \\ &= w_{y_i}^{(k)} Y_i - \delta_i\frac{1 - w_{y_i}^{(k)}}{2}(s_i \tilde{Y}_i^{(k)} - \frac{1}{2}\delta_i Y_i) + O(\delta_i^3) \end{align}\] where \[\begin{align} \tilde{X}_i^{(k)} = \left(\frac{1 - w_{y_i}^{(k)}}{2} Y_i - \frac{1 + w_{y_i}^{(k)}}{2}Z_i \right),\;\tilde{Y}_i^{(k)} = \left(\frac{1 + w_{x_i}^{(k)}}{2} Z_i - \frac{1 - w_{x_i}^{(k)}}{2}X_i \right) \;. \end{align}\] The resulting toggling-frame Hamiltonian reads \[\begin{align} \label{eq:app95RA95tfH} H^{(k)} (\pmb{\delta}) = \sum_{i<j}\left(J^X_{ij}w^{(k)}_{x_i}w^{(k)}_{x_j}X_i X_j + J^Y_{ij}w^{(k)}_{y_i}w^{(k)}_{y_j} Y_i Y_j \right) + H^{(k)}_{\text{RA}} = H^{(k)}(\pmb{\delta} = 0) + H^{(k)}_{\text{RA}} \end{align}\tag{39}\] where \[\begin{align} H^{(k)}_{\text{RA}} = \sum_{O = X, Y}\sum_{i<j} J_{ij}^O \bigg( \left(s_i\delta_i \frac{1 - w^{(k)}_{o_i}}{2} w^{(k)}_{o_j} \tilde{O}_i^{(k)} O_j + s_j \delta_j w^{(k)}_{o_i}\frac{1 - w^{(k)}_{o_j}}{2} O_i \tilde{O}_j^{(k)} \right) + \\ + s_i s_j \delta_i \delta_j \frac{1 - w^{(k)}_{o_i}}{2}\frac{1 - w^{(k)}_{o_j}}{2} \tilde{O}_i^{(k)}\tilde{O}_j^{(k)} +\frac{1}{2}\left(\delta_i^2 \frac{1 - w^{(k)}_{o_i}}{2} w^{(k)}_{o_j} + \delta_j^2 w^{(k)}_{o_i}\frac{1 - w^{(k)}_{o_j}}{2} \right)O_i O_j \bigg) + O(\delta^3) \;. \end{align}\] This results tells us that the rotation angle errors will ultimately lead to errors in the average Hamiltonian of order \(O(\delta)\). In the next subsection we derive a condition to remove the \(O(\delta)\) error from the average Hamiltonian.

13.2 Robustness conditions↩︎

We now show how it is possible to use the sign functions \(s_i\) to remove the order \(O(\delta)\) in the rotation angle error term of the average Hamiltonian. We start by labelling a single Walsh sequence (or a cycle of multiple Walsh sequences) with an index \(l\). This means that the cycles \(\tilde{U}(\tau)\) (which can be made by multiple Walsh sequence) will be now labelled as \(\tilde{U}^{(l)}(\tau)\), such that \[(\tilde{U}(\tau))^{T/\tau} \rightarrow \prod_{l=1}^{T/\tau}\tilde{U}^{(l)}(\tau) \;.\] For later convenience, we choose the index \(l\) to have a periodicity \(L\), i.e. \[\tilde{U}^{(l+L)}(\tau) = \tilde{U}^{(l)}(\tau) \;.\] This is also known as second averaging [27]. We now consider the sign functions \(s_i\) to be dependent on this new index \(l\), i.e. \(s_i \rightarrow s_i^l\). Note that this does not change the results for Walsh sequences without pulse errors. Let us further choose \(s_i^l = w^{(l)}_{e_i}\). We then introduce the second index in the toggling-frame Hamiltonians \(H^{(k, l)}\), and the doubly averaged Hamiltonian \[\tilde{H}'_{\text{RA}} = \frac{1}{L}\sum_{l=1}^L\tilde{H}^{(l)} = \frac{1}{nL}\sum_{l=1}^L \sum_{k=1}^n H^{(k, l)}\] where \(L\) is the period of the longest Walsh function \(w_{e_i}\). Assuming \(x_i, y_i \neq 0\), performing the double average over the toggling-frame Hamiltonians in Eq. 39 , we obtain \[\begin{align} \tilde{H}'_{\text{RA}} = \sum_{O = X, Y}\sum_{i<j} J_{ij}^O \Bigg(\delta_{o_i, o_j}\left( 1 - \frac{\delta_i^2 + \delta_j^2}{4} \right)O_i O_j + \left(\delta_{0, e_i}\delta_i \mathcal{O}_i^j O_j + \delta_{0, e_j} \delta_j O_i \mathcal{O}_j^i \right) + \delta_{e_i, e_j} \delta_i \delta_j \mathcal{C}_{ij} \Bigg) + O(\delta^3) \end{align}\] where the sum over \(O = X, Y\) runs also over \(o = x, y\) for the Walsh indices, and \[\begin{align} \mathcal{O}_i^j =& \frac{1}{n}\sum_{k=1}^n \frac{1 - w^{(k)}_{o_i}}{2}w^{(k)}_{o_j}\tilde{O}^{(k)}_i \\ \mathcal{C}_{ij} =& \frac{1}{n}\sum_{k=1}^n \frac{1 - w^{(k)}_{o_i}}{2}\frac{1 - w^{(k)}_{o_j}}{2}\tilde{O}^{(k)}_i\tilde{O}^{(k)}_j \;. \end{align}\] From the equations above, it can be noticed that choosing \(e_i \neq 0\) \(\forall i\) removes exactly the undesired term of \(O(\delta)\). We will refer to this condition as the robustness condition against rotation angle errors. The easiest choice to implement the robustness condition is \(e_i = 1\) \(\forall i\), which amounts to performing alternatively \(\pi-\) and \(-\pi-\) pulses in consecutive Walsh sequences. Another choice is instead \(e_i = i\), which furthermore removes the term \(\sim \delta_i \delta_j\) from the doubly averaged Hamiltonian. Upon choosing \(e_i = i\), the doubly averaged Hamiltonian reads \[\begin{align} \tilde{H}'_{\text{RA}} = \sum_{i<j} \left( 1 - \frac{\delta_i^2 + \delta_j^2}{4} \right)\left( J_{ij}^X\delta_{x_i, x_j}X_i X_j + J_{ij}^Y\delta_{y_i, y_j}Y_i Y_j \right) + O(\delta^3) \end{align}\] which amounts to the target Hamiltonian, plus a \(O(\delta^2)\) inhomogeneous renormalization of the interactions.

A numerical example of the implementation of this robustness condition is shown in Fig. 9. In particular, we consider the cluster state preparation protocol via Ising dynamics (presented in Sec. 4) for \(N=6\) qubits, and \(\alpha = 1.2\). We consider a uniform sampling of each \(\delta_i\) from the interval \([-\epsilon_{\text{RA}}J, \epsilon_{\text{RA}}J]\), for different values of \(\epsilon_{RA}\). We quantify the errors by plotting the fidelity error with respect to the cluster state, averaged over \(64\) \(\delta_i\) samples, both for the case in which the robustness condition is implemented (blue lines) and not implemented (red lines), and compare it with the pure Trotter errors. It can be observed that implementing the robustness condition reduces the fidelity error by orders of magnitude. While the fidelity error for robust Walsh sequences still saturates for small values of \(\tau\) (due to the residual \(O(\delta^2)\) term), this saturation happens for much smaller values of \(\tau\), leading to a relevant improvement of the protocol.

14 Finite pulse duration errors↩︎

In this section we consider the effect of pulses of finite duration with respect to the ideal case, in which the pulses are instantaneous, and their dynamics is decoupled from the Hamiltonian dynamics given by interactions. We first derive the effect on the average Hamiltonian, which are of order \(O(t_p/\tau)\), and then we show a robustness condition which completely removes the errors in AHT, upon a slight modification of the simulation procedure.

14.1 Average Hamiltonian theory for finite pulse duration↩︎

Let consider an interval \(k\) of a Walsh sequence, where we want to apply a pulse \(P^{(k)}\) at the beginning, and another pulse \((P^{(k)})^{-1}\) at the end. We assume that the pulses have a square-wave form in time. In order to apply a \(s_i\pi-\)pulse in the \([O]_i\) direction, the physical implementation is done by means of the Hamiltonian \[H_p = \frac{\pi}{2t_p}\sum_i s_i [O]_i\] where \[P^{(k)} = e^{-it_p H_p},\;(P^{(k)})^{-1} = e^{it_p H_p}\] and \(t_p\) is the time duration of the pulse. In general, during the application of a pulse, the resource Hamiltonian will still make the qubits of the underlying system interact. Therefore, the evolution during an interval \(k\) of a Walsh sequence will be \[U_{\frac{\tau}{n}}^{(k)} = e^{-it_p(H_R - H_p)}e^{-i(\tau/n - 2t_p)H_R}e^{-it_p(H_R + H_p)} \;.\] In the case of instantaneous pulses (\(t_p=0\)), considered throughout the text, the evolution above becomes \[U_{\frac{\tau}{n}}^{(k)} = (P^{(k)})^{-1}e^{-i(\tau/n)H_R}P^{(k)} = e^{-i(\tau/n)H^{(k)}}\] i.e. the pulses and the interacting evolution are decoupled, as expected. To analyze the relation between the instantaneous pulses case and the finite pulses case, we put ourselves in the interaction picture with respect to the pulsing Hamiltonians \(H_p\), obtaining [55] \[\begin{align} e^{-it_p(H_R \pm H_p)} = e^{\mp it_p H_p}\mathcal{U}_{\pm}(t_p) &= (P^{(k)})^{\pm 1}\mathcal{U}_{\pm}(t_p) \\ \frac{d\mathcal{U}_{\pm}(t)}{dt} = -i(e^{\pm it H_{p}} H_R e^{ \mp it H_{p}})\mathcal{U}_{\pm}(t) &= -i((Q^{(k)}(t))^{\mp 1} H_R (Q^{(k)}(t))^{\pm 1}\mathcal{U}_{\pm}(t) \end{align}\] where we have defined the partially applied pulse \[\begin{align} Q^{(k)}(t) = e^{-itH_p} = \bigotimes_i q^{(k)}_i(t) \;. \end{align}\] Hence, we have \[\begin{align} U^{(k)}_{\frac{\tau}{n}} = (P^{(k)})^{-1}\mathcal{U}_-(t_p) e^{-i(\tau/n - 2t_p))H_R} P^{(k)} \mathcal{U}_+(t_p) \;. \end{align}\] Defining the toggling-frame Hamiltonian \(H^{(k)} = (P^{(k)})^{-1}H_R P^{(k)}\), we have \[U^{(k)}_{\frac{\tau}{n}} = \mathcal{U}'_-(t_p) e^{-i(\tau/n - 2t_p))H^{(k)}}\mathcal{U}_+(t_p)\] where \[\frac{d\mathcal{U}'_-(t)}{dt} = (P^{(k)})^{-1} \frac{d\mathcal{U}_-(t)}{dt} P^{(k)}= -i(Q^{(k)}(t) H^{(k)} (Q^{(k)}(t))^{-1})\mathcal{U}_-'(t) \;.\] We can then define the instantaneous toggling-frame Hamiltonian for the interval \(k\) as follows: \[\begin{align} H^{(k)}(t) = (Q^{(k)}(t))^{-1} H_R Q^{(k)}(t)&,\;t<t_p \\ H^{(k)}(t) = (P^{(k, 1)})^{-1} H_R P^{(k, 1)} = H^{(k)}&,\;t_p<t<\frac{\tau}{n} - t_p \\ H^{(k)}(t) = Q^{(k)}(t) H^{(k)} (Q^{(k)}(t))^{-1} &,\;\frac{\tau}{n} - t_p < t < \frac{\tau}{n} \;. \end{align}\] The average Hamiltonian over the interval \(k\) is then [45] \[H^{(k)}_{\text{av}} = \frac{n}{\tau}\int_0^{\tau/n} dt\;H^{(k)}(t)\] and over the sequence is \[\begin{align} \tilde{H} = \frac{1}{n} \sum_{k=1}^n H^{(k)}_{\text{av}} = \frac{1}{n} \sum_{k=1}^n \frac{n}{\tau} \int_0^{\tau/n} dt\;H^{(k)}(t) \;. \end{align}\] It is already possible to note that the average Hamiltonian over the interval \(k\) reads \[\begin{align} H^{(k)}_{\text{av}} =& \left(1 - \frac{2nt_p}{\tau} \right)H^{(k)} + \frac{n}{\tau} \int_0^{t_p}dt\;\left((Q^{(k)}(t))^{-1} H_R Q^{(k)}(t) + Q^{(k)}(t) H^{(k)} (Q^{(k)}(t))^{-1} \right) = \\ =& \left(1 - \frac{2nt_p}{\tau} \right)H^{(k)} + \frac{2n}{\tau}\int_0^{t_p}dt\;(Q^{(k)}(t))^{-1} H_R Q^{(k)}(t) \end{align}\] from which we can define the dimensionless error parameter \[\epsilon_{\text{FP}} := \frac{2nt_p}{\tau}\] which amounts to the fraction of time of the interval \(k\) during which pulses are applied. Finally, we define the average Hamiltonian over the whole sequence for finite pulses \[\label{eq:app95avg95H95FP95general} \tilde{H}_{\text{FP}} = \left(1 - \epsilon_{\text{FP}} \right)\tilde{H} + \epsilon_{\text{FP}}\frac{1}{n}\sum_{k=1}^n\frac{1}{t_p}\int_0^{t_p}dt\;(Q^{(k)}(t))^{-1} H_R Q^{(k)}(t)\tag{40}\] where \(\tilde{H}\) is the target Hamiltonian. We conclude then that we only need to compute the second term on the right-hand side to estimate the finite pulse duration effects.

14.2 Average Hamiltonian with finite pulse duration effects↩︎

We now apply the result above to the Walsh sequence protocol considered in the text. We consider the long-range \(XY\) model as a resource Hamiltonian \[\begin{align} H_R = \sum_{i<j} (J_{ij}^X X_i X_j + J_{ij}^Y Y_i Y_j) \;. \end{align}\] In order to compute the instantaneous toggling-frame Hamiltonian in an interval \(k\), we need to keep track of the evolution of the instantaneous toggling-frame \(X\) and \(Y\) operators, which evolve during the first pulse (\(0<t<t_p\)) as \[X^{(k)}_i (t) :=(q_i^{(k)}(t))^{-1} X_i q_i^{(k)}(t) = \frac{1 + w^{(k)}_{x_i}}{2} X_i + \frac{1 - w^{(k)}_{x_i}}{2}\left( X_i \cos(\pi (t / t_p ) ) + s_i\tilde{X}^{(k)}_i \sin(\pi (t / t_p )) \right)\] and \[\begin{align} Y^{(k)}_i (t) :=(q_i^{(k)}(t))^{-1} Y_i q_i^{(k)}(t) = \frac{1 + w^{(k)}_{y_i}}{2} Y_i + \frac{1 - w^{(k)}_{y_i}}{2}\left( Y_i \cos(\pi (t / t_p ) ) + s_i\tilde{Y}^{(k)}_i \sin(\pi (t / t_p )) \right) \end{align}\] where \[\begin{align} \tilde{X}_i^{(k)} = \left(\frac{1 - w^{(k)}_{y_i}}{2} Y_i - \frac{1 + w^{(k)}_{y_i}}{2}Z_i \right),\;\tilde{Y}_i^{(k)} = \left(\frac{1 + w^{(k)}_{x_i}}{2} Z_i - \frac{1 - w^{(k)}_{x_i}}{2}X_i \right) \;. \end{align}\] We can now compute the integral in Eq. 40 as \[\begin{align} \frac{1}{t_p}\int_0^{t_p} dt\;(Q^{(k)}(t))^{-1} H_R Q^{(k)}(t) = \frac{1}{\pi}\int_0^{\pi} d\theta\;H^{k}(\theta) = \frac{1}{\pi} \sum_{O = X, Y} \sum_{i<j} J_{ij}^O \int_0^{\pi} d\theta\;O_i^{(k)}(\theta)O_j^{(k)}(\theta) \end{align}\] where \[\begin{align} O^{(k)}_i (\theta) = \frac{1 + w^{(k)}_{o_i}}{2} O_i + \frac{1 - w^{(k)}_{o_i}}{2}\left( O_i \cos\theta + s_i\tilde{O}^{(k)}_i \sin\theta \right) \;. \end{align}\] Making use of the following trigonometric integrals \[\begin{align} \int_0^{\pi} d\theta = \pi,\; \int_0^{\pi} d\theta \cos\theta = 0,\; \int_0^{\pi} d\theta \sin\theta = 2, \\ \int_0^{\pi} d\theta \sin\theta \cos\theta = 0,\; \int_0^{\pi} d\theta \cos^2\theta = \frac{\pi}{2},\; \int_0^{\pi} d\theta \sin^2\theta = \frac{\pi}{2} \;. \end{align}\] we obtain \[\begin{align} \frac{1}{\pi}\int_0^{\pi} d\theta\;O^{(k)}_i(\theta) O^{(k)}_j(\theta) = \frac{1}{8}\left( 3 - w^{(k)}_{o_i} - w^{(k)}_{o_j} + 3 w^{(k)}_{o_i}w^{(k)}_{o_j} \right) O_i O_j + \frac{1}{2}\frac{1 - w^{(k)}_{o_i}}{2} \frac{1 - w^{(k)}_{o_j}}{2} s_i s_j \tilde{O}_i^{(k)} \tilde{O}_j^{(k)} + \\ + \frac{2}{\pi}\left(s_i \frac{1 - w^{(k)}_{o_i}}{2}\frac{1 + w^{(k)}_{o_j}}{2}\tilde{O}_i^{(k)} O_j + s_j \frac{1 + w^{(k)}_{o_i}}{2}\frac{1 - w^{(k)}_{o_j}}{2} O_i \tilde{O}_j^{(k)} \right) \;. \end{align}\]

Performing the sum over the intervals \(k\), we obtain that finite pulse duration effects contribute to the average Hamiltonian in the following way \[\begin{align} \frac{1}{n}\sum_{k=1}^n \frac{1}{t_p} \int_0^{t_p}dt\;(Q^{(k)}(t))^{-1} H_R Q^{(k)}(t) &= \sum_{O = X, Y}\sum_{i<j}J^O_{ij}\left( \frac{3}{8}(1 + \delta_{o_i, o_j})O_i O_j + \frac{s_i s_j}{2}\mathcal{C}_{ij} + \frac{2}{\pi}\left( s_i {\mathcal{O}'}_i^j O_j + s_j O_i {\mathcal{O}'}_j^i \right)\right) = \\ &= \frac{3}{8}\tilde{H} + \frac{3}{8}H_R + \tilde{H}_{\text{FP,err}}^{(1)} + \tilde{H}_{\text{FP,err}}^{(2)} \end{align}\] where \[\begin{align} \mathcal{C}_{ij} =& \frac{1}{n}\sum_{k=1}^n \frac{1 - w^{(k)}_{o_i}}{2}\frac{1 - w^{(k)}_{o_j}}{2}\tilde{O}^{(k)}_i\tilde{O}^{(k)}_j \\ {\mathcal{O}'}_i^j =& \frac{1}{n}\sum_{k=1}^n \frac{1 - w^{(k)}_{o_i}}{2}\frac{1 + w^{(k)}_{o_j}}{2}\tilde{O}^{(k)}_i \end{align}\] and \[\begin{align} \tilde{H}^{(1)}_{\text{FP, err}} =& \frac{1}{2}\sum_{O = X, Y}\sum_{i<j}J_{ij}^O s_i s_j\mathcal{C}_{ij} \\ \tilde{H}^{(2)}_{\text{FP, err}} =& \frac{2}{\pi}\sum_{O = X, Y}\sum_{i<j}J_{ij}^O\left( s_i {\mathcal{O}'}_i^j O_j + s_j O_i {\mathcal{O}'}_j^i \right) \;. \end{align}\] Summing all the contribution to the average Hamiltonian, we obtain \[\begin{align} \tilde{H}_{\text{FP}} =\left(1 - \frac{5 \epsilon_{\text{FP}}}{8} \right)\tilde{H} + \frac{3\epsilon_{\text{FP}}}{8}H_R + \epsilon_{\text{FP}}(\tilde{H}^{(1)}_{\text{FP, err}} + \tilde{H}^{(2)}_{\text{FP, err}}) \;. \end{align}\] In the next subsection, we show how to remove the last two terms in AHT using the same second averaging as done with rotation angle errors, and to engineer \(\tilde{H}\) up to Trotter errors even with finite pulse duration effects.

14.3 Robustness conditions↩︎

As already done with rotation angle errors, we now show that it is possible to exploit the sign functions \(s_i\) to mitigate the finite pulse duration errors. We choose the sign functions to change from a cycle to the other as \(s_i^l = w^{(l)}_{e_i}\). Averaging the Hamiltonian errors over \(L\) cycles, \(L\) being the periodicity of the Walsh functions \(w_{e_i}\) with largest Walsh index, gives \[\begin{align} \frac{1}{L}\sum_{l = 1}^L \tilde{H}^{(1, l)}_{\text{FP, err}} =& \frac{1}{2}\sum_{O = X, Y}\sum_{i<j}J_{ij}^O \delta_{e_i, e_j}\mathcal{C}_{ij} \\ \frac{1}{L}\sum_{l = 1}^L \tilde{H}^{(2, l)}_{\text{FP, err}} =& \frac{2}{\pi}\sum_{O = X, Y}\sum_{i<j}J_{ij}^O\left( \delta_{0, e_i} {\mathcal{O}'}_i^j O_j + \delta_{0, e_j} O_i {\mathcal{O}'}_j^i \right) \end{align}\] Then, choosing \(e_i \neq e_j\) \(\forall i, j\) and \(e_i\neq 0\) \(\forall i\), the doubly averaged Hamiltonian becomes \[\label{eq:app95corr95avg95H95FP} \begin{align} \tilde{H}'_{\text{FP}} =\left(1 - \frac{5 \epsilon_{\text{FP}}}{8} \right)\tilde{H} + \frac{3\epsilon_{\text{FP}}}{8}H_R \;. \end{align}\tag{41}\] To recover \(\tilde{H}\) in AHT, we first take care of the second term, which it is possible to get rid of by simulating the interval \(k=1\) (for which \(H^{(1)} = H_R\)) for a reduced time. This can be seen by expanding in Eq. 41 the average over pulses with respect to the first term \[\begin{align} \left(\tilde{H}'_{\text{FP}} - \frac{1}{n}H_R\right) + \frac{1}{\tau}\left(\frac{\tau}{n} - \frac{3nt_p}{4} \right) H_R =\left(1 - \frac{5 \epsilon_{\text{FP}}}{8} \right)\tilde{H} \end{align}\] resulting in the reduced time \(\frac{\tau}{n} - \frac{3nt_p}{4}\) for the interval \(k=1\). Finally, the factor \((1 - \frac{5 \epsilon_{\text{FP}}}{8})\) in front of \(\tilde{H}\) just leads to the renormalization of the total simulation time \(T\), which can be taken care of by increasing it as \(T \rightarrow T\left(1 - \frac{5 \epsilon_{\text{FP}}}{8}\right)^{-1}\). We note that these robustness conditions exactly remove finite pulse duration errors in AHT.

We benchmark numerically the effect of these robustness conditions with the cluster state preparation protocol via Ising dynamics for \(N=6\) qubits. The numerical results for the fidelity error are plotted in Fig. 9, where we observe that the fidelity error, being in AHT when the robustness conditions are not respected, collapse to the Trotter error when they are met instead.

As a last note, we highlight how the robustness conditions to mitigate finite pulse duration errors (\(e_i \neq e_j\), \(e_i\neq 0\)) also imply the robustness against rotation angle errors. This means that it is possible to mitigate both errors at the same time, without additional requirements than those asked to mitigate finite pulse duration errors. While this mitigates exactly the individual contributions of these two sources of errors, there might be additional errors due to their mutual contribution, which are not considered here. A numerical example for the combined mitigation of finite pulse duration and rotation angle errors in plot in Figs. 4, 9, which shows a reduction of the fidelity error in the cluster state preparation protocol of orders of magnitude.

15 Dynamical decoupling with Walsh functions↩︎

We derive the conditions for dynamical decoupling with Walsh functions from external disorder fields. As done in the main text, we consider a resource Hamiltonian which also has space-dependent single-qubit fields \[H_R' = H_R + H_{\text{ext}}\] with \[H_{\text{ext}} = \sum_i (h^x_i X_i + h^y_i Y_i + h^z_i Z_i) \;.\] We now compute the effect of a Walsh sequence on the Hamiltonian term \(H_{\text{ext}}\). Applying a Walsh sequence, we note that the \(X\) disorder term would become \[\tilde{H}_{\text{ext}}^X = \frac{1}{n}\sum_{k=1}^n \sum_i h^x_i w^{(k)}_{x_i} X_i = \sum_i \left(\frac{1}{n}\sum_{k=1}^n w^{(k)}_{x_i} \right) h^X_i X_i = \sum_i \left(\frac{1}{n}\sum_{k=1}^n w_0^{(k)} w^{(k)}_{x_i} \right) h^X_i X_i = \sum_i \delta_{0, x_i} h^X_i X_i\] since \(w^0_k = 1\;\forall k\), and applying the orthonormality condition for Walsh functions. The same reasoning can be applied to the \(Y\) part. Finally, we note that the sign of the \(Z\) operator is constrained by the signs of \(X\) and \(Y\), since \[p^{(k)}_i Z (p^{(k)}_i)^{-1} = i p^{(k)}_i X (p^{(k)}_i)^{-1} p^{(k)}_i Y (p^{(k)}_i)^{-1} = i (w^{(k)}_{x_i}X_i) (w^{(k)}_{y_i}Y_i) = w^{(k)}_{x_i} w^{(k)}_{y_i} Z_i\] which implies \[\tilde{H}_{\text{ext}}^Z = \sum_i \left(\frac{1}{n}\sum_{k=1}^n w^{(k)}_{x_i} w^{(k)}_{y_i} \right) h^Z_i Z_i = \sum_i \delta_{x_i, y_i} h^Z_i Z_i\] applying the orthonormality conditions for Walsh functions.

16 Derivation of the quantum circuit for the surface code↩︎

We show how to express the quantum circuit for the stabilizer measurement Fig. 5a) in terms of Ising interactions and single-qubit rotations. Our starting point is the expression of the gate \(C_{XX}=HC_{NOT}H\) \[C_{XX} = H( \ket{0}\bra{0}\otimes \mathbf{1} +\ket{1}\bra{1}\otimes X)H= \ket{0}_X\bra{0}\otimes \mathbf{1} +\ket{1}_X\bra{1}\otimes X\] with \(\ket{a}_X=H\ket{a}\) the eigenstates of the \(X\) operator. This can be mapped to Ising interactions as follows \[\begin{align} C_{XX} &=\exp(i \pi (\ket{1}_X\bra{1}\otimes \ket{1}_X\bra{1})) \nonumber \\ &= \exp(i (\pi/4) (1-X)\otimes (1-X)) \nonumber \\ &= R_{XX}^{\pi/2}(R_X^{\pi/2}\otimes R_X^{\pi/2}) \end{align}\] with \(R_{XX}^{\pi/2}=\exp(i X\otimes X \pi / 4)\), \(R_X^{\pi/2}=\exp(-i X \pi / 4)\). Note that the last equality holds up to an irrelevant global phase. This leads directly to the circuits shown in Fig. 5.

References↩︎