Robust Oscillator-Mediated Phase Gates Driven by Low-Intensity Pulses

Robust qubit-qubit interactions mediated by bosonic modes are central to many quantum technologies. Existing proposals combining fast oscillator-mediated gates with dynamical decoupling require strong pulses or fast control over the qubit-boson coupling. Here, we present a method based on dynamical decoupling techniques that leads to faster-than-dispersive entanglement gates with low-intensity pulses. Our method is general, i.e., it is applicable to any quantum platform that has qubits interacting with bosonic mediators via longitudinal coupling. Moreover, the protocol provides robustness to fluctuations in qubit frequencies and control fields, while also being resistant to common errors such as frequency shifts and heating in the mediator as well as crosstalk effects. We illustrate our method with an implementation for trapped ions coupled via magnetic field gradients. With detailed numerical simulations, we show that entanglement gates with infidelities of $10^{-3}$ or $10^{-4}$ are possible with current or near-future experimental setups, respectively.

In this regard, Manovitz et al. [59] have experimentally shown that the MS gate can be combined with pulse sequences given the ability to tune and turn on-and-off the qubit-boson coupling as many times as the number of DD pulses introduced.Another possibility explored theoretically is to combine an always-on qubit-boson coupling with strong π pulses [51,60,61].Although this is possible in certain trapped-ion architectures, turning on-and-off the qubit-boson coupling may be not practical in other platforms.On the other hand, the use of strong π pulses is experimentally challenging since high-power controls are needed, while these induce crosstalk and hinder the applicability in multimode scenarios.
In this article, we design a DD sequence with lowintensity π pulses -named TQXY16-that achieves fasterthan-dispersive entangling gates using static (i.e.non-tunable) qubit-oscillator coupling.Importantly, our gates decouple from dephasing, pulse imperfections, and unwanted finitepulse effects, leading to high-fidelity.Furthermore, we demonstrate the versatility of our protocol by incorporating techniques that lead to additional resilience to decoherence on the bosonic mediator and potential crosstalk effects.Although our method is general, we exemplify its performance in radiofrequency controlled trapped ions demonstrating infidelities within the 10 −3 threshold at state-of-the-art experimental conditions, and of 10 −4 in near-future setups.

A. Gate with instantaneous pulses
We consider a system that comprises two qubits and a bosonic mode -with frequencies ω 1 , ω 2 and νcoupled via longitudinal coupling [2][3][4][5][6][7] (here, and throughout the paper, H is H/ℏ, meaning all Hamiltonians are given in units of angular frequency), Here, a † (a) is the creation (annihilation) operator of the bosonic mode, S µ ≡ σ µ 1 + σ µ 2 with µ ∈ x, y, z are collective qubit operators, and ην is the coupling strength.Also, note that H 0 is written in a rotating frame with respect to (w.r.t) the qubit free-energy Hamiltonian H q = µ ω µ σ z µ /2.We assume the usual experimental scenario η ≪ 1, thus we stay arXiv:2209.14817v2[quant-ph] 30 May 2023 away from other paradigms that require stronger qubit-boson couplings [62][63][64][65][66]. H 0 contains no driving fields, while in our method we drive the qubits for two main reasons: (i) Accelerate the gate by making the qubits rotate at a frequency close to the bosonic frequency ν, and (ii) Protection of the gate from qubit noise of the form ϵ j (t)σ z j /2 leading to dephasing.When driving the qubits, H 0 is completed with the term H d (t) = µ=x,y Ω µ (t)S µ /2.In an interaction picture w.r.t H d + νa † a we get where µ=x,y,z f µ (t being the time-ordered propagator.See supplementary note 1 for additional details. If driving fields are delivered as instantaneous π pulses (note this requires Ω x,y ≫ ν during the application of the pulse) spaced τ/2 apart, f x,y (t) can be neglected and f z (t) = 1(−1) if the number of applied pulses is even (odd), see the grey solid line in Fig. 1(a).For the moment we consider instantaneous pulses, while later we treat the realistic case of non-instantaneous ones.As π pulses are applied periodically, f z (t) takes the form of a function with period τ such that f z (t) = ∞ n=1 f n cos (nωt), where ω = 2π/τ and f n = 2 τ τ 0 dt ′ f z (t ′ ) cos (nωt ′ ).Hence, under the assumption of instantaneous pulses we get whilst setting an interpulse spacing τ/2 = τ k /2 such that ω = ω k ≈ ν/k leads to a resonant qubit-boson interaction via the kth harmonic (from now on τ → τ k and ω → ω k , where the subscript k refers to the kth harmonic).As η ≪ 1, the terms in Eq. ( 3) that rotate with frequencies ±|ν − nω k | (where n k) and ±|ν + nω k | can be substituted, using the rotating-wave approximation, by their second-order contribution (here, and in the rest of the paper, second-order stands for second order in η) leading to where ξ k = ν − kω k is the detuning w.r.t. the kth harmonic and ) is an effective spin-spin coupling constant that contains contributions from all harmonics.Note that, as η ≪ 1 contributions of higher order in η can be neglected.See supplementary note 2 for additional details.The propagator associated to Hamiltonian (3) is where α(t) = −iην t 0 dt ′ f z (t ′ )e iνt ′ ≈ −ην f k /(2ξ k )(e iξ k t − 1) and 1 9 e T 9 A 4 w F B T D k q 1 X S f S n Q S k Z p T j J O / F C i O g I x h g 2 9 B 0 n + o k 0 x c m 9 o l R e n Z f S F O h t q f q 7 4 k E A q X G g W 8 6 A 9 B D N e + l 4 n 9 e O 9 b 9 y 0 7 C w i j W G N L Z o n 7 M b S 3 s N A + 7 x y R S z c e G A J X M 3 G r T I U i g 2 q S W N y G 4 8 y 8 v k k a l 7 J 6 X K 7 d n x e p V F k e O H J F j U i I u u S B V c k N q p E 4 o e S T P 5 J W 8 W U / W i / V u f c x a l 6 x s 5 o D 8 g f X 5 A 3 Z w l 3 c = < / l a t e x i t >

Re(↵)
< l a t e x i t s h a 1 _ b a s e 6 4 = " 6 s Z a B 4 X y G r e o c h R V 6 S c T 8 R k A b J g = " > A A A C A n i c b V B N S 8 N A E N 3 4 W e t X 1 J N 4 C R a h X k p S R D 0 W v e i t g v 2 A p p T J d Entangling gates with instantaneous π pulses: (a) Rabi frequency Ω(t) and modulation function f z (t) during a period τ k for k = 5.For comparison, we plot cos (kω k t) in green.(b) Phasespace trajectory of α(t) during the application of a pulse sequence with k = 1 and k = 5 in blue and green, respectively.
where C is the phase-space trajectory followed by α(t).Note that, if the gate time is chosen as t g = 2π/|ξ k |, α(t g ) ≈ 0 at the end of the gate, making the gate insensitive to the bosonic state.To satisfy condition θ(t g ) = π/8, we choose After a time t g the propagator U(t) approximates to exp (i π 8 S 2 z ).For two qubits, this is equivalent (up to a global qubit rotation) to the CPHASE gate, and transforms the state | ++⟩ into the Bell state It is noteworthy that the choice of ξ k (thus τ k ) is, in general, not trivial, as both f k and J k depend on τ k .However, in the cases discussed here this dependance does not hold, making the choice of ξ k direct.See supplementary note 3 for analytic expressions for f k and J k .
For instantaneous π pulses one finds f k = f ins k = 4 kπ sin ( kπ 2 ).Thus, if resonance is achieved via a low harmonic, e.g.k = 1, the gate time is t k=1 g ≈ π 2 /4ην, a factor π/2 longer than the original MS gate.On the other hand, for sufficiently large harmonics the gate is mostly governed by the dispersive term 1 2 η 2 νJ k→∞ S 2 z ≈ η 2 νS 2 z in Eq. ( 4), leading to t d g = π/8η 2 ν.We define faster-than-dispersive gates as those that satisfy t k g /t d g < 1.For example, in the case η = 0.01 and k = 1 we find a faster-than-dispersive gate with t k=1 g /t d g ≈ 1/16.This is, the gate is 16 times faster than the dispersive one.In Fig. 1(b), we show α(t) for η = 0.03, with k = 1 and 5.
It is noteworthy that the analysis conducted above is valid for N q qubits homogeneously coupled to the same bosonic mode, i.e. S µ → N q j=1 σ µ j .For two qubits with inhomogeneous coupling, i.e. S z → S z = (η 1 σ z 1 + η 2 σ z 2 )/η where η j ≪ 1, the method also yields to a CPHASE gate, however, in this case, the correct expression for the detuning ξ k is that in which every η is substituted by √ η 1 η 2 .For larger η, the rotating-wave approximation is not justified and terms neglected from Eq. (3) to Eq. (4) will lead to significant residual qubit-boson entanglement at the end of the gate.See supplementary note 4 for additional details.

B. Gate with low-intensity pulses
In what follows, we discuss the realistic case of noninstantaneous pulses.For standard top-hat pulses the Fourier coefficient f k that quantifies the strength of the qubit-boson interaction reads (see supplementary note 3 for the derivation) Notice that for low-intensity pulses -defined as those holding Ω < νthe value of f k decays with (Ω/ν) 2 .As a result, achieving faster-than-dispersive gates is no longer possible.Note that f k directly relates to the non-dispersive contribution in θ(t), and, through condition θ(t g ) = π/8, to the gate time t g .
To solve this problem and optimize the strength of the qubit-boson interaction, we propose to modulate the Rabi frequency during the execution of each π pulse.Specifically, we pose the following ansatz for f z (t) where t π is the π pulse duration, and t i and t m = t i + t π /2 are the initial and central points of the pulse.Note that the Rabi frequency is then given by Ω(t . For the envelope function β(t), we propose where t r = t m + bt π and t l = t m − bt π .The free parameters b and c serve to control the width of the envelope function β(t), while d is proportional to its amplitude.Suitable values for b, c and d for the first harmonics are shown in table II (Methods).
From now on, we assume t π = τ k /2, i.e., the pulse extends over a whole period τ k /2, leading to solutions with the lowest intensities.As a result of our pulse design with suitable b and c, the value for the Fourier coefficient f k is given by f m k = −4d/πk where |d| can take values from 0 to |d max | > 1. See supplementary note 3 for the derivation.Since now f k depends on d, this serves to control the strength of the interaction, thus the duration of the gate t g .Also, d relates to the amplitude of the pulse, thus to the maximum value of the Rabi frequency Ω pp .Typically, we look for large values of d, bounded by d max or by the experimentally available Ω pp .Now we describe the recipe to design faster-than-dispersive gates using low-intensity pulses.First we choose a value for the harmonic k.Larger k allow for lower pulse intensities at the price of longer gates.Second, we use Eq. ( 8) to generate the modulation function f z (t) and the Rabi frequency Ω(t) for different values of d, and calculate both the gate time t g and Ω pp = max[|Ω(t)|].We note that the obtained Ω(t) can lead to pulses along arbitrary axes (e.g.X or Y).In particular, for reasons described later, we target gates formed by concatenating blocks of 16 pulses.For that, the gate time t g must be 8Nτ k , where N is an integer number.This translates into the condition (ν < l a t e x i t s h a 1 _ b a s e 6 4 = " q A Z 2 O K N 6 W e h d g J Z 4 < l a t e x i t s h a 1 _ b a s e 6 4 = " J x 2 1 c 4 D i h 7 i R 3 h P 5 2 8 S x 3 z q < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 / E p D G J / 5 L h 1 s 3    < l a t e x i t s h a _ b a s e = " J x c D i h i R h P S x z q  and . The final step is to select the values of d for which this last condition is satisfied.As a result, we obtain all possible gates within the harmonic k as well as the corresponding values for Ω pp .Figure 2 (b) shows values of t g and Ω pp obtained following the previous prescription for η = 0.005 and k = 7, 9. Notice that there are plenty of solutions giving faster-than-dispersive gates, i.e. t g /t d g < 1, using low-intensity pulses with values of Ω pp well below the frequency ν.
As an example, we choose two solutions within the 9th harmonic, where τ k extends over approximately 9 oscillator periods.In Fig. 2 (c) the shapes of Ω(t) and f z (t) are displayed for cases N = 5 and 10.Notice that Ω(t) achieves a larger amplitude when N = 5.As a consequence, it generates a faster gate.This is shown in Fig. 2 (d), where the two-qubit gate phase θ(t) related to the N = 5 gate reaches the target value π/8 faster than the N = 10 gate or the dispersive gate.
The reason for choosing the gate time as an integer multiple of 8τ k has to do with an efficient decoupling from finite-pulse effects produced by the terms f x,y (t) neglected in Eq. ( 3).In the same way, the XY8≡XYXYYXYX pulse structure assures cancelation of σ z type noise, as well as of Rabi frequency fluctuations.In Fig. 2(a) the Rabi frequency is plotted (Ω x (t) when "X"; Ω y (t) when "Y") for an XY8 block.
To understand the elimination of finite-pulse effects, we calculate the second-order Hamiltonian of Eq. ( 2) after a XYXY block leading to (see supplementary note 5 for the derivation) If our gate contains only XYXY blocks, H XYXY adds to Hamiltonian (4) spoiling a high-fidelity performance.To overcome this problem, we use a two-step strategy.Firstly, we concatenate XYXY and YXYX blocks (which form a XY8 block) such that the term B k (a † a)S z gets refocused.Note that in the presence of bosonic decoherence, this term will induce qubit dephasing.Secondly, we cancel the remaining term J ⊥ k (S 2 x + S 2 y ) by driving the two qubits with opposite phases every second XY8 block.This is, when rotating by an angle π the phase of the second qubit's driving This changes the sign of the σ u 1 σ µ 2 terms with u, µ ∈ {x, y}, leading to refocusing of the term J ⊥ k (S 2 x + S 2 y ) after every pair of XY8 blocks.Note that the second step requires the ability to address each qubit individually, and assumes that N q = 2, i.e. S µ is given by the sum of two qubit operators.In the absence of individual addressing, one can incorporate the term J ⊥ k (S 2 x + S 2 y ) into the gate, but then the operation applied is not equivalent to the CPHASE gate.For a discussion regarding this alternative gate, as well as the extension to the multiqubit case, see supplementary note 6.
Summarizing, our two-qubit gates are generated by nesting TQXY16≡XY8 (+) XY8 (−) blocks, where TQXY16 stands for "two-qubit" XY16, while XY8 (±) imply qubits driven in phase or in anti-phase as discussed in the previous paragraph, while, importantly, each π pulse is implemented according to the designs for f z (t) and β(t) presented in Eqs.(8,9).

C. Trapped-ion implementation & Numerical results
We benchmark our method by simulating its performance in a pair of trapped ions in a static magnetic field gradient [7].In this scenario, qubit frequencies ω µ take values around (2π) × 10 GHz, ν = (2π) × 220 kHz is the frequency of the centre-of-mass vibrational mode, η = γ e g B /8ν √ ℏ/Mν is an effective Lamb-Dicke factor where γ e = (2π) × 2.8 MHz/G, g B is the magnetic field gradient, and M is the ion mass.The two-ion system has a second vibrational mode "b" with its corresponding qubit-boson coupling.Thus, Hamiltonian (1) is replaced by H 0 + H 2M , where z .The addition of H 2M changes the dispersive coupling in Eq. ( 4) as 2 ), which must be taken into account when following the prescription to calculate the valid gates.This step can be done for an arbitrary amount of spectator modes, given that the mode frequencies ν m fulfil the condition η f n ≪ |ν m − nω k | for all odd n.
Although we simulate the performance of the gate with the two-mode Hamiltonian H f = H (±)  d + H 0 + H 2M (see column ∆I 2M in table I), due to computational limitations we use the single-mode Hamiltonian H s = H (±)  d +H 0 +H eff 2M instead, where H eff 2M = 1 3 νη 2 rS 2 z is the second-order contribution of H 2M .See supplementary note 7 for additional details.Here, H (±)  d stands for H (+)  d (H (−) d ) every first (second) half of a TQXY16 block.We investigate two regimes: (i) η = 0.005 (g B = 19.16T/m), which is the state-of-the-art of current experiments [18,54], and (ii) η = 0.04 (g B = 153.2T/m), which can be reached in near future setups [39].
In regime (i), we consider three different gates, all within the 9th harmonic.The first gate (G1), with a duration t g = 1.64 ms, appears after five TQXY16 blocks with pulse length t π = 20.5 µs reaching Ω pp = (2π) × 124 kHz.The second gate (G2) with gate time t g = 3.28 ms uses ten TQXY16 blocks with pulse length t π = 20.5 µs reaching Ω pp = (2π)×77.8kHz.The third gate (G3), with the gate-time t g = 3.94 ms and Ω pp = (2π) × 78.69 kHz, uses twelve blocks, each with a different pulse length and detuning, while it incorporates a technique to mitigate errors due to mode decoherence, see supplementary note 8.In regime (ii) we consider a gate within the 5th harmonic (G4).This gate occurs after two TQXY16 blocks where t g = 368 µs, t π = 11.5 µs, and Ω pp = (2π) × 80.9 kHz.For further details regarding pulse parameters, see supplementary note 7.
The performance of the four gates in the presence of distinct error sources is shown in table I.Each simulated experiment starts from the state |+ x + y ⟩ and targets the Bell-state , while in all cases we consider an initial motional thermal state with n = 1 [54].Other initial states result in similar values for the fidelity.
In the 2nd and 3rd columns of table I we show the gate error I = 1 − F obtained by concatenating XY8 or TQXY16 blocks, respectively.Here, F = ⟨ Φ+ |ρ| Φ+ ⟩/ Tr(ρ 2 ) [67], where ρ is the final state after tracing out the bosonic states.Notice that TQXY16 blocks achieve a clearly superior performance due to efficient decoupling from finite pulse effects.For these, I TQXY16 ≤ 10 −6 for all gates except G4, where finite the residual qubit-boson entanglement limits the error to approximately 10 −5 .In the fourth column we evaluate the effect of the second mode b by numerically simulating the two-mode Hamiltonian H f (initialising the second mode b in a thermal state with n = 1), which results in I 2M .The infidelities relative to the previous case (i.e.∆I 2M = I 2M − I TQXY16 ) are given in the "∆I 2M " column of table I. Again, the effect of the second mode is relevant only for G4, which contributes 3.8×10 −5 to the total error.Importantly, this demonstrates that our gate is compatible with the presence of spectator modes.
To investigate the effect of crosstalk, we add the term 54 and 20.34 MHz for regimes (i) and (ii), respectively.The results are given in the "∆I CT " column of table I.In contrast to the effect of the spectator mode, crosstalk is most harmless with the G4 gate.This is expected, as G4 operates with a larger qubit detuning ∆ω than the rest, while using a similar Rabi frequency.To reduce the impact of crosstalk, we combine our pulses with sin 2 -shaped ramps at the beginning and end of each pulse, see supplementary note 7, and optimize the length of the ramp using numerical simulations.The resulting infidelities are shown in the "∆I CT * " column.Note that the sin 2 ramp reduces the value of I by at least an order of magnitude in most cases.

Re(↵)
< l a t e x i t s h a 1 _ b a s e 6 4 = " 6 s < l a t e x i t s h a 1 _ b a s e 6 4 = " P j 7 e v u y J A W K S F 7 J i J d e l 8 q + E a l 8 = " > l a t e x i t > I < l a t e x i t s h a 1 _ b a s e 6 4 = " P j 7 e v u y J A W K S F 7 J i J d e l 8 q + E a l 8 = " > < l a t e x i t s h a 1 _ b a s e 6 4 = " 8 P A 5 W 8 5 y 9 o 0 c g 0 k 9 a a K B y o T 6 H g g FIG. 3. Sensitivity to errors: (a) Gate error I versus 1/T * 2 for the gates G1 (blue), G2 (red), G3 (purple) and G4 (green).In (b) and (c), the error is shown for Rabi frequency shifts Ω(t) → (1 ± δ Ω )Ω(t) and in the mode frequency ν → (1 ± δ ν )ν.(d) Error under heating for different rates ṅ.(e) Trajectory of α(t) for gates G1 (circular, blue), G2 (circular, red), G3 (non-circular, purple), and G4 (circular, green).
Table I shows that mode heating is is the main source of error for gates in regime (i).This, along with dephasing and crosstalk, limits the fidelity of these gates to the 10 −3 regime.Despite its longer gate duration, G3 achieves better perfor-mance in terms of motion-induced errors than G1 and G2, proving the validity of the mode decoherence protecting technique described in the supplementary note 8. Finally, table I shows that G4 is the most robust w.r.t.experimental imperfections.This is reasonable since it uses a larger η and is an order of magnitude faster than the other gates.In particular, we find that G4 achieves infidelities on the 10 −4 regime, mainly limited by residual qubit-boson entanglement caused by off-resonant harmonics and the spectator mode.Note that the influence of this error has been taken into account in the supplementary note 3, where we also discuss potential effects of micromotion.

III. DISCUSSION
We have presented a DD sequence (TQXY16) based on the delivery of low-intensity π pulses that achieve faster-thandispersive two-qubit gates.Without the need of any numerical optimisation, we have designed entangling gates which are robust to fluctuations in qubit frequencies and control fields, as well as to finite-pulse effects hindering a high-fidelity performance.In addition, we have demonstrated the versatility of our protocol to adopt forms that provide an increased robustness against crosstalk and mode decoherence.
Our scheme is best suited for systems i) using longitudinal qubit-boson coupling with η ≪ 1, ii) where dephasing is the main source of qubit decoherence, and iii) where the Rabi frequencies Ω(t) are of the order (or far below) the mode frequencies ν.This is the case, e.g., for spin qubits coupled to microwave cavities [5,6].In Bosco et al. [6], η ∼ 10 −2 and Ω/ν ∼ 0.1.Superconducting qubit architectures exploiting longitudinal qubit-boson coupling have also been proposed [3,4].Our method is also well suited for these systems when working with small η.
Finally, we tested the performance of our protocol in trapped ions coupled via static magnetic field gradients, where conditions (i), (ii), and (iii) are perfectly satisfied.Compared to existing multi-level schemes [18,39], our method has the advantage of using only two levels, which lowers the experimental requirements.Compared to previous pulsed DD methods [51], our method has the advantage of using realistic pulse intensities.Using detailed numerical simulations, we have obtained infidelities within the 10 −3 threshold at state-of-the-art conditions, and in the 10 −4 regime in near-future setups.
Our numerical simulations for dephasing consider an additional ±δω/2S z term in H s , where δω = √ 2/T * 2 .In all three cases, each point is the average error obtained by a positive (e.g.ω j → ω j + δω) and a negative (e.g.ω j → ω j − δω) displacement.If n is even, value of the integral is zero.In n is odd, f n is and its value will depend on the form of f z (x).For convenience, we make a last change of variable x → 1/2 − x leading to where, now, x = 0 is the central point of the pulse.In the following we calculate f n for the cases of instantaneous, top-hat, and modulated π pulses.
Top-Hat Pulses: For top-hat pulses, f z (x) = sin (πx/τ π ) from x = 0 to x = τ π /2, while f z (x) = 1 for the rest.Here, τ π = 2t π /τ = 2π/Ωτ.We divide the integral accordingly, dx sin (nπx). (S19) For the first part, we use that 2 sin A sin B = cos (A − B) − cos (A + B), and, then, and this first part gives The second part gives 1 nπ cos (nπτ π /2).The whole integral gives or Modulated Pulses: The ansatz proposed for the modulation function is where t π is the duration of the π pulse, t i and t m = t i + t π /2 are the initial and central points of the pulse, t r = t m + bt π and t l = t m − bt π .Also, 0 < b < 0. In the appropriate regime, these integrals result in [S2] and respectively.

SUPPLEMENTARY NOTE 4: RESIDUAL SPIN-PHONON ENTANGLEMENT AND MICROMOTION
We assume a generic multiqubit-boson Hamiltonian of the form In this case, the first-order unitary operator is We assume the ideal pure state ρ i is given after applying the second-order evolution operator U g = e iθ(t g )S2 z .If we assume [S 2 z , H(t)] = 0, the reduced final state is approximately given (considering only the first and second-order evolution operators) ρ = Tr a {U(t)U g ρ 0 ⊗ ρ a U † g U † (t)} = Tr a {U(t) ρ i ⊗ ρ a U † (t)}, where t is now the final time t g .If η ≪ 1, we can expand U(t) in η, and we get where we use Tr a (a n ρ a ) = Tr a ((a † ) n ρ a ) = 0, and where S (t) = iν t 0 S (t ′ )dt ′ .Expanding the fidelity, defined here as and f a e S n j W K x S 9 s j e D u 0 z 8 j B Q h Q 6 1 X + O r 0 F U 0 i J p E K Y k z b 9 2 L s p k Q j p 4 J N 8 p 3 E s J j Q E R m w t q W S R M x 0 0 9 m 1 E / f U K n 0 3 V N q W R H e m / p 5 I S W T M O A p s Z 0 R w a B a 9 q f i f 1 0 4 w v O q m X M Y J M k n n i 8 J E u K j c 6 e t u n 2 t G U Y w t I V R z e 6 t L h 0 Q T i j a g v A 3 B X 3 x 5 m T T O y / 5 F u X J X K V a v s z h y c A w n U A I f L q E K t f a e S n j W K x S 9 s j e D u 0 z 8 j B Q h Q 6 1 X + O r 0 F U 0 i J p E K Y k z b 9 2 L s p k Q j p 4 J N 8 p 3 E s J j Q E R m w t q W S R M x 0 0 9 m 1 E / f U K n 0 3 V N q W R H e m / p 5 I S W T M O A p s Z 0 R w a B a 9 q f i f 1 0 4 w v O q m X M Y J M k n n i 8 J E u K j c 6 e t u n 2 t G U Y w t I V R z e 6 t L h 0 Q T i j a g v A 3 B X 3 x 5 m T T O y / 5 F u X J X K V a v s z h y c A w n U A I f L q E K t < l a t e x i t s h a 1 _ b a s e 6 4 = " A X Z E f w F 1 z 9 u T s p e m H / w g + Form of the modulation functions f x,y,z (t) for a single XY4 block with top hat pulses and with t π = τ/4 = π/Ω th .In the XY4 block, the first and third pulse rotate the qubit along the x axis, that is, Ω x (t) = Ω(t), Ω y (t) = 0.For the second and fourth pulses, the same is true for the y axis.
Here, the sine components of f z (t) are zero because the integral τ 0 f z (t) sin (nωt) is zero for all n ∈ N.This can be proven by dividing the integral in two parts ({0, t 1 } and {t 1 , t 2 }) and using the symmetry property f z (t 1 + t) = f z (t 1 − t).In a similar fashion, the symmetry properties f y (t f ⊥ (t + t 1 /2, 0) cos (nωt/2)dt.Note that e n is zero if n is even.In the following, we will calculate the second-order Hamiltonian associated to Eq. (2) of the main text after a XY4 pulse sequence.After the Fourier expansion, this looks like where c x,y n = a x,y n + ib x,y n , and cx,y n is the complex conjugate.As e n is zero for even n, nω and nω/2 never coincide, and, if ην ≪ ω, the second order Hamiltonian is (see Eq. ( S6)) where As we indicated above, e n is zero for even n.Thus, J xy k = 0 for XY4 sequences.The rest of the terms will have a non-zero value and ought to be removed by applying refocusing techniques.The second-order Hamiltonian for the inverse XY4 sequence, that is, YXYX, is equivalent except for the term 1  2 η 2 ν[B ′ k (a † a + 1) + B ′′ k ]S z , which has the opposite sign.Concatenating the XYXY and the YXYX blocks in what is called the XY8 sequence, we will achieve a partial refocusing of this term.Refocusing the − 1 2 η 2 νJ ⊥ k (S 2 x + S 2 y ) term is possible by driving the two qubits with opposite phases in a second XY8 block.The effective Hamiltonian will be equivalent except that, now, the operators S x,y will become S (−)  x,y = σ x,y 1 − σ x,y 2 , changing the sign of terms σ u 1 σ µ 2 where u, µ ∈ {x, y}.
where a † (a) and b † (b) are the creation (annihilation) operators associated to the longitudinal center-of-mass and breathing modes of the two-ion crystal.To simulate the performance of the gate under the influence of two modes, we evolve the state according to Hamiltonian H f = H (±) d + H 0 + H 2M , starting from motional thermal states with n = 1, and truncating the Hilbert space of the first (a) and second (b) mode to 16 and 12 states, respectively.

E. Effective Hamiltonian of the driven system
Here we justify the use of the term H eff 2M as a right quantity to account for the second-order effect of the breathing mode, without explicitly taking it into account in our simulations.For that, we first calculate the effective Hamiltonian of the driven system for the two-mode case.
After moving into an interaction picture w.r.t.H d + νa † a + √ 3νb † b, neglecting the terms going with f x,y (t) and expanding f z (t) as a Fourier series, the Hamiltonian of the driven, two-mode system is f n ηνS (+) z ae −iνt (e inω k t + e −inω k t ) − 3 −1/4 2 f n ηνS (−) z be −i √ 3νt (e inω k t + e −inω k t ) + H.c.. (S69) In Hamiltonian (S69), the only resonant term is the one going with e −i(ν−kω k ) and its complex conjugate, i.e. 1 2 f k ηνS z ae −i(ν−kω k )t + H.c.. Using Eq. (S6), the influence of the rest of the terms is well approximated by the second order Hamiltonian where . Because (S (−) z ) 2 = 4I − S 2 z , Hamiltonian (S70) is, up to a global phase, The effective Hamiltonian describing the driven system is then f k ηνS z a(e −i(ν−kω k )t + H.c.) − 1 2 η 2 ν(J a k − J b k /3)S 2 z , (S72) which is equivalent to Hamiltonian (4) in the main text but with J k → J k − J b k /3.In our numerical simulations, we do not want to include the dynamics of the breathing bosonic mode.Yet, we want to capture its second-order effect correctly.This is achieved by the following Hamiltonian where r = J b k / ∞ n=1 f 2 n .Following the procedure described from Eq. (S68) to Eq. (S72), but now with H, one realises that the effective Hamiltonian corresponding Eq. (S73) is equivalent to H ′ (t).

F. Crosstalk
In the two-ion system, the frequency of qubit j is related with the intensity of the magnetic field in the ion's equilibrium position z 0 j , and both frequencies differ by ∆ω = ω 2 − ω 1 = γ e g B d. Here, the distance between the ions is d = (e 2 /2πε 0 Mν 2 ) 1/3 where e is the electric charge of the electron and ε 0 is the vacuum permittivity.For two 171 Yb + ions each with mass M = 171 amu and trapped with longitudinal frequency ν = (2π) × 220 kHz, this value gives ∆ω = 2.54 and 20.34 MHz for g B = 19.57and 153.2 T/m, respectively.
The results of the simulations show that, for the Rabi frequencies we use, crosstalk terms have a non-negligible effect.To reduce their impact we combine each pulse with a sin 2 -shaped ramp at the beginning and end of each pulse.Each pulse is then constructed following where f z (t) is defined in Eqs.(8,9) of the main text and the pulse parameters can be found in Appendix A. Also, we freely change a factor δΩ π multiplying the Rabi frequency, i.e.Ω(t) → (1 + δΩ π )Ω(t), to ensure that t π 0 dt ′ Ω(t ′ ) = π.We simulate the performance of the gate for different durations of the sin 2 -shaped ramp t ramp .The obtained infidelities are shown in Fig. S2(a).As it is shown, the achieved infidelity is in general better for the ramped case.The infidelities presented in table I of the main text are calculated using ramps of size 149, 295, 1260 and 49 ns, for gates G1-4, respectively.Because of the shift δΩ π , the maximum values of the Rabi frequencies Ω pp /2π change to 125.2, 79.05, 82.09, and 81.14 kHz, respectively.The Rabi frequencies Ω x (t) and Ω y (t) are shown in Fig. S2(b), for gates G1-4, and for the first eight pulses.

G. Heating of the bosonic mode
To include the effect of motional heating we solve the master equation where ρ is the density matrix and the Lindblad superoperator is with N = [exp ( ℏν k B T ) − 1] −1 and T = 300 K.For N ≫ 1, Γ relates to the heating rate as ṅ = Γ N. For our simulations, we use a heating rate of ṅ = 35 ph/s for regime (i).We derived this from the rate ṅref = 120 ph/s given in Ref. [S7] for ν ref = (2π) × 120 kHz, and using the relation ṅ = ṅref × (ν ref /ν) 2 [S8].For regime (ii), we use a heating rate of ṅ = 100 ph/s [S9].
t e x i t s h a 1 _ b a s e 6 4 = " kP K l d X v b J H p m i 4 k B F B a + D w z A f Y U = " > A A A C A n i c b V B N S 8 N A E N 3 4 W e t X 1 Z N 4 C R a h X k p S R D 0 W v X i s Y j + g K W W y n b Z L N 9 m w u x F K K F 7 8 K 1 4 8 K O L V X + H N f + O m z U Fb H w w 8 3 p t h Z p 4 f c a a 0 4 3 x b S 8 s r q 2 v r u Y 3 8 5 t b 2 z m 5 h b 7 + h R C w p 1 q n 1 b L 2 w U N 7 e 2 d 3 b t v f 2 m E o n E p I E F E 7 I d g C K M R q S h q W a k H U s C P G C k F Y y u M 7 / 1 Q K S i I r r X 4 5 h 0 O Q w i G l I M 2 k g 9 + 9 A X M Z G g h Y y A k / S W T 8 o + s H g I p z 2 7 5 F b c K Z x F 4 u W k h H L U e / a X 3 x c 4 4 S T S m I F S H c + N d T c F q S l m Z F L 0 E 0 V i w C M Y k I 6 h 2 T 7 V T a c v T J w T o / S d U E h T k X a m 6 u + J F L h S Yx 6 Y T g 5 6 q O a 9 T P z P 6 y Q 6 v O y m N I o T T S I 8

5 <
r H e r Y 9 Z 6 5 K V z x y g P 7 A + f w B 0 2 p d 2 < / l a t e x i t > l a t e x i t s h a 1 _ b a s e 6 4 = " d K 5 5 j g O 9 x k 5 S Z B v v 3 S F o Y / S m d C x 6 4 e 4 8 2 9 M 2 1 l o 6 4 E L h 3 P u 5 d 5 7 Y s m o N p 7 3 7 a y s r q 1 v b J a 2 y t s 7 u 3 v 7 7 s F h W y e p w q S F E 5 a o b o w 0 r / 4 b b / 4 b 0 3 Y P 2 v p g 4 P H e D D P z wo Q z b V z 3 2 y m s r K 6 t b x Q 3 S 1 v b O 7 t 7 5 f 2 D l o 5 T R W i T x D x W n R B r y p m k T c M M p 5 1 E U S x C T t v h 6 G b q t 5 + o 0 i y W D 2 a c 0 E D g g W Q R I 9 h Y 6 d G / E 3 S A z 3 y Z 9 s o V t + rO g J a J l 5 M K 5 G j 0 y l 9 + P y a p o N I Q j r X u e m 5 i g g w r w w i n k 5 K f a p p g M s I

5 w 9 H
F J S I < / l a t e x i t > cos (k! k t) < l a t e x i t s h a 1 _ b a s e 6 4 = " c 2 P 4 + U l 2 5 z a FV O T w 3 r C d 2 K D I x 1 4 = " > A A A B / H i c b V B N S 8 N A E N 3 U r 1 q / o j 1 6 W S y C p 5 o U U Y 9 F L x 4 r 2 A 9 o Q t l s N + 3 S z S b s T o Q Q 6 l / x 4 k E R r / 4 Q b / 4 b t 2 0 O 2 v p g 4 P H e D D P z g k R w D Y 7 z b Z X W 1 j c 2 t 8 r b l Z 3 d v f 0 D + / C o o + N U U d a m s Y h V L y C a C S 5 Z G z g I 1 k s U I 1 E g W D e Y 3 M 7 8 7 i N T m s f y A b K E + R E Z S R 5 y S s B I A 7 v q O o 4 H P G L a k y m G 8 4 a X 8 I F d c + r O H H i V u A W p o Q K t g f 3 l D W O a R k w C F U T r v u s k 4 O d E Aa e C T S t e q l l C 6 I S M W N 9 Q S c w 6 P 5 8 f P 8 W n R h n i M F a m J O C 5 + n s i J 5 H W W R S Y z o j A W C 9 7 M / E / r 5 9 C e O 3 n X 4 g G d 4 h T d H O S / O u / M x b 1 1 x s p k j + A P n 8 w f / u I 7 D < / l a t e x i t > f z (t) < l a t e x i t s h a 1 _ b a s e 6 4 = " W P e e b j c P U 4 o c I x 5 T b k D 6 8 E r G t k 8 = " > A A A B 6 X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y G o B 6 D X j x J F P O A Z A m z k 9 l k y

4 <
5 e p 9 t V S 7 z u L I w w m c w j l 4 c A k 1 u I U 6 N I B B C M / w C m / O y H l x 3 p 2 P R W v O y W a O 4 Q + c z x 8 a b I 0 X < / l a t e x i t > l a t e x i t s h a 1 _ b a s e 6 4 = " P l 1 / R 8 l P

5 < l a t e x i t s h a 1 _ b a s e 6 4 =
n / T b p c I T N i a A l l i t t b C e t T R Z m x 6 R R s C N 7 8 y 4 u k f l L 2 z s u n 9 6 e l y v U s j j w c w C E c g w c X U I F b q E I N G P T g G V 7 h z R H O i / P u f E x b c 8 5 s Z h / + w P n 8 A a D 5 j W E = < / l a t e x i t >N = " l / T M O R K M f s c r B y j Q a a q I W s 6 z i 4 s = " > A A A B 6 3 i c b V B N S w M x E J 2 t X 7 V + V T 1 6 C R b B U 9 m V U r 0 I R S + e p I L 9 g H Y p 2 T T b h i b Z J c k K Z e l f 8 O J B E a / + I W / + G 7 P t H r T 1 w c D j v R l m 5 g U x Z 9 q 4 7 r d T W F v f 2 N w q b p d 2 d v f 2 D 8 q H R 2 0 d J Y r Q F o l 4 p L o B 1 p Q z S V u G G U 6 7s a J Y B J x 2 g s l t 5 n e e q N I s k o 9 m G l N f 4 J F k I S P Y Z N L 9 t e c O y h W 3 6 s 6 7 y O I p w A q d w D h 5 c Q g P u o A k t I D C G Z 3 i F N 0 c 4 L 8 6 7 8 7 F o L T j 5 z D H 8 g f P 5 A w l c j Z c = < / l a t e x i t > FIG. 2. Amplitude modulated π pulses: (a) Ω(t) for a single XY8 block.The whole sequence here is a concatenation of 2N blocks.(b) Gate time t g as a function of Ω pp , for values of d between 0 and d max .Values satisfying t g = 8Nτ k are represented by square and round markers for k = 9 and k = 7, respectively.(c) Ω(t) and f z (t) for cases N = 5 with d = 1.888 (blue) and N = 10 with d = 0.908 (red) of the 9th harmonic.(d) Gate phase θ(t) for N = 5 (dashed blue line), N = 10 (solid red line), and the dispersive case (dotted line).

2 <
l a t e x i t s h a 1 _ b a s e 6 4 = " kP K l d X v b J H p m i 4 k B F B a + D w z A f Y U = " > A A A C A n i c b V B N S 8 N A E N 3 4 W e t X 1 Z N 4 C R a h X k p S R D 0 W v X i s Y j + g K W W y n b Z L N 9 m w u x F K K F 7 8 K 1 4 8 K O L V X + H N f + O m z U Fb H w w 8 3 p t h Z p 4 f c a a 0 4 3 x b S 8 s r q 2 v r u Y 3 8 5 t b 2 z m 5 h b 7 8 6 A 9 B D N e + l 4 n 9 e O 9 b 9 y 0 7 C w i j W G N L Z o n 7 M b S 3 s N A + 7 x y R S z c e G A J X M 3 G r T I U i g 2 q S W N y G 4 8 y 8 v k k a l 7 J 6 X K 7 d n x e p V F k e O H J F j U i I u u S B V c k N q p E 4 o e S T P 5 J W 8 W U / W i / V u f c x a l 6 x s 5 o D 8 g f X 5 A 3 Z w l 3 c = < / l a t e x i t >

1 √
S37) where ⟨•⟩ = Tr(•ρ i ).For final state |Φ + ⟩ = t e x i t s h a 1 _ b a s e 6 4 = " U E p 9 S M 9 / b 8 h M u + 2 I W E L V 9 N I b b S o = " > A A A B 7 X i c b V D L S g N B E O z 1 G e M r 6 t H L Y h D i J e x K U I 9 B L x 4 j m A c k I c x O Z p M x s z P L T K 8 Y l v y D F w + K e P V / v P k 3 T p I 9 a G 1 C D O l B 4 g G d 4 h T d H O S / O u / Mx b 1 1 x s p k j + A P n 8 w f 8 q o 7 B < / l a t e x i t > f x (t) < l a t e x i t s h a 1 _ b a s e 6 4 = " S r 6d D P s Q / Z E k 9 U L Q 8 3 W A I w c h P i E = " > A A A B 7 X i c b V B N S w M x E M 3 W r 1 q / q h 6 9 B I t Q L 2 V X i n o s e v F Y w X 5 A u 5 R s m m 1 j s 8 m S z A r L 0 v / g x Y M i X v 0 / 3 v w 3 p u 0 e t P X B w O O 9 G W b m B b H g B l z 3 2 y m s r W 9 s b h W 3 S z u 7 e / s H 5 c O j t l G J p q x F l V C 6 G x D D B J e s B R w E 6 8 a a k S g Q r B N M b m d + 5 4 l p w 5 V 8 g D R m f k R G k o e c E r B S O x y k V T g f l C t u z Z 0 D r x I v J x W U o z k o f / W H i i Y R k 0 A F M a b n u T H 4 G d H A q W D T U j 8 x L C Z 0 Q k a s Z 6 k k E T N + N r 9 2 i s + s M s S h 0 r Y k 4 L n 6 e y I j k T F p F N j O i M D Y L H s z 8 T + v l 0 B 4 7 W d c x g k w S R e L w k R g U H j 2 O h 5 y z S i I 1 B J C N b e 3 Y j o m m l C w A Z V s C N 7 y y 6 u k f V H z L m v 1 + 3 q l c Z P H U U Q n 6 B R V k Y e uU A P d o S Z q I Y o e 0 T N 6 R W + O c l 6 c d + d j 0 V p w 8 p l j 9 A f O 5 w / + M Y 7 C < / l a t e x i t > f y (t) < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 / E p D G J / 5 L h 1 s 3 x A Y x i 6 6 R w 2 B 3 E = " > A A A B 7 X i c b V D L S g N B E O z 1 G e M r 6 t H L Y h D i J e x K U I 9 B L x 4 j m A c k I c x O Z p M x s z P L T K 8 Q l / y D F w + K e P V / v P k 3 T p I 9 a G 1 C D O l B 4 g G d 4 h T d H O S / O u / M x b 1 1 x s p k j + A P n 8 w f / u I 7 D < / l a t e x i t > f z (t) < l a t e x i t s h a _ b a s e = " N b 7 k a o o b / E i w d F v P p T v P l v 3 L Y 5 a O u D g c d 7 M 8 z M C 1 P O l H a c b 6 u 0 t r 6 x u V X e r u z s 7 u 1 X 7 Y P D j k o y S W i b J D y R v R A r y l l M 2 5 p p T n u p p F i E n H b D 8 c 3 M 7 z 5 S q V g S P + h J S n 2 B h z G L G M H a S I F d 9 e 4 E H e I g 9 6 R A e j Q N 7 J p T along with trigonometric identities sin (a ± b) = sin a cos b ± cos a sin b and cos (a ± b) = cos a cos b ∓ sin a sin b, can be used to prove that a x n = − cos (3πn/4)e n , b x n = − sin (3πn/4)e n , (S56) a y n = cos (πn/4)e n , b y n = sin (πn/4)e n , (S57) where e n = [cos (nπ) − 1] 2 τ t π /2 0 e inωt + e −inωt )S z + (c x n S x + cy n S y )e inωt/2 + (c x n S x + c y n S y )e −inωt/2 + H.c., (S58)

TABLE I .
(10or budget: Column XY8 and TQXY16 show the infidelities after an evolution with Hamiltonians H (+) d + H 0 + H eff 2M and H s , respectively.The remaining columns show infidelities relative to the TQXY16 case (e.g.∆I 2M = I 2M − I TQXY16 ), taking into account various experimental imperfections.In columns ∆I 2M , ∆I CT , and ∆I CT * , infidelities obtained considering a second mode, crosstalk, and crosstalk with the sin 2 ramp are shown.In ∆I δΩ , ∆I δν , and ∆I T 2 we show relative infidelities considering static shifts of δΩ = 5 × 10 −3 , δν = 10 −5 , and δω = (2π) × 2 √ 2 kHz.ṅ shows the error considering heating with rates ṅ1 = 35 ph/s and ṅ2 = 100 ph/s for regimes (i) and (ii), respectively.The last column shows the overall error obtained by summing the values of all columns except those in I XY8 and ∆I CT .Gate I XY8 I TQXY16 ∆I 2M ∆I CT ∆I CT * ∆I T 2 ∆I δΩ ∆I δν ∆Iṅ I total(10

TABLE II .
Suitable pulse parameters.Suitable values the for b, c and |d max |. |d max | corresponds to the maximum value of |d| for which a physical pulse (i.e, | f z |(t) ≤ 1) can still be generated.