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

Iñigo Arrazola<EMAIL_ADDRESS>Vienna Center for Quantum Science
and Technology, Atominstitut, TU Wien, 1040 Vienna, Austria  Jorge Casanova
Department of Physical Chemistry, University of the Basque Country UPV/EHU,
Apartado 644, 48080 Bilbao, Spain  EHU Quantum Center, University of the
Basque Country UPV/EHU, Leioa, Spain  IKERBASQUE, Basque Foundation for
Science, Plaza Euskadi 5, 48009 Bilbao, Spain

###### Abstract

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. This represents an improvement
of more than an order of magnitude over the state-of-the-art.

_Introduction.-_ High-fidelity entanglement generation among qubits is crucial
for quantum information processing Nielsen . In most platforms, entangling
gates come via direct interactions (e.g. hyperfine fields among nuclear spins)
or via a bosonic mediator. Examples of the latter include solid-state qubits
coupled to microwave resonators Jin12 ; Billangeon15 ; Richer16 ; Beaudoin16 ,
or trapped ions sharing vibrational modes Mintert01 . In this scenario, the
paradigmatic Mølmer-Sørensen (MS) gate Sorensen99 ; Solano99 ; Sorensen00 ;
Sackett00  and related schemes Milburn00 ; Leibfried03 ; Zheng04  reach
entanglement operations with inherent robustness to uncertainties in the
bosonic state.

In recent years, several gate schemes have been developed that operate on the
same MS principle but are also robust against other sources of error
Valahu22_2 . Some examples are frequency or amplitude-modulated gates
decoupling from mode decoherence Cross15 ; Haddadfarshi16 ; Webb18 ; Shapira18
; Zarantonello19 ; Sutherland20 , from spectator modes Hayes12 ; Choi14 ;
Green15 ; Leung18 ; Schafer18 ; Lu19 ; Milne20  or from deviations in the
qubit-boson coupling strength Royer17 ; Ge19 ; Burd21 ; Shapira22 . In another
vein, dynamical decoupling (DD) is a well established paradigm to protect
qubits against decoherence Viola98 ; Ban98 . Some continuous DD techniques
have been demonstrated to be suitable for quantum gate implementations
Timoney11 ; Tan13 ; Harty16 ; Weidt16 ; Guo18 , while pulsed DD methods
achieve increased robustness employing suited sequences such as XY8 Carr54 ;
Meiboom58 ; Gullion69 ; Souza12 ; Kabytayev14 ; Munuera20 ; Dong21 ; Ezzell22
or AXY Casanova15 ; Casanova17 ; Arrazola18 . However, the use of pulsed DD to
protect oscillator-mediated gates has been mostly limited to dispersive
regimes Piltz13 ; Qiu21 ; Barthel22 ; Morong22 , and to few spin-echo Paik16 ;
Ballance16 ; Bazavan22  or rotary-echo Tan13  pulses. Note that the
application of several $\pi$ pulses is desirable for efficient elimination of
time-varying noise.

Interestingly, Manovitz et al. Manovitz17  have experimentally shown that the
Mølmer-Sørensen gate can be combined with pulse sequences given the ability to
tune and turn on-and-off the qubit-boson coupling. Another possibility
explored theoretically is to combine an always-on qubit-boson coupling with
strong $\pi$ pulses Rabl10 ; Arrazola18 ; Rosenfeld21 . 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 $\pi$ pulses is experimentally challenging since high-power
controls are needed, while these induce crosstalk and hinder the applicability
in multimode scenarios.

In this Letter, we design a DD sequence with finite-width $\pi$ pulses –named
TQXY16– that achieves faster-than-dispersive entangling gates using static
(i.e. non-tunable) qubit-oscillator coupling. Importantly, our gates decouple
from dephasing, pulse imperfections, and unwanted finite-pulse 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 radio frequency
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.

_The method.-_ We consider a system that comprises two qubits and a bosonic
mode –with frequencies $\omega_{1},\omega_{2}$ and $\nu$– coupled via
longitudinal coupling Mintert01 ; Jin12 ; Billangeon15 ; Richer16 ; Beaudoin16
,

$H_{0}=\nu a^{\dagger}a+\eta\nu(a+a^{\dagger})S_{z}.$ (1)

Here, $a^{\dagger}(a)$ is the creation (annihilation) operator of the bosonic
mode, $S_{\mu}\equiv\sigma_{1}^{\mu}+\sigma_{2}^{\mu}$ with $\mu\in{x,y,z}$
are collective qubit operators, and $\eta\nu$ is the coupling strength. We
assume the usual experimental scenario $\eta\ll 1$, thus we stay away from
other paradigms that require stronger qubit-boson couplings GarciaRipoll03 ;
Duan04 ; Steane14 ; Bentley13 ; Sameti21 . $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 $\nu$, and (ii) Protection of the gate from qubit noise of the form
$\epsilon_{j}(t)\sigma_{j}^{z}/2$ leading to dephasing. When driving the
qubits, $H_{0}$ is completed with the term
$H_{d}(t)=\sum_{\mu=x,y}\Omega_{\mu}(t)S_{\mu}/2$. In an interaction picture
with respect to (w.r.t) $H_{d}+\nu a^{\dagger}a$ we get

$H(t)=\eta\nu(ae^{-i\nu t}+a^{\dagger}e^{i\nu
t})\sum_{\mu=x,y,z}f_{\mu}(t)S_{\mu},$ (2)

where $\sum_{\mu=x,y,z}f_{\mu}(t)S_{\mu}=U^{\dagger}_{d}(t)S_{z}U_{d}(t)$ with
$U_{d}(t)=\mathcal{T}\exp{[-i\int_{0}^{t}H_{d}(t^{\prime})dt^{\prime}]}$ being
the time-ordered propagator.

If driving fields are delivered as instantaneous $\pi$ pulses (note this
requires $\Omega_{x,y}\gg\nu$ during the application of the pulse) spaced
$\tau_{k}/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(b).
For the moment we consider instantaneous pulses, while later we treat the
realistic case of non-instantaneous ones. As $\pi$ pulses are applied
periodically, $f_{z}(t)$ takes the form of a function with period $\tau_{k}$
such that $f_{z}(t)=\sum_{n=1}^{\infty}f_{n}\cos{(n\omega_{k}t)}$, where
$\omega_{k}=2\pi/\tau_{k}$ and
$f_{n}=\frac{2}{\tau_{k}}\int_{0}^{\tau_{k}}dt^{\prime}f_{z}(t^{\prime})\cos{(n\omega_{k}t^{\prime})}$.
Hence, under the assumption of instantaneous pulses we get

$H(t)=\eta\nu\sum_{n=1}^{\infty}f_{n}\cos{(n\omega_{k}t)}(ae^{-i\nu
t}+a^{\dagger}e^{i\nu t})S_{z},$ (3)

whilst setting an interpulse spacing $\tau_{k}/2$ such that
$\omega_{k}\approx\nu/k$ leads to a resonant qubit-boson interaction via the
$k$th harmonic. As $\eta\ll 1$, the terms in Eq. (3) that rotate with
frequencies $\pm|\nu-n\omega_{k}|$ (where $n\neq k$) and
$\pm|\nu+n\omega_{k}|$ can be substituted by their second-order contribution
leading to

$H(t)\approx\frac{1}{2}\eta\nu f_{k}(ae^{-i\xi_{k}t}+{\rm
H.c.})S_{z}-\frac{1}{2}\eta^{2}\nu J_{k}S_{z}^{2},$ (4)

where $\xi_{k}=\nu-k\omega_{k}$ is the detuning w.r.t. the $k$th harmonic and
$J_{k}=f_{k}^{2}/4+\sum_{n\neq k}f_{n}^{2}/(1-n^{2}/k^{2})$ is an effective
spin-spin coupling constant that contains contributions from all harmonics.
The propagator associated to Hamiltonian (3) is

$U(t)=\exp{\\{[\alpha(t)a^{\dagger}-\alpha^{*}(t)a]S_{z}\\}}\times\exp{[i\theta(t)S_{z}^{2}]}$
(5)

where $\alpha(t)=-i\eta\nu\int_{0}^{t}dt^{\prime}f_{z}(t^{\prime})e^{i\nu
t^{\prime}}\approx-\eta\nu f_{k}/(2\xi_{k})(e^{i\xi_{k}t}-1)$ and

$\theta(t)={\rm
Im}\int_{\mathcal{C}}\alpha\,d\alpha\approx\frac{\eta^{2}\nu^{2}f_{k}^{2}}{4\xi_{k}}\Big{[}t-\frac{\sin{(\xi_{k}t)}}{\xi_{k}}\Big{]}+\frac{1}{2}\nu\eta^{2}J_{k}t,$
(6)

where $\mathcal{C}$ is the phase-space trajectory followed by $\alpha(t)$. We
choose $\tau_{k}$ such that
$\xi_{k}=2\eta\nu\Big{\\{}\sqrt{f_{k}^{2}+4\eta^{2}J_{k}^{2}}+2\eta
J_{k}\Big{\\}}$ for $J_{k}>0$. After a time $t_{\rm g}=2\pi/|\xi_{k}|$ the
propagator $U(t)$ approximates to $\exp{(i\frac{\pi}{8}S_{z}^{2})}$, which
leads to
$|\\!++\rangle\rightarrow|\Phi^{+}\rangle=\frac{1}{\sqrt{2}}(|\\!++\rangle+i|\\!--\rangle)$.

Figure 1: (a) Bell state fidelity vs final gate time $t_{\rm g}^{k}$ for
sequences where $\tau_{k}=2\pi k/(\nu-\xi_{k})$. The dispersive case
corresponds to a gate generated by Eq. (1) without $\pi$ pulses. (b) Rabi
frequency $\Omega(t)$ and modulation function $f_{z}(t)$ during a period
$\tau_{k}$ for $k=5$. For comparison, we plot $\cos{(k\omega_{k}t)}$ in green.
(c) Phase-space trajectory of $\alpha(t)$ during the application of a pulse
sequence with $k=1$ and $k=5$ in blue and green, respectively.

For instantaneous $\pi$ pulses one finds
$f_{k}=\frac{4}{k\pi}\sin{(\frac{k\pi}{2})}$. Thus, if resonance is achieved
via a low harmonic, e.g. $k=1$, the gate time is $t_{\rm
g}^{k=1}\approx\pi^{2}/4\eta\nu$. On the other hand, for sufficiently large
harmonics the gate is mostly governed by the dispersive term
$\frac{1}{2}\eta^{2}\nu J_{k\rightarrow\infty}S_{z}^{2}\approx\eta^{2}\nu
S_{z}^{2}$ in Eq. (4), leading to $t_{\rm g}^{\rm d}=\pi/8\eta^{2}\nu$. We
define faster-than-dispersive gates as those that satisfy $t_{\rm g}/t_{\rm
g}^{\rm d}<1$. For example, in the case $\eta=0.01$ and $k=1$ we find a
faster-than-dispersive gate with $t_{\rm g}^{k=1}/t_{\rm g}^{\rm d}\approx
1/16$. This is, the gate is $16$ times faster than the dispersive one.

Figure 1(a) shows the Bell-state fidelity reached by gates following Eq. (5)
for the first harmonics. More specifically, we compare the final state
$\sigma={\rm Tr}_{a}[U(t_{\rm
g})|\psi_{0}\rangle\langle\psi_{0}|U^{\dagger}(t_{\rm g})]$, (${\rm Tr}_{a}$
stands for the trace over the bosonic states and
$|\psi_{0}\rangle=|\\!++\rangle\otimes|0\rangle$) with the Bell state
$|\Phi^{+}\rangle$ using $\mathcal{F}=|{\rm
Tr}(|\Phi^{+}\rangle\langle\Phi^{+}|\sigma)|/\sqrt{{\rm
Tr}(\sigma\sigma^{\dagger})}$ Wang08 . All cases use $\eta=0.03$. In grey we
show the fidelity for the dispersive gate obtained using Eq. (1), i.e.,
without adding pulses. As expected, the $k=1$ gate (in blue) is the fastest,
and the $k=5$ gate (green) stands between the $k=1$ and the dispersive gate.

In what follows, we discuss the realistic case of non-instantaneous pulses.
For standard top-hat pulses the Fourier coefficient $f_{k}$ that quantifies
the strength of the qubit-boson interaction reads (see supplemental material
SupMat  (SM) for the derivation)

$f^{\rm
th}_{k}\approx\frac{f_{k}}{1-\nu^{2}/\Omega^{2}}\cos{\bigg{(}\frac{\pi\nu}{2\Omega}\bigg{)}}.$
(7)

Notice that for low-intensity pulses –defined as those holding $\Omega<\nu$–
the value of $f^{\rm th}_{k}$ decays with $(\Omega/\nu)^{2}$. As a result,
achieving faster-than-dispersive gates is no longer possible.

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 $\pi$ pulse. Specifically, we pose the following _ansatz_ for $f_{z}(t)$

$f_{z}(t)=\cos[{\pi(t-t_{i})/t_{\pi}]}+\beta(t)\sin{[k\omega_{k}(t-t_{m})]},$
(8)

where $t_{\pi}$ is the $\pi$ pulse duration, and $t_{i}$ and
$t_{m}=t_{i}+t_{\pi}/2$ are the initial and central points of the pulse. For
the envelope function $\beta(t)$, we propose

$\beta(t)=\frac{d}{\pi kb}\sin(\pi k/2)\bigg{[}{\rm
erf}\bigg{(}\frac{t-t_{l}}{ct_{\pi}}\bigg{)}-{\rm
erf}\bigg{(}\frac{t-t_{r}}{ct_{\pi}}\bigg{)}\bigg{]},$ (9)

where $t_{r}=t_{m}+bt_{\pi}$ and $t_{l}=t_{m}-bt_{\pi}$. The free parameters
$b$ and $c$ serve to control the width of the envelope function $\beta(t)$,
while $d$ is proportional to its amplitude. Suitable values for $b,c$ and $d$
for the first harmonics are shown in Appendix A. From now on, we assume
$t_{\pi}=\tau_{k}/2$, i.e., the pulse extends over a whole period
$\tau_{k}/2$, leading to solutions with the lowest intensities. As a result of
our pulse design, the value for the Fourier coefficient is given by $f^{\rm
m}_{k}=-{4d}/{\pi k}$ where $|d|$ can take values from $0$ to $|d_{\rm max}|$.
See SM SupMat  for the derivation. Since $f_{k}^{\rm m}$ depends on $d$, this
serves to control the strength of the interaction, thus the duration of the
gate $t_{\rm g}$. Also, $d$ relates to the amplitude of the pulse, thus to the
maximum value of the Rabi frequency $\Omega_{\rm pp}$. Typically, we look for
large values of $d$, bounded by $d_{\rm max}$ or by the experimentally
available $\Omega_{\rm 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$. The larger $k$
is, the lower is the pulse intensity. Second, we use Eq. (8) to generate the
modulation function $f_{z}(t)$ and calculate the gate time $t_{\rm g}$ and
$\Omega_{\rm pp}$ for different values of $d$. Note that $\Omega_{\rm pp}$ is
determined by calculating the maximum of $\Omega(t)=-\frac{\partial
f_{z}(t)}{\partial t}\times[1-f^{2}_{z}(t)]^{-1/2}$, see SM SupMat . We note
that the obtained $\Omega(t)$ leads 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_{\rm g}$ must
be $8N\tau_{k}$, where $N$ is an integer number. This translates into the
condition $(\nu-\xi_{k})/8k|\xi_{k}|\in\mathbb{N}$ (note $t_{\rm
g}=8N\tau_{k}$, while $t_{\rm g}=2\pi/|\xi_{k}|$ and $\tau_{k}=2\pi
k/(\nu-\xi_{k})$). 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 $\Omega_{\rm
pp}$.

Figure 2:  (a) $\Omega(t)$ for a single XY8 block. The whole sequence here is
a concatenation of $2N$ blocks. (b) Gate time $t_{\rm g}$ as a function of
$\Omega_{\rm pp}$, for values of $d$ between $0$ and $d_{\rm max}$. Values
satisfying $t_{\rm g}=8N\tau_{k}$ are represented by square and round panels
for $k=9$ and $k=7$, respectively. (c) $\Omega(t)$ and $f_{z}(t)$ for cases
$N=5$ with $d=1.888$ (blue) and $N=10$ with $d=0.908$ (red) of the $9$th
harmonic. (d) Gate phase $\theta(t)$ for $N=5$ (blue), $N=10$ (red), and the
dispersive case (dashed).

Figure 2 (b) shows values of $t_{\rm g}$ and $\Omega_{\rm pp}$ obtained
following the previous prescription. In particular, these have been calculated
for a system with $\eta=0.005$ and for two different harmonics $k=9$ (square
panels) and $k=7$ (round panels). Notice that, for $\eta=0.005$, there are
plenty of solutions giving faster-than-dispersive gates, i.e. $t_{\rm
g}/t_{\rm g}^{\rm d}<1$, using low-intensity pulses with values of
$\Omega_{\rm pp}$ well below the mode frequency $\nu$.

As an example, we choose two solutions within the $9$th harmonic, where
$\tau_{k}$ extends over approximately $9$ oscillator periods. In Fig. 2 (c)
the shapes of $\Omega(t)$ and $f_{z}(t)$ are displayed for cases $N=5$ and
$N=10$. Notice that $\Omega(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 $\theta(t)$ related to the $N=5$ gate (dashed blue)
reaches the target value $\pi/8$ faster than the $N=10$ gate (solid red) or
the dispersive gate (dotted line).

The reason for choosing the gate time as an integer multiple of $8\tau_{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$\equiv$XYXYYXYX pulse structure assures cancelation of $\sigma^{z}$ type
noise, as well as of Rabi frequency fluctuations. In Fig. 2(a) the Rabi
frequency is plotted ($\Omega_{x}(t)$ when “X”; $\Omega_{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 $H_{\rm
XYXY}=\frac{1}{2}\eta^{2}\nu\Big{\\{}-J_{k}^{\perp}(S_{x}^{2}+S_{y}^{2})+[B_{k}^{\prime}(a^{\dagger}a+1)+B_{k}^{\prime\prime}]S_{z}\Big{\\}}$,
see SM SupMat . Thus, if our gate contains only XYXY blocks, $H_{\rm 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}^{\prime}(a^{\dagger}a+1)+B_{k}^{\prime\prime}]S_{z}$ gets refocused.
Secondly, we cancel the remaining term $J_{k}^{\perp}(S_{x}^{2}+S_{y}^{2})$ by
driving the two qubits with opposite phases every second XY8 block. When
rotating by an angle $\pi$ the phase of the second qubit’s driving
$H_{d}^{(+)}=\sum_{\mu=x,y}\Omega_{\mu}(t)S^{(+)}_{\mu}/2$ becomes
$H_{d}^{(-)}=\sum_{\mu=x,y}\Omega_{\mu}(t)S^{(-)}_{\mu}/2$ instead, where
$S_{\mu}^{(\pm)}\equiv\sigma_{1}^{\mu}\pm\sigma_{2}^{\mu}$. This changes the
sign of the $\sigma_{1}^{u}\sigma_{2}^{\mu}$ terms with $u,\mu\in\\{x,y\\}$,
leading to refocusing of the term $J_{k}^{\perp}(S_{x}^{2}+S_{y}^{2})$ after
every pair of XY8 blocks.

Summarizing, our two-qubit gates are generated by nesting
TQXY16$\equiv$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 $\pi$ pulse is implemented
according to the designs for $f_{z}(t)$ and $\beta(t)$ presented in Eqs. (8,
9).

_Numerical results.-_ We benchmark our method by simulating its performance in
a pair of trapped ions in a static magnetic field gradient Mintert01 . In this
scenario, $\nu=(2\pi)\times 220$ kHz is the frequency of the centre-of-mass
vibrational mode, $\eta=\gamma_{e}g_{B}/8\nu\sqrt{\hbar/M\nu}$ is an effective
Lamb-Dicke factor where $\gamma_{e}=(2\pi)\times 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_{\rm 2M}$, where $H_{\rm
2M}=\sqrt{3}\nu b^{\dagger}b-3^{-1/4}\eta\nu(b+b^{\dagger})S_{z}^{(-)}$. The
addition of $H_{\rm 2M}$ changes the dispersive coupling in Eq. (4) as
$J_{k}\rightarrow J_{k}-1/3\sum_{n=1}^{\infty}f_{n}^{2}/(1-n^{2}/3k^{2})$,
which must be taken into account when following the prescription to calculate
the valid gates. Although we simulate the performance of the gate with the
full Hamiltonian $H_{f}=H^{(\pm)}_{d}+H_{0}+H_{\rm 2M}$, due to computational
limitations we use (unless otherwise stated $H_{s}=H_{d}^{(\pm)}+H_{0}+H_{\rm
2M}^{\rm eff}$ instead, where $H_{\rm 2M}^{\rm
eff}=\frac{1}{3}\nu\eta^{2}rS_{z}^{2}$ is the second-order contribution of
$H_{\rm 2M}$ SupMat . Here, $H_{d}^{(\pm)}$ stands for $H_{d}^{(+)}$
($H_{d}^{(-)}$) every first (second) half of a TQXY16 block.

We investigate two regimes: (i) $\eta=0.005$ ($g_{B}=19.16\,\rm T/m$), which
is the state-of-the-art of current experiments Webb18 ; Barthel22 , and (ii)
$\eta=0.04$ ($g_{B}=153.2\,\rm T/m$), which can be reached in near future
setups Weidt16 .

In regime (i), we consider three different gates, all within the $9$th
harmonic. The first gate (G1), with a duration $t_{\rm g}=1.641$ ms, appears
after five TQXY16 blocks with pulse length $t_{\pi}=20.51\ \mu$s reaching
$\Omega_{\rm pp}=(2\pi)\times 124.5$ kHz. The second gate (G2) with gate time
$t_{\rm g}=3.277$ ms uses ten TQXY16 blocks with pulse length $t_{\pi}=20.48\
\mu$s reaching $\Omega_{\rm pp}=(2\pi)\times 77.85$ kHz. The third gate (G3)
uses twelve TQXY16 blocks, each with a different pulse length and detuning,
while it incorporates a technique to mitigate errors due to mode decoherence,
see SM SupMat . The list of detunings is $\vec{\xi}_{k}=(2\pi)\times[$1.24,
0.31, 0.64, 0.09, 0.55, 0.05, 0.54, 0.06, 0.57, 0.14, 0.73, 0.80$]$ kHz, the
gate-time is $t_{\rm g}=3.936$ ms, and $\Omega_{\rm pp}=(2\pi)\times 78.69$
kHz. In regime (ii) we consider a gate within the $5$th harmonic. This gate
occurs after two TQXY16 blocks where $t_{\rm g}=368~{}\mu$s,
$t_{\pi}=11.51~{}\mu$s, and $\Omega_{\rm pp}=(2\pi)\times 80.88$ kHz. For
further details regarding pulse parameters, see Appendix A.

The performance of the four gates in the presence of distinct error sources is
shown in table 1. Each simulated experiment starts from the state
$|\\!+_{x}+_{y}\rangle$ and targets the Bell-state
$|\tilde{\Phi}^{+}\rangle=\frac{1}{\sqrt{2}}(|\\!+_{x}+_{y}\rangle+i|\\!-_{x}-_{y}\rangle)$,
while in all cases we consider an initial motional thermal state with
$\bar{n}=1$.

In the 2nd and 3rd columns of table 1 we show the Bell-state infidelity
$\mathcal{I}=1-\mathcal{F}$ obtained by concatenating XY8 or TQXY16 blocks,
respectively. Notice that TQXY16 blocks achieve a clearly superior performance
due to efficient decoupling from finite pulse effects. For these,
$\mathcal{I}_{\rm TQXY16}\leq 10^{-6}$ for all gates except G4, where finite
pulse effects limit the infidelity to approximatelly $10^{-5}$. In the fourth
column we evaluate the effect of the second mode $b$ by numerically simulating
$H_{f}$ (initialising the second mode $b$ in a thermal state with
$\bar{n}=1$), which results in $\mathcal{I}_{\rm 2M}$. The infidelities
relative to the previous case (i.e. $\mathcal{I}_{\rm 2M}-\mathcal{I}_{\rm
TQXY16}$) are given in the “2M” column of table 1. Again, the effect of the
second mode is relevant only for G4, which contributes $3.8\times 10^{-5}$ to
the total infidelity. Importantly, this demonstrates that our gate is
compatible with the presence of spectator modes.

To investigate the effect of crosstalk, we add the term
$H^{(\pm)}_{c}=\sum_{\mu=x,y}\frac{\Omega_{\mu}(t)}{2}(\sigma_{2}^{-}e^{-i\Delta\omega
t}\pm\sigma_{1}^{-}e^{i\Delta\omega t}+{\rm H.c.})$ to $H_{s}$, where
$\Delta\omega/(2\pi)=2.54$ and $20.34$ MHz for regimes (i) and (ii),
respectively. The results are given in the “CT” column of table 1. To reduce
the impact of crosstalk, we combine our pulses with sin2-shaped ramps at the
beginning and end of each pulse, see SM SupMat , and optimize the length of
the ramp using numerical simulations. The resulting infidelities are shown in
the “CT*” column. Note that the sin2 ramp reduces the value of $\mathcal{I}$
by at least an order of magnitude in most cases.

Robustness w.r.t. common errors such as dephasing over qubits due to
$\omega_{j}\rightarrow(1\pm\delta\omega)\,\omega_{j}$ (with $T_{2}^{*}\approx
500\mu$s Barthel22 ), Rabi-frequency fluctuations (i.e.
$\Omega(t)\rightarrow(1\pm\delta_{\Omega})\,\Omega(t)$ with
$\delta_{\Omega}=5\times 10^{-3}$, and shifts on the mode frequency,
$\nu\rightarrow(1\pm\delta_{\nu})\,\nu$, with $\delta_{\nu}=10^{-5}$ Johnson16
, is shown in columns 7-9 of table 1. Furthermore, in column $\dot{\bar{n}}$
we give the relative infidelities for G1-3 and G4 for mode heating rates
$\dot{\bar{n}}=35$ and $100$ ph/s, respectively. For further details, see
Appendix B and SM SupMat .

Table 1 shows that mode heating is is the main source of infidelity for gates
in regime (i). This, along with dephasing and crosstalk, limits the fidelity
of these gates to the $10^{-3}$ regime. Interestingly, despite its longer gate
duration, G3 achieves better performance in terms of motion-induced errors
than G1 and G2, proving the validity of the mode decoherence protecting
technique described in the SM SupMat . In addition, G3 also shows the best
overall performance, achieving $99.91\%$ fidelity within current experimental
limits. We note that this represents an improvement of more than an order of
magnitude w.r.t the current record Weidt16 . Finally, table 1 shows that G4 is
the most robust w.r.t. experimental imperfections. This is reasonable since it
uses a larger $\eta$ and is an order of magnitude faster than the other gates,
leading to a total infidelity of $1.14\times 10^{-4}$ which represents an
improvement of two orders of magnitude w.r.t. the state of the art.

Table 1: Error budget: Column XY8 and TQXY16 show the infidelities after an evolution with Hamiltonians $H_{s}=H_{d}^{(+)}+H_{0}+H_{\rm 2M}^{\rm eff}$ and $H_{s}$, respectively. The remaining columns show infidelities relative to the TQXY16 case (e.g. $\mathcal{I}_{\rm 2M}-\mathcal{I}_{\rm TQXY16}$), taking into account various experimental imperfections. In columns 2M, CT, and CT∗, infidelities obtained considering a second mode, crosstalk, and crosstalk with the sin2 ramp are shown. In $\delta\Omega$, $\delta\nu$, and $T_{2}^{*}$ we show relative infidelities considering static shifts of $\delta\Omega=5\times 10^{-3}$, $\delta\nu=10^{-5}$, and $\delta\omega=2\sqrt{2}$ kHz. $\dot{\bar{n}}$ shows the infidelity considering heating with rates $\dot{\bar{n}}_{1}=35$ ph/s and $\dot{\bar{n}}_{2}=100$ ph/s for regimes (i) and (ii), respectively. The last column shows the overall infidelity obtained by summing the values of all columns except those in XY8 and CT.  $\mathcal{I}\,{(10^{-4})}$ | XY8 | TQXY16 | 2M | CT | CT* | $T_{2}^{*}$ | $\delta\Omega$ | $\delta\nu$ | $\dot{\bar{n}}$ | Total
---|---|---|---|---|---|---|---|---|---|---
G1 | $5.50$ | $0.01$ | $0.04$ | $24.8$ | $2.26$ | $2.34$ | $0.21$ | 0.28 | $5.65$ | $10.8$
G2 | $28.7$ | $<\\!\\!10^{-2}$ | $<\\!0.01$ | $3.20$ | $0.95$ | $2.45$ | $0.32$ | $1.01$ | $19$ | $23.7$
G3 | $41.3$ | $<\\!\\!10^{-2}$ | $<\\!0.01$ | $134$ | $1.82$ | $2.35$ | $0.31$ | $0.26$ | $4.71$ | $9.45$
G4 | $>\\!\\!10^{3}$ | $0.12$ | $0.38$ | $0.23$ | $0.01$ | $0.43$ | $0.09$ | $<\\!\\!10^{-2}$ | $0.11$ | $1.14$

_Conclusions.-_ We present a DD sequence (TQXY16) based on the delivery of
low-intensity $\pi$ pulses that achieve faster-than-dispersive two-qubit
gates. Without the need of any numerical optimisation, we desing 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 demonstrate the versatility of our protocol to
adopt forms that provide an increased robustness against crosstalk and mode
decoherence. We exemplified the method performance in trapped ions setups
obtaining infidelities within the $10^{-3}$ threshold at state-of-the-art
conditions, and in the $10^{-4}$ regime in near-future setups.

###### Acknowledgements.

_Acknowledgements.-_ I. A. would like to thank P. Rabl and J. S. Pedernales
for useful discussions. I. A. acknowledges support from the European Union’s
Horizon2020 research and innovation programme under Grant Agreement No. 899354
(SuperQuLAN). J. C. acknowledges the Ramón y Cajal (RYC2018-025197-I) research
fellowship, the financial support from Spanish Government via EUR2020-112117
and Nanoscale NMR and complex systems (PID2021-126694NB-C21) projects, the EU
FET Open Grant Quromorphic (828826), the ELKARTEK project Dispositivos en
Tecnologías Cuánticas (KK-2022/00062), and the Basque Government grant
IT1470-22.

## Appendix A: Values for pulse parameters

Table 2: Suitable values the for $b$, $c$ and $|d_{\rm max}|$. $k$ | $3$ | $5$ | $7$ | $9$ | $11$ | $13$ | $15$
---|---|---|---|---|---|---|---
$b$ | $0.33$ | $0.30$ | $0.29$ | $0.33$ | $0.34$ | $0.35$ | $0.30$
$c$ | $0.035$ | $0.04$ | $0.05$ | $0.042$ | $0.03$ | $0.035$ | $0.035$
$d_{\rm max}$ | $-2.3$ | $-1.5$ | $2.4$ | $2.3$ | $1.9$ | $1.7$ | $2.3$

In table 2 we show suitable values of $b$ and $c$ given the harmonic $k$.
Also, we show the maximum value of $|d|$ for which a physical pulse (i.e,
$|f_{z}|(t)\leq 1$) can still be generated with the ansatz given in Eqs. (8,
9). For gates G1, G2 and G4, the selected values of $d$ are $1.915$, $0.933$
and $-0.321$, respectively.

## Appendix B: Behaviour of G1-4 regarding distinct error rates

Fig. A1:  (a) Infidelity $\mathcal{I}$ versus $1/T_{2}^{*}$ for the gates G1
(blue), G2 (red), G3 (purple) and G4 (green). In (b) and (c), the infidelity
is shown for Rabi frequency shifts
$\Omega(t)\rightarrow(1\pm\delta_{\Omega})\Omega(t)$ and in the mode frequency
$\nu\rightarrow(1\pm\delta_{\nu})\nu$. (d) Infidelity under heating for
different rates $\dot{\bar{n}}$. (e) Trajectory of $\alpha(t)$ for gates G1
(circular, blue), G2 (circular, red), G3 (non-circular, purple), and G4
(circular, green).

Here, we show the performance of our gate w.r.t. different error rates. More
specifically, in Fig. A1(a), we plot $\mathcal{I}$ w.r.t. the dephasing time
$T_{2}^{*}$. In particular, our simulations consider an additional
$\pm\delta\omega/2S_{z}$ term in $H_{s}$, where
$\delta\omega=\sqrt{2}/T_{2}^{*}$ and each point is the average infidelity
obtained by a positive and a negative displacement. Here, we assume that qubit
frequency fluctuations are static, which requires that the correlation-time of
the noise $\tau_{c}$ is much larger than the period $\tau_{k}$. In Fig. A1(b)
and (c) we show $\mathcal{I}$ w.r.t. static shifts in the Rabi frequency
$\Omega(t)\rightarrow(1\pm\delta_{\Omega})\Omega(t)$ and the mode frequency,
$\nu\rightarrow(1\pm\delta_{\nu})\nu$. In Fig. A1 (d), the effect of motional
heating is shown. For this, we calculate the final state solving a master
equation for different heating rates $\dot{\bar{n}}$, see SM SupMat . In Fig.
A1 (e), the phase-space trajectory of $\alpha(t)$ is shown for gates G1-4.

## References

  * (1) M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information (Cambridge University press, Cambridge, 2000).
  * (2) P.-Q. Jin, M. Marthaler, A. Shnirman, and G. Schön, Strong Coupling of Spin Qubits to a Transmission Line Resonator, Phys. Rev. Lett. 108, 190506 (2012).
  * (3) P.-M. Billangeon, J. S. Tsai, and Y. Nakamura, Circuit-QED-based scalable architectures for quantum information processing with superconducting qubits, Phys. Rev. B 91, 094517 (2015).
  * (4) S. Richer and D. DiVincenzo, Circuit design implementing longitudinal coupling: A scalable scheme for superconducting qubits, Phys. Rev. B 93, 134501 (2016).
  * (5) F. Beaudoin, D. Lachance-Quirion, W. A. Coish, and M. Pioro-Ladrière, Coupling a single electron spin to a microwave resonator: controlling transverse and longitudinal couplings, Nanotechnology 27, 464003 (2016).
  * (6) F. Mintert, and C. Wunderlich, Ion-trap quantum logic using long-wavelength radiation, Phys. Rev. Lett. 87, 257904 (2001).
  * (7) A. Sørensen and K. Mølmer, Quantum computation with ions in thermal motion, Phys. Rev. Lett. 82, 1971 (1999).
  * (8) E. Solano, R. L. de Matos Filho, and N. Zagury, Deterministic Bell states and measurement of the motional state of two trapped ions, Phys. Rev. A 59, R2539(R) (1999).
  * (9) A. Sørensen and K. Mølmer, Entanglement and quantum computation with ions in thermal motion, Phys. Rev. A 62, 022311 (2000).
  * (10) C. A. Sackett, D. Kielpinski, B. E. King, C. Langer, V. Meyer, C. J. Myatt, M. Rowe, Q. A. Turchette, W. M. Itano, D. J. Wineland, and C. Monroe, Experimental entanglement of four particles, Nature 404, 256 (2000).
  * (11) G. J. Milburn, S. Schneider, and D. F. V. James, Ion Trap Quantum Computing with Warm Ions, Fortschr. Phys. 48, 801 (2000).
  * (12) D. Leibfried, B. DeMarco, V. Meyer, D. Lucas, M. Barrett, J. Britton, W. M. Itano, B. Jelenkovic, C. Langer, T. Rosenband, and D. J. Wineland, Experimental demonstration of a robust, high-fidelity geometric two ion-qubit phase gate, Nature 422, 412 (2003).
  * (13) S.-B. Zheng, Unconventional geometric quantum phase gates with a cavity QED system, Phys. Rev. A 70, 052320 (2004).
  * (14) C. H. Valahu, I. Apostolatos, S. Weidt, and W. K. Hensinger, Quantum control methods for robust entanglement of trapped ions, J. Phys. B: At. Mol. Opt. Phys. 55 204003 (2022).
  * (15) A. W. Cross and J. M. Gambetta, Optimized pulse shapes for a resonator-induced phase gate, Phys. Rev. A 91, 032325 (2015).
  * (16) F. Haddadfarshi and F. Mintert, High fidelity quantum gates of trapped ions in the presence of motional heating, New J. Phys. 18, 123007 (2016).
  * (17) A. E. Webb, S. C. Webster, S. Collingbourne, D. Bretaud, A. M. Lawrence, S. Weidt, F. Mintert, and W. K. Hensinger, Resilient Entangling Gates for Trapped Ions, Phys. Rev. Lett. 121, 180501 (2018).
  * (18) Y. Shapira, R. Shaniv, T. Manovitz, N. Akerman, and R. Ozeri, Robust Entanglement Gates for Trapped-Ion Qubits, Phys. Rev. Lett. 121, 180502 (2018).
  * (19) G. Zarantonello, H. Hahn, J. Morgner, M. Schulte, A. Bautista-Salvador, R. F. Werner, K. Hammerer, and C. Ospelkaus, Robust and resource-efficient microwave near-field entangling ${}^{9}{\rm Be}^{+}$ gate, Phys. Rev. Lett. 123, 260503 (2019).
  * (20) R. T. Sutherland, R. Srinivas, S. C. Burd, H. M. Knaack, A. C. Wilson, D. J. Wineland, D. Leibfried, D. T. C. Allcock, D. H. Slichter, and S. B. Libby, Laser-free trapped-ion entangling gates with simultaneous insensitivity to qubit and motional decoherence, Phys. Rev. A 101, 042334 (2020).
  * (21) D. Hayes, S. M. Clark, S. Debnath, D. Hucul, I. V. Inlek, K. W. Lee, Q. Quraishi, and C. Monroe, Coherent Error Suppression in Multiqubit Entangling Gates, Phys. Rev. Lett. 109, 020503 (2012).
  * (22) T. Choi, S. Debnath, T. A. Manning, C. Figgatt, Z.-X. Gong, L.-M. Duan, and C. Monroe, Optimal Quantum Control of Multimode Couplings between Trapped Ion Qubits for Scalable Entanglement, Phys. Rev. Lett. 112, 190502 (2014).
  * (23) Todd J. Green and Michael J. Biercuk, Phase-Modulated Decoupling and Error Suppression in Qubit-Oscillator Systems, Phys. Rev. Lett. 114, 120502 (2015).
  * (24) P. H. Leung, K. A. Landsman, C. Figgatt, N. M. Linke, C. Monroe, and K. R. Brown, Robust 2-Qubit Gates in a Linear Ion Crystal Using a Frequency-Modulated Driving Force, Phys. Rev. Lett. 120, 020501 (2018).
  * (25) V. M. Schäfer, C. J. Ballance, K. Thirumalai, L. J. Stephenson, T. G. Ballance, A. M. Steane, and D. M. Lucas, Fast quantum logic gates with trapped-ion qubits, Nature 555, 75 (2018).
  * (26) Y. Lu, S. Zhang, K. Zhang, W. Chen, Y. Shen, J. Zhang, J.-N. Zhang, and K. Kim, Global entangling gates on arbitrary ion qubits, Nature 572, 363 (2019).
  * (27) A. R. Milne, C, L. Edmunds, C. Hempel, F. Roy, S. Mavadia, and M. J. Biercuk, Phase-Modulated Entangling Gates Robust to Static and Time-Varying Errors, Phys. Rev. Applied 13, 024022 (2020).
  * (28) S. Puri and A. Blais, High-Fidelity Resonator-Induced Phase Gate with Single-Mode Squeezing, Phys. Rev. Lett. 116, 180501 (2016).
  * (29) B. Royer, A. L. Grimsmo, N. Didier, and A. Blais, Fast and High-Fidelity Entangling Gate through Parametrically Modulated Longitudinal Coupling, Quantum 1, 11 (2017).
  * (30) W. Ge, B. C. Sawyer, J. W. Britton, K. Jacobs, J. J. Bollinger, and M. Foss-Feig, Trapped Ion Quantum Information Processing with Squeezed Phonons, Phys. Rev. Lett. 122, 030501 (2019).
  * (31) S. C. Burd, R. Srinivas, H. M. Knaack, W. Ge, A. C. Wilson, D. J. Wineland, D. Leibfried, J. J. Bollinger, D. T. C. Allcock, and D. H. Slichter, Quantum amplification of boson-mediated interactions, Nat. Phys., 17, 898 (2021).
  * (32) Y. Shapira, S. Cohen, N. Akerman, A. Stern, and R. Ozeri, Robust two-qubit trapped ions gates using spin-dependent squeezing, arXiv:2207.01660 (2022).
  * (33) L. Viola and S. Lloyd, Dynamical suppression of decoherence in two-state quantum systems, Phys. Rev. A 58, 2733 (1998).
  * (34) M. Ban, Photon-echo technique for reducing the decoherence of a quantum bit, J. Mod. Opt. 45 2315 (1998).
  * (35) N. Timoney, I. Baumgart, M. Johanning, A. F. Varón, M. B. Plenio, A. Retzker, and C. Wunderlich, Quantum gates and memory using microwave-dressed states, Nature 476, 185 (2011).
  * (36) T. R. Tan, J. P. Gaebler, R. Bowler, Y. Lin, J. D. Jost, D. Leibfried, and D. J. Wineland, Demonstration of a Dressed-State Phase Gate for Trapped Ions, Phys. Rev. Lett. 110, 263002 (2013).
  * (37) T. P. Harty, M. A. Sepiol, D. T. C. Allcock, C. J. Ballance, J. E. Tarlton, and D. M. Lucas, High-fidelity trapped-ion quantum logic using near-field microwaves, Phys. Rev. Lett. 117, 140501 (2016).
  * (38) S. Weidt, J. Randall, S. C. Webster, K. Lake, A. E. Webb, I. Cohen, T. Navickas, B. Lekitsch, A. Retzker, and W. K. Hensinger,Trapped-ion quantum logic with global radiation fields, Phys. Rev. Lett. 117, 220501 (2016).
  * (39) Q. Guo, S.-B. Zheng, J. Wang, C. Song, P. Zhang, K. Li, W. Liu, H. Deng, K. Huang, D. Zheng, X. Zhu, H. Wang, C.-Y. Lu, and J.-W. Pan, Dephasing-Insensitive Quantum Information Storage and Processing with Superconducting Qubits, Phys. Rev. Lett. 121, 130501 (2018).
  * (40) H. Y. Carr and E. M. Purcell, Effects of diffusion on free precession in nuclear magnetic resonance experiments, Phys. Rev. 94, 630 (1954).
  * (41) S. Meiboom and D. Gill, Modified spin-echo method for measuring nuclear relaxation times, Rev. Sci. Instrum. 29, 688 (1958).
  * (42) T. Gullion, D. B. Baker, and M. S. Conradi, New, compensated carr-purcell sequences, J. Magn. Reson. (1969) 89, 479 (1990).
  * (43) A. M. Souza, G. A. Álvarez, and D. Suter, Robust dynamical decoupling, Phil. Trans. R. Soc. A. 370, 4748 (2012).
  * (44) C. Kabytayev, T. J. Green, K. Khodjasteh, M. J. Biercuk, L. Viola, and K. R. Brown, Robustness of composite pulses to time-dependent control noise, Phys. Rev. A 90, 012316 (2014).
  * (45) C. Munuera-Javaloy, I. Arrazola, E. Solano, and J. Casanova, Double quantum magnetometry at large static magnetic fields, Phys. Rev. B 101, 104411 (2020).
  * (46) L. Dong, I. Arrazola, X. Chen, and J. Casanova, Phase-Adaptive Dynamical Decoupling Methods for Robust Spin-Spin Dynamics in Trapped Ions, Phys. Rev. Applied 15, 034055 (2021).
  * (47) N. Ezzell, B. Pokharel, L. Tewala, G. Quiroz, and D. A. Lidar, Dynamical decoupling for superconducting qubits: a performance survey, arXiv:2207.03670 (2022).
  * (48) J. Casanova, Z.-Y. Wang, J. F. Haase, and M. B. Plenio, Robust dynamical decoupling sequences for individual-nuclear-spin addressing, Phys. Rev. A 92, 042304 (2015).
  * (49) J. Casanova, Z.-Y. Wang, and M. B. Plenio, Arbitrary nuclear-spin gates in diamond mediated by a nitrogen-vacancy-center electron spin, Phys. Rev. A 96, 032314 (2017).
  * (50) I. Arrazola, J. Casanova, J. S. Pedernales, Z.-Y. Wang, E. Solano, and M. B. Plenio, Pulsed dynamical decoupling for fast and robust two-qubit gates on trapped ions, Phys. Rev. A 97, 052312 (2018).
  * (51) Ch. Piltz, B. Scharfenberger, A. Khromova, A.F. Varón, and C. Wunderlich, Protecting conditional quantum gates by robust dynamical decoupling, Phys. Rev. Lett. 110, 200501 (2013).
  * (52) J. Qiu, Y. Zhou, C.-K. Hu, J. Yuan, L. Zhang, J. Chu, W. Huang, W. Liu, K. Luo, Z. Ni, X. Pan, Z. Yang, Y. Zhang, Y. Chen, X.-H. Deng, L. Hu, J. Li, J. Niu, Y. Xu, T. Yan, Y. Zhong, S. Liu, F. Yan, and D. Yu, Suppressing Coherent Two-Qubit Errors via Dynamical Decoupling, Phys. Rev. Applied 16, 054047 (2021).
  * (53) P. Barthel, P. H. Huber, J. Casanova, I. Arrazola, D. Niroomand, T. Sriarunothai, M. B. Plenio, and C. Wunderlich, Robust Two-Qubit Gates Using Pulsed Dynamical Decoupling, arXiv:2208.00187 (2022).
  * (54) W. Morong, K. S. Collins, A. De, E. Stavropoulos, T. You, and C. Monroe, Engineering dynamically decoupled quantum simulations with trapped ions, arXiv:2209.05509 (2022).
  * (55) H. Paik, A. Mezzacapo, M. Sandberg, D. T. McClure, B. Abdo, A. D. Córcoles, O. Dial, D. F. Bogorin, B. L. T. Plourde, M. Steffen, A. W. Cross, J. M. Gambetta, and J. M. Chow, Experimental Demonstration of a Resonator-Induced Phase Gate in a Multiqubit Circuit-QED System, Phys. Rev. Lett. 117, 250502 (2016).
  * (56) C. J. Ballance, T. P. Harty, N. M. Linke, M. A. Sepiol, and D. M. Lucas, High-fidelity quantum logic gates using trapped-ion hyperfine qubits, Phys. Rev. Lett. 117, 060504 (2016).
  * (57) O. Băzăvan, S. Saner, M. Minder, A. C. Hughes, R. T. Sutherland, D. M. Lucas, R. Srinivas, and C. J. Ballance, Synthesizing a $\hat{\sigma}_{z}$ spin-dependent force for optical, metastable, and ground state trapped-ion qubits, arXiv:2207.11193 (2022).
  * (58) T. Manovitz, A. Rotem, R. Shaniv, I. Cohen, Y. Shapira, N. Akerman, A. Retzker, and R. Ozeri, Fast Dynamical Decoupling of the Mølmer-Sørensen Entangling Gate, Phys. Rev. Lett. 119, 220505 (2017).
  * (59) P. Rabl, S. J. Kolkowitz, F. H. L. Koppens, J. G. E. Harris, P. Zoller, and M. D. Lukin, A quantum spin transducer based on nanoelectromechanical resonator arrays, Nat. Phys. 6, 602 (2010).
  * (60) E. Rosenfeld, R. Riedinger, J. Gieseler, M. Schuetz, and M. D. Lukin, Efficient Entanglement of Spin Qubits Mediated by a Hot Mechanical Oscillator, Phys. Rev. Lett. 126, 250505 (2021).
  * (61) B. Lekitsch, S. Weidt, A. G. Fowler, K. Mølmer, S. J. Devitt, C. Wunderlich, and W. K. Hensinger, Blueprint for a microwave trapped ion quantum computer, Sci. Adv. 3, e1601540 (2017).
  * (62) J. J. García-Ripoll, P. Zoller, and J. I. Cirac, Speed Optimized Two-Qubit Gates with Laser Coherent Control Techniques for Ion Trap Quantum Computing, Phys. Rev. Lett. 91, 157901 (2003).
  * (63) L.-M. Duan, Scaling Ion Trap Quantum Computation through Fast Quantum Gates, Phys. Rev. Lett. 93, 100502 (2004).
  * (64) A. M. Steane, G. Imreh, J. P. Home, and D. Leibfried, Pulsed force sequences for fast phase-insensitive quantum gates in trapped ions, New J. Phys. 16, 053049 (2014).
  * (65) C. D. B. Bentley, A. R. R. Carvalho, D. Kielpinski, and J. J. Hope, Fast gates for ion traps by splitting laser pulses, New J. Phys. 15, 043006 (2013).
  * (66) Mahdi Sameti, Jake Lishman, and Florian Mintert, Strong-coupling quantum logic of trapped ions, Phys. Rev. A 103, 052603 (2021).
  * (67) X. Wang, C.-S. Yu, and X. Yi, An alternative quantum fidelity for mixed states of qubits, Phys. Lett. A 373, 58 (2008).
  * (68) Supplemental Material
  * (69) K. G. Johnson, J. D. Wong-Campos, A. Restelli, K. A. Landsman, B. Neyenhuis, J. Mizrahi, and C. Monroe, Active stabilization of ion trap radiofrequency potentials, Review of Scientific Instruments 87, 053110 (2016).

Supplemental Material: Robust Oscillator-Mediated Phase Gates Driven by Low-
Intensity Pulses

##  I Fourier coefficients $f_{n}$ and dispersive coupling $J_{k}$

In this section, we give analytical expressions for the Fourier coefficients
$f_{n}$ for the cases of instantaneous, top-hat, and modulated pulses. For the
case of instantaneous pulses, $f^{\rm\,i}_{n}=\frac{4}{n\pi}\sin{(n\pi/2)}$.
For top-hat pulses,

$f^{\rm th}_{n}=\frac{4\sin{(n\pi/2)}\cos{(n\pi
t_{\pi}/\tau)}}{n\pi(1-n^{2}t_{\pi}^{2}/4\tau^{2})}$ (S1)

where $t_{\pi}$ is the duration of the $\pi$-pulse, and $\tau$ is the period
of the modulation function $f_{z}(t)$. For the modulated case, the expressions
of the Fourier coefficients are

$f^{\rm m}_{1}=1-\frac{4d\sin(\pi k/2)}{\pi kb}\Bigg{\\{}\frac{e^{-[(k-1)\pi
c]^{2}/4}}{(k-1)\pi}\sin{[(k-1)\pi b]}-\frac{e^{-[(k+1)\pi
c]^{2}/4}}{(k+1)\pi}\sin{[(k+1)\pi b]}\Bigg{\\}}$ (S2)

for $n=1$, $f^{\rm m}_{k}=-4d/\pi k$ for $n=k$, and

$f^{\rm m}_{n}=-\frac{4d\sin{(n\pi/2)}\sin(\pi k/2)}{\pi
kb}\Bigg{\\{}\frac{e^{-[(k-n)\pi c]^{2}/4}}{(k-n)\pi}\sin{[(k-n)\pi
b]}-\frac{e^{-[(k+n)\pi c]^{2}/4}}{(k+n)\pi}\sin{[(k+n)\pi b]}\Bigg{\\}}.$
(S3)

for the rest. Here, $k$ is the value of the selected harmonic, and $b,c$ and
$d$ are pulse parameters. Suitable values of these can be found in Appendix A.
For the single mode case, the dispersive coupling is

$J_{k}=f_{k}^{2}/4+\sum_{n\neq k}^{\infty}f_{n}^{2}/(1-n^{2}/k^{2}),$ (S4)

and we find that its value is sufficiently accurate truncating the sum at
$2k$.

###  I.1 Derivations

A periodic modulation function $f_{z}(t)$ can be written as a sum of infinite
Fourier harmonics

$f_{z}(t)=\sum_{n=1}^{\infty}f_{n}\cos{(n\omega t)},$ (S5)

where $\omega=2\pi/\tau$ and
$f_{n}=\frac{2}{\tau}\int_{0}^{\tau}dt^{\prime}f_{z}(t^{\prime})\cos{(n\omega
t^{\prime})}$. To calculate the coefficients $f_{n}$, the change of variable
$t^{\prime}=x\tau/2$ is appropiate. The expression then changes to

$f_{n}=\int_{0}^{2}dx\ f_{z}(x)\cos{(n\pi x)},$ (S6)

where $f_{z}(x)$ has now periodicity of $x_{\tau}=2$. Because $f_{z}(x)$ is
always symmetric with respect to (w.r.t.) the $x=1$ point, i.e.
$f_{z}(2-x)=f_{z}(x)$, we divide the integral in two parts

$f_{n}=\int_{0}^{1}dx\ f_{z}(x)\cos{(n\pi x)}+\int_{1}^{2}dx\
f_{z}(x)\cos{(n\pi x)}.$ (S7)

With a change of variable $x\rightarrow 2-x$ in the second integral, and using
that $\cos{(\pi n(2-x))}=\cos{\pi nx}$, this results in

$f_{n}=2\int_{0}^{1}f(x)\cos{(\pi nx)}\ dx.$ (S8)

Moreover, $f_{z}(x)$ is also antisymmetric w.r.t. the point $x=1/2$, i.e.
$f_{z}(1-x)=-f_{z}(x)$. We further divide the integral,

$f_{n}=2\int_{0}^{1/2}dx\ f_{z}(x)\cos{(n\pi x)}+2\int_{1/2}^{1}dx\
f_{z}(x)\cos{(n\pi x)},$ (S9)

and, with the change of variable $x\rightarrow 1-x$ in the second integral, we
get

$f_{n}=2[1-\cos{(n\pi)}]\int_{0}^{1/2}dx\ f_{z}(x)\cos{(n\pi x)}.$ (S10)

If $n$ is even, value of the integral is zero. In $n$ is odd, $f_{n}$ is

$f_{n}=4\int_{0}^{1/2}dx\ f_{z}(x)\cos{(n\pi x)},$ (S11)

and its value will depend on the form of $f_{z}(x)$. For convenience, we make
a last change of variable $x\rightarrow 1/2-x$ leading to

$f_{n}=4\sin{(n\pi/2)}\int_{0}^{1/2}dx\ f_{z}(x)\sin{(n\pi x)},$ (S12)

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 $\pi$
pulses.

Instantaneous Pulses: For instantaneous pulses, $f_{z}(x)=1$, thus,
$f_{n}=\frac{4}{n\pi}\sin{(n\pi/2)}$.

Top-Hat Pulses: For top-hat pulses, $f_{z}(x)=\sin{(\pi x/\tau_{\pi})}$ from
$x=0$ to $x=\tau_{\pi}/2$, while $f_{z}(x)=1$ for the rest. Here,
$\tau_{\pi}=2t_{\pi}/\tau=2\pi/\Omega\tau$. We divide the integral
accordingly,

$f_{n}/4\sin{(n\pi/2)}=\int_{0}^{\tau_{\pi}/2}dx\ \sin{(\pi
x/\tau_{\pi})}\sin{(n\pi x)}+\int_{\tau_{\pi}/2}^{1/2}dx\sin{(n\pi x)}.$ (S13)

For the first part, we use that $2\sin A\sin B=\cos{(A-B)}-\cos{(A+B)}$, and,
then,

$\int_{0}^{\tau_{\pi}/2}dx\ \sin{(\pi x/\tau_{\pi})}\sin{(n\pi
x)}=\frac{\sin{[\pi(\tau_{\pi}^{-1}-n)\tau_{\pi}/2]}}{2\pi(\tau_{\pi}^{-1}-n)}-\frac{\sin{[\pi(\tau_{\pi}^{-1}+n)\tau_{\pi}/2]}}{2\pi(\tau_{\pi}^{-1}+n)},$
(S14)

and this first part gives

$\frac{\tau_{\pi}}{2\pi}\Big{(}\frac{1}{1-n\tau_{\pi}}-\frac{1}{1+n\tau_{\pi}}\Big{)}\cos{(n\pi\tau_{\pi}/2)}.$
(S15)

The second part gives $\frac{1}{n\pi}\cos{(n\pi\tau_{\pi}/2)}$. The whole
integral gives

$f_{n}=\frac{4}{n\pi}\sin{(n\pi/2)}\cos{(n\pi\tau_{\pi}/2)}\Big{[}\frac{n\tau_{\pi}}{2}\frac{2n\tau_{\pi}}{1-n^{2}\tau_{\pi}^{2}}+1\Big{]}=\frac{4\sin{(n\pi/2)}\cos{(n\pi\tau_{\pi}/2)}}{n\pi(1-n^{2}\tau_{\pi}^{2})},$
(S16)

or

$f_{n}=\frac{4\sin{(n\pi/2)}\cos{(n\pi
t_{\pi}/\tau)}}{n\pi(1-n^{2}t_{\pi}^{2}/4\tau^{2})}.$ (S17)

Modulated Pulses: The ansatz proposed for the modulation function is

$f(t)=\cos[{\pi(t-t_{i})/t_{\pi}]}+\frac{d}{\pi kb}\sin(\pi k/2)\bigg{\\{}{\rm
erf}\bigg{(}\frac{t-t_{l}}{ct_{\pi}}\bigg{)}-{\rm
erf}\bigg{(}\frac{t-t_{r}}{ct_{\pi}}\bigg{)}\bigg{\\}}\sin{[k\omega_{k}(t-t_{m})]},$
(S18)

where $t_{\pi}$ is the duration of the $\pi$ pulse, $t_{i}$ and
$t_{m}=t_{i}+t_{\pi}/2$ are the initial and central points of the pulse,
$t_{r}=t_{m}+bt_{\pi}$ and $t_{l}=t_{m}-bt_{\pi}$. Also, $0<b<0.5$ and $c\ll
1$. Moreover, we extend the pulse over $\tau/2$, i.e. $t_{\pi}=\tau/2$. The
modulation function with the introduced change of variables is then

$f_{z}(x)=\sin{(\pi x)}-\frac{d}{\pi kb}\sin(\pi k/2)\sin{(k\pi x)}\,{\rm
erfc}\bigg{(}\frac{x-b}{c}\bigg{)},$ (S19)

where ${\rm erfc}(x)=1-{\rm erf}(x)$ and the Fourier coefficient $f_{n}$ is
given by

$f^{\rm m}_{n}=4\sin{(n\pi/2)}\Bigg{\\{}\int_{0}^{1/2}dx\ \sin{(\pi
x)}\sin{(n\pi x)}-\frac{d}{\pi kb}\sin(\pi k/2)\int_{0}^{1/2}dx\ \sin{(k\pi
x)}\,{\rm erfc}\bigg{(}\frac{x-b}{c}\bigg{)}\sin{(n\pi x)}\Bigg{\\}}.$ (S20)

If $n=1$, the first part gives $1/4$. If $n\neq 1$ this is $0$. The second
part can be rewritten as

$-\frac{a}{2\pi kb}\sin(\pi k/2)\int_{0}^{1/2}dx\ [\cos{[(k-n)\pi
x]}-\cos{[(k+n)\pi x]}]\,{\rm erfc}\bigg{(}\frac{x-b}{c}\bigg{)},$ (S21)

and, for $n=k$, this simplifies to

$-\frac{d}{2\pi kb}\sin(\pi k/2)\int_{0}^{1/2}dx\ [1-\cos{(2k\pi x)}]\,{\rm
erfc}\bigg{(}\frac{x-b}{c}\bigg{)}.$ (S22)

If $1/2-b\gg c$, then ${\rm erfc}[(1/2-b)/c]\approx 0$ and the upper limit of
the integral can taken to be $\infty$. With a last change of variables
$x\rightarrow x+b$, the second part is

$-\frac{d}{2\pi kb}\sin(\pi k/2)\int_{-b}^{\infty}dx\
\\{\cos{[\pi(k-n)(x+b)]}-\cos{[\pi(k+n)(x+b)]}\\}\,{\rm erfc}(x/c).$ (S23)

and, for $n=k$,

$-\frac{d}{2\pi kb}\sin(\pi k/2)\int_{-b}^{\infty}dx\
\\{1-\cos{[2k\pi(x+b)]}\\}\,{\rm erfc}(x/c).$ (S24)

In the appropriate regime, these integrals result in Ng68

$-\frac{d\sin(\pi k/2)}{\pi kb}\Bigg{\\{}\frac{e^{-[(k-n)\pi
c]^{2}/4}}{(k-n)\pi}\sin{[(k-n)\pi b]}-\frac{e^{-[(k+n)\pi
c]^{2}/4}}{(k+n)\pi}\sin{[(k+n)\pi b]}\Bigg{\\}}.$ (S25)

and

$f^{\rm m}_{k}\approx-\frac{4d}{\pi k},$ (S26)

respectively.

##  II Derivation of the stroboscopic second-order Hamiltonian

Here we derive the second-order Hamiltonian associated to a sequence
concatenating $N_{XY4}$ XY4 blocks, where each block is composed by four
pulses with alternating phases, i.e. XYXY. The evolution of each block can be
divided in four parts. During the first part,
$U_{d}(t)=\exp{[-i/2\int_{t_{0}}^{t}\Omega(t)S_{x}]}$, and, thus,

$H(t)=\eta\nu(ae^{-i\nu t}+a^{\dagger}e^{i\nu
t})[f(t,t_{0})S_{z}+f_{\perp}(t,t_{0})S_{y}],$ (S27)

where $f(t,t_{0})=\cos{[\int_{t_{0}}^{t}\Omega(t^{\prime},t_{0})dt^{\prime}]}$
and
$f_{\perp}(t,t_{0})=\sin{[\int_{t_{0}}^{t}\Omega(t^{\prime},t_{0})dt^{\prime}]}$.
During the second part we have that
$U_{d}(t)=\exp{[-i/2\int_{t_{1}}^{t}\Omega(t,t_{1})S_{y}]}(-\sigma_{1}^{x}\sigma_{2}^{x})$
and

$H(t)=\eta\nu(ae^{-i\nu t}+a^{\dagger}e^{i\nu
t})[-f(t,t_{1})S_{z}-f_{\perp}(t,t_{1})S_{x}].$ (S28)

For the third and fourth parts, we have that
$U_{d}(t)=\exp{[-i/2\int_{t_{2}}^{t}\Omega(t,t_{2})S_{x}]}(-\sigma_{1}^{z}\sigma_{2}^{z})$
and
$U_{d}(t)=\exp{[-i/2\int_{t_{3}}^{t}\Omega(t,t_{3})S_{y}]}(-\sigma_{1}^{y}\sigma_{2}^{y})$,
respectively. The corresponding Hamiltonians are

$H(t)=\eta\nu(ae^{-i\nu t}+a^{\dagger}e^{i\nu
t})[-f(t,t_{2})S_{z}-f_{\perp}(t,t_{2})S_{y}]\,\,\,\,{\rm
and}\,\,\,\,H(t)=\eta\nu(ae^{-i\nu t}+a^{\dagger}e^{i\nu
t})[-f(t,t_{3})S_{z}+f_{\perp}(t,t_{3})S_{x}],$ (S29)

respectively. When concatenating $N_{XY4}$ blocks where
$t_{0}=0,t_{1}=\tau/2,t_{2}=\tau,t_{3}=3\tau/2$ and $t_{4}=2\tau$, one can use
the following Hamiltonian to describe the dynamics

$H(t)=\eta\nu(ae^{-i\nu t}+a^{\dagger}e^{i\nu
t})[f_{x}(t)S_{x}+f_{y}(t)S_{y}+f_{z}(t)S_{z}],$ (S30)

where

$\displaystyle f_{x}(t)$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}-f_{\perp}(t,t_{1}+2n\tau)\Big{\\{}\Theta[t-t_{1}-2n\tau]-\Theta[t-t_{2}-2n\tau]\Big{\\}}+f_{\perp}(t,t_{3}+2n\tau)\Big{\\{}\Theta[t-t_{3}-2n\tau]-\Theta[t-t_{4}-2n\tau]\Big{\\}},$
(S31) $\displaystyle f_{y}(t)$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}f_{\perp}(t,2n\tau)\Big{\\{}\Theta[t-2n\tau]-\Theta[t-t_{1}-2n\tau]\Big{\\}}-f_{\perp}(t,t_{2}+2n\tau)\Big{\\{}\Theta[t-t_{2}-2n\tau]-\Theta[t-t_{3}-2n\tau]\Big{\\}},$

and

$\displaystyle f_{z}(t)$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{2N}f(t,n\tau)\Big{\\{}\Theta[t-n\tau]-\Theta[t-t_{1}-n\tau]\Big{\\}}-f(t,t_{1}+n\tau)\Big{\\{}\Theta[t-t_{1}-n\tau]-\Theta[t-t_{2}-n\tau]\Big{\\}},$

where $\Theta(t-t_{i})$ is the Heaviside step function centred in $t_{i}$. The
form of the modulation functions defined in Eqs. (S31)-(II) is shown in Fig.
S1 for the case of top-hat $\pi$ pulses with length $t_{\pi}=\tau/4$.

Fig. S1:  Form of the modulation functions $f_{x,y,z}(t)$ for a single XY4
block with top hat pulses and with $t_{\pi}=\tau/4=\pi/\Omega_{\rm th}$.

For sequences built by concatenating XY4 blocks, the periods of $f_{z}(t)$ and
$f_{x,y}(t)$ differ by a factor of two. After applying the Fourier
decomposition to Eqs. (S31)-(II), these become

$f_{x}(t)=\sum_{n=1}^{\infty}a^{x}_{n}\cos{(n\omega
t/2)}+b^{x}_{n}\sin{(n\omega
t/2)},\,\,\,\,\,\,\,\,\,\,\,\,f_{y}(t)=\sum_{n=1}^{\infty}a^{y}_{n}\cos{(n\omega
t/2)}+b^{y}_{n}\sin{(n\omega t/2)},$ (S32)

and

$f_{z}(t)=\sum_{n=1}^{\infty}f_{n}\cos{(n\omega t)},$ (S33)

where now $n$ stands for the number of the Fourier harmonic,
$\omega=2\pi/\tau$, and the Fourier coefficients are given by

$a^{x,y}_{n}=\frac{1}{\tau}\int_{0}^{2\tau}f_{x,y}(t)\cos{(n\omega
t/2)},\,\,\,\,\,\,\,\,\,\,\,\,b^{x,y}_{n}=\frac{1}{\tau}\int_{0}^{2\tau}f_{x,y}(t)\sin{(n\omega
t/2)},$ (S34)

and

$f_{n}=\frac{2}{\tau}\int_{0}^{\tau}f_{z}(t)\cos{(n\omega t)}.$ (S35)

Here, the sine components of $f_{z}(t)$ are zero because the integral
$\int_{0}^{\tau}f_{z}(t)\sin{(n\omega t)}$ is zero for all $n\in\mathbb{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+t_{2})=-f_{y}(t)$, $f_{x}(t+t_{2})=-f_{x}(t)$,
$f_{y}(t_{1}/2-t)=f_{y}(t_{1}/2+t)$, and $f_{x}(3t_{1}/2-t)=f_{x}(3t_{1}/2+t)$
along with trigonometric identities $\sin{(a\pm
b)}=\sin{a}\cos{b}\pm\cos{a}\sin{b}$ and $\cos{(a\pm
b)}=\cos{a}\cos{b}\mp\sin{a}\sin{b}$, can be used to prove that

$\displaystyle a^{x}_{n}$ $\displaystyle=$ $\displaystyle-\cos{(3\pi
n/4)}e_{n},\ \ b^{x}_{n}=-\sin{(3\pi n/4)}e_{n},$ (S36) $\displaystyle
a^{y}_{n}$ $\displaystyle=$ $\displaystyle\cos{(\pi n/4)}e_{n},\ \ \ \
b^{y}_{n}=\sin{(\pi n/4)}e_{n},$ (S37)

where
$e_{n}=[\cos{(n\pi)}-1]\frac{2}{\tau}\int_{0}^{t_{\pi}/2}f_{\perp}(t+t_{1}/2,0)\cos{(n\omega
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

$H=\frac{1}{2}\eta\nu a^{\dagger}e^{i\nu
t}\sum_{n=1}^{\infty}\Big{\\{}f_{n}(e^{in\omega t}+e^{-in\omega
t})S_{z}+(\bar{c}_{n}^{x}S_{x}+\bar{c}_{n}^{y}S_{y})e^{in\omega
t/2}+({c}_{n}^{x}S_{x}+{c}_{n}^{y}S_{y})e^{-in\omega t/2}\Big{\\}}+{\rm
H.c.},$ (S38)

where $c_{n}^{x,y}=a_{n}^{x,y}+ib_{n}^{x,y}$, and $\bar{c}_{n}^{x,y}$ is the
complex conjugate. As $e_{n}$ is zero for even $n$, $n\omega$ and $n\omega/2$
never coincide, and, if $\eta\nu\ll\omega$, the second order Hamiltonian is

$H_{\rm
XYXY}^{(2)}=\frac{1}{2}\eta^{2}\nu\Big{\\{}-J_{k}^{\parallel}S_{z}^{2}-J_{k}^{\perp}(S_{x}^{2}+S_{y}^{2})+[B_{k}^{\prime}(a^{\dagger}a+1)+B_{k}^{\prime\prime}]S_{z}+J_{k}^{xy}(\sigma_{1}^{x}\sigma_{2}^{y}+\sigma_{1}^{y}\sigma_{2}^{x})\Big{\\}},$
(S39)

where

$\displaystyle J_{k}^{\parallel}$ $\displaystyle=$ $\displaystyle
f_{k}^{2}/2+\sum_{n\neq
k}^{\infty}\frac{f_{n}^{2}}{1-n^{2}\omega^{2}/\nu^{2}},$ (S40) $\displaystyle
J_{k}^{\perp}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{n=1}^{\infty}\frac{e_{n}^{2}}{1-n^{2}\omega^{2}/4\nu^{2}},$
(S41) $\displaystyle B_{k}^{\prime}$ $\displaystyle=$ $\displaystyle
2\sum_{n=1}^{\infty}\frac{n\omega}{\nu}\frac{e_{n}^{2}}{1-n^{2}\omega^{2}/4\nu^{2}}\sin{(n\pi/2)},$
(S42) $\displaystyle B_{k}^{\prime\prime}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\frac{e_{n}^{2}}{1-n^{2}\omega^{2}/4\nu^{2}}\sin{(n\pi/2)},$
(S43) $\displaystyle J_{k}^{xy}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\frac{e_{n}^{2}}{1-n^{2}\omega^{2}/4\nu^{2}}\cos{(n\pi/2)}.$
(S44)

For this, we used that

$H=\sum_{j}A_{j}^{\dagger}e^{i\Delta_{j}t}+A_{j}e^{-i\Delta_{j}t}\approx\sum_{j}\frac{[A_{j}^{\dagger},A_{j}]}{\Delta_{j}}$
(S45)

for $||A_{j}||\ll|\Delta_{j}|$ and
$||A_{j}||\ll|\Delta_{j}-\Delta_{j^{\prime}}|$ ($j\neq j^{\prime}$).

As we indicated above, $e_{n}$ is zero for even $n$. Thus, $J_{k}^{xy}=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
$\frac{1}{2}\eta^{2}\nu[B_{k}^{\prime}(a^{\dagger}a+1)+B_{k}^{\prime\prime}]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 $-\frac{1}{2}\eta^{2}\nu
J_{k}^{\perp}(S_{x}^{2}+S_{y}^{2})$ 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}^{(-)}=\sigma_{1}^{x,y}-\sigma_{2}^{x,y}$, changing the sign of terms
$\sigma_{1}^{u}\sigma_{2}^{\mu}$ where $u,\mu\in\\{x,y\\}$.

##  III Robustness to mode decoherence

Here, we explain how our gates can follow non-circular phase-space
trajectories and obtain enhanced robustness against the mode decoherence.
First, let us recall the propagator associated to an oscillator-mediated two-
qubit phase gate,

$U(s)=\exp{[(\alpha(s)a^{\dagger}-\alpha^{*}(s)a)S_{z}]}\times\exp{[i\theta(s)S_{z}^{2}]}.$
(S46)

where $\alpha(s)$ describes the phase-space trajectory, $s$ is proportional to
time and

$\theta(s)={\rm
Im}\int_{0}^{s}ds^{\prime}\alpha^{*}(s^{\prime})\dot{\alpha}(s^{\prime})$
(S47)

is the two-qubit phase. In the usual case, $\alpha(s)$ follows a circular
trajectory $\alpha(s)\approx\alpha_{\rm circ}(s)=R_{\rm circ}(1-e^{is})$,
where $R_{\rm circ}$ is determined by the target two-qubit phase. At $s_{\rm
g}=2\pi$, the mode-displacement is zero, $\alpha_{\rm circ}(s_{\rm g})=0$, and
the phase is $\theta(s_{\rm g})=2\pi R_{\rm circ}^{2}$. Typically,
$\theta(s_{\rm g})=\pi/8$ is targeted, leading to $R_{\rm circ}=1/4$. That is,
an operation where $\alpha(s)$ follows the phase-space path of a $1/4$ radius
circumference equals a $\exp{[i\pi S_{z}^{2}/8]}$ phase-gate at $s_{\rm
g}=2\pi$.

To become more robust against mode decoherence, Haddadfarshi and Mintert
Haddadfarshi16_2  propose to follow trajectories related to cardioids instead.
The simplest cardioid is parametrized by $\alpha_{\rm card}(s)=R_{\rm
card}(e^{is}-e^{is/2})$, where $s=\\{0,4\pi\\}$. At $s_{\rm g}=4\pi$,
$\alpha_{\rm card}(s_{\rm g})=0$, and the phase is $\theta(s_{\rm g})=6\pi
R_{\rm card}^{2}$, which, with the usual condition $\theta(s_{\rm g})=\pi/8$
leads to a cardioid with $R_{\rm card}={(4\sqrt{3})}^{-1}$.

In our case, the two-qubit phase has a contribution from the circular
trajectory, but also from the dispersive terms,

$\theta(t_{\rm g})\approx\theta_{\rm circ}(t_{\rm g})+\theta_{\rm disp}(t_{\rm
g})=2\pi R_{\rm circ}^{2}+\frac{1}{2}\nu\eta^{2}J_{k}t_{\rm g},$ (S48)

see Eq. (6) of the main text. Because of that, the radius of the circular
trajectory is smaller than $1/4$, and it is given by

$R_{\rm circ}=\frac{1}{4}\Big{\\{}\sqrt{1+4\eta^{2}J_{k}^{2}/f_{k}^{2}}+2\eta
J_{k}/f_{k}\Big{\\}}^{-1}.$ (S49)

For the circular case, $R_{\rm circ}$ and $\xi_{k}$ can be determined exactly.
However, we will describe a method to approximate its value since it will be
useful for the non-circular path. For that, we approximate the final time of
the gate as $t_{\rm g}=8N\tau_{k}\approx 16\pi kN/\nu$, where $\nu$ is the
mode frequency, $k$ is the selected harmonic and $N$ is the number of TQXY16
blocks in the gate. This approximation leads to the following equation

$2\pi R_{\rm circ}^{2}+8\pi k\eta^{2}J_{k}N\approx\pi/8,$ (S50)

or

$R_{\rm circ}\approx\frac{1}{4}\sqrt{1-64k\eta^{2}J_{k}N}.$ (S51)

In the case of the cardioid, we have that

$\theta(t_{\rm g})\approx\theta_{\rm card}(t_{\rm g})+\theta_{\rm disp}(t_{\rm
g})=6\pi R_{\rm card}^{2}+\frac{1}{2}\nu\eta^{2}J_{k}t_{\rm g},$ (S52)

and, using the same approximation, we get

$R_{\rm card}\approx\frac{1}{4\sqrt{3}}\sqrt{1-64k\eta^{2}J_{k}N}.$ (S53)

Notice that, once $\eta$, $k$ and $N$ are defined, Eq. (S53) gives us a good
guess about the value of the parameter $R_{\rm card}$ characterising the
cardioid resulting in $\theta(t_{\rm g})=\pi/8$.

For our $\alpha(t)$ to follow a trajectory different from that of a circle, we
allow the variation of the pulse length $t_{\pi}=\tau_{k}/2$ (thus $\xi_{k}$)
after every TQXY16 block. After the $j$-th TQXY16 block, at $t_{j}$, the value
of $\alpha_{j}\equiv\alpha(t_{j})$ is

$\alpha_{j}=-i\eta\nu\int_{0}^{t_{j}}dt^{\prime}f_{z}(t^{\prime})e^{i\nu
t^{\prime}}=-i\eta\nu\int_{t_{j-1}}^{t_{j}}dt^{\prime}f_{z}(t^{\prime})e^{i\nu
t^{\prime}}+\alpha_{j-1}\approx-i\frac{\eta\nu
f_{k}}{2}\int_{t_{j-1}}^{t_{j}}dt^{\prime}e^{i\xi_{k}^{j}t^{\prime}}+\alpha_{j-1}=\frac{\eta\nu
f_{k}}{2\xi_{k}^{j}}(1-e^{i\xi_{k}^{j}(t_{j}-t_{j-1})})e^{i\xi_{k}^{j}t_{j-1}}+\alpha_{j-1}.$
(S54)

where $t_{0}=0$ and $\alpha_{0}=0$. Using that the length of each block is
$t_{j}-t_{j-1}=8\tau_{k}^{j}=16\pi k/(\nu-\xi_{k})$, this can be rewritten as

$\alpha_{j}=R_{j}\,\big{(}1-e^{i\varphi_{j}}\big{)}e^{i\phi_{j}}+\alpha_{j-1},$
(S55)

where $R_{j}=\eta\nu f_{k}/2\xi_{j}$, $\varphi_{j}=16\pi
k\xi_{j}/(\nu-\xi_{j})$, and
$\phi_{j}=\xi_{j}\times\sum_{j^{\prime}=1}^{j-1}\varphi_{j^{\prime}}/\xi_{j^{\prime}}$.
Also, $\phi_{1}=0$ and, for clarity, $\xi_{j}\equiv\xi_{k,j}$. Notice that
varying $\xi_{j}$ changes the radius and the direction of the $j$-th
trajectory and that concatenating trajectories with different radiuses $R_{j}$
will lead to a non-circular displacement in the phase space.

Fig. S2: Determination of $\xi_{j}$ for every TQXY16 block for G3. (a) From
all $\xi_{1}=\\{0,{\nu}/{8k}\\}$, we find the value of $\xi_{1}$ for which
$\alpha_{1}$ crosses $\alpha_{\rm card}(s)$ for $s=\\{0,2\pi\\}$. In (b) and
(c), the same is plotted but for $\xi_{2}$ and $\xi_{6}$, respectively. (d)
shows the determination of the last six values for $\xi_{j}$, by matching the
respective $\alpha_{j}$ with $\alpha_{\rm card}(s)$, for $s=\\{2\pi,4\pi\\}$.

To find the values of $\xi_{j}$ leading to a path that approximates that
defined by the cardioid characterised by $\alpha_{\rm card}(s)=R_{\rm
card}(e^{is}-e^{is/2})$, where $R_{\rm card}$ is given in Eq. (S53), we follow
this procedure: (i) We choose $\eta$ and $k$ and, following the recipe in the
main text, select a gate composed by $N$ TQXY16 blocks. Note that, when
selecting $N$, the pulse parameters $b,c,$ and $d$ have also been determined.
To building the cardioid trajectory, all blocks will stick to the same values
of $b,c$ and $d$, while we will allow the variation of $t_{\pi}$ for every
block. (ii) We propose a target cardioid with $R_{\rm
card}=\frac{1}{4\sqrt{3}}\sqrt{1-64k\eta^{2}J_{k}\tilde{N}}$ where $\tilde{N}$
is an integer number approximately $122\%$ larger than than $N$. The reason
behind this is that the perimeter of the cardioid is $122\%$ larger than that
of the circle. (iii) Now, ee find all $\xi_{j}$ numerically. Once $\tilde{N}$
is defined, we obtain $\tilde{N}$ equations like Eq. (S55), the first one with
two unknowns, i.e. $\alpha_{1}$ and $\xi_{1}$, and the next ones depending
also on $\alpha_{j-1}$ and previous values of $\xi_{j}$. We find the values of
$\alpha_{1}$ and $\xi_{1}$ by numerically minimising
$|\alpha_{j}(\xi^{\prime}_{j})-\alpha_{\rm card}(s)|$ for variables
$\xi^{\prime}_{j}=\\{0,{\nu}/{8k}\\}$ and $s=\\{0,2\pi\\}$. In Fig. S2(a), we
illustrate this by showing $\alpha_{\rm card}(s)$ for $s=\\{0,2\pi\\}$ in
blue, $\alpha_{1}(\xi^{\prime}_{1})$ for $\xi^{\prime}_{1}=\\{0,{\nu}/{8k}\\}$
in black. The value of $\alpha_{1}$ that crosses the cardioid $\alpha_{\rm
card}(s)$ for $s=\\{0,2\pi\\}$ is shown with a red panel. After the value of
$\alpha_{1}$ (thus $\xi_{1}$) has been numerically determined, we follow the
same procedure for $\alpha_{2}$ and $\xi_{2}$, see Fig. S2(b). We repeat this
procedure until we find all $\alpha_{j}$ and $\xi_{j}$ for the first and
second half of the cardioid with $s=\\{0,2\pi\\}$ and $s=\\{2\pi,4\pi\\}$
respectively, see Figs. S2(c) and (d). The last step (iv) consist in checking
the values of $\alpha_{N}$ and $\theta(t_{\rm g})$, which should give $0$ and
$\pi/8$, respectively. If this is not case, we slightly change the values of
$d$, $R_{\rm card}$ or $\tilde{N}$, and run the protocol described in (iii)
until $\alpha_{N}\approx 0$ and $\theta(t_{\rm g})\approx\pi/8$. Note that
$\alpha_{N}$ and $\theta(t_{\rm g})$ can be calculated exactly using
$\alpha(t_{\rm g})=-i\eta\nu\int_{0}^{t_{\rm
g}}dt^{\prime}f_{z}(t^{\prime})\exp{(i\nu t^{\prime})}$ and $\theta(t_{\rm
g})=\int_{C}\alpha\,d\alpha$, where $f_{z}(t)$ is fully determined given all
$\xi_{k}$.

##  IV Details of numerical simulations

###  IV.1 Two-mode Hamiltonian

The Hamiltonian describing two trapped ions under a static magnetic field
gradient $\partial B/\partial z=g_{B}$, in an interaction picture w.r.t. free
energy term of the qubits
$\frac{\omega_{1}}{2}\sigma_{1}^{z}+\frac{\omega_{2}}{2}\sigma_{2}^{z}$ is
Arrazola18_2

$H_{0}+H_{\rm 2M}=\nu a^{\dagger}a+\sqrt{3}\nu
b^{\dagger}b+\eta\nu(a+a^{\dagger})S_{z}^{(+)}-3^{-1/4}\eta\nu(b+b^{\dagger})S_{z}^{(-)},$
(S56)

where $a^{\dagger}$($a$) and $b^{\dagger}$ ($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}^{(\pm)}+H_{0}+H_{\rm 2M}$, starting from motional
thermal states with $\bar{n}=1$, and truncating the Hilbert space of the first
($a$) and second ($b$) mode to $16$ and $12$ states, respectively.

###  IV.2 Effective Hamiltonian of the driven system

Here we justify the use of the term $H_{\rm 2M}^{\rm eff}$ 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}+\nu
a^{\dagger}a+\sqrt{3}\nu b^{\dagger}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

$H^{\prime}(t)=\sum_{n=1}^{\infty}\frac{1}{2}f_{n}\eta\nu S_{z}^{(+)}ae^{-i\nu
t}(e^{in\omega_{k}t}+e^{-in\omega_{k}t})-\frac{3^{-1/4}}{2}f_{n}\eta\nu
S_{z}^{(-)}be^{-i\sqrt{3}\nu t}(e^{in\omega_{k}t}+e^{-in\omega_{k}t})+{\rm
H.c.}.$ (S57)

In Hamiltonian (S57), the only resonant term is the one going with
$e^{-i(\nu-k\omega_{k})}$ and its complex conjugate, i.e.
$\frac{1}{2}f_{k}\eta\nu S_{z}ae^{-i(\nu-k\omega_{k})t}+{\rm H.c.}$. Using Eq.
(S45), the influence of the rest of the terms is well approximated by the
second order Hamiltonian

$H_{\rm eff}=-\frac{1}{2}\eta^{2}\nu
J^{a}_{k}(S_{z}^{(+)})^{2}-\frac{1}{6}\eta^{2}\nu J^{b}_{k}(S_{z}^{(-)})^{2},$
(S58)

where $J^{a}_{k}=f_{k}^{2}/4+\sum_{n\neq k}^{\infty}f_{n}^{2}/(1-n^{2}/k^{2})$
and $J^{b}_{k}=\sum_{n=1}^{\infty}f_{n}^{2}/(1-n^{2}/3k^{2})$. Because
$(S_{z}^{(-)})^{2}=4\mathbb{I}-S_{z}^{2}$, Hamiltonian (S58) is, up to a
global phase,

$H_{\rm eff}=-\frac{1}{2}\eta^{2}\nu(J^{a}_{k}-J^{b}_{k}/3)S_{z}^{2},$ (S59)

The effective Hamiltonian describing the driven system is then

$H^{\prime}(t)\approx\frac{1}{2}f_{k}\eta\nu
S_{z}a(e^{-i(\nu-k\omega_{k})t}+{\rm
H.c.})-\frac{1}{2}\eta^{2}\nu(J^{a}_{k}-J^{b}_{k}/3)S_{z}^{2},$ (S60)

which is equivalent to Hamiltonian (4) in the main text but with
$J_{k}\rightarrow 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

$H=H_{d}+\nu
a^{\dagger}a+\eta\nu(a+a^{\dagger})S_{z}+\frac{1}{3}\nu\eta^{2}rS_{z}^{2},$
(S61)

where $r=J_{k}^{b}/\sum_{n=1}^{\infty}f_{n}^{2}$. Following the procedure
described from Eq. (S56) to Eq. (S60), but now with $H$, one realises that the
effective Hamiltonian corresponding Eq. (S61) is equivalent to
$H^{\prime}(t)$.

###  IV.3 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_{j}^{0}$,
and both frequencies differ by
$\Delta\omega=\omega_{2}-\omega_{1}=\gamma_{e}g_{B}d$. Here, the distance
between the ions is $d=(e^{2}/2\pi\varepsilon_{0}M\nu^{2})^{1/3}$ where $e$ is
the electric charge of the electron and $\varepsilon_{0}$ is the vacuum
permittivity. For two 171Yr+ ions each with mass $M=171$ amu and trapped with
longitudinal frequency $\nu=(2\pi)\times 220$ kHz, this value gives
$\Delta\omega=2.54$ and $20.34$ MHz for $g_{B}=19.57$ and $153.2\,\rm T/m$,
respectively.

Assuming we address only two levels of each ion, the Hamiltonian of the driven
qubits is

$H_{dq}=\frac{\omega_{1}}{2}\sigma_{1}^{z}+\frac{\omega_{2}}{2}\sigma_{2}^{z}+\Omega(t)(\sigma_{1}^{x}+\sigma_{2}^{x})\cos{(\omega_{1}t-\phi_{1})}+\Omega(t)(\sigma_{1}^{x}+\sigma_{2}^{x})\cos{(\omega_{2}t-\phi_{2})},$
(S62)

where $\Omega_{1}=\Omega_{2}$ for simplicity, and $\phi_{j}$ is the phase of
driving $j$. In an interaction picture w.r.t.
$\omega_{1}/2\sigma_{1}^{z}+\omega_{2}/2\sigma_{2}^{z}$, and after neglecting
fast oscillating terms, the Hamiltonian reads

${H}_{dq}^{I}\approx\frac{\Omega(t)}{2}(\sigma_{1}^{+}e^{i\phi_{1}}+\sigma_{2}^{+}e^{i\phi_{1}}e^{-i\Delta\omega
t}+{\rm
H.c.})+\frac{\Omega(t)}{2}(\sigma_{2}^{+}e^{i\phi_{2}}+\sigma_{1}^{+}e^{i\phi_{2}}e^{i\Delta\omega
t}+{\rm H.c.})$ (S63)

Notice that, if $\Omega(t)\ll\Delta\omega$ , the terms oscillating with
$\Delta\omega$ can also be neglected. This means that, in this regime, the
ions can be then selectively addressed with global fields. For the case
$\phi_{1}=\phi_{2}=\phi$, Hamiltonian (S63) becomes

${H}_{dq}^{I}\approx\frac{\Omega_{x}(t)}{2}S^{(+)}_{x}+\frac{\Omega_{x}(t)}{2}(\sigma_{2}^{-}e^{-i\Delta\omega
t}+\sigma_{1}^{-}e^{i\Delta\omega t}+{\rm H.c.})$ (S64)

or

${H}_{dq}^{I}\approx\frac{\Omega_{y}(t)}{2}S_{y}^{(+)}+\frac{\Omega_{y}(t)}{2}(i\sigma_{2}^{-}e^{-i\Delta\omega
t}+i\sigma_{1}^{-}e^{i\Delta\omega t}+{\rm H.c.}),$ (S65)

depending $\phi=0$ or $\pi/2$. In the main text, we describe both cases using
the more general Hamiltonian

$H_{d}^{(+)}+H_{c}^{(+)}=\frac{\Omega_{x}(t)}{2}S^{(+)}_{x}+\frac{\Omega_{y}(t)}{2}S_{y}^{(+)}+\frac{\Omega_{x}(t)}{2}(\sigma_{2}^{-}e^{-i\Delta\omega
t}+\sigma_{1}^{-}e^{i\Delta\omega t}+{\rm
H.c.})+\frac{\Omega_{y}(t)}{2}(i\sigma_{2}^{-}e^{-i\Delta\omega
t}+i\sigma_{1}^{-}e^{i\Delta\omega t}+{\rm H.c.}),$ (S66)

where $\Omega_{x}(t)=0$ ($\Omega_{y}(t)=0$) when $\phi=\pi/2$ ($\phi=0$). For
the case where the qubits are driven with opposite phase, i.e.
$\phi_{2}=\phi_{1}+\pi=\phi+\pi$, Hamiltonian (S66) is

$H_{d}^{(-)}+H_{c}^{(-)}=\frac{\Omega_{x}(t)}{2}S^{(-)}_{x}+\frac{\Omega_{y}(t)}{2}S_{y}^{(-)}+\frac{\Omega_{x}(t)}{2}(\sigma_{2}^{-}e^{-i\Delta\omega
t}-\sigma_{1}^{-}e^{i\Delta\omega t}+{\rm
H.c.})+\frac{\Omega_{y}(t)}{2}(i\sigma_{2}^{-}e^{-i\Delta\omega
t}-i\sigma_{1}^{-}e^{i\Delta\omega t}+{\rm H.c.})$ (S67)

instead. The Hamiltonian used for the simulations with crosstalk are then
$H_{s}+H^{(+)}_{d}+H^{(+)}_{c}$ and $H_{s}+H^{(-)}_{d}+H^{(-)}_{c}$, for every
first and second half of a TQXY16 block, respectively. Doing so, we assume
that the phase of the drivings can be changed instantaneously.

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 sin2-shaped ramp at the beginning and end of each
pulse. Each pulse is then constructed following

$\Omega(t)=\begin{cases}-\frac{\partial f_{z}(0)}{\partial
t}\times[1-f^{2}_{z}(0)]^{-1/2}\times\sin^{2}({\pi t/2t_{\rm ramp}})\,\,\,{\rm
if}\,\,\,0<t\leq t_{\rm ramp}\\\ -\frac{\partial f_{z}(t)}{\partial
t}\times[1-f^{2}_{z}(t)]^{-1/2}\,\,\,{\rm if}\,\,\,t_{\rm
ramp}<t\leq(t_{\pi}-t_{\rm ramp})\\\ -\frac{\partial f_{z}(t_{\pi})}{\partial
t}\times[1-f^{2}_{z}(t_{\pi})]^{-1/2}\times\sin^{2}({\pi(t-t_{\pi})/2t_{\rm
ramp}})\,\,\,{\rm if}\,\,\,(t_{\pi}-t_{\rm ramp})<t\leq t_{\pi}\end{cases},$
(S68)

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
$\delta\Omega_{\pi}$ multiplying the Rabi frequency, i.e.
$\Omega(t)\rightarrow(1+\delta\Omega_{\pi})\Omega(t)$, to ensure that
$\int_{0}^{t_{\pi}}dt^{\prime}\Omega(t^{\prime})=\pi$.

Fig. S3:  Optimisation of ramp duration (a) Bell state infidelity with respect
to different ramp durations $t_{\rm ramp}$. For gates G1-4, we choose ramps
with length $149,295,1260$ and $49$ ns, respectively. (b) The form of the Rabi
frequencies $\Omega_{x}(t)$ (light line) and $\Omega_{y}(t)$ (dark line) is
shown for gates from G1 (top figure) to G4 (bottom figure).

We simulate the performance of the gate for different durations of the
sin2-shaped ramp $t_{\rm ramp}$. The obtained infidelities are shown in Fig.
S3(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 $\delta\Omega_{\pi}$, the maximum values of
the Rabi frequencies $\Omega_{\rm pp}/2\pi$ change to $125.2,79.05,82.09$, and
$81.14$ kHz, respectively. The Rabi frequencies $\Omega_{x}(t)$ and
$\Omega_{y}(t)$ are shown in Fig. S3(b), for gates G1-4, and for the first
eight pulses.

###  IV.4 Heating of the bosonic mode

To include the effect of motional heating we solve the master equation

$\dot{\rho}=-i[H_{s},\rho]+\mathcal{L}(\rho),$ (S69)

where $\rho$ is the density matrix and the Lindblad superoperator is

$\displaystyle\mathcal{L}(\rho)=\frac{\Gamma}{2}(\bar{N}+1)(2a\rho
a^{\dagger}-a^{\dagger}a\rho-\rho
a^{\dagger}a)+\frac{\Gamma}{2}\bar{N}(2a^{\dagger}\rho a-aa^{\dagger}\rho-\rho
aa^{\dagger})$

with $\bar{N}=[\exp{(\frac{\hbar\nu}{k_{\rm B}T})}-1]^{-1}$ and $T=300$ K. For
$\bar{N}\gg 1$, $\Gamma$ relates to the heating rate as
$\dot{\bar{n}}=\Gamma\bar{N}$. For our simulations, we use a heating rate of
$\dot{\bar{n}}=35$ ph/s for regime (i). We derived this from the rate
$\dot{\bar{n}}_{\rm ref}=120$ ph/s given in Ref. Barthel22_2  for $\nu_{\rm
ref}=(2\pi)\times 120$ kHz, and using the relation
$\dot{\bar{n}}=\dot{\bar{n}}_{\rm ref}\times(\nu_{\rm ref}/\nu)^{2}$
Brownnutt15 . For regime (ii), we use a heating rate of $\dot{\bar{n}}=100$
ph/s Private .

## References

  * (1) E. W. Ng and M. Geller, A Table of Integrals of the Error Functions, J. Res. Natl. Inst. Stand. Technol. 73B 1 (1968).
  * (2) F. Haddadfarshi and F. Mintert, High fidelity quantum gates of trapped ions in the presence of motional heating, New J. Phys. 18, 123007 (2016).
  * (3) I. Arrazola, J. Casanova, J. S. Pedernales, Z.-Y. Wang, E. Solano, and M. B. Plenio, Pulsed dynamical decoupling for fast and robust two-qubit gates on trapped ions, Phys. Rev. A 97, 052312 (2018).
  * (4) P. Barthel, P. H. Huber, J. Casanova, I. Arrazola, D. Niroomand, T. Sriarunothai, M. B. Plenio, and C. Wunderlich, Robust Two-Qubit Gates Using Pulsed Dynamical Decoupling, arXiv:2208.00187 (2022).
  * (5) M. Brownnutt, M. Kumph, P. Rabl, and R. Blatt, Ion-trap measurements of electric-field noise near surfaces, Rev. Mod. Phys. 87, 1419 (2015).
  * (6) Private communication with the W. K. Hensinger group.