Bespoke Pulse Design for Robust Rapid Two-Qubit Gates with Trapped Ions

Two-qubit gate performance is vital for scaling up ion-trap quantum computing. Optimized quantum control is needed to achieve reductions in gate-time and gate error-rate. We describe two-qubit gates with addressed Raman beams within a linear trapped-ion chain by a quantum master equation (QME). The QME incorporates the single-ion two-photon effective Rabi frequency, Autler-Townes and vibrational Bloch-Siegert energy shifts, off-resonant transitions, Raman and Rayleigh scattering, laser-power fluctuations, motional heating, cross-Kerr phonon coupling, laser spillover, asymmetric addressing beams and an imperfect initial motional ground state, with no fitting parameters. Whereas state-of-the-art methods are oblivious to these effects in the gate design procedure. We employ global optimization to design pulse sequences for achieving a robust rapid two-qubit gate for seven trapped $^{171}$Yb$^{+}$ ions by optimizing over numerically integrated QME solutions. Here, robust means resilient against slow drift of motional frequencies, and rapid means gate execution where the effective Rabi frequency is comparable to the detuning of the laser from the ion's bare electronic transition. Our robust quantum control delivers rapid high-quality two-qubit gates in long ion chains, enabling scalable quantum computing with trapped ions.


INTRODUCTION
Quantum computers are implemented on various platforms-trapped ions [1,2], superconducting circuits [3,4], photonic systems [5] and neutral atoms [6]. Although logic cycles are faster for superconducting systems compared to ion traps, ion-trap systems currently feature better operation fidelity and qubits with longer coherence time and higher connectivity, making them one of the leading platforms for achieving noisy intermediatescale quantum (NISQ) and post-NISQ scalable computing [7][8][9]. Key performance indicators for entangling gates in trapped ions include peak power, logic-cycle time τ , connectivity, and Bell-state preparation infidelity I [2,7], estimated from measurements of the output populations and parity-oscillation amplitudes [10,11].
High-quality rapid two-qubit gate (2QG) control policies for preparing Bell states are vital for successful universal quantum computation, yet these gates have been achieved only in few-ion systems [12][13][14][15]. 2QGs for long ion chains operate on adiabatic time scales, where slow quantum information processing suffers from significant amounts of decoherence and fluctuations [16][17][18][19][20].
Current control approaches for 2QGs employ certain simplifications, such as the presumption of unitary evolution, to obtain a closed mathematical expression for trapped-ion dynamics and employ greedy optimization to devise control policies [21][22][23][24][25][26][27][28][29][30][31][32][33][34]. Therefore, such control approaches neither scale properly to fast 2QGs for chains of more than two trapped ions, nor deliver feasible control sequences. Thus, we develop and validate a comprehensive model of controlled open-system dynam-ics described by a quantum master equation (QME). We then numerically integrate the QME and use global optimization to devise quantum control policies for robust rapid 2QGs. Our approach incorporates the detrimental effects of noise and decoherence directly into the gate design procedure. We provide feasible control sequences that deliver 2QG infidelity below 0.01 for a chain of seven trapped 171 Yb + ions in this work, limited by the specific noise environment in the setup considered. We apply Fourier analysis to motional mode dynamics and to control sequences in order to identify and interpret the underlying physics. Our QME method surpasses the standard approaches and delivers rapid, robust 2QG control policies that maintain high fidelity for long ion chains under realistic conditions, thereby removing the major obstacle to implementing scalable quantum computation in trapped ions [2].
Our aim is to reduce the gate duration τ significantly while respecting an in-principle lower bound due to the quantum speed limit [35]. Reducing τ has been explored experimentally for the 'light-shift' gate employing two trapped 43 Ca + ions, achieving τ = 1.6 µs, I = 0.002, and ∼200 mW peak power [36]. Scaling this experimental technique to several ions is very difficult due to the utilization of standing waves with precise ion indexing, and the restriction to qubit states that exhibit linear Zeeman splitting [26].
Our 2QG design method involves four steps. First, we develop a comprehensive model of trapped-ion dynamics represented by a QME. Second, we validate our algorithm for computing the 4×4 reduced density matrix for the two target ions based on empirical data. Third, we cast the arXiv:2212.00702v2 [quant-ph] 25 May 2023 control problem of 2QG design as a feasibility problem. Finally, we employ an off-the-shelf global optimization algorithm to search for feasible robust, rapid 2QGs.

Model.
A linear Paul trap with N alkali-like Λ-type [38] ions (each labelled ȷ ∈ [N ] := {1, . . . , N }) are prepared in an electronic (meta)stable stationary state |0⟩; another (meta)stable stationary state is |1⟩. For a third (meta)stable stationary state |2⟩, three driving fields are applied to each ion: a beam is red-detuned by ∆ from the |0⟩ ↔ |2⟩ transition, and the other two beams are detuned by ∆ ± µ(t) from the |1⟩ ↔ |2⟩ transition, which we call chromatic components for 'red' and 'blue' detuning. Collectively, these three beams drive a bichromatic stimulated Raman transitions for |0⟩ ↔ |1⟩, with effective Rabi frequency Ω ȷ (t) for the ȷ th ion. The phase difference 2ϕ ȷ between the two chromatic components is effectively constant [39]. The three beams are arranged in a counterpropagating geometry to increase net momentum transfer along one of the principal axes of the trap with a wave-vector difference ∆k and to cancel the common-mode phase fluctuation of the bichromatic fields on the spin-dependent force [39]. The time-dependent chromatic beams interact with the ions for t ∈ [0, τ ].
Ion motion is given by N collective oscillatory modes. The rotating quadrature operator is x l (t) := a l exp{−i ((ν l + δ)t + φ l )} + hc, l ∈ [N ], (1) with a l the phonon-annihilation operator and hc the Hermitian conjugate. Here ν = (ν l ) ∈ R N represents motional angular frequencies, with l = 1 for centre-of-mass (CoM) motion of the Wigner crystal, next, l = 2 is for the rocking (tilt) mode and so on up to l = N . The offset in motional-mode frequencies, caused by slow drift in the overall trapping strength, is represented by δ. Phase offsets are incorporated into φ = (φ l ) ∈ [0, 2π] N . The coupling parameter between ion ȷ and mode l in the interaction is described by the Lamb-Dicke parameters η ȷl := b ȷl ∆k ℏ /2m(ν l +δ) with b ȷl the normal-mode coupling parameters and m the ion's mass [40]. The ȷ th ion's position operator is β ȷ (t) : The internal states of each ion are treated as a twolevel system (2LS) with |2⟩ adiabatically decoupled from the Raman transition with Ω ȷ (t) ≤ Ω max . The laser detuning µ(t) ≡ µ is constant, and we restrict ourselves to identical pulses for each of the two (r and s) target ions for the 2QG and zero driving for the other ions. The ion motion is not perfectly cooled to the ground state: η 2 ȷl (2n l + 1) ≪ 1∀ȷ, l withn l the mean phonon occupation number of mode l. Note that our condition differs from the original definition of a 'warm' ion, which presupposes that all but the motional mode used for logical operations have zero phonons [37]. All motional modes are in thermal equilibrium. The temperature T for all motional modes is specified by mean phonon numbern l for each mode l [40,41].

QUANTUM MASTER EQUATION
The interaction Hamiltonian between the ȷ th ion and the electromagnetic field, in the interaction picture with the rotating-wave approximation with respect to qubit resonance frequency, is (2) with ladder operators σ + = |1⟩ ⟨0| = σ − † . We use the interaction picture to avoid numerical instabilities while integrating the QME. The system Hamiltonian for all ions is The Autler-Townes shift is calculated and included as a pulse-dependent drift of laser detuning. The cross-Kerr coupling is taken into account by adding the term χn l n l ′ [42][43][44] to H(t), where n l is the phonon-number operator of motional mode l and the coupling constant χ is taken from the experimental data reported in Ref. [44]. The laser spillover onto other ions is 2% of the Rabi frequency on neighbouring ions.
Problem 1 (Find a physical solution). Input: number N of trapped ions, labels (r, s) of any pair of ions, motional angular frequencies ν, temperature T for the CoM mode, physical bounds for Ω(t) and µ, δ tol and I ⊚ . The ions' internal states are initialized in |0⟩. Output: Shortest (r, s)-2QG duration τ such that (Ω(t), µ) exists and satisfies and, for |Φ ± ⟩ := |00⟩ ± i |11⟩, with ρ the QME solution. The problem is thus to reduce τ , which is a functional of Ω(t) and µ, while ensuring that I does not increase and the 2QG is robust. 2QG speed-up typically requires increased peak Rabi frequency [55]. Formally, we adapt Problem 1 mathematically as a feasibility problem. First, we discretize time by introducing a time mesh comprising m equal segments so t ∈ (τ /m)Z m+1 with m a hyperparameter, meaning a parameter that is 'tuned' for optimization. As m = 2N + 1 suffices [21,22], we replace Ω(t) by Ω ∈ [0, Ω max ] m . The functional τ [Ω(t), µ] is thus replaced by the function τ (Ω, µ).

APPROACH
We develop a code to simulate an ion-trap system, which involves the time-dependent Hamiltonian, jump operators with time-dependent rates, input state descriptions, thermal preparation of motional modes and Bellstate preparation fidelity estimator. Subsequently, we unravel the QME employing the quantum trajectories theory [56] and use the C++ 'Quantum trajectory class library' [57] to integrate the resulting stochastic quantum trajectories. We also develop our own code for the state-of-the-art (SotA) method for pulse shaping [21,22].
We solve the feasibility problem by employing an off-the-shelf global-optimization (GO) algorithm, specifically differential evolution [58,59], to search for a feasible Ω and µ in the region [0, Ω max ] m × [µ min , µ max ] over the uniform measure. We then compare our method to SotA. Whereas our GO method uses the numerically integrated solutions of the QME to compute I, the SotA method relies on certain simplifications to obtain a closed mathematical expression for trapped-ion dynamics. SotA  (3) TABLE II. We use existing data together with independently measured device noise characteristics to validate our method. A successful validation means reproducing the accessible empirical data within experimental error. The table presents empirical and estimated results for 2QGs designed using the state-of-the-art amplitude and frequency modulation method [60]. Here τ corresponds to 2QG duration, Ω0(error) to Rabi frequency, P (error) to even-parity population andP to estimated even-parity population.
and other standard approaches then obtain a pulse sequence by minimizing the phase-space closure criterion without regard to the geometric phase condition required for the entangling 2QGs. The pulse sequence is later scaled to meet the geometric phase condition, assuming that the dynamics of motional phase-space trajectories are linear in Ω(t). Finally, we gain physical insight into pulse sequences and motional mode dynamics through Fourier analysis and by investigating the infidelity contributions by different noise sources.
We validate our algorithm for computing the 4 × 4 reduced density matrix for the two target ions as follows. The predicted I and even-parity population P must agree with empirical results within experimental error. We choose these two quantities to validate as the experimental results along with all requisite experimental parameters and device noise characteristics are available to us [55,60].
We simulate a chain of seven 171 Yb + -ions in a linear Paul trap with motional frequencies on a transverse axis spanning from 3.07 MHz (CoM) to 2.96 MHz (7 th mode) and η ȷ1 ≡ 0.065 ∀ȷ based on the experimental parameters in Ref. [55]. The qubit states, |F, m F ⟩, are encoded into 171 Yb + hyperfine states |0, 0⟩ and |1, 0⟩ with separation ω 0 = 2π × 12.642821 GHz. Stimulated Raman transfer occurs by virtually exciting 6 2 P1 /2 and 6 2 P3 /2 . One Raman beam is global, it illuminates the entire crystal, whereas the others are a pair of tightly focused beams on the target ions.
In order to validate our algorithm, we first estimate I for the (3,4)-2QG with τ = 254 µs [60]. Our resultant I = 0.005(3) agrees, within experimental error, with the empirical result I = 0.007 (4). Second, Table II summarizes predicted and empirical values of P for the (4,5)-2QG over five values of τ [55]. The data set shows good theoretical-empirical agreement, considering the uncertainties in η ȷl and Ω 0 [55].

RESULTS
Now we show that our GO method delivers superior pulse sequences compared to SotA in the sense that our GO method always works when SotA works, and our GO method succeeds even when SotA fails. In Figs. 1(a) and 1(b) we present two SotA pulse sequences and three GO pulse sequences for a fast 2QG. In this case, all GO pulse sequences are feasible, but neither SotA pulse sequence is feasible unless I ⊚ is adequately increased.
We now discuss GO and SotA performance for an increasingly fine time mesh hyperparamaterized by m. We observe that, for SotA, increasing m beyond 2N + 1 increases both I and peak Rabi frequency Ω peak := max Ω for fixed τ and µ. These increases are evident in Fig. 1(b). Specifically, no feasible SotA pulse sequence for m = 30 exists unless we double the Ω max constraint and increase I ⊚ to 0.50. On the contrary, the m = 30 GO pulse sequence is feasible under the same constraint as is the m = 15 GO pulse sequence.
Computing I with the jump terms for individual noise sources removed, reveals the contribution of different sources of decoherence. In our case, the most detrimental are motional dephasing and laser power fluctuation, which account for 23.4% and 22.7% of the total I, respectively. This kind of analysis can help identify the most profitable avenues for hardware improvements. Now we discuss pulse-sequence robustness against long-term drift δ of motional frequencies. Figure 1(c) presents a plot of I versus δ for the (3,4)-2QG in a seven-ion chain. The GO pulse sequences yield superior I and meet the feasibility condition I ⊚ = 0.01 for the experimentally relevant domain |δ| ≤ δ tol = 1.5 kHz. The SotA pulse sequence does not achieve I ⊚ = 0.01 and is minimized outside of the δ tol ≤ 1.5 kHz domain. Figure 1(c) also indicates the greater robustness that is achieved with respect to δ by increasing, either Ω max or m for our method. The displacement in the minimum of the robustness curve for the SotA pulse sequence is due to not accounting for the vibrational Bloch-Siegert shift in the Hamiltonian [61]. Now we Fourier analyze Ω to gain physical insight. We sample the SotA and GO pulse sequences with Ωmax /2π = 1129 kHz and Ωmax /2π = 840 kHz, respectively, at f s = 10 MHz to minimize aliasing. The discrete Fourier trans-formsΩ of the sampled SotA and GO pulse sequences are shown in Fig. 1(d), where f ∈ (1/τ )Z ( fs τ /2)+1 . The spectral width of both pulse sequences is about 28.57 kHz. The feature around 428 kHz is an artifact due to discreteness imposed by the segment size of 2.3 µs. The broad bump in both spectra starting from about 142 kHz to 285 kHz corresponds to contributions about two times the pulse segment size of 2.3 µs. The SotA pulse sequence exhibits a greater contribution in this region, which corresponds to the square-wave-like pattern of the pulse se-quence with a period of twice the pulse segment size. The GO pulse sequence does not manifest such pronounced periodicity.
Phase profiles for SotA and GO are shown in Fig. 1(e). We note that significant differences in argΩ can indicate a deficiency in the SotA pulse-sequence design method. SotA and GO phases are similar for all frequencies (mod 2π) except between 150 kHz and 250 kHz, where the phase of the GO pulse is antisymmetric about the midpoint 200 kHz. This antisymmetry signifies weak coupling to the CoM mode in comparison to the SotA solution, showing that the GO pulse succeeds because GO avoids channelling energy to the CoM and instead concentrates the force onto the lower-frequency modes. Now we analyze the spectral properties of phase-space trajectories for the sixth and seventh motional modes. These trajectories are defined by the complex-valued a l (t) := tr(|++⟩⟨++| ⊗ a l ρ(t)), where |+⟩ := |0⟩ + |1⟩, and are solved for the SotA pulse sequence with Ωmax /2π = 970 kHz and µ /2π = 2.88 MHz, and plotted in Figs. 2(a) and 2(b). We observe that the phase-space trajectories are jagged when solved using our QME, which is due to off-resonant coupling of the driving field to the bare electronic transition. The SotA model is oblivious to this effect. Therefore, SotA pulse sequences fail to close the phase-space trajectories into orbits at the end of 2QG implementation, which leads to residual entanglement between the electrons and the motional modes.
The Fourier transform of the real component of the seventh motional mode is shown in Fig. 2(c). The QME approach shows frequency contributions around multiples of µ, which correspond to the higher-order effects of the off-resonant bare electronic transitions.

CONCLUSION
Although we provide a study of our 2QG design method on a simulated chain of seven 171 Yb + -ions in a linear Paul trap, our method is readily extensible to longer chains of ions, gate duration and physical constraints. We validate our model by computing I for two of the state-of-the-art 2QG methods within experimental error. Thus, we show that the physics we extract applies across other control scenarios, error environments, and system sizes considered within the scope of this work. The maximum dimension of the Hilbert space we use for simulating the chain of seven 171 Yb + -ions is 320000, which includes the Hilbert space for two ions plus two more neighbouring ions due to laser spill-over. Additionally, our optimized implementation of the QME integrator dynamically lowers or increases the size of Fock space for motional modes based on their occupation probability and adjacency to µ. The parallelized implementation of our GO method takes about 45 min on a computer with 40 cores to obtain a feasible 2QG solution, which is well within the time frame required for an experimental test. We expect that with further optimization this time can be considerably reduced as well.
Our work uncovers principles for 2QG design that are informed by a comprehensive model of trapped ions. To minimize off-resonant carrier modulation, the pulse sequence should begin and end smoothly at zero strength since the amplitude of the modulation depends linearly on the Rabi frequency. Moreover, the pulse symmetry, together with the appropriate detuning, contribute to close the phase-space trajectories and improve the robustness.