Power-optimal, stabilized entangling gate between trapped-ion qubits

To achieve scalable quantum computing, improving entangling-gate fidelity and its implementation-efficiency are of utmost importance. We present here a linear method to construct provably power-optimal entangling gates on an arbitrary pair of qubits on a trapped-ion quantum computer. This method leverages simultaneous modulation of amplitude, frequency, and phase of the beams that illuminate the ions and, unlike the state of the art, does not require any search in the parameter space. The linear method is extensible, enabling stabilization against external parameter fluctuations to an arbitrary order at a cost linear in the order. We implement and demonstrate the power-optimal, stabilized gate on a trapped-ion quantum computer.

Representing and processing information according to the laws of quantum physics, a quantum computer may surpass the computational power of a classical computer by many orders of magnitude, and is expected to transform areas such as machine learning [1,2], cryptosystems [3], materials science [4,5], and finance [6,7], to name only a few. Improving the reliability of quantum computation beyond the level of today's machines [8][9][10] is therefore critical to promote the quantum computer from a subject of academic interest to a powerful tool for solving problems of practical importance and utility.
The trapped-ion quantum information processor (TIQIP) is one of the most promising architectures for achieving a universal, programmable quantum computer, operating according to the gate model of quantum computing. Apart from a set of single-qubit gates, only a single entangling, two-qubit gate is necessary for achieving this goal. Today's TIQIPs [8,9] typically use an Ising xx gate, following the Mølmer-Sørensen protocol [11][12][13], as the two-qubit native gate. Its best reported fidelity is 99.9% [14,15], which may be compared with the best reported fidelity 99.9999% of single-qubit gates [16]. A host of pulse-shaping techniques have been devised [9,12,13,17,18] to better control the underlying trapped-ion quantum systems for more efficient xx gate implementation, while reducing errors.
Highlighting the importance of efficient and robust implementation of xx gates, Fig. 1 shows the resource requirements for various quantum computations. For this figure and for near-term, pre-fault-tolerant (FT) quantum computers, we considered variational quantum eigensolvers that compute the ground state of the water molecule [19], a material spin-dynamics undergoing state-evolution according to the Heisenberg Hamiltonian [5], a quantum approximate optimization algorithm addressing a maximum-cut problem relevant for various optimization problems [20], the widely-employed quantum * Electronic address: blumel@ionq.co † Electronic address: grzesiak@ionq.co ‡ Electronic address: nam@ionq.co Fourier transform subroutine [21], quantum factoring of a 1024-bit integer [22], which is meaningful for cybersecurity, and data-driven quantum-circuit learning for certain visual patterns [2]. The resource-cost metric used in the pre-FT regime considered here is the gate count for two-qubit xx gates, since these are the gates that limit algorithm performance. For the FT regime, in addition to the FT-regime-optimized versions of Heisenberg-Hamiltonian simulations, the quantum Fourier transform, and integer factoring, we considered Jellium-and Hubbard-model simulations for condensed-matter systems [23], the Femoco simulation [4], relevant for a certain nitrogen-fixation process that can help make fertilizer production more efficient, and solving difficult instances of satisfiability problems [24,25]. The resourcecost metric used for the FT regime is the number of t gates. Note that each t gate in FT quantum computing requires, e.g., a distillation process, typically referred to as a magic-state factory [26,27]. Each distillation process for the t-gate implementation in the FT regime requires at least a few tens of two-qubit gates, such as xx gates, at the native, hardware-implementation-level [26,27]. Optimizing the xx-gate implementation on a TIQIP is therefore critical for both pre-FT and FT quantum computing, and the construction and experimental demonstration of robust, power-optimal pulses for xx-gate implementation on a TIQIP is the focus of this paper.

I. POWER-OPTIMAL TWO-QUBIT ENTANGLING GATE
There is only so much power that optical components can withstand. But more importantly, increased power leads to reduced gate fidelity due to, e.g., carrier coupling [28], ion-ion crosstalk [9,29], and spontaneous emission from intermediate Raman levels [30]. Therefore, it is important to construct gates that, in addition to being stabilized against control-parameter fluctuations, require the least amount of power possible, i.e., they need to be power-optimal. In this paper we present a comprehensive, scalable approach to the construction of stabilized, power-optimal xx gates that is based on the Mølmer-FIG. 1: Resource requirement of various quantum circuits as a function of system size. For the pre-FT regime, the resource cost is measured in terms of required number of xx gates. For the FT regime, the resource cost is measured in t gates, where each FT t gate requires tens of hardwareimplementation-level xx gates [26,27]. Shown are the watermolecule ground-state computation (Water) [19], Heisenberg-Hamiltonian simulation (Heisenberg) [5], maximum-cut optimization (QAOA) [20], the quantum Fourier transform (QFT) [21], integer factoring (Shor) [22], generative-model quantum machine learning (QML) [2], Jellium-and Hubbardmodel simulations (Jellium, Hubbard) [23], and the Femoco simulation (Femoco) [4]. Grover's algorithm [24] implementation (not shown) that solves known difficult satisfiability problems [25] requires 2000 qubits and 2·10 27 t gates. See Supplementary Information (SI) section S1 for details.
Sørensen Hamiltonian of an ion chain interacting with laser pulses: Here i and p label the ions and the motional modes, respectively, η i p is the Lamb-Dicke parameter, σ i x and g i (t) are the Pauli-x operator and the pulse function acting on ion i, and ω p and a p are the frequency and mode operator of motional mode p, respectively. A judicious choice of pulse functions generates an xx gate that induces entanglement between two trapped-ion qubits, defined by the unitary operator where θ ij = 4χ ij denotes the degree of entanglement between ions i and j. To induce the desired xx gate in practice, all motional modes of the ion chain need to be decoupled from the computational states of the qubits at the end of the gate operation [8,9,11], leaving only the spins entangled. For an N -qubit system, and assuming that the same pulses are directed at ions i and j, these constraints are of the form where τ , a free parameter, is the pulse duration. The degree of entanglement between qubits i and j is obtained To find a power-optimal pulse, we require that the norm of g is minimized, while g still satisfies (3) and (4) exactly. This can be achieved by expanding g in a complete basis, e.g., the Fourier-sine basis according to g(t) = n=1 A n sin(2πnt/τ ), which spans the entire desired function space over the gate-time interval τ . Restricted to a finite sub-space with basis-function amplitudes A n , n = 1, 2, .., N A , with sufficiently many N A basis functions, the constraint (3) can be written in linear algebraic form as shown on the right-hand side of (3), where M pn is the time integral of the product between the nth basis function and e iωpt . Therefore, to satisfy (3), all that is required is to draw amplitudes A n from the null-space of M [31], where the null space is defined as the vector space that is mapped to zero under the action of M . Similarly, the constraints (4) can be denoted in linear algebraic form, as stated in the second line of (4), where the matrix D has elements D nm , defined as the psum of the double integrals in (4) of the product between sin[ω p (t 2 − t 1 )] and the nth and mth basis functions that stem from expanding the two g functions. Thus, defining the symmetric matrix S = (D + D T )/2, (4) can be satisfied, including the requirement of minimal norm of g, by finding the appropriate linear combination of the null-space vectors of M that combine to the eigenvector of R with maximal absolute eigenvalue, where R is the null-space projected matrix S.
Our approach is linear and satisfies the two constraints (3) and (4) exactly. Since (3) can be split into even and odd symmetry components, it is possible to consider only N out of 2N real constraints of (3) and induce, at our discretion, a pulse that is symmetric or anti-symmetric about its center. Additionally, because, e.g., the Fourier basis is complete in its respective symmetry class, the resulting g(t) is provably optimal in minimizing the norm of g(t), which corresponds to minimizing the average power required to induce a xx gate. There is no iteration of any kind necessary. For instance, searching for an optimal solution in the parameter space, such as in [12,13,17,18,29,32,33], is not necessary in our approach. Since only matrix operations are required to arrive at the optimal pulse solution, the optimal pulse is obtained in time O(N 3 A ) [34]. Figures 2a and b show the amplitudes A n for a sample pulse function of the form g(t) = N A n A n sin(2πnt/τ ) for N A = 10000 and τ = 300µs. As expected, the |A n | are large for 2πn/τ ≈ ω p , showing that, to induce the desired xx interaction between two qubits via motional modes, the frequency components of the pulse function g(t) need to be reasonably close to the motional-mode frequencies. We confirmed that a N A = 1000 basis-function solution essentially results in the same A n spectrum, vi-sually indistinguishable from that with N A = 10, 000, when overplotted. This demonstrates the robustness of our method with respect to the basis size.
The pulse function g(t) corresponding to the amplitudes A n shown in Figures 2a and b is shown as the green line in Fig. 3(a). Since g(t) is a fast-oscillating function, it is instructive to write it in the form where Ω(t) is the envelope function of the pulse (orange line in Fig. 3a and ψ(t) = t 0 µ(t )dt is the phase function, where µ(t) is the detuning function [8,9]. We see that the amplitude of the pulse function is relatively flat, which implies that the average power minimization obtained by the g-norm minimization is essentially as good as minimizing the peak power of the pulse. Figure 3b shows the detuning function µ(t). Consistent with the large Fourier amplitudes near the motional-mode frequencies in Fig. 2a, the demodulated frequency hovers around the motional-mode frequencies.
Because the minimal-power pulse function can be determined efficiently, it is straightforward to investigate the power requirement of the optimal pulse as a function of system size. Figure 2c shows the maximal power of the optimal pulse max t g(t), obtained with N A = 1000, for system sizes ranging from 5 to 100 ions. The power is consistent with our analytical bounds (see SI section S5). Additionally, since according to the analytical results the power requirement is inversely proportional to the gate duration, the power optimum implies gate-time optimum for a given power budget. Thus, for a given amount of maximally available power, the power-optimal pulse is the fastest possible for xx gate execution. The ion displacement in position-momentum phase space for each mode ω p entering into the computation of our sample pulse function g(t) shown in Figs. 2a and b, is shown in Fig. 2d.

II. CONTROL-PULSE STABILIZATION
Because the pulse function is constructed using a completely linear method, any additional linear constraints may be added, which still results in a power-optimal pulse when generated according to the steps discussed in the previous section. As an example, we show here how to stabilize the pulse against errors in external parameters, such as mode-frequency fluctuations.
To stabilize against fluctuations of ω p , we start by expanding the number of constraints (3). Explicitly, we add where k denotes the order of stabilization. Since the additional constraints in (6) are linear, all we need to do to stabilize the pulse up to Kth order is to include the additional linear equations (6) in the coefficients matrix M . The decoupling between the computational states of the qubits and the motional modes is thus achieved exactly, and the pulse is stabilized by using up N (K + 1) degrees of freedom. Figure 2d shows the phase-space trajectories for the stabilized pulse with K = 3. Compared to the K = 0 case, the general structure of the pulse with K = 3 remains the same -the frequency components are centered around the motional-mode frequencies and the phasespace closure is guaranteed. In Fig. 4a, the infidelity of stabilized pulses K = 0, 1, .., 8 is shown as a function of the extent of the ω p fluctuations. Considered are pulses with duration τ = 300µs over the five-ion chain considered in the previous section. The widths of the infidelity curves, extracted at infidelity of 0.001, increase from ∼0.1kHz to ∼13kHz as K is increased from 0 to 8 (see Fig. 4b). The power requirement of the stabilized pulses for each K is shown in Fig. 4c; the requirement scales linearly in K. Figure 4d shows the width-scaling for each K as a function of different choices of gate duration τ . The effect of the stabilization increases inversely proportional to the gate duration.
The most straightforward way of experimentally implementing our AMFM pulses is via an arbitrary waveform generator (AWG) [35,36]. However, if an AWG is not available, implementation via the decomposition (5) (demodulation) is also possible (see Supplementary Information, section S9, for more detail).
In various contexts our protocol has already been implemented, tested, and verified experimentally. In [8] it was used in its simpler amplitude-modulated version to benchmark one of the IonQ quantum computers. In [36] it was used as the basis for demonstrating a new fidelity trade-off scheme called extended null-space (ENS) [36]. By sacrificing negligible amounts of fidelity, via ENS, an "add-on" to our basic AMFM scheme, additional power savings of up to an order of magnitude can be realized. This demonstrates that our protocol is extensible and adaptable.
In both applications our power-optimal protocol has proved its experimental utility.
While our protocol has already found experimental applications, its stabilizing effect in its exact [see (3)] version, employing simultaneous amplitude and frequency modulation, has not yet been demonstrated. This is done in the following section.

III. EXPERIMENTAL DEMONSTRATION
We implement the power-optimal, stabilized entangling gate on a trapped-ion quantum computer that has been described elsewhere in detail [8]. A chain of seven ions is sideband cooled according to the protocol detailed in [37]. The two end ions are used to obtain equal spacing between the middle five individually addressed ions, which comprise the qubit register, and single-and twoqubit gate operations are driven by a two-photon Raman FIG. 2: Power-optimal pulse function g and spin-dependent force applied to an ion qubit. See SI section S4 for the relevant parameters of the sample case of five qubits considered here. a. Fourier-sine coefficients |An| of the pulse function g(t) = n An sin(2πnt/τ ), τ = 300µs. The tails of |An| decay according to ∼ 1/n. The arrows indicate the locations of the mode frequencies. The signs of An are color coded, i.e., negative An are depicted with orange and positive An are depicted with gray. The top scale shows the basis frequencies n/τ , which are in resonance with the mode frequencies at the locations of the arrows. b. Same as a but for a basis of 10,000 states. The main features occur at the same positions in n and in frequency, which shows convergence and basis independence. c. Scaling of the maximal power, maxt g(t) as a function of the number of qubits N . Gate time τ = 500µs. Orange circles: numerical results, heuristic bound. Gray curve: analytic bound derived in SI section S5. d. The time-dependent displacement in the phase space of ion number 1, where K = 0 (pink) and K = 3 (blue) are shown. The trajectories start and end at the origin.  Fig. 2a. a. The optimal pulse functionĝ(t) (thin green line) with its amplitude function (thick orange line), obtained by demodulatingĝ(t) as detailed in SI section S9. b. Detuning function µ(t) obtained by frequency demodulating the pulse function, using the method described in SI section S9. The frequencies of the motional modes are shown as the five horizontal lines. The sample motional-mode frequencies used to generate this pulse are listed in Table S1 and the set of η parameters used are listed in Table S2 in SI section S4.
transition at 355 nm. As described in [8] the coherence of our quantum computer has been completely characterized. In particular, the fidelity of two-qubit gates was determined to be ≈ 96%, measured via parity contrast and partial tomography as described in [8,14,15]. The propagator of (1) can be written in the form [12] where and xx(θ ij ) is defined in (2). For α p = 0, which is guaranteed exactly according to our protocol, an AMFM gate solution, computed from the measured motional-mode spectrum, implements the unitary in (2) over two qubits with ≈ 96% fidelity. Since α p = 0 corresponds to zero mode-frequency drift, the ≈ 4% loss in fidelity are not due to mode-frequency drift but are due to other processes, such as beam-steering errors, laser-power fluctuations, fluctuating ambient electric and magnetic fields, etc.. Starting in the intial state |00 , the ideal gate operator (2) strictly preserves the even-parity population P even = P 00 + P 11 , where P mn is the population in |mn after the gate pulse is over. When motional-mode frequencies start to drift, α p becomes nonzero and the propagator V in (7) is "switched on". This causes population transfer into the odd-parity states |01 and |10 , which, in turn, causes a reduction of P even . Thus, P even is a sensitive probe of stabilization of an AMFM gate against mode-frequency drift. We strongly emphasize that mode frequency drift causes a reduction of fidelity that is completely independent of the other sources of infidelity mentioned above. Thus, if the mode frequency drift is compensated, i.e., V is "switched off" by constructing a steering pulse g(t) that assures α p ≈ 0 to high order, as accomplished by our protocol, the gate fidelity is unchanged from its value at zero mode frequency offset. This was verified experimentally in [36] in the way of spot-checks using contrast measurements for selected values of mode frequency drift. Figure 5 shows the even-parity population P even as a function of mode-frequency offset, measured after applying an xx(π/2) gate designed for zero mode-frequency offset, to the initially prepared two-qubit state |00 . Three different pulses were used to implement xx(π/2), with moment stabilization orders K = 1, 4, 7. The grey line in each frame of Fig. 5 shows the analytical infidelity (4/5) p (η i p ) 2 +(η j p ) 2 |α p | 2 [30], with a constant 4% offset to account for the other errors that are independent of the gate-frequency offset. The shaded region indicates the range of mode-frequency offsets for each moment stabilization where the theoretical fidelity F > 0.99. The width of the shaded region can be seen to increase for increasing moment. In fact, for K = 1 stabilization (top frame of Fig. 5) acceptable fidelity is achieved only over a range of less than 1 kHz, outside of which P even drops precipitously to around 65%. K = 4 stabilization (middle frame of Fig. 5) achieves stability over a mode-frequency range of about ±2 kHz, a factor-2 improvement of allowed mode-frequency fluctuations accompanied by a restoration of the gate fidelity back up to ≈ 96%. In the interval between 1 kHz and 2 kHz we see an improvement of P even of about 30%. This is not a small effect. This is a sizable stabilization effect adding to the proven utility of our protocol. K = 7 stabilization (bottom frame of Fig. 5) shows a similar ≈ 30% improvement of P even over an even larger frequency interval ≈ ±3 kHz), a factor-3 improvement of the stable frequency range. More details concerning our The infidelity (see SI section S8 for detail) as a function of the motional-mode frequency drift ∆f . All mode frequencies were drifted according to ωp → ωp + 2π∆f . b. The width of the infidelity curves in a for various error tolerances = 10 −3 , 10 −5 , and 10 −7 , as a function of the highest moment K of stabilization. c. The maximal power requirement maxt |gK (t)| of the control pulses, as a function of the order or moment of stabilization K. The power requirement suggests a linear scaling in the moment of stabilization. d. The width of the infidelity curves for various different orders of stabilization K = 0, 2, 4, and 6, as a function of the gate duration τ for a fixed error tolerance level = 10 −3 . The data suggests ∼ 1/τ scaling of the width. experiment can be found in SI section S11.

IV. DISCUSSION
Application of our protocol is not limited to the case of the trapped ion-chain quantum computer architecture [8,9]. It can be applied to all quantum computing architectures that rely on phase-space closure, such as the QCCD architecture [38]. It is even applicable to the superconducting quantum computer architecture that also relies on pulse shaping [39].
Our work is most closely related to [28,40], who use a multi-tone approach to construct stabilized two-qubit gates. However, there are substantial differences. For instance: (i) The authors of [28,40] formulate and demonstrate their protocol for the special case of trapped two-ion crystals, and since it relies on an analytical solution for the gate pulses, it is not scalable to N -ion, manymode systems. (ii) The multi-tone approach in [28,40] requires the different tones (which can be interpreted as Fourier components) of g(t) to be equi-spaced. While this is very closely related to our expansion of g(t) in a Fourier-sine series, equi-spacing in the frequency of basis functions is not necessary in our approach. Any basis can be used to express g(t). For instance, we found that short pulses are better represented in a polynomial basis, whose Fourier components are not equi-spaced. Thus, any kind of chirped pulse can be used in our method. (iii) While, indeed, stability against several types of errors was demonstrated in [28,40], these conditions, to much higher order than demonstrated in [28,40], are more easily implemented according to our protocol, which requires only added linear equations for each stability condition FIG. 5: Experimental demonstration of K-moment stabilized two-qubit gates. A single xx(π/2) gate is applied to the fiducial state |00 , and the even-parity population is plotted as a function of the gate frequency offset. For increasing moment stabilization (shown are K = 1, 4, 7), the region of high even-parity population increases, as indicated by the detuning-robust region shaded in blue. The gray line shows the analytical expression (S60), valid in the low-error limit, with a 4% offset to account for other sources of infidelity that are independent of the gate frequency. Notice that our experiment reproduces the oscillating fine structure of the fidelity expression (S60). Error bars on the experimental data are 1σ confidence intervals, sampled from a binomial distribution, and each point represents 300 realizations of the experiment.
or constraint. (iv) As discussed in SI sections S15, S16, and S17, our protocol can be extended to include stability against drifts in the degree of entanglement. (v) Power optimization was not considered in [28,40]. Thus our protocol, presenting a unified, scalable approach to arbitrarily stabilized N -ion gate construction, goes far beyond the methods presented in [28,40] and is thus a substantial advance over existing methods. The ability to symmetrize the pulse solution gives rise to potential additional room for robustness with respect to errors. Since, e.g., the inner products between a symmetric pulse function g(t) and the antisymmetric part of the constraint e iωpt in (3) are zero (see SI section S3 for details), akin to echos, as long as the pulse function is modified symmetrically due, e.g., to implementation defects, half of the null-space conditions in (3) is still exactly satisfied, while leaving the error entirely in the part of the constraint with the opposite symmetry. The knowledge that the error lies in the oppositely symmetric part leaves room for a secondary echo, wherein the sign of the errors may be flipped. In case the errors cannot be manipulated to be echoed out in this particular symmetry, other symmetries may be considered at the pulse-construction level, rendering our approach an integrated protocol that can be designed to be robust against different symmetry classes of errors.
The methodology and paradigm used to construct our power-optimal pulses is general and can readily accept incomplete bases to result in pulses that are subjected to additional constraints. Beyond the AM, FM, PM pulses that can be obtained by appropriately demodulating the pulse function according to (5), a step-pulse approach that has been used in the literature [12,13,29] can also be derived. In fact, as shown in SI section S12, by carefully tuning the gate duration time, even the beneficial symmetric pulse structure can be preserved.
Adding to the generality is the application to the Efficient Arbitrary Simultaneously Entangling (EASE) gates [29], where any combination of quadratically many pairs of qubits can be entangled to any degree of entanglement. As detailed in SI section S13, because our approach is linear and the EASE-gate approach is amenable to any linear approach, it is straightforward to adapt the pulseconstruction method presented here to the EASE protocol. Together with the power-or time-optimality guaranteed in our pulse construction by design, then, the EASE gate equipped with our method enables one of the fastest ways to implement as many entangling gates as possible in a TIQIP.
The moment-stabilization adds robustness against those errors induced by not properly decoupling the computational states from the motional modes. However, unitary errors in the computational space may still linger, since the entanglement degree is sensitive to, e.g., ω p fluctuations (see SI section S14). In practice, this may be fended off by calibrations, i.e., by monitoring how the degree of entanglement changes over time and adjusting the amplitude of the illuminating beams to compensate for this change. Note that the shape of the pulse does not change; amplitude scaling suffices. If frequent calibrations are impractical, active stabilization in the entanglement degree χ ij may be implemented, which involves finding pulses that, apart from satisfying (3), (4), and (6), also satisfy where K χ is the maximal order of desired χ stabilization. In SI sections S15, S16, and S17, we offer three methods which either achieve (9) approximately (see S15) or completely (see S16 and S17). While the projection method presented in SI section S15 achieves (9) only approximately, and only for K χ = 1, it does reduce the errors originating from fluctuating mode frequencies substantially. A qualitative improvement over the projection method are the two-pulse moments method presented in S16 and the hybrid method, presented in S17. These two methods fulfill (9) exactly, and consequently actively stabilize the degree of entanglement over a large interval of frequency drift [see SI section S16, Fig. S6a and SI section S17, Fig. S7]. However, constructing the two pulses characterizing these methods is computationally more expensive and the pulses require more power. Therefore, whether to use the simple projection method described in S15 or the two-pulse methods described in S16 and S17, depends on experimental conditions and available computational resources. Stabilization against other parameters, such as the Lamb-Dicke parameters η i p or the amplitude of the laser beam, encoded in the norm of g(t), is also possible. No- and ∆ g are small constants, does not affect the decoupling condition (3). The effect of the errors is confined to the degree of entanglement (4), i.e., χ ij →χ ij +∆χ ij , where ∆χ ij is the error in χ ij that arises from ∆ η i p and ∆ g . As discussed in SI section S18, this can be adequately compensated by, e.g., a broadband compensation sequence applicable for the two-qubit case.
We note that, while technically challenging, in principle, it is possible to directly implement our Fourier-basis pulse solution g(t) using a multi-tone laser. As shown in SI section S19, by implanting N A different colors with different amplitudes to the beams that address ions, then locking the phases of them, we can induce the desired evolution of the xx on a TIQIP. The technique here is similar to the discrete multi-tone, widely used in communication lines. The development of the technology in the optical regime remains as a promising avenue for future research.

V. CONCLUSION
Formally speaking, there are infinitely many smooth solutions that qualify as adequate pulse functions. Out of these infinitely many possible solutions, our protocol extracts the power-optimal solution without any iterations or parameter scans. Including symmetry and stabilization, the solution is also robust against errors. With an AWG, as already demonstrated [36], our AMFM pulses can be implemented directly experimentally to produce an optimal xx gate. If AWG technology is not available, as in our experiments reported in this paper, implementation via demodulation (see SI section S9) is also an option. Demodulation of our AMFM pulses also shows that their amplitude function is nearly flat, i.e., average power minimization is an excellent approximation of maximal power minimization. Indeed, an exact analytical bound on power and its comparison to the demodulation results shows that the optimal solution is close to the bound. With the xx gate implemented using the pulse constructed according to our protocol, just about any quantum algorithms can now be implemented with minimal power requirement, or in the shortest possible time for a given power budget, at the two-qubit gate, physical implementation level. This provides decisive advantages in improving both noisy, near-term TIQIPs and fault-tolerant TIQIPs to come in the future.
For case (i), we considered pre-FT HF+7 and HF+21 cases, where HF denotes the Hartree-Fock method detailed in [19], and 7 and 21 denote different approximation qualities. The xx gate counts for the two cases are available in Fig. 2b of [19].
For case (ii), we considered the Heisenberg Hamiltonian applied to spins with their connectivity specified by (k−d−n) graphs, where k denotes the degree, d denotes the distance, and n denotes the number of vertices of the graph. Specifically, the graphs considered are (3−5−70), (4−4−98), and (5−3−72). For the pre-FT cases we used cnot gate counts reported in the pre-FT part of Table I of [5]. For the FT cases, we used t-gate counts reported in the FT part, specifically the RUS part, of the same table.
For case (iii), we considered the quantum approximate optimization algorithm in the pre-FT regime with eight stages, based on its performance compared to the well-known instance of semidefinite programming called Goemans-Williamson approximation algorithm [41]. The graphical representation of how the quantum algorithm solving the maximum cut problem performs with stage numbers 2 0 , 2 1 , .., 2 5 may be found in Fig. 2 of [20]. Each stage requires n(n − 1)/2 xx gates, as can be seen from Eq. (7) of [20].
For case (iv), we considered the approximate quantum Fourier transform [21], where all controlled-rotation gates with rotation angles less than π/2 b , b = log 2 (n), where n is the number of qubits, are removed. For the pre-FT regime, one xx gate was expended per controlled-rotation gate. For the FT regime, see Table 1 of [21].
For case (v), we used the implementation presented in [22]. While an explicit resource cost is not available, an estimate is available in section A of the appendix of [42]. The implementation in [22] uses 4n 3 + O(n 2 log(n)) gates and 3n + 6 log(n) + O(1) qubits, assuming an arbitrary two-qubit gate may be implemented. For the pre-FT regime, each arbitrary two-qubit gate costs three cnot or xx gates, as per [43]. For the FT regime, see the discussion section A of the appendix of [42], which results in 16n 3 t gates.
For case (vi), we largely base the resource counts on Table 1 of [2], where several sample instances of bars-and-stripes patterns are explicitly considered for n ranging from 4 to 100. The expected xx gate counts are computed assuming the all-to-all connectivity available in the trapped-ion quantum information processor (TIQIP), and we used four layers in the training circuit (see Fig. 1 of [2] for further information) that worked well for a small system with n = 4.
For case (vii), Tables 3 and 4 of [23] detail the FT resource-cost for several different cases.
For case (viii), see Table 1 of [4] for the Femoco simulation. We used a serial version of structure 1 with accuracy of simulation of 10 −3 Hartree. Note this case appears to demand a fairly large amount of resources compared to other examples considered herein. This is due in part to the other examples being tailor-developed and co-designed for minimal resource requirements. Subsequent development since [4] has led to a reduction in the resource requirements by several orders of magnitude, inching closer to the cluster of points that appear in Fig. 1 of the main text. See [44] for details.
We also considered Grover's algorithm solving certain difficult instances of a Boolean satisfiability problem [24] with n variables and m clauses. Specifically, we considered "hole12", "Urq7 5", "chnl11x20', and "fpga13.12" problems, where the names were taken verbatim from Table 1 of [25]. To construct the FT circuit, we used k-control Toffoli gates to implement the Grover oracle [24], where k is the length of a clause. Specifically, we used m clean ancilla qubits to compute the satisfiability of m clauses individually, and used a m-control Toffoli gate with an additional ancilla qubit to implement the oracle. Whenever possible, we used relative Toffoli gates in [45] to reduce the t counts, while keeping track of the number of recyclable ancilla qubits in implementing the multi-control Toffoli gates. Together with a ncontrol Toffoli gate for the Grover diffusion operator, we obtain for "hole12" 2053 qubits and 2.094 · 10 27 t gates, for "Urq7 5" 4627 qubits and 2.031 · 10 40 t gates, for "chnl11x20" 8879 qubits and 4.931 · 10 70 t gates, and for "fpga13.12" 2717 qubits and 1.538 · 10 39 t gates, where we used π/4 · √ 2 n iterations for near-optimal results. Of course it is challenging to realize on the order of 10 27 , . . . , 10 70 quantum gates. However the presented scaling corroborates the need for efficient implementations of quantum gates in less demanding circumstances.

S2. Ising gate on a trapped-ion quantum information processor
The participating ions of an Ising xx gate couple to all motional modes [11], and have to be decoupled from the motional modes at the end of the gate. The relevant equations are [12]  where τ is the length of the pulse, i is the ion number, N is the total number of ions, p is the mode number, P is the total number of modes, Φ i is the initial phase, ω p are the motional-mode frequencies, and Ω i (t) is the amplitude function, i.e., the time-dependent Rabi frequency. The time-dependent phases ψ i (t) in (S1) are defined as where µ i (t) is the detuning function. In order not to start the pulse abruptly, we require Φ i = 0. For ease of presentation, we also assume, from now on, that the same pulse shape acts on all N ions, such that, together with the assumption of vanishing initial phase, (S1) acquires the simplified form where If the pulse acts simultaneously on ions i and j, the gate angle ϕ ij of the xx gate is given by [12] where and η i p is the Lamb-Dicke parameter [46], which describes the coupling strength of ion number i to motional-mode number p. A maximally entangling gate is achieved for ϕ ij = ±π/4. According to (S6), χ ij = χ ji , i.e., ϕ ij = 2χ ij , so that a maximally entangling gate requires Since both Ω(t) and sin[ψ(t)] are unknown, we combine them into one single pulse function Thus, for given motional-mode frequencies ω p and Lamb-Dicke parameters η i p , our task is to find a pulse g(t), which solves (S3) and produces |χ ij | = π/8 with minimal power requirement. Known solution methods include amplitude-modulation techniques [12,13], which require fixed detuning frequency µ 0 , frequency-modulation techniques [17], which require a given shape of the pulseenvelope function Ω(t), and phase modulation [18].
Our approach goes beyond previously demonstrated approaches in that we modulate amplitude, frequency, and phase simultaneously. In addition, we use a linear method, which yields the optimal pulse shape directly, without any iterations or parameter searches, using exclusively linear-algebra techniques. Figure S1 shows a comparison between RMS Rabi frequencies, [ Ω 2 (t) ] 1/2 /(2π) in MHz, required for fully entangling gates using AM and AMFM pulses. The data shows the lowest RMS Rabi frequencies obtained when scanning the gate time in steps of 1µs. For every gate time, the detuning of AM pulses is scanned in steps of 10 Hz, and the lowest RMS Rabi frequency obtained is presented on the plot. The data points were computed without any added stabilization to external sources of errors. While the required RMS Rabi frequency is lower for AMFM pulses, the stability to external errors is comparable. The execution times of the classical computations required to find the power-optimal pulses for every gate time are comparable, with 0.73 seconds for an AMFM pulse and 0.93 seconds for an AM pulse. This, however, changes drastically if higher-order moment stabilization is included, which requires the AM pulse to be represented by a large number of segments. In this case, the number of segments will approach the number of AMFM basis states and the execution time of the AMFM method will be shorter by a large factor.

S3. Symmetry classes
Since g(t) is a real function, the P complex equations (S3) for p = 1, . . . , P are equivalent to 2P real equations It follows that if (S9) is satisfied, any linear combination h p (t) = A p cos(ω p t) + B p sin(ω p t) (S10) We define two special linear combinations and which satisfy i.e., h p (t) are even and odd functions with respect to τ /2. We also define i.e., the even and odd components of the pulse g(t).
We call g (+) (t) the positive-parity pulse and g (−) (t) the negative-parity pulse. Then the P equations are satisfied automatically, which implies that for given parity, we have to satisfy only P real, nontrivial equations In analogy to the definition of the two parities for the pulse function g(t), we may also define even and odd pulse envelope functions, Ω (±) (t), and even and odd detuning functions, µ (±) (t), which are even and odd functions with respect to τ /2 according to respectively. For the examples presented in this paper, we choose pulses where both the pulse-envelope function Ω(t) and the pulse-detuning function µ(t) are of positive parity. This entails that ψ(t), according to (S4), has odd parity with respect to τ /2, so that sin[ψ(t)] is also of odd parity, resulting in a pulse function g (−) (t) of odd parity.
Thus, to illustrate our pulse-generation method, we will in the following focus on negative-parity pulses, g (−) (t), constructed from a positive-parity pulse-envelope function Ω (+) (t), negative-parity sin[ψ (−) (t)], and positiveparity pulse-detuning function, µ (+) (t). Since the pulse function is of negative parity, we expand the pulse into a Fourier-sine series according to where A n , n = 1, . . . , N A , are real expansion amplitudes and N A is chosen large enough to achieve convergence. The expansion (S19) provides the additional benefit of switching g (−) (t) off continuously at t = τ without a discontinuous jump to g = 0 at t = τ . It is straightforward to show that the expansion (S19) is indeed odd with respect to τ /2. The expansion (S19) is complete, i.e., any pulse function g (−) (t) with g (−) (t = 0) = g (−) (t = τ ) = 0 can be represented this way. Expanding the entire pulse g (−) (t) as a whole, and not Ω(t) and µ(t) separately, is natural, since neither Ω(t) nor µ(t) are known. In fact, expansion of the entire pulse function g(t) is the key idea that motivated our method of AMFM pulse construction.

S4. Pulse construction
We focus in this section on computing the poweroptimized pulse function g (−) (t) for a given set of motional-mode frequencies ω p and Lamb-Dicke parameters η i p , i = 1, . . . , N , p = 1, . . . , P . Since in this case, according to (S16), the P equations Using the expansion (S19) and the explicit form (S13) of h (−) p (t), we obtain the following set of real, linear equations where In matrix notation we may write (S21) in the form where M is the P × N A coefficient matrix of (S21) and A is the amplitude vector of length N A . In order for (S23) to have non-trivial solutions, we require N A > P . In general, then, M in (S23) will have rank P , and there exist N 0 = N A − P non-trivial solutions A (α) of (S23), α = 1, . . . , N 0 . Since N A > P , the matrix M is a rectangular matrix. This suggests to multiply (S23) from the left with the transpose, M T , of M , which turns (S23) into the eigenvalue problem, where Γ = M T M is a symmetric matrix, and we are looking for the N 0 eigenvectors A (α) of Γ with eigenvalues 0. The N 0 nontrivial vectors A (α) with eigenvalues 0 span the kernel of the matrix Γ, also known as the null space of Γ. Numerically diagonalizing Γ, its eigenvalues typically are of the order of 10 −12 in the null space, and several orders of magnitude larger in the complementary space. Thus, the transition from the null space to the complementary space is sharp, with eigenvalues jumping many orders of magnitude at the transition point. Therefore, the null space can be identified clearly and unambiguously. Without restriction of generality we may also assume that the null-space vectors are normalized. Since all null-space vectors A (α) have the common eigenvalue 0, the null space is degenerate. Thus, any linear combination of the N 0 null-space vectors A (α) are also null-space vectors, and we may assume that the A (α) form an orthonormal basis of the null space according to where δ αβ is the Kronecker symbol. Our goal now is to linearly combine the orthonormal null-space vectors A (α) with real expansion amplitudes Λ α to find the optimal null-space vectorˆ is optimal in the sense that it produces |χ ij | = π/8, according to (S7), and has the smallest possible norm which entails the smallest possible average power needed to execute a maximally entangling xx gate. Using (S27) with (S8) and (S7) in (S6), we obtain where D is a real N A × N A matrix with matrix elements Sinceˆ A T Dˆ A is a scalar, we can also writê is a symmetric matrix. Using (S32) and (S31) in (S29) we now obtain where Λ is the vector of expansion amplitudes Λ α , α = 1, . . . , N 0 , and R is the symmetric, reduced N 0 × N 0 matrix with matrix elements Since R is symmetric, it can be diagonalized, where, since R is real and symmetric, the eigenvectors V (k) can be assumed orthonormal. We now linearly combine the vector of expansion amplitudes Λ from the set of vectors V (k) according to According to (S28), we now have to determine the expansion amplitudes v k such that under the condition   Geometrically, (S35) is a principal-axis transformation, V (k) are the N 0 principal directions of R in the null space, (S37) is a N 0 -dimensional sphere of radius γ, and (S38) is a N 0 -dimensional conic section with principal axes |λ k | −1/2 . Thus, geometrically speaking, we are looking for the smallest sphere that touches the conic section. This is obviously achieved if the sphere is inscribed in the conic section and just touches the conic section along the principal axis with the smallest length, i.e., the largest |λ k | Thus, our optimization problem is solved: The optimal pulse (S27) is constructed with the help of the am-plitudesˆ where k max is the index of the eigenvalue λ k of (S35) with the largest modulus |λ k |, and where To illustrate the method discussed in this section, we show in the main text in Figs. 4a and b the optimal pulseĝ(t) obtained for N = 5 ions and P = 5 motional modes for mode frequencies and Lamb-Dicke parameters as shown in Tables S1 and S2, respectively. The pulse has a symmetric envelope function and is amplitude as well as frequency modulated.

S5. Analytical lower bound of required peak pulse power
In this section we derive an exact, closed-form, integral-free, analytical expression for the lower bound of the minimally required pulse power needed to operate an i ↔ j XX gate. To be specific, throughout this section we choose P = N , which is the mode in which our quantum computer is operated [8]. Generalizing the lower bound to the case P = N is straightforward.
We define We also define where we defined which, for fixed N , is essentially a constant, which depends only weakly on τ , i.e., For instance, for an 80 µs pulse, and the mode frequencies and η values listed in Tables S1 and S2, respectively, the first term in (S45) is 2 × 10 −5 while the second term is 5 × 10 −8 . Therefore, in practice, the second term in (S45) may be neglected. Since the Lamb-Dicke parameters η j p are proportional to the jth component of a unit vector [46], we have, on average, η j p ∼ 1/ √ N , which then, because of (S46), implies With these definitions, and using the Cauchy-Schwarz inequality for integrals, we obtain: Using sin 2 [ψ(t)] ≤ 1, which is valid for all arguments ψ(t), the most straightforward, exact estimate for σ is Using this in the inequality (S48) and solving for Ω max , we obtain or, transitioning to lab frequency, This is the formula used to compute the analytical lower bounds of minimally required power to operate an XX gate, stated in the lower half of implies that no pulse exists, even in principle, that would require lower power than indicated by (S50) [(S51), respectively] to operate an XX gate. We also see that, because of (S47), Ω max (f max , respectively) scales like ∼ N 1/4 . In many cases (S50) [(S51), respectively] may be sharpened if lower (µ min ) and upper (µ max ) bounds for the detuning function µ(t) are available (see, e.g., the main text Fig. 3b ), i.e., We define Since ψ(t), according to (S2), is defined via an integral, and since µ(t) > 0 for all t, ψ(t) is a monotonically : Analytical lower bounds of minimally required analytically computed peak power (lower triangle) and numerically computed peak power of optimal pulses (upper triangle) for gate combinations i ↔ j, mode frequencies as listed in Table S1, Lamb-Dicke parameters η as listed in Table S2, and τ = 300 µs. Powers quoted are in kHz. Basis size: NA = 1000; nmin = 1. increasing function of t. Therefore, in (S43), we may change variables from t to ψ to obtain [ψτ − cos(2ψ0 + ψτ ) sin(ψτ )] ≤ 1 2µmin (ψτ + 1) where, in the last inequality, we used (S53). With (S54), the inequality (S48) can now be stated in the form or, solved for Ω max , Transitioning from angular frequency to lab frequency in Hz, we obtain This is our central result. No pulse exists with a power lower than stated in (S57) if χ ij is determined by (4).
To illustrate our analytical result, we show in Table S3 a comparison between our analytical lower limit of peak pulse power and numerically obtained peak pulse powers for our sample case of N = 5 ions and P = 5 motionalmode frequencies as listed in Tables S1 and S2. We see that our analytical result is indeed lower than all numerically obtained peak pulse powers, but that both are qualitatively close.

S6. Power and execution time scaling
The execution time of our linear pulse-construction algorithm is dominated by two diagonalizations, i.e., the diagonalization of the matrix Γ = M T M [see (S24)] and the reduced matrix R [see (S34)]. The dimension of Γ is N A ×N A , and the dimension of R is (N A −P )×(N A −P ). Therefore, the execution time of our algorithm scales like ∼ N 3 A . Since, in general, N A P , the execution time is dominated by N A and depends on P only via N A > P , which is needed for a nontrivial null space. Therefore, the overall scaling is dominated by N A and the algorithm scales like ∼ N 3 A . We confirmed the ∼ N 3 A scaling of our algorithm in numerous pulse-generation runs.
We also investigated the scaling of pulse power in N with up to N = 50 ions. For our investigation of power scaling we generated motional-mode frequencies and Lamb-Dicke parameters according to the procedure outlined in [47]. We used simulated ion positions, approximately equi-spaced with a spacing of about 5 µm and a frequency ratio of axial to radial trap frequencies of ω x /ω r = 0.088. We focused on operating an XX gate between ions 1 and 3. For these parameters and for N = 50 particles we obtained an average motionalmode frequency spacing of ∆f = 1.46 kHz. We found that our algorithm is stable only if τ ∆f ≈ 1. Therefore, for our power-scaling simulations, we chose τ = 500 µs. The result of our power-scaling simulations in a basis of N A = 1000 states is shown in the main text Fig. 2c. We see that the power scales approximately like N 1/4 , which is consistent with the analytical power scaling (S50) with (S47) (gray full line in Fig. 2c). As pointed out in the main text, knowledge of power scaling is important since, apart from possibly damaging optical components when applying too much power, increasing power also enhances important sources of errors.

S7. Power optimality requires identical pulses
In this section we show that if we do not actively stabilize the degree of entanglement χ, a gate is power optimal if ions i and j participating in a two-qubit gate are illuminated with identical laser pulses, i.e., g i (t) = g j (t) = g(t). To show this, let A and B be the expansion amplitudes of g i and g j , respectively. Then, the degree of entanglement is χ = A T R B, where the symmetric matrix R is defined in (S34). Define P 2 A = A T A and P 2 b = B T B. Then, the task is to minimize P 2 = P 2 A + P 2 B under the constraint χ. Thus, the target function to be minimized is F ( A, B) = P 2 − λχ, where λ is a Lagrangian parameter. This yields two equations: From (S58) and (S59) we obtain immediately P 2 A = λ A T R B/2 = λχ/2 and P 2 B = λ B T R A/2 = λ A T R B/2 = λχ/2.
Thus, P A = P B , i.e., for power optimality the same power must be directed at both ions. From From (S58) and (S59) we further obtain A = (λR/2)(λR A/2) = λ 2 R 2 A/4 and B = (λR/2)(λR B/2) = λ 2 R 2 B/4, i.e., A and B satisfy the same eigenvalue equation. This means that, up to normalization, A and B are the same. Together with P A = P B , we now have g i (t) = g j (t) = g(t).

S8. Stabilization against mode-frequency fluctuations
In this section we show that our linear approach lends itself naturally to a method of constructing pulses that stabilize the fidelity of the xx gate against mode drifts and mode fluctuations. Due to uncontrollable effects, such as stray electromagnetic fields, build-up of charge in the trap due to photoionization or temperature fluctuations, the frequencies of the motional modes, ω p , will drift or fluctuate in time. Therefore, in a typical quantumcomputer run, one would determine the current values of ω p and the associated pulseĝ(t). However, typically over a timespan of minutes, the motional-mode frequencies ω p will drift with typical excursions of ∆ω p /(2π) ≈ 1 kHz. If we now useĝ(t), determined on the basis of the original mode frequencies ω p , in the situation of the drifted modes, ω p + ∆ω p , the set of equations (S1) are no longer fulfilled, resulting in a reduction of the fidelity of the xx gate. A simple estimate for the infidelity increase due to the now non-zero α's in (S1) is presented in [30]. According to [30], at zero temperature of the motional-mode phonons, the infidelity,F , is approximately given bŷ This suggests stabilizing the fidelity of the quantum computer against mode drifts and fluctuations by requiring that α ip be stationary up to nth order with respect to variations in ω p . This is easily accomplished by adding the following set of equations to the set of equations (S1): Because of the presence of the factor t k in the integrand of (S61), we call this extension of our linear approach the moments approach. Adding the moments equations (S61) to the set (S1) does not change the linearity of our method. The same techniques can be applied in solving this extended system of linear equations as was described in the main text.

S9. Demodulation of pulses
The optimal pulse functionsĝ(t) are simultaneously amplitude-, frequency-, and phase-modulated pulses. In this section we show how to demodulate the pulseĝ(t), i.e., how to separateĝ(t) into its amplitude function Ω(t) and its detuning function µ(t).
The first step of our demodulation procedure is to find the zeros ζ j ofĝ(t). This is numerically unproblematic, since the detuning function µ(t) is bounded away from zero, which means that degeneracies of nontrivial zeros (ζ j > 0) do not occur. In addition, in numerous simulation runs, we observed that the envelope function Ω(t) was always bounded away from zero. Therefore, in order not to complicate the discussion, we may also assume that Ω(t) does not have any zeros. Thus, all the zeros inĝ(t) are caused by zeros of sin[ψ(t)], i.e., ψ(ζ j ) is a multiple of π. Since no degenerate zeros occur, we have even more, namely where N z is the total number of zeros ofĝ(t), including the zero ζ 0 = 0 at t = 0 and ζ Nz−1 = τ at t = τ . We now approximate the detuning function µ(t) as a constant between zeros ofĝ(t), i.e., With (S2) and (S62) this entails As an example of frequency demodulation, the main text Fig. 3b shows the result of the detuning function µ(t) for the pulse shown in the main text Fig. 3a. We see that µ(t) hovers about the middle motional mode, staying away from the strongly heating mode with the highest motional frequency. Since g(t), for the example shown in Fig. 3a, has a dense set of N z = 387 zeros, µ(t) approximated by as many piece-wise constant plateaus appears as a smooth function on the scale of Fig. 3b. We now turn to extracting the pulse envelope function Ω(t) fromĝ(t). Differentiating (5) and evaluating the result at the zeros ζ j ofĝ(t) yieldŝ where we used (S4) and (S62). This equation can be solved for Ω(ζ j ) with the result where we inserted the factor σ = −ĝ (ζ 1 )/|ĝ (ζ 1 )|, which ensures that Ω(t) is "right-side up", i.e., if it does not change sign, Ω(t) > 0 for all t. Sinceĝ(t), according to (S27), is represented by a Fourier series, it is trivial to obtainĝ nÂ n cos 2πn t τ (S67) and thusĝ (−) (ζ j ). The values of the detuning function µ(ζ j ) may be obtained in several ways. We may use spline interpolation of the data set of values µ j as defined in (S64), or, as we found, with sufficient accuracy, simply use (i) µ(ζ j ) = µ j , (ii) µ(ζ j ) = µ j+1 , or (iii) µ(ζ j ) = (µ j+1 + µ j )/2. We used method (i) to obtain the pulse envelope function Ω(t)/(2π) (heavy orange line in the main text Fig. 3a) of the pulseĝ (−) (t), shown as the thin green line in the main text Fig. 3a. The main text Fig. 3a shows that our amplitude demodulation technique presented above works very well and accurately extracts the envelope function. At this point we may wonder how well the exact pulsê g (−) (t) is approximated by the pulseg (−) (t), i.e., the pulse reconstructed via (S8) from the amplitude and detuning functions obtained by demodulatingĝ (−) (t) according to the above procedures. Therefore, to get a first impression of the accuracy of our pulse demodulation method, we compute where, for j = 1, 2, . . . , N z − 1, and Notice that Ω j in (S70) does not contain the factor σ as in (S66), since this time we do not need the "right-side up" pulse, but the pulse that has the same sign of the amplitude asĝ (−) (t). For the example shown in the main text Fig. 3a, we obtain ∆g 2 = 1.3 × 10 −5 . Hence, the pulseg (−) (t) reconstructed from the demodulated pulsê g (−) (t) is sufficiently accurate to guarantee high-fidelity gates.

S10. Stabilization against pulse-timing errors
Once τ and the mode frequencies ω p are given, our method determines the amplitudes A n , which are then fixed when implementing the pulse g(t). However, the experimental clock may run fast or slow, which results in pulse-timing errors. To stabilize against pulse-timing errors of this nature, we require, for all l = 0, . . . , L and p = 1, . . . , N : Since L = 0, according to (3) of the main text, is already satisfied, (S72) represents LN linear equations that may be added to the N (K + 1) linear equation of motionalmode stabilization. Apart from providing explicit formulas for stabilization against pulse-timing errors, this also provides a template for stabilizing against other types of errors and parameter fluctuations, and shows that it is straightforward to extend our method by any number of such linear constraints, as long as the dimension of the available null-space is not exhausted.

S11. Implementation details
In this section, we present the pulse-level implementation details. Figure S2 shows the amplitude and frequency profiles of pulses, represented according to equation (5) in the main text, with moment stabilization orders K = 1, 4, 7, used to implement the xx gates on our 5-qubit, 7-ion TIQIP. The motional-mode frequencies for each of the experiments are reported in Tables S4 and S5. The two-qubit gates were performed on qubits 3 and 4, with indexing starting at 0.
We note that the even-parity population, used in Fig. 5 to demonstrate stabilization against mode-frequency drift, is not a complete characterization of the resulting quantum gate in terms of fidelity. However, the evenparity population is the most relevant metric to use when evaluating the general experimental performance of our methods. The formalism for stabilizing phase space closure at the end of the gate with respect to gate frequency offset, relevant to Fig. 5, aims to minimize |α p | in (3), which is most closely related to the odd-parity population, i.e., 1 minus the even-parity population. Moreover, fidelity contains within itself other sources of error accumulated during the parity contrast measurement (e.g., intensity noise differential, phase stability, laser-beam steering, overlap errors of the lasers with the positions of the trapped ions, etc.), which would mask the effect of stabilization against motional-mode drift. Thus, we decided to present the even-parity population as the most relevant quantity of interest to our work.    Fig. S2.

S12. Fixed-detuning step pulses
Possibly the most widely studied type of fixed-detuning pulses are segmented step pulses [12]. According to this method, the detuning function µ(t) is set to a constant, i.e., µ(t) = µ 0 = const for 0 ≤ t ≤ τ , and the pulse interval [0, τ ] is broken up into N seg > P equi-spaced intervals [t j−1 , t j ], t 0 = 0, t j = j∆t, ∆t = τ /N seg , j = 1, . . . , N seg , in which the pulse amplitude is set to a constant, i.e., For this type of pulses our methods are directly applicable with only two minor modifications. (i) For given µ 0 , we choose the gate length τ such that J = µ 0 τ /π is an integer. This way, still requiring that Ω(t) is an even function with respect to τ /2, we obtain even-or oddparity pulses,ĝ (±) (t) = Ω(t) sin(µ 0 t), for J odd or even, respectively. Since J needs to be an integer to obtain the desired symmetry classes, τ can take only discrete values. However, since for quantum computer hardware of practical interest (for instance, Yb-ion quantum computers [9]), the detuning µ 0 is such that J is a large integer (of the order of 1000), the discretization of τ is of no consequence in practice. (ii) The second modification concerns the computation of the matrix M defined in (3). For step pulses, we let M pn → M pj , where, including both negative-and positive-parity pulses, we have Defining A = (Ω 1 , Ω 2 , . . . , Ω Nseg ), all the procedures outlined in Section S4 can now be applied to construct step pulses.   Figure S3 shows an example of a negative-parity step pulse, generated for the same set of motional-mode frequencies and Lamb-Dicke parameters as in the main text Fig. 3a. Although Fig. S3 shows the negative-parity pulse with the lowest peak-power requirement that we found in the detuning interval from µ 0 /(2π) = 2.2 MHz to µ 0 /(2π) = 2.6 MHz, we see that this pulse is about 10% higher in peak power than the pulse shown in the main text Fig. 3a. This is expected, since fixed-detuning pulses lack the additional degrees of freedom that are associated with being able to modulate the detuning. However, in analogy to the main text Fig. 3a, we see that the pulse tends to be relatively flat in amplitude, a feature we observed in all power-optimized pulses we generated in the course of numerous simulations.
There are several reasons why step pulses should be replaced with AMFM pulses. In our opinion the two leading reasons are (i) amplitude-, frequency-, and phasemodulated pulses have lower power requirement as seen when comparing the main text Fig. 3a and Fig. S3 and (ii) in contrast to the sharp transitions in power levels characteristic for step pulses, amplitude-, frequency-, and phase-modulated pulses have a smooth pulse envelope, which eliminates ringing and the Gibbs phenomenon [48] that accompanies sudden changes in power levels.
In contrast to the straightforward construction of our amplitude-, frequency-, and phase-modulated pulses, finding the optimal pulse for step pulses requires a search in the 4D parameter space consisting of the number of segments, N seg , the detuning µ 0 , the integer J, and the parity (±) of the pulse. While N seg is discrete, and we found that good convergence is already achieved with relatively few segments, in terms of parity there are only two cases to check, and, if τ is pre-specified to a certain value, say, τ = 300 µs ± 1 µs, the range of J that falls into this interval is not large, and, moreover, J is discrete, searching for the optimal detuning µ 0 requires considerable computational overhead that is avoided using our "single-shot" AMFM approach.
Concluding this subsection, we can say that our new linear algorithm is certainly general enough to encom-pass the important class of step pulses. Thus, if such pulses are required to run, e.g., existing quantum computers with existing controller hardware which requires step pulses as input, our method can be used to generate these pulses efficiently and directly.

S13. Efficient arbitrary simultaneously entangling gates
In this section, we show how to use our method in conjunction with the Efficient Arbitrary Simultaneously Entangling (EASE) gate protocol detailed in [29]. To see how this may be achieved, the only thing that is required is to show that the equations to be solved are isomorphic. In particular, the null-space condition (S23) is of the same structure as Eq. (2) of [29] and the degree-ofentanglement condition (S33) is of the same structure as Eq. (3) of [29], which fully specify the problem of solving for the EASE-gate pulse shapes. The rest of the EASEgate protocol follows immediately. The resulting pulse shapes can implement up to N (N − 1)/2 xx gates simultaneously in a short time for a given power budget.

S14. Sensitivity of the degree of entanglement
In this section we now explore the effects of motionalmode drifts on the gate angle χ. For N = 5 ions, two cases are investigated. (i) All modes ω p drift in unison from 0 Hz to +2π × 500 Hz and (ii) individual modes drift independently. For case (ii), we simulated a case in which, chosen randomly, and with random signs of the drift direction, ω 1 drifts from 0 Hz to +2π × 500 Hz, ω 2 drifts from 0 Hz to −2π × 400 Hz, ω 3 drifts from 0 Hz to +2π×300 Hz, ω 4 drifts from 0 Hz to −2π×500 Hz, and ω 5 drifts from 0 Hz to +2π × 400 Hz. Figure S4 shows that although all drift amplitudes are substantially smaller than 1 kHz, the effect on the gate angle χ is substantial.
The strong sensitivity of χ with respect to drifts in ω p is due to the amplification effect of the relatively long pulse duration. In order to compute χ, we have to evaluate the double integral (S6). Under the integral we have the term sin[ω p (t 2 −t 1 )], and if we replace ω p by ω p +∆ω p , then the sin[ω p (t 2 − t 1 )] term becomes, in linear order, sin[ω p (t 2 − t 1 )] + cos[ω p (t 2 − t 1 )]∆ω p (t 2 − t 1 ). Now, while |∆ω p | is at most 2π × 500 Hz, which looks small, and indicates that we might be able to neglect the second term, when we multiply the second term with 300 µs, which is the maximum of t 2 − t 1 , we get 2π × 0.0005 MHz × 300 µs = 0.94, which is large. In fact, this term is so large that the linearization approximation breaks down. Therefore, the pulse length is the amplification mechanism and explains the strong sensitivity of χ to relatively small drifts in ω p . It also underpins the observed sensitivity (see Fig. S4) with a detailed qualitative analytical understanding.
In order to counteract drifts in χ, we suggest to monitor the value of χ continuously and readjust the laser power that drives the xx gate if χ drifts away. This is a valid correction mechanism since the set of equations (S1) depends only on the shape of the pulse, but not on the pulse amplitude. Therefore, without compromising the validity of (S1), the power can be continuously adjusted to keep χ within tolerable bounds. Of course, it may be difficult in practice to continuously monitor and readjust χ. Nevertheless, at least in principle, this is a possible correction and stabilization mechanism. In analogy to our moments approach for active stabilization of the α conditions (S1), it is also possible to encode active stabilization of χ in the pulse shape itself. S15. Single-pulse active stabilization of the degree of entanglement: Projection method Ideally, to actively stabilize χ ij against ω p fluctuations, integrated in the pulse-shape construction, we should require where K χ is the maximal desired degree of χ stabilization. Since all pulse shapes, regardless of their maximal degree of stabilization K χ , need to satisfy both the decoupling conditions (S1) between the motional modes and the computational states and the degree-of-entanglement condition (S7) (where "π/8" may be replaced by the actual desired degree of entanglement), we may write where and To understand the consequences of (S75) [(S76), respectively], we spectrally decompose R (k) p according to where λ ν,p have the same sign. However, we can prove analytically (the proof is lengthy and not shown here) and confirmed numerically, that, for instance, R (k=1) p is a definite matrix for all p, i.e., the eigenvalues of R (k=1) p are all non-zero and have the same sign, which makes it impossible to satisfy (S81) for k = 1. We did, however, notice that only a few of the eigenvalues of R (k=1) p are particularly large in absolute magnitude, which may be the ultimate reason for the strong sensitivity of χ ij in linear order (see Fig. S4). This observation suggests a strategy for actively stabilizing χ ij against ω p -fluctuations: Projecting those components of the spectra of R (k) p out of the null-space of M (which can be assumed to already include stabilization of (S1) against ω p -fluctuations) that correspond to the eigenvalues with the largest absolute values. If we project out L such components from each of the R (k) p matrices, this leaves us with a null space of N 0 = N 0 −P K χ L dimensions that now, to a large degree, actively stabilizes χ ij against ω p -fluctuations. Following this projection step, we now use the techniques presented in S4, applied to the reduced null space of N 0 dimensions, to satisfy the degree-of-entanglement condition with the smallest possible average power.
To illustrate this technique, and focusing on the case of uniform drift of the motional modes from 0 to 500 Hz as defined above, we present in Fig. S4a the result of projecting 1,2,3,4, and 12 states from the null space that correspond to the eigenvalues with largest absolute values of R (1) p , p = 1, . . . , 5. We see that already for a single projected state we achieve noticeable stabilization that improves further for 2, 3, and 4 projected states. This improvement continues if more states are projected, reaching an optimum ("sweet spot") for 12 projected states. This is illustrated in Fig. S4b, which shows the normalized χ for 11, 12, and 13 projected states. Therefore, while we found that projecting relatively few states always results in improved active stabilization, "overprojection" should be avoided, since it is both costly in power and does not improve χ stabilization any further. In fact, as expected, active χ stabilization, in analogy with stabilizing α, requires increased levels of power. For example, for the case shown in Fig. S4, the projection of 1, 2, 3, 4, and 12 states requires power levels of 1.5, 1.7, 2.0, 2.1, and 3.3 with respect to the power level without projection. But we also see that projection is relatively inexpensive compared with the significant amount of stabilization gained.
The projection technique works for all orders k ≥ 1 of R (k) p . However, not much is gained by continuing the projection beyond k = 1. The reason is the following. In our example, the best result, obtained by projecting the first 12 states, brings down the variation in the relative χ from 35% to just 1.5%. These 1.5%, however, are mostly due to the residual slope of the first-order stabilization, so that second-order stabilization would not contribute much, other than computational effort and power expended. We see this in the following way. The slope of the first-order stabilization for 12 projected states at 0 motional-mode drift is 3×10 −5 /Hz. Therefore, at 500 Hz motional-mode drift, the variation in the relative value of χ is 0.015, i.e. 1.5%. This is exactly the amount we read off in Fig. S4b for the case of 12 projected states. Therefore, the residual variation is mostly due to the first order, and stabilizing the second order will have a negligible effect. Nevertheless, as shown in Fig. S4a, and given an unstabilized variation of 35%, active first-order stabilization via projection, which brings this variation down to about 1.5% (see Fig. S4b), is already significant for the stabilization of quantum-computer operation. While, as discussed above, it is not possible to implement the moments strategy (S75) with a single pulse, working with two different pulses, each directed at a different one of the two ions participating in the gate, (S75) can in fact be realized and results in the moments methods described in the following two sections. While it is impossible to satisfy (S75) using identical pulses directed at both ions i and j, the condition (S75) can be satisfied using different pulses directed at ions i and j. Let g (i) (t) and g (j) (t) be the pulse functions directed at ions i and j, respectively, and letˆ F andˆ G be the null-space expansion amplitudes of g (i) (t) and g (j) (t), respectively. Then, condition (4) in the main text implies which has to be solved under the conditions [see (S75)] where F and G are any un-normalized versions ofˆ F andˆ G. Unlike in the single-pulse case discussed in S15, where, due to the observed definiteness of R p , it was impossible to satisfy (S75), working with two different pulses, F and G, (S83) can be satisfied as soon as F and G are orthogonal to each other with respect to R (k) p . An explicit solution of (S82) and (S83) can be constructed in the following way. We start by choosing G to be the power-optimal pulse as computed in the single-pulse case described in S4. Then, we construct the space P, which is spanned by the vectors v p G, p = 1, . . . , P , k = 1, . . . , K χ . We also define the space Q, which is the orthogonal complement of P with respect to the null space. With these definitions, taking F out of Q, we obtain the most power-optimal solution for F by choosing F =Q G, whereQ projects into the Q space. At this point the normalizations of F and G are still two free parameters. We use the first scaling freedom to obtain |ˆ F | = |ˆ G| (symmetric pulse power) and use the second scaling freedom to satisfy (S82).
For the same case of uni-directional ω p drift as used in Fig. S4, the performance of active two-pulse χ stabilization is illustrated in Fig. S5 and can be compared with the performance of the single-pulse projection method illustrated in Fig. S4. We see that while the projection method at best cuts the relative χ drift error down to 1.5% at 500 Hz drift, the relative error in the case of K χ = 5 stabilization orders is smaller than 1 permille even at 500 Hz drift. In addition, while in the projection method the stabilization is of linear order around zero motional-mode drift, the moments method produces favorable non-linear stabilization of order K χ + 1 around zero motional-mode drift (see Fig. S5). The power of the moments method is best appreciated in Fig. S5 by comparing the behavior of the unstabilized χ (the nearly vertical, purple line in Fig. S5 around zero motionalmode drift; K χ = 0) with the shape of χ as a function of motional-mode drift for different K χ ≥ 1. Even for small K χ ≥ 1, substantial stabilization is observed.
In Fig. S6 we show a summary of various aspects of the two-pulse active χ stabilization method. Figure S6a shows the infidelity of the two-pulse active-stabilization method as a function of mode-frequency drift. We see that infidelities 10 −3 can be achieved for motionalmode frequency drifts in the range ±200 Hz for all K χ ≥ 1 and and an infidelity of 10 −4 can be achieved over a motional-mode frequency drift in the range ±500 Hz for K χ = 5. Figure S6b shows a summary of frequency widths of χ stabilization for infidelity cut-offs of 10 −3 , 10 −5 , and 10 −7 , respectively. On average we observe a linear increase of frequency width with the χ stabilization order. While these results are promising, they do come with a price. Figure S6c shows the power requirement as a function of stabilization order K χ . In constrast with the linear power scaling of α stabilization shown in Fig. 4c,  Fig. S6c indicates that the required pulse power in the case of χ stabilization increases exponentially with the stabilization order. This, however, should not discourage us from using the moments method for active two-pulse stabilization of χ, since, as shown in Figs. S5 and S6, significant χ stabilization is already achieved for relatively small K χ . Figure S6d shows the scaling of frequency width for several stabilization orders K χ as a function of gate duration τ . Similar to the results shown in Fig. 4d, we observe that the frequency width is inversely proportional to the gate duration.
FIG. S6: Stabilization of the control pulses. a. Infidelity (see SI section S8 for detail) as a function of the motional-mode frequency drift ∆f . All mode frequencies were drifted according to ωp → ωp + 2π∆f . b. The width of the infidelity curves in a for various error tolerances = 10 −3 , 10 −5 , and 10 −7 , as a function of the highest moment Kχ of stabilization. c. The maximal power requirement maxt |gK (t)| of the control pulses as a function of the highest moment Kχ of stabilization. The power requirement suggests an exponential scaling of power in the order of stabilization Kχ. d. Width of the infidelity curves for various different orders of stabilization Kχ = 0, 2, and 4, as a function of the gate duration τ for a fixed error tolerance level = 10 −3 . The data suggests ∼ 1/τ scaling of the width. S17. Two-pulse active stabilization of the degree of entanglement: Hybrid method Combining the most advantageous features of both the moments and projection methods, we arrive at the hybrid method. In this method, for desired stabilization order K χ , we first construct the pulseˆ G as described in section S15, i.e., as a single pulse in a null space of dimension N 0 = N 0 − P L, where we projected out those L eigenvectors from each matrix R (1) p , p = 1, . . . , P , that correspond to the L eigenvalues with the largest absolute values. We then proceed with the construction of the pulseˆ F as described in section S16. Thus, even before construction of the pulseˆ F , we eliminate those null-space components from the pulseˆ G that potentially produce the most sen-sitivity in χ. As documented in Fig. S7, and compared with the single-pulse projection method (see S15) and the bare two-pulse moments method (see S16), the hybrid method yields the best results in terms of active χ stabilization. Figure S7a, for the same case of uniform modefrequency drift as used in Figs. S4, S5, and S6, shows the result of χ stabilization for K χ = 1 and L vectors with the largest absolute values of their eigenvalues projected out. We see that already for the case K χ = 1, compared with the bare moment method for K χ = 1 (see Fig. S5), the gain in stabilization is substantial. This observation is important since, because of the exponential power cost of the two-pulse method in K χ , it is advantageous, for a given stabilization target, to stabilize with the smallest possible K χ . Figure S7b shows that with a slightly larger K χ [K χ = 3 in Fig. S7b], a substantial broadening shows that compared with a. larger Kχ results in a broader stabilization region, which is widened even further by projection.
of the stabilization region can be achieved. Figure S7b also shows that by projecting out L = 10 vectors, χ stabilization on the level of 10 −5 can be achieved for 500 Hz motional-mode drift.

S18. Broadband sequence
A host of compensation pulse sequences that mitigate the errors in the single-qubit gate are known (see [51] and the references therein) and similar techniques can indeed be used to mitigate the errors in χ ij that arise from, e.g., relative offsets in η p i or g(t). We refer interested readers to [49,50], where a variety of broadband behaviors have been explored. Below, we show an example for completeness.

S19. Direct implementation of Fourier-basis pulse function
According to [12], in the Lamb-Dicke regime, the interaction Hamiltonian for the ion-chain system, subjected to a dual-tone, symmetric blue-and red-sideband beam with detuning ±µ, in the x basis is where a p and a † p are the annihilation and creation operators of the pth motional mode, respectively. Consider now a multi-tone beam with amplitudes Ω i,n (t) and detuning frequencies ±µ n , where n = 1, 2, .., N A . This results in the Hamiltonian where N A is chosen sufficiently large to achieve convergence. Inserting (S92) in (S91), we obtain (S94) which induces the system evolution over the gate time τ described by as shown in [12] using Magnus' formula. Inserting (S94) in (S95), together with (S92) and (S93), we obtain, up to a global phase, where α * ip denotes its complex conjugate, and (S98) Comparing (S97) and (S98) with (3) and (4), respectively, together with (S93), and assuming g i (t) = g j (t) = g(t) = N A n=1 A n sin(2πnt/τ ), as we did in the main text, we see that Ω i,n (t) = A n and µ n = 2πn/τ implements the pulse function that implements the desired xx gate.