ERROR LIMITING PROTOCOL FOR THE CONSTRUCTION OF TWO-QUBIT GATES IN AN ION-TRAP QUANTUM COMPUTING SYSTEM

Abstract

A method of performing a two-qubit gate operation includes computing, by a classical computer, a control pulse to be applied to a pair of trapped ions in a plurality of trapped ions in a quantum processor, each of the plurality of trapped ions having two frequency-separated states defining a qubit, wherein computing the control pulse comprises: computing a pulse function of the control pulse based on a phase-space closure condition and an auxiliary condition, and computing the pulse function of the control pulse further based on a gate angle condition, and applying, by a system controller, the control pulse, having the computed pulse function, to the pair of trapped ions.

Description

BACKGROUND

Field

The present disclosure generally relates to a method of performing computations in a quantum computing system, and more specifically, to a method of constructing two-qubit gate operations in a quantum computing system that includes a chain of trapped ions.

Description of the Related Art

Among physical systems upon which it is proposed to build large-scale quantum computers, is a chain of ions (i.e., charged atoms), which are trapped and suspended in vacuum by electromagnetic fields. The ions have internal hyperfine states which are separated by frequencies in the several GHz range and can be used as the computational states of a qubit (referred to as “qubit states”). These hyperfine states can be controlled using radiation provided from a laser, or sometimes referred to herein as the interaction with laser beams. The ions can be cooled to near their motional ground states using such laser interactions. The ions can also be optically pumped to one of the two hyperfine states with high accuracy (preparation of qubits), manipulated between the two hyperfine states (single-qubit gate operations) by laser beams, and their internal hyperfine states detected by fluorescence upon application of a resonant laser beam (read-out of qubits). A pair of ions can be controllably entangled (two-qubit gate operations) by a qubit-state dependent force using laser pulses that couple the ions to the collective motional modes of a chain of trapped ions, which arise from the Coulombic interaction between the ions. In general, entanglement occurs when pairs or chains of ions (or particles) are generated, interact, or share spatial proximity in ways such that the quantum state of each ion cannot be described independently of the quantum state of the others, even when the ions are separated by a large distance.

Realistic fault-tolerant quantum computing at reasonable overhead requires two-qubit gates with the highest possible fidelity, typically an infidelity of ≤10⁻⁴. However, in a commonly used phase-sensitive architecture, even under noise-free ideal conditions, control pulses constructed by conventional methods do not provide Mølmer-Sørensen XX gates based on the Raman scheme with an infidelity of ≤10⁻⁴.

Therefore, there is a need for a method of constructing gate operations with respect to the parameter drifts or fluctuations within an acceptable error in quantum computation.

SUMMARY

Embodiments of the present disclosure provide a method of performing a two-qubit gate operation. The method includes computing, by a classical computer, a control pulse to be applied to a pair of trapped ions in a plurality of trapped ions in a quantum processor, each of the plurality of trapped ions having two frequency-separated states defining a qubit, wherein computing the control pulse comprises computing a pulse function of the control pulse based on a phase-space closure condition and an auxiliary condition, and computing the pulse function of the control pulse further based on a gate angle condition, and applying, by a system controller, the control pulse, having the computed pulse function, to the pair of trapped ions.

Embodiments of the present disclosure also provide a quantum computing system. The quantum computing system includes a quantum processor comprising a plurality of trapped ions, wherein each of the trapped ions having two frequency-separated states defining a qubit, a classical computer configured to compute a control pulse to be applied to a pair of trapped ions in the quantum processor, wherein computing the control pulse comprises computing a pulse function of the control pulse based on a phase-space closure condition and an auxiliary condition, and computing the pulse function of the control pulse further based on a gate angle condition, and a system controller configured to apply the control pulse, having the computed pulse function, to the pair of trapped ions.

Embodiments of the present disclosure further provide a quantum computing system. The quantum computing system includes non-volatile memory having a number of instructions stored therein which, when executed by one or more processors, causes the quantum computing system to perform operations comprising computing, by a classical computer, a control pulse to be applied to a pair of trapped ions in a plurality of trapped ions in a quantum processor, each of the plurality of trapped ions having two frequency-separated states defining a qubit, wherein computing the control pulse comprises computing a pulse function of the control pulse based on a phase-space closure condition and an auxiliary condition, and computing the pulse function of the control pulse further based on a gate angle condition, and applying, by a system controller, the control pulse, having the computed pulse function, to the pair of trapped ions.

BRIEF DESCRIPTION OF THE DRAWINGS

So that the manner in which the above-recited features of the present disclosure can be understood in detail, a more particular description of the disclosure, briefly summarized above, may be had by reference to embodiments, some of which are illustrated in the appended drawings. It is to be noted, however, that the appended drawings illustrate only typical embodiments of this disclosure and are therefore not to be considered limiting of its scope, for the disclosure may admit to other equally effective embodiments.

FIG. 1 is a schematic partial view of an ion trap quantum computing system according to aspects of this disclosure.

FIG. 2 depicts a schematic view of an ion trap for confining ions in a chain according to aspects of this disclosure.

FIG. 3 depicts a schematic energy diagram of each ion in a chain of trapped ions according to aspects of this disclosure.

FIG. 4 depicts a qubit state of an ion represented as a point on the surface of the Bloch sphere.

FIGS. 5A, 5B, and 5C depict a few schematic collective transverse motional mode structures of a chain of five trapped ions.

FIGS. 6A and 6B depict schematic views of a motional sideband spectrum of each ion and a motional mode according to aspects of this disclosure.

FIG. 7 depicts the relationship between a hardware-implemented quantum computer (REALITY, governed by the Hamiltonian Ĥ_R), the corresponding model quantum computer (MODEL, governed by the model Hamiltonian Ĥ_M) and the control pulses (CONTROL, constructed on the basis of the control Hamiltonian Ĥ_C) that control both the actual hardware-implemented quantum computer (REALITY) and the model quantum computer (MODEL).

FIG. 8 depicts a number of gates vs. Φ-infidelity for N=36 ions, uncalibrated, τ=300 μs AMFM pulses, not requiring the Φ condition. The infidelities are in units of 10⁻¹.

To facilitate understanding, identical reference numerals have been used, where possible, to designate identical elements that are common to the figures. In the figures and the following description, an orthogonal coordinate system including an X-axis, a Y-axis, and a Z-axis is used. The directions represented by the arrows in the drawing are assumed to be positive directions for convenience. It is contemplated that elements disclosed in some embodiments may be beneficially utilized on other implementations without specific recitation.

DETAILED DESCRIPTION

Embodiments described herein are generally related to a method of performing a computation in a quantum computing system, and more specifically, to a method of constructing control pulses that implement entangling gates in a quantum computing system that includes a chain of trapped ions. The method can include a process of constructing the control pulses that implement stabilized entangling gate operations used in the computational process performed by a quantum computing system.

Embodiments of the disclosure include a quantum computing system that is able to perform a quantum computation process by use of a classical computer, a system controller, and a quantum processor. The classical computer performs supporting tasks including selecting a quantum algorithm to be used, computing quantum circuits to run the quantum algorithm, and outputting results of the execution of the quantum circuits by use of a user interface. A software program for performing the tasks is stored in a non-volatile memory within the classical computer. The quantum processor includes trapped ions that are coupled with various hardware, including lasers to manipulate internal hyperfine states (qubit states) of the trapped ions and photomultiplier tubes (PMTs) to read-out the internal hyperfine states (qubit states) of the trapped ions. The system controller receives from the classical computer instructions for controlling the quantum processor, and controls various hardware associated with controlling any and all aspects used to run the instructions for controlling the quantum processor, and transmits a read-out of the quantum processor and thus output of results of the read-out to the classical computer. In some embodiments, the classical computer will then utilize the computational results based on the output of results of the read-out to form a results-set that is then provided to a user in the form of results displayed on a user interface, stored in a memory and/or transferred to another computational device for solving technical problems.

In Section III, various Hamiltonians used in gate construction and error analysis are shown. In Section IV, construction methods, focusing on the AMFM pulse-construction method that are used for generating test pulses are discussed. In Section V, a gate-simulation method is shown, which propagates an initial state |ψ(t=0)> into a final state |ψ(τ)>, where τ is the gate time, including the 2-qubit computational space and a portion of the phonon space large enough to obtain converged results. In Section VI, various fidelity measures are defined that are used to assess the quality of various Hamiltonians used to construct XX gates. In Section VII, numerical simulation results are shown. In Section VIII, a simple pulse-scaling method (calibration) is shown that may be used to eliminate coherent errors that result in XX gates that deviate from the target degree of entanglement. In Section IX, sources of errors that are incurred by a general pulse function g(t) are analyzed. This Method may be used in particular for a pulse function g(t) constructed on the basis of the standard Hamiltonian to control a quantum computer assumed to be governed by the full Mølmer-Sørensen Hamiltonian. The analysis is conducted analytically, and error integrals are evaluated numerically to assess error magnitudes. In Section X, a linear method of eliminating an important class of coherent errors are presented. It is this method, together with pulse calibration that allows us to suppress coherent XX gate errors to the level of ≤10⁻¹. In Section XI, the scaling of the results to the case of chains consisting of 32 ions is discussed.

I. General Hardware Configurations

FIG. 1 is a schematic partial view of an ion trap quantum computing system according to one embodiment. The ion trap quantum computing system 100 includes a classical (digital) computer 102, a system controller 104 and a quantum processor that is a chain 106 of trapped ions (i.e., five shown) that extend along the Z-axis. Each ion in the chain 106 of trapped ions is an ion having a nuclear spin I and an electron spin S such that a difference between the nuclear spin I and the electron spin S is zero, such as a positive ytterbium ion, ¹⁷¹Yb⁺, a positive barium ion ¹³³Ba⁺, a positive cadmium ion ¹¹¹Cd⁺ or ¹¹³Cd⁺, which all have a nuclear spin

$I=\frac{1}{2}$

and the ²S_1/2 hyperfine states. In some embodiments, all ions in the chain 106 of trapped ions are the same species and isotope (e.g., ¹⁷¹Yb⁺). In some other embodiments, the chain 106 of trapped ions includes one or more species or isotopes (e.g., some ions are ¹⁷¹Yb⁺ and some other ions are ¹³³Ba⁺). In yet additional embodiments, the chain 106 of trapped ions may include various isotopes of the same species (e.g., different isotopes of Yb, different isotopes of Ba). The ions in the chain 106 of trapped ions are individually addressed with separate laser beams. The classical computer 102 includes a central processing unit (CPU), memory, and support circuits (or I/O). The memory is connected to the CPU, and may be one or more of a readily available memory, such as a read-only memory (ROM), a random-access memory (RAM), floppy disk, hard disk, or any other form of digital storage, local or remote. Software instructions, algorithms and data can be coded and stored within the memory for instructing the CPU. The support circuits (not shown) are also connected to the CPU for supporting the processor in a conventional manner. The support circuits may include conventional cache, power supplies, clock circuits, input/output circuitry, subsystems, and the like.

An imaging objective 108, such as an objective lens with a numerical aperture (NA), for example, of 0.37, collects fluorescence along the Y-axis from the ions and maps each ion onto a multi-channel photo-multiplier tube (PMT) 110 for measurement of individual ions. Non-copropagating Raman laser beams from a laser 112, which are provided along the X-axis, perform operations on the ions. A diffractive beam splitter 114 creates an array of static Raman beams 116 that are individually switched using a multi-channel acousto-optic modulator (AOM) 118 and is configured to selectively act on individual ions. A global Raman laser beam 120 is configured to illuminate all ions simultaneously. In some embodiments, individual Raman laser beams (not shown) each illuminate individual ions. The system controller (also referred to as a “RF controller”) 104 controls the AOM 118 and thus controls laser pulses to be applied to trapped ions in the chain 106 of trapped ions. The system controller 104 includes a central processing unit (CPU) 122, a read-only memory (ROM) 124, a random-access memory (RAM) 126, a storage unit 128, and the like. The CPU 122 is a processor of the system controller 104. The ROM 124 stores various programs and the RAM 126 is the working memory for various programs and data. The storage unit 128 includes a nonvolatile memory, such as a hard disk drive (HDD) or a flash memory, and stores various programs even if power is turned off. The CPU 122, the ROM 124, the RAM 126, and the storage unit 128 are interconnected via a bus 130. The system controller 104 executes a control program which is stored in the ROM 124 or the storage unit 128 and uses the RAM 126 as a working area. The control program will include software applications that include program code that may be executed by a processor in order to perform various functionalities associated with receiving and analyzing data and controlling any and all aspects of the methods and hardware used to create the ion trap quantum computer system 100 discussed herein.

FIG. 2 depicts a schematic view of an ion trap 200 (also referred to as a Paul trap) for confining ions in the chain 106 according to one embodiment. The confining potential is exerted by both a static (DC) voltage and a radio frequency (RF) voltage. A static (DC) voltage is applied to end-cap electrodes 210 and 212 to confine the ions along the Z-axis (also referred to as an “axial direction” or a “longitudinal direction”). The ions in the chain 106 are nearly evenly spaced in the axial direction due to the Coulomb interaction between the ions. In some embodiments, the ion trap 200 includes four hyperbolically shaped electrodes 202, 204, 206, and 208 extending along the Z-axis.

During operation, a sinusoidal voltage V₁ (with an amplitude V_RF/2) is applied to an opposing pair of electrodes 202, 204 and a sinusoidal voltage V₂ with a phase shift of 180° from the sinusoidal voltage V₁ (and the amplitude V_RF/2) is applied to the other opposing pair of electrodes 206, 208 at a driving frequency ω_RF, generating a quadrupole potential. In some embodiments, a sinusoidal voltage is only applied to one opposing pair of electrodes 202, 204, and the other opposing pair 206, 208 is grounded. The quadrupole potential creates an effective confining force in the X-Y plane perpendicular to the Z-axis (also referred to as a “radial direction” or “transverse direction”) for each of the trapped ions, which is proportional to the distance from a saddle point (i.e., a position in the axial direction (Z-direction)) at which the RF electric field vanishes. The motion in the radial direction (i.e., direction in the X-Y plane) of each ion is approximated as a harmonic oscillation (referred to as secular motion) with a restoring force towards the saddle point in the radial direction and can be modeled by spring constants k_x and k_y, respectively. In some embodiments, the spring constants in the radial direction are modeled as equal when the quadrupole potential is symmetric in the radial direction. However, undesirably in some cases, the motion of the ions in the radial direction may be distorted due to some asymmetry in the physical trap configuration, a small DC patch potential due to inhomogeneity of a surface of the electrodes, or the like and due to these and other external sources of distortion the ions may lie off-center from the saddle points. The Paul trap described herein is just one example of the types of traps that can be used as the ion trap 200. Other types of traps, including surface traps, can also be used for this purpose although their operation may be somewhat different.

FIG. 3 depicts a schematic energy diagram 300 of each ion in the chain 106 of trapped ions according to one embodiment. Each ion in the chain 106 of trapped ions is an ion having a nuclear spin I and an electron spin S such that a difference between the nuclear spin I and the electron spin is zero. In one example, each ion may be a positive Ytterbium ion ¹⁷¹Yb⁺, which has a nuclear spin

$I=\frac{1}{2}$

and the ²S_1/2 hyperfine states (i.e., a multiplet of electronic states of which two are used as computational states) with an energy split corresponding to a frequency difference (referred to as a “carrier frequency”) of ω₀₁/2π=12.642812 GHz. In other examples, each ion may be a positive barium ion ¹³³Ba⁺, a positive cadmium ion ¹¹¹Cd⁺ or ¹¹³Cd⁺, which all have a nuclear spin

$I=\frac{1}{2}$

and the ²S_1/2 hyperfine states. A qubit is formed with two hyperfine states, denoted as |d and |1, where the hyperfine ground state (i.e., the lowest energy state of the ²S_1/2 hyperfine states) is chosen to represent |0. Hereinafter, the terms “hyperfine states”, “internal hyperfine states”, and “qubit states” may be interchangeably used to represent |0 and |1. Further, the terms “trapped ions,” “ions,” and “qubits” may be interchangeably used. Each ion may be cooled (i.e., the kinetic energy of the ion may be reduced) to near the motional ground state |0_p for any motional mode with no phonon excitation (i.e., n_ph=0) by known laser cooling methods, such as Doppler cooling or resolved sideband cooling, and then the qubit state prepared in the hyperfine ground state |0 by optical pumping. Here, |0 represents the individual qubit state of a trapped ion whereas |0_p with the subscript denotes the motional ground state for a motional mode p of a chain 106 of trapped ions.

An individual qubit state of each trapped ion may be manipulated by, for example, a mode-locked laser at 355 nanometers (nm) via the excited ²P_1/2 level (denoted as |e). As shown in FIG. 3, a laser beam from the laser may be split into a pair of non-copropagating laser beams (a first laser beam with frequency ω₁ and a second laser beam with frequency ω₂) in the Raman configuration, and detuned by a one-photon transition detuning frequency Δ=ω₁−ω_0e with respect to the transition frequency ω_0e between |0 and |e, as illustrated in FIG. 3. A two-photon transition detuning frequency δ includes adjusting the amount of energy that is provided to the trapped ion by the first and second laser beams, which, when combined, is used to cause the trapped ion to transfer between the hyperfine states |0 and |1. When the one-photon transition detuning frequency Δ is much larger than the two-photon transition detuning frequency (also referred to simply as “detuning frequency”) δ=ω₁−ω₂−ω₀₁ (hereinafter denoted as ±μ, μ being a positive value), the single-photon Rabi frequencies Ω_0e(t) and Ω_1e(t) (which are time-dependent, and are determined by amplitudes and phases of the first and second laser beams), at which Rabi flopping between states |0 and |e and between states |1 and |e respectively occur, and the spontaneous emission rate from the excited state |e, Rabi flopping between the two hyperfine states |0 and |1 (referred to as a “carrier transition”) is induced at the two-photon Rabi frequency Ω(t). The two-photon Rabi frequency Ω(t) has an intensity (i.e., absolute value of amplitude) that is proportional to Ω_0eΩ_1e/2Δ, where Ω_0e and Ω_1e are the single-photon Rabi frequencies due to the first and second laser beams, respectively. Hereinafter, this set of non-copropagating laser beams in the Raman A configuration to manipulate internal hyperfine states of qubits (qubit states) may be referred to as a “composite pulse” or simply as a “pulse,” described by a pulse function g(t) and the resulting time-dependent pattern of the two-photon Rabi frequency Ω(t) may be referred to as an “amplitude” of a pulse. The detuning frequency δ=ω₁−ω₂−ω₀₁ may be referred to as the detuning frequency of the composite pulse or the detuning frequency of the pulse. The amplitude of the two-photon Rabi frequency Ω(t), which is determined by amplitudes of the first and second laser beams, may be referred to as an “amplitude” of the composite pulse, such that the pulse function g(t) is now represented as

$g\left(t\right)=\Omega \left(t\right)\sin \left\{{\int }_0^t\mu \left(t^{\prime }\right){\mathrm{dt}}^{\prime }+\phi \left(t\right)\right\},$

where φ(t) is a “phase” of the composite pulse that may be time dependent.

It should be noted that the particular atomic species used in the discussion provided herein is just one example of atomic species which have stable and well-defined two-level energy structures when ionized and an excited state that is optically accessible, and thus is not intended to limit the possible configurations, specifications, or the like of an ion trap quantum computer according to the present disclosure. For example, other ion species include alkaline earth metal ions (Be⁺, Ca⁺, Sr⁺, Mg⁺, and Ba⁺) or transition metal ions (Zn⁺, Hg⁺, Cd⁺).

FIG. 4 is provided to help visualize a qubit state of an ion represented as a point on the surface of the Bloch sphere 400 with an azimuthal angle ϕ and a polar angle θ. Application of the composite pulse as described above, causes Rabi flopping between the qubit state |0 (represented as the north pole of the Bloch sphere) and |1 (the south pole of the Bloch sphere) to occur. Adjusting time duration and the composite pulse flips the qubit state from |0 to |1 (i.e., from the north pole to the south pole of the Bloch sphere), or the qubit state from |1 to |0 (i.e., from the south pole to the north pole of the Bloch sphere). This application of the composite pulse is referred to as a “π-pulse”, Further, by adjusting time duration and the composite pulse, the qubit state |0 may be transformed to a superposition state |0+1, where the two qubit-states |0 and |1 are added and equally-weighted in-phase (a normalization factor of the superposition state is omitted hereinafter for convenience) and the qubit state |1 to a superposition state |0-|1, where the two qubit-states |0 and |1 are added equally-weighted but out of phase. This application of the composite pulse is referred to as a “π/2-pulse”. More generally, a superposition of the two qubit-states |0 and |1 that are added and equally weighted is represented by a point that lies on the equator of the Bloch sphere. For example, the superposition states |0±|1 correspond to points on the equator with the azimuthal angle ϕ being zero and π, respectively. The superposition states that correspond to points on the equator with the azimuthal angle ϕ are denoted as |0+e^iϕ|1. Transformation between two points on the equator (i.e. a rotation about the Z-axis on the Bloch sphere) can be implemented by shifting phases of the composite pulse.

II. Entanglement Formation

FIGS. 5A, 5B, and 5C depict a few schematic structures of collective transverse motional modes (also referred to simply as “motional mode structures”) of a chain 106 of five trapped ions, for example. Here, the confining potential due to a static voltage V_S applied to the end-cap electrodes 210 and 212 is weaker compared to the confining potential in the radial direction. The collective motional modes of the chain 106 of trapped ions in the transverse direction are determined by the Coulomb interaction between the trapped ions combined with the confining potentials generated by the ion trap 200. The trapped ions undergo collective transversal motions (referred to as “collective transverse motional modes,” “collective motional modes,” or simply “motional modes”), where each mode has a distinct energy (or equivalently, a frequency) associated with it. A motional mode having the p-th lowest frequency is hereinafter referred to as |n_php, where n_ph denotes the number of motional quanta (in units of energy excitation, referred to as phonons) in the motional mode, and the number of motional modes P in a given transverse direction is equal to the number of trapped ions N in the chain 106. FIGS. 5A-5C schematically illustrate examples of different types of collective transverse motional modes that may be experienced by five trapped ions that are positioned in a chain 106. FIG. 5A is a schematic view of the common motional mode P having the highest energy where P is the number of motional modes. In the common motional mode |n_phP, all ions oscillate in phase in the transverse direction. FIG. 5B is a schematic view of the tilt motional mode |n_phP-1 which has the second highest energy. In the tilt motional mode, ions on opposite ends move out of phase in the transverse direction (i.e., in opposite directions). FIG. 5C is a schematic view of the higher-order motional mode |n_phP-3 which has a lower energy than that of the tilt motional mode |n_phP-1, and in which the ions move in a more complicated mode pattern.

It should be noted that the particular configuration described above is just one among several possible examples of a trap for confining ions according to the present disclosure and does not limit the possible configurations, specifications, or the like of traps according to the present disclosure. For example, the geometry of the electrodes is not limited to the hyperbolic electrodes described above. In other examples, a trap that generates an effective electric field causing the motion of the ions in the radial direction as harmonic oscillations may be a multi-layer trap in which several electrode layers are stacked, and an RF voltage is applied to two diagonally opposite electrodes, or a surface trap in which all electrodes are located in a single plane on a chip. Furthermore, a trap may be divided into multiple segments, adjacent pairs of which may be linked by shuttling one or more ions or coupled by photon interconnects. A trap may also be an array of individual trapping regions arranged closely to each other on a micro-fabricated ion trap chip. In some embodiments, the quadrupole potential has a spatially varying DC component in addition to the RF component described above.

In an ion trap quantum computer, the motional modes may act as a data bus to mediate entanglement between two qubits and this entanglement is used to perform an XX gate operation. That is, each of the two qubits is entangled with the motional modes, and then the entanglement is transferred to an entanglement between the two qubits by using motional sideband excitations, as described below. FIGS. 6A and 6B schematically depict views of a motional sideband spectrum for an ion in the chain 106 in a motional mode |n_php having frequency ω_p according to one embodiment. As illustrated in FIG. 6B, when the detuning frequency of the composite pulse is zero (i.e., the frequency difference between the first and second laser beams is tuned to the carrier frequency, δ=0), simple Rabi flopping between the qubit states |0 and |1 (carrier transition) occurs. When the detuning frequency of the composite pulse is positive (i.e., the frequency difference between the first and second laser beams is tuned higher than the carrier frequency, δ=μ>0, referred to as a blue sideband, by cop), Rabi flopping between combined qubit-motional states |0|n_php and |1|n_ph+1_p occurs (i.e., a transition from the p-th motional mode with n_ph-phonon excitations denoted by |n_php to the p-th motional mode with n_ph+1-phonon excitations denoted by |n_ph+1_p occurs while the qubit state |0 flips to |1). When the detuning frequency of the composite pulse is negative (i.e., the frequency difference between the first and second laser beams is tuned lower than the carrier frequency by δ=−μ<0, referred to as a red sideband, by ω_p), Rabi flopping between combined qubit-motional states |0|n_php and |1|n_ph−1_p occurs (i.e., a transition from |n_php to |n_ph−1 with one less phonon excitations occurs while the qubit state |0 flips to |1). A π/2-pulse on the blue sideband applied to a qubit transforms the combined qubit-motional state |0|n_ph)_p into a superposition of |0|n_ph)_p and |1|n_ph+1)_p. A π/2-pulse on the red sideband applied to a qubit transforms the combined qubit-motional state |0|n_php into a superposition of |0|n_php and |1|n_ph−1)_p. When the two-photon Rabi frequency Ω(t) is smaller as compared to the detuning frequency δ=±μ, the blue sideband transition or the red sideband transition may be selectively driven. Thus, a qubit can be entangled with a desired motional mode by applying the right type of pulse, such as a π/2-pulse, which can be subsequently entangled with another qubit, leading to an entanglement between the two qubits that is needed to perform an XX-gate operation in an ion trap quantum computer.

By controlling and/or directing transformations of the combined qubit-motional states as described above, an XX-gate operation may be performed on two qubits (i-th and j-th qubits). In general, the XX-gate operation (with maximal entanglement) respectively transforms two-qubit states |0_i|0_j, |0_i|1_j, |1_i|0_j, and |1_i|1_j as follows:

$\begin{array}{c}{{{{{{❘0\rangle }_i❘0\rangle }_j\to ❘0\rangle }_i❘0\rangle }_j-i❘1\rangle }_i❘1\rangle }_j\\ {{{{{{❘0\rangle }_i❘1\rangle }_j\to ❘0\rangle }_i❘1\rangle }_j-i❘1\rangle }_i❘0\rangle }_j\\ {{{{{{❘1\rangle }_i❘0\rangle }_j\to ❘1\rangle }_i❘0\rangle }_j-i❘0\rangle }_i❘1\rangle }_j\\ {{{{{{❘1\rangle }_i❘1\rangle }_j\to ❘1\rangle }_i❘1\rangle }_j-i❘0\rangle }_i❘0\rangle }_j\end{array}$

For example, when the two qubits (i-th and j-th qubits) are both initially in the hyperfine ground state |0 (denoted as |0_i|0_j) and subsequently a π/2-pulse on the blue sideband is applied to the i-th qubit, the combined state of the i-th qubit and the motional mode |0_i|n_php is transformed into a superposition of |0_i|_php and |1_i|_ph+1, and thus the combined state of the two qubits and the motional mode is transformed into a superposition of |0_i|0_j|n_php and |1_i|0_j|n_ph+1_p. When a π/2-pulse on the red sideband is applied to the j-th qubit, the combined state of the j-th qubit and the motional mode |0_j|n_php is transformed to a superposition of |0_j|n_php and |1_j|n_ph−1_p and the combined state |0_i|0_j|n_php is transformed into a superposition of |0_i|0_j|n_php and |0_i|1_j|n_ph−1_p.

Thus, applications of a π/2-pulse on the blue sideband on the i-th qubit and a π/2-pulse on the red sideband on the j-th qubit may transform the combined state of the two qubits and the motional mode |0_i|_0j|n_php into a superposition of |0_i|0_j|n_php and |1_i|1_j|n_phP, the two qubits now being in an entangled state. For those of ordinary skill in the art, it should be clear that two-qubit states that are entangled with a motional mode having a different number of phonon excitations from the initial number of phonon excitations (i.e., |1_i|0_j|n_ph+1_p and |0_i|0_j|_ph−1) can be removed by a sufficiently complex pulse sequence, and thus the combined state of the two qubits and the motional mode after the XX-gate operation may be considered disentangled as the initial number of phonon excitations n_ph in the p-th motional mode stays unchanged at the end of the XX-gate operation. Thus, qubit states before and after the XX-gate operation will be described below generally without including the motional modes.

More generally, the combined state of i-th and j-th qubits, transformed by the application of control pulses described by pulse functions g_i(t) and g_j(t), respectively, for duration τ (referred to as the “gate time”), can be described in terms of a degree of entanglement (also referred to as a “gate angle”) χ_ij as follows:

$\begin{array}{c}{{{{{{❘0\rangle }_i❘0\rangle }_j\to \cos \left(2{\chi }_{\mathrm{ij}}\right)❘0\rangle }_i❘0\rangle }_j-i\sin \left(2{\chi }_{\mathrm{ij}}\right)❘1\rangle }_i❘1\rangle }_j\\ {{{{{{❘0\rangle }_i❘1\rangle }_j\to \cos \left(2{\chi }_{\mathrm{ij}}\right)❘0\rangle }_i❘1\rangle }_j-i\sin \left(2{\chi }_{\mathrm{ij}}\right)❘1\rangle }_i❘0\rangle }_j\\ {{{{{{❘1\rangle }_i❘0\rangle }_j\to \cos \left(2{\chi }_{\mathrm{ij}}\right)❘1\rangle }_i❘0\rangle }_j-i\sin \left(2{\chi }_{\mathrm{ij}}\right)❘0\rangle }_i❘1\rangle }_j\\ {{{{{{❘1\rangle }_i❘1\rangle }_j\to \cos \left(2{\chi }_{\mathrm{ij}}\right)❘1\rangle }_i❘1\rangle }_j-i\sin \left(2{\chi }_{\mathrm{ij}}\right)❘0\rangle }_i❘0\rangle }_j\end{array}$

where

${\chi }_{\mathrm{ij}}={\sum }_{p=1}^P{\eta }_p^i{\eta }_p^j{\int }_0^{\tau }d{t}_2{\int }_0^{t_2}d{t}_1{g}_i\left(t_1\right)g_i\left(t_2\right)\sin \left[{\omega }_p\left(t_2-t_1\right)\right],$

and η_pⁱ is the Lamb-Dicke parameter that quantifies the coupling strength between the i-th qubit and the p-th motional mode having the frequency ω_p, and P is the number of the motional modes (equal to the number N of ions in the chain 106).

The entanglement interaction between two qubits described above can be used to perform an XX-gate operation. The XX-gate operation (XX gate) along with single-qubit gate operations (R gates) forms a set of gates {R, XX} that can be used to build a quantum computer that is configured to perform desired computational processes. Among several known sets of logic gates by which any quantum algorithm can be decomposed, a set of logic gates, commonly denoted as {R, XX}, is native to a quantum computing system of trapped ions described herein. Here, the R gate corresponds to manipulation of individual qubit states of trapped ions, and the XX gate (also referred to as an “entangling gate”) corresponds to manipulation of the entanglement of two trapped ions.

III. Trapped-ion Architecture and Hamiltonian

The trapped-ion architecture, i.e., chains of trapped ions, coherently controlled via Raman hyperfine transitions, is one of the most promising routes to scalable quantum computing. This quantum computer architecture is used both in laboratory experiments and in the emerging quantum computing industry. For both the current era of noisy intermediate-scale quantum computing and the anticipated era of fault-tolerant, error-corrected quantum computing, two-qubit gates of the highest possible fidelity are essential. While fault-tolerant quantum computing and quantum error-correction may, in principle, be achieved with two-qubit gates of modest fidelity, the overhead, i.e., the number of physical qubits required for one error-corrected logical qubit depends on the native fidelity of the physical gates and may be enormous for physical two-qubit gates of only modest fidelity. A reasonable amount of overhead in fault-tolerant quantum computing can be achieved only if the physical two-qubit gates themselves have a high native fidelity. Typically, for realistic, tolerable overhead, a physical two-qubit infidelity of ≤10⁻⁴ is recommended. Two-qubit gate infidelities close to this target have indeed already been achieved. However, the experimental demonstrations of high-fidelity two-qubit gates are restricted to two-qubit gates in very short ion chains, i.e., chains consisting of up to four ions. Moreover, to date, even in these cases two-qubit gate infidelities of ≤10⁻⁴ have not yet been achieved experimentally. Two adversaries stand in the way of achieving two-qubit gate infidelities ≤10⁻⁴: Random noise and deterministic, coherent control errors. Even in the shortest chains (two to four ions stored simultaneously), the target of ≤10⁻⁴ infidelity may not be achieved if the Hamiltonian used to design the control pulses for two-qubit gate implementation does not accurately enough reflect the reality of the quantum computer's hardware implementation. What this means is illustrated in FIG. 7. The actual quantum computer, i.e., the reality, is governed by a Hamiltonian Ĥ_R. Reality can never be captured exactly. It can only be modeled approximately. Consequently a model of the quantum computer is constructed, replacing the unknown Hamiltonian Ĥ_R with Ĥ_M, where it is hoped that Ĥ_M˜Ĥ_R to a high accuracy. Both the quantum computer (Ĥ_R) and its model (Ĥ_M) are controlled by control pulses that are constructed on the basis of a Hamiltonian Ĥ_C. Ideally, Ĥ_C=Ĥ_R. However, since Ĥ_R is unknown, the best possible control Hamiltonian is Ĥ_C=Ĥ_M. However, in most cases Ĥ_M is too complicated to use for efficiently constructing control pulses that frequently also have to be computed in real time. Therefore, Ĥ_C is chosen as a compromise, close enough to Ĥ_M to ensure acceptable control of the quantum computer, but simple enough to ensure efficient control-pulse construction.

One of the most basic tasks of a quantum computer is to construct two-qubit gates. In the embodiments described herein, two-qubit XX gates constructed according to the Mølmer-Sørensen scheme are considered. Ideally, given an input state, |ψ_in, the quantum computer, governed by Ĥ_R, is expected to turn |ψ_in into |ψ_out=XX|ψ_in, where XX is the ideal XX gate. However, since Ĥ_R≠Ĥ_C, this will not happen in practice. Therefore, |ψ_out^R, i.e., the output state produced by Ĥ_R, is different from XX|ψ_in. However, assuming that Ĥ_M is very close to Ĥ_R, it is possible, at least approximately, to assess the quality of the control pulses by computing the model output state |ψ_out^M and comparing it with the ideal output state XX|ψ_in, for instance, by computing the overlap |ψ_out^M|XX|ψ_in|². Only for H_C=Ĥ_M is it expected that this overlap equals 1. However, if Ĥ_M is close to Ĥ_R, but Ĥ_C is chosen as the manageable standard Hamiltonian Ĥ_S, it is expected that this overlap and other fidelity measures (see Section VI) differ from 1. In this case, because of the insufficient quality of the control pulses, the resulting two-qubit gates may not be accurate enough for the target infidelity (for instance, <10⁻⁴). To this end, constructing a realistic model Hamiltonian Ĥ_M for Raman-controlled trapped ions, it is shown that currently employed control-pulse construction techniques based on the standard Hamiltonian Ĥ_S, i.e., Ĥ_C=Ĥ_S are not accurate enough to achieve the ≤10⁻⁴ infidelity target for realistic fault-tolerant quantum computing. This is shown by constructing control pulses on the basis of Ĥ_S, which are then used to compute two-qubit XX gates in a 7-ion chain governed by Ĥ_M, assuming that Ĥ_M is sufficiently close to Ĥ_R. The XX gate simulations of the 7-ion chain are performed with a gate-simulator code that is accurate on the 10⁻⁷ level and takes both computational levels and phonon states explicitly into account. This way it is shown that even assuming ideal conditions, i.e., neglecting all incoherent noise sources, control-pulse construction needs to be improved to reach the ≤10⁻⁴ infidelity goal.

There are two principal methods of implementing a trapped-ion chain quantum computer using the stimulated Raman scheme: Phase sensitive and phase insensitive. Both schemes have advantages and disadvantages. The phase-sensitive scheme is focused since it is technically more straightforward to implement and is currently used by commercial quantum computing companies such as IonQ, Inc. Therefore, the focus is to isolate the leading coherent error sources in the phase-sensitive architecture and to construct methods that allow us to eliminate these error sources. It is shown that on the basis of a new linear scheme of pulse construction, and in the absence of all incoherent error sources, XX gate infidelities ≤10⁻⁴ can be reached.

The entangling gate can be constructed from the model Hamiltonian Ĥ_M (referred to as the full Mølmer-Sørensen Hamiltonian) that governs a chain of ions in a linear Paul trap (e.g., the chain 106 of trapped ions, illuminated simultaneously by two laser beams, one red detuned and one blue detuned). The effective two-level Hamiltonian:

$\begin{array}{cc}\widehat{H}={\widehat{H}}_0+{\sum }_{j=1}^N\left(\frac{\hslash {\Omega }_{\mathrm{eff}}\left(t\right)}{2}\right)e^{-i\left[\left(\Delta \overrightarrow{k}\right)·{\overrightarrow{r}}_j-\mu t-\mathrm{\Delta \phi }\right]}{\hat{\sigma }}_{-}^{\left(j\right)}+h.c.,& \left(1\right)\end{array}$

can be obtained via adiabatic elimination according to the Raman Λ scheme. Here, N is the number of ions in the chain,

$\begin{array}{cc}{\widehat{H}}_0={\sum }_{p=1}^N{\omega }_p{\hat{a}}_p^{†}{\hat{a}}_p& \left(2\right)\end{array}$

is the phonon Hamiltonian with the phonon frequencies ω_p and associated phonon creation and destruction operators {circumflex over (α)}_p^† and â_p, respectively, the j sum is over the ions in the chain, Ω_eff(t) is the time-dependent effective Rabi frequency, Δ{right arrow over (k)} is the wavenumber difference between the two Raman lasers, {right arrow over (r)}_j is the position operator of ion number j in the chain, and Δφ is the phase difference between the two Raman lasers.

To construct a Mølmer-Sørensen (MS) gate, the ions are illuminated simultaneously with blue (+μ) and red (−μ) shifted light. In the interaction representation with respect to Ĥ₀, the MS Hamiltonian is:

$\begin{array}{cc}{\widehat{H}}_{MS}\left(t\right)={\sum }_j{\Omega }_j\cos \left(\mu t-{\varphi }_j^m\right)\left\{\cos \left[{\sum }_p{\eta }_p^j\left({\hat{a}}_p^{†}·e^{i{\omega }_{p}t}+{\hat{a}}_p{e}^{-i{\omega }_{p}t}\right)\right]{\hat{\sigma }}_y^{\left(j\right)}+\sin \left[{\sum }_p{\eta }_p^j\left({\hat{a}}_p^{†}{e}^{i{\omega }_{p}t}+{\hat{a}}_p{e}^{-i{\omega }_{p}t}\right)\right]{\hat{\sigma }}_x^{\left(j\right)}\right\},& \left(3\right)\end{array}$

where μ is the detuning frequency, η_p^j are the Lamb-Dicke parameters, and {circumflex over (σ)}_x, {circumflex over (σ)}_y, {circumflex over (σ)}_z are the Pauli operators. The Hamiltonian (3) may be generalized according to:

$\begin{array}{cc}{\Omega }_j\cos \left(\mu t-{\varphi }_j^m\right)\to {\Omega }_j\left(t\right)\cos \left[{\int }_0^t\mu \left(t^{\prime }\right){\mathrm{dt}}^{\prime }+{\varphi }_j^{\left(0\right)}\right]\to g_j\left(t\right),& \left(4\right)\end{array}$

where g_j(t) may be any time-dependent pulse function, i.e., it includes amplitude-modulated (AM) pulses, frequency-modulated (FM) pulses, phase-modulated (PM) pulses, and simultaneously amplitude- and frequency-modulated (AMFM) pulses. Thus, written in terms of the most general pulse function g_j(t), the MS Hamiltonian (3) becomes:

$\begin{array}{cc}{\widehat{H}}_{MS}={\sum }_j{g}_j\left(t\right)\left\{\cos \left[\widehat{V_j}\left(t\right)\right]{\hat{\sigma }}_y^{0)}+\sin \left[\widehat{V_j}\left(t\right)\right]{\hat{\sigma }}_x^{\left(j\right)}\right\},& \left(5\right)\end{array}$

where, for later convenience, the Lamb-Dicke operators are defined as:

$\begin{array}{cc}{\hat{V}}_j\left(\tau \right)={\sum }_p{\eta }_p^j\left({\hat{a}}_p^{†}{e}^{i{\omega }_{p}t}+{\hat{a}}_p{e}^{-i{\omega }_{p}t}\right),& \left(6\right)\end{array}$

which satisfy

$\begin{array}{cc}\left[{\hat{V}}_j\left(t_1\right),{\hat{V}}_l\left(t_2\right)\right]=-2i{\sum }_p{\eta }_p^j{\eta }_p^l\sin \left[{\omega }_p\left(t_1-t_2\right)\right].& \left(7\right)\end{array}$

Expansion of the MS Hamiltonian Ĥ_MS in powers of {circumflex over (V)} leads to a family of model Hamiltonians:

$\begin{array}{cc}{\widehat{H}}_M^{\left(N_c,N_s\right)}={\sum }_j{g}_j\left(t\right)\left\{{\sum }_{\begin{array}{c}n=0\\ \mathrm{neven}\end{array}}^{N_c}{\left(-1\right)}^{n/2}\frac{{\hat{\sigma }}_y^{\left(j\right)}}{n!}{\hat{V}}_j^n\left(t\right)+{\sum }_{\begin{array}{c}m=1\\ \mathrm{modd}\end{array}}^{N_s}{\left(-1\right)}^{\left(m-1\right)/2}\frac{{\hat{\sigma }}_x^{\left(j\right)}}{m!}{\hat{V}}_j^m\left(t\right)\right\},& \left(8\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}{\widehat{H}}_M^{\left(\infty ,\infty \right)}={\widehat{H}}_{MS}.& \left(9\right)\end{array}$

For later convenience, the standard Hamiltonian is also defined as:

$\begin{array}{cc}{\widehat{H}}_S={\widehat{H}}_M^{\left(-2,1\right)}={\sum }_j\hslash g_j\left(t\right)\widehat{V_j}\left(t\right){\hat{\sigma }}_x^{\left(j\right)},& \left(10\right)\end{array}$

which is frequently used as the basis for control-pulse construction. Here, the superscript N_c=−2 codes for the complete suppression (omission) of the cosine term in (8). The test pulses used in the example are also all constructed on the basis of the standard Hamiltonian Ĥ_S (as discussed in detail below in Section IV). The expanded model Hamiltonians (8) are use in Section VII to assess the accuracy of expansions of the MS Hamiltonian Ĥ_MS in powers of the Lamb-Dicke operators {circumflex over (V)}_j(t).

IV. Construction of Control Pulses for Entangling Gate Operations

Quantum computation can be performed in a quantum computing system, such as the ion trap quantum computing system 100, using a set of quantum gate operations including single-qubit gate operations (R gates) and two-qubit gate operations, such as XX-gate operations (XX gates). Although the methods for applying such basic building blocks of quantum computation have been established, there are errors, which result from experimental parameter drift or fluctuations, such as the vibrational mode frequencies of an ion chain, in the hardware of the quantum computing system. These errors are mainly due to the lack of knowledge about the changes in the computing environment and the properties of quantum computing hardware within the quantum computing system. Thus, stabilizing gate operations, i.e., the task of computing and implementing control pulses in the quantum computing system to be robust against the errors is needed to provide scalable and reliable quantum computation results.

While the pulse functions g(t) may be any pulse functions constructed by any of the known pulse construction schemes (e.g., AM, PM, FM, AMFM, etc.), to be specific, and for the purpose of illustrating the methods, but not restricting the generality of the methods in any way, in the embodiments described herein, the pulse functions g(t) of a control pulse are taken to be sine-AMFM pulses:

$\begin{array}{cc}g\left(t\right)={\sum }_{n_{\min }}^{n_{\max }}{B}_n\sin \left({\omega }_{n}t\right),& \left(11\right)\end{array}$

where the basis spans the states n∈{n_min, n_max}, the Fourier coefficients B_n are real amplitudes,

$\begin{array}{cc}{\omega }_n=\frac{2\pi n}{\tau }& \left(12\right)\end{array}$

are the basis frequencies, and r is the gate time. The pulse functions (11) fulfill the “carrier condition”

$\begin{array}{cc}{\int }_0^{\tau }g\left(t\right)\mathrm{dt}=0,& \left(13\right)\end{array}$

which, as shown in Section IX, ensures that significant, first-order contributions to the gate infidelity due to the “carrier term”, i.e., the n=0 term in the expansion (8) of the MS Hamiltonian, are eliminated, and higher-order contributions to the gate infidelity are much reduced.

The phase-space closure condition requires:

$\begin{array}{cc}{\alpha }_p^j=-{\eta }_p^j{\int }_0^{\tau }{g}_j\left(t\right)e^{i{\omega }_{p}t}\mathrm{dt}=0\left(p,j=1,\dots ,N\right),& \left(14\right)\end{array}$

such that the trapped ions that are displaced from their initial positions as the motional modes are excited by the delivery of the control pulse return to their initial motional states. This condition must be imposed for all motional modes (p=1,2, . . . , N) and all ions (j=1, . . . ,N) in the chain of trapped ions. The sine-AMFM pulses (11) are known to fulfill (14) exactly. The number of ions in the chain is chosen as N=7, as an example, because it is large enough to be realistic, yet small enough to enable direct numerical simulations of the quantum dynamics of the chain, explicitly including a large phonon space. In addition, in the embodiments described herein, g_j(t)=g(t) (=1, . . . , N) is chosen for illustrative purposes but is not intended to be a restriction of generality. That is, here, whenever a two-qubit gate is constructed between two ions, it is assumed that each of the two ions is irradiated with laser light controlled by the same pulse function g(t). This is a common choice, used in laboratories and in commercial quantum computers.

Among the pulse functions g(t) that fulfill the phase-space closure condition and the carrier condition, the pulse function g(t) that satisfies the gate angle condition, which requires the pulse g(t) to produce the desired gate angle χ according to:

$\begin{array}{cc}\chi =2{\sum }_{p=1}^N{\eta }_p^{j_1}{\eta }_p^{j_2}{\int }_0^{\tau }\mathrm{dt}{\int }_0^t{\mathrm{dt}}^{\prime }g\left(t\right)g\left(t^{\prime }\right)\sin \left[{\omega }_p\left(t-t^{\prime }\right)\right],& \left(15\right)\end{array}$

is selected. A maximally entangling gate is obtained for χ=π/4. The mode frequencies ω_p, p=1, . . . ,7, and the Lamb-Dicke parameters η_p^j, p,j=1, . . . ,7, are listed in Tables I and II, respectively.

TABLE I
Motional-mode frequencies ωp/(2π) (in MHz), p = 1, . . . ,
7, used to construct the 7-ion AMFM test pulses:
mode number p ωp/(2π) (MHz)
1             2.953
2             2.966
3             2.984
4             3.006
5             3.029
6             3.049
7             3.060

TABLE II
Lamb-Dicke parameters ηpj, p, j = 1, . . . ,
7, used to construct the 7-ion AMFM test pulses:
p	j = 1	j = 2	j = 3	j = 4	j = 5	j = 6	j = 7
1	0.01033	−0.03355	0.05478	−0.06313	0.05478	−0.03355	0.01033
2	0.02226	−0.05627	0.05036	0.00000	−0.05036	0.05627	−0.02226
3	−0.03494	0.05644	0.00760	−0.05823	0.00760	0.05644	−0.03494
4	0.04644	−0.03074	−0.05488	0.00000	0.05488	0.03074	−0.04644
5	0.05503	0.00994	−0.03678	−0.05637	−0.03678	0.00994	0.05503
6	−0.05848	−0.04491	−0.02433	0.00000	0.02433	0.04491	0.05848
7	0.04143	0.04143	0.04143	0.04143	0.04143	0.04143	0.04143

For computing analytical infidelities the following identities that follow directly from (14) with the expansion (11) and hold for all p=1, . . . , N are used:

$\begin{array}{cc}{\int }_0^{\tau }g\left(t\right)e^{\pm i{\omega }_{p}t}\mathrm{dt}=0,& \left(16\right)\end{array}$ $\begin{array}{cc}{\int }_0^{\tau }g\left(t\right)\sin \left({\omega }_{p}t\right)\mathrm{dt}=0,& \left(17\right)\end{array}$ $\begin{array}{cc}{\int }_0^{\tau }g\left(t\right)\cos \left({\omega }_{p}t\right)\mathrm{dt}=0,& \left(18\right)\end{array}$ $\begin{array}{cc}{\sum }_{n_{\min }}^{n_{\max }}{B}_n\left(\frac{{\omega }_n}{{\omega }_n^2-{\omega }_p^2}\right)=0.& \left(19\right)\end{array}$

For later use in the following sections, the following functions are defined:

$\begin{array}{cc}\tilde{\chi }={\sum }_{\mathrm{jp}}{\left({\eta }_p^j\right)}^2{\int }_0^{\tau }\mathrm{dt}{\int }_0^t{\mathrm{dt}}^{\prime }g\left(t\right)g\left(t^{\prime }\right)\sin \left[{\omega }_p\left(t-t^{\prime }\right)\right],& \left(20\right)\end{array}$ $\begin{array}{cc}G\left(t\right)={\int }_0^{t}g\left(t^{\prime }\right){\mathrm{dt}}^{\prime }={\sum }_n\left(\frac{B_n}{{\omega }_n}\right)\left[1-\cos \left({\omega }_{n}t\right)\right]=2{\sum }_n\left(\frac{B_n}{{\omega }_n}\right){\sin }^2\left({\omega }_{n}t/2\right),& \left(21\right)\end{array}$ $\begin{array}{cc}Q\left(w\right)={\int }_0^{\tau }g\left(t\right)e^{\mathrm{iwt}}\mathrm{dt}=\left[e^{iw\tau }-1\right]{\sum }_n{B}_n\left(\frac{{\omega }_n}{w^2-{\omega }_n^2}\right),& \left(22\right)\end{array}$ $\begin{array}{cc}f\left(w\right)={\int }_0^{\tau }g\left(t\right)G\left(t\right)e^{\mathrm{iwt}}\mathrm{dt}=\left(\frac{1}{2}\right)\left[e^{lw\tau }-1\right]{\sum }_{nm}\left(\frac{B_n{B}_m}{{\omega }_n}\right)\left\{\frac{{\omega }_m-{\omega }_n}{{\left({\omega }_m-{\omega }_n\right)}^2-w^2}+\frac{{\omega }_m+{\omega }_n}{{\left({\omega }_m+{\omega }_n\right)}^2-w^2}\right\},& \left(23\right)\end{array}$ $\begin{array}{cc}{S}_p\left(w\right)={\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)\sin \left[{\omega }_p\left(t_1-t_2\right)\right]e^{{\mathrm{iwt}}_1}={\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)\sin \left[{\omega }_p\left(t_1-t_2\right)\right]e^{{\mathrm{iwt}}_2}={\sum }_{\mathrm{nm}}\frac{B_n{B}_m}{2\left(w^2-{\omega }_m^2\right)}\left\{\left(\frac{2{\omega }_n{\omega }_m}{{\omega }_n^2-4w^2}\right)e^{\mathrm{iw}\tau }\sin \left(w\tau \right)+{\mathrm{iw}}^2\left(e^{\mathrm{iw}\tau }-1\right)\left[\frac{1}{{\left({\omega }_n-{\omega }_m\right)}^2-w^2}-\frac{1}{{\left({\omega }_n-{\omega }_m\right)}^2-w^2}\right]\right\},& \left(24\right)\end{array}$ $\begin{array}{cc}{J}_p={\int }_0^{\tau }g\left(t\right)G^2\left(t\right)e^{i{\omega }_{p}t}\mathrm{dt},& \left(25\right)\end{array}$ $\begin{array}{cc}Z\left(w_1,w_2\right)={\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)e^{{\mathrm{iw}}_1{t}_1}{e}^{{\mathrm{iw}}_2{t}_2}={\sum }_{\mathrm{nm}}{B}_n{B}_m\left(\frac{1}{{\omega }_m^2-w_2^2}\right)\left\{\left(\frac{{\omega }_n{\omega }_m}{w_1^2-{\omega }_n^2}\right)\left[e^{{\mathrm{iw}}_1\tau }-1\right]-\left(\frac{{\omega }_n{\omega }_m\left({\omega }_m^2-{\omega }_n^2+w_1^2+4w_1{w}_2+3w_2^2\right)}{\left[{\left({\omega }_n+{\omega }_m\right)}^2-{\left(w_1+w_2\right)}^2\right]\left[{\left({\omega }_n-{\omega }_m\right)}^2-{\left(w_1+w_2\right)}^2\right]}\right)\left[e^{i\left(w_1+w_2\right)\tau }-1\right]\right\},& \left(26\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}\Phi \left[\eta ,g\right]={\sum }_{p=1}^N{\eta }_p^1{\eta }_p^2{\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)G\left(t_2\right)\sin \left[{\omega }_p\left(t_1-t_2\right)\right]=\frac{\chi \tau }{4\pi }{\sum }_{n=n_{\min }}^{n_{\max }}\frac{B_n}{n},& \left(27\right)\end{array}$

where (15) is used and

$\begin{array}{cc}{\sum }_n\left(\frac{B_n}{n}\right)=\left(\frac{2\pi }{{\tau }^2}\right){\int }_0^{\tau }\mathrm{dt}{\int }_0^t{\mathrm{dt}}^{\prime }g\left(t^{\prime }\right){\mathrm{dt}}^{\prime }.& \left(28\right)\end{array}$

V. Gate Simulator

Given a Hamiltonian and enough computer power, the time evolution of the computational states including phonon excitations can always be solved to any desired accuracy. For a 7-ion chain, the complete state |ψ(t) can be expanded in the combined Hilbert space of computational states |a, b, a, b∈{0,1}, and phonon states |m₁m₂ . . . m₇, m_j=0,1, . . . , m_j^max, j=1, . . . ,7, according to:

$\begin{array}{cc}|\psi \left(t\right)\rangle ={\sum }_{\underset{m_1{m}_2\dots m_7}{a,b}}{A}_{m_1{m}_2\dots m_7}^{\left(a,b\right)}\left(t\right)|a,b\rangle |m_1{m}_2\dots m_7\rangle ,& \left(29\right)\end{array}$

where m_j^max≥0 is the maximal phonon occupation number included in the basis. The set {m₁^max, . . . , m₇^max} is referred to as the phonon scheme. With (29), the time-dependent Schrödinger equation results in the amplitude equations:

$\begin{array}{cc}i\hslash {\stackrel{.}{A}}_{n_1{n}_2\dots n_7}^{\left(a,b\right)}\left(t\right)={\sum }_{\underset{m_1{m}_2\dots m_7}{a^{\prime },b^{\prime }}}\langle a,b|\langle n_1\dots n_7|\hat{H}\left(t\right)|m_1\dots m_7\rangle |a^{\prime },b^{\prime }\rangle A_{m_1\dots m_7}^{\left(a^{\prime },b^{\prime }\right)}\left(t\right),& \left(30\right)\end{array}$

where Ĥ(t) may be any of the model Hamiltonians defined in Section III. The set of equations (30) are ordinary first-order differential equations that may be solved with any standard numerical differential equations solver. For simplicity and straightforward numerical error control, an elementary fourth-order Runge-Kutta solver with constant step size is chosen. With this simple integrator, the inventors have been able to obtain a relative accuracy of the numerical solution of better than 10⁻⁷. It should be noted that the gate simulator is not limited to 7 ions. It scales trivially to any number of ions.

Convergence in the phonon scheme is assessed by using the gate simulator to compute |ψ(τ) for the Hamiltonian Ĥ(t)=Ĥ_S(t)=Ĥ^−2,1)(t). In this case, full phase-space closure is required, because the control pulses are constructed on the basis of the same Hamiltonian that governs the time evolution. Only if enough phonon states are included in the basis used by the gate-simulator, will the phase space be closed and thus the phonon-space truncation according to the phonon scheme will guarantee an accurate result.

VI. Fidelity Measures

The ideal XX gate propagator is:

$\begin{array}{cc}{\hat{U}}_{\mathrm{ideal}}=e^{i\chi {\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}},& \left(31\right)\end{array}$

where χ is the gate angle. For starting state |ψ₀ this produces the ideal final state:

$\begin{array}{cc}|{\psi }_{\mathrm{ideal}}\rangle ={\hat{U}}_{\mathrm{ideal}}|{\psi }_0\rangle .& \left(32\right)\end{array}$

However, a given gate pulse g(t), in general, produces a gate propagator of the form:

$\begin{array}{cc}{\hat{U}}_{\mathrm{actual}}=e^{i\left[{\mathrm{\chi \sigma }}_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}+\lambda Ê\right]},& \left(33\right)\end{array}$

where Ê is a Hermitian error operator and λ is its error strength. In this case, and for given initial state |ψ₀, the state fidelity F_S is defined as:

$\begin{array}{cc}{F}_S=|\langle {\psi }_{\mathrm{ideal}}|{\psi }_{\mathrm{actual}}\rangle {|}^2=|\langle {\psi }_0|{\hat{U}}_{\mathrm{ideal}}^{†}{\hat{U}}_{\mathrm{actual}}|{\psi }_0\rangle {|}^2.& \left(34\right)\end{array}$

There are two types of error operators, i.e., those that commute with the product of the Pauli operators σ_x⁽¹⁾σ_x⁽²⁾ and those that do not. If only a single, commuting error operator Ê is present, the state fidelity F_S, up to second order in the error strength λ, is:

$\begin{array}{cc}{F}_S=|\langle {\psi }_0❘e^{-i{\mathrm{\chi \sigma }}_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}+i\lambda Ê}{e}^{i{\mathrm{\chi \sigma }}_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}+i\lambda Ê}|{\psi }_0\rangle {|}^2=|\langle {\psi }_0|e^{i\lambda Ê}|{\psi }_0\rangle {|}^2=1-{\lambda }^2{\sigma }_{Ê}^2,& \left(35\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}{\sigma }_{Ê}^2=\langle {\psi }_0\left|{Ê}^2\right|{\psi }_0\rangle -\langle {\psi }_0\left|Ê\right|{\psi }_0{\rangle }^2.& \left(36\right)\end{array}$

This means that the state infidelity F_S:

$\begin{array}{cc}{\overline{F}}_S=1-F_S,& \left(37\right)\end{array}$

up to second order in the error strength λ, is:

$\begin{array}{cc}{\stackrel{\_}{F}}_S={\lambda }^2{\sigma }_{\hat{E}}^2.& \left(38\right)\end{array}$

The state infidelity F_S is proportional to the square of the error strength λ.

In case the error operator Ê does not commute with the product of the Pauli operators σ_x¹α_x⁽²⁾, the Baker-Hausdorff-Campbell formula can be used to obtain the state fidelity F_S, up to second order in the error strength λ,

$\begin{array}{cc}{{F_s=❘{\psi }_0❘{\hat{U}}_{\mathrm{ideal}}^{†}{\hat{U}}_{\mathrm{actual}}❘{\psi }_0\rangle {❘}^2=❘\langle {\psi }_0❘e^{-i\chi {\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}}{e}^{-i\chi {\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}+i\lambda \hat{E}}❘{\psi }_0\rangle {❘}^2=❘\langle {\psi }_0❘e^{i\lambda \hat{C}+i{\lambda }^2\hat{D}+\cdots }❘{\psi }_0\rangle ❘}^2=❘\langle {\psi }_0❘1+i\lambda \left(\hat{C}+\lambda \hat{D}\right)-\frac{1}{2}{{\lambda }^2\left(\hat{C}+\lambda \hat{D}\right)}^2+\cdots ❘{\psi }_0\rangle ❘}^2=1-{\lambda }^2{\sigma }_{\hat{C}}^2,& \left(39\right)\end{array}$

where Ĉ and {circumflex over (D)} are Hermitian operators that can be expressed as linear combinations of nested multi-commutators of the product of the Pauli operators σ_x⁽¹⁾σ_x⁽²⁾ and the error operator Ê. This implies that in this case, too, the state infidelity F_S is proportional to the square of the error strength λ.

Most of the error operators Ê that appear in the present context do not commute with the product of the Pauli operators σ_x⁽¹⁾σ_x². However, if the error operator Ê can be written in the form Ê=Â{circumflex over (Ω)}, where  is a Hermitian error operator in the computational space and {circumflex over (Ω)} is a Hermitian error operator in the phonon space, and if the error operator  in the computational space fulfills the anti-commutation relation:

$\begin{array}{cc}\left\{{\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)},\hat{A}\right\}={\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)},\hat{A}+A{\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}=0,& \left(40\right)\end{array}$

the state infidelity F_S can be computed explicitly, which also provides an explicit expression for the operator Ĉ in (39). Examples of error operators  that fulfill (40) are:

$\begin{array}{cc}A\in \left\{{\sigma }_y^{\left(1\right)},{\sigma }_y^{\left(2\right)},{\sigma }_z^{\left(1\right)},{\sigma }_z^{\left(2\right)}{\sigma }_x^{\left(1\right)}{\sigma }_y^{\left(2\right)},{\sigma }_x^{\left(1\right)}{\sigma }_z^{\left(2\right)}{\sigma }_x^{\left(2\right)}{\sigma }_y^{\left(1\right)},{\sigma }_x^{\left(2\right)}{\sigma }_z^{\left(1\right)}\right\}.& \left(41\right)\end{array}$

Since they act in different spaces, the following conditions are satisfied:

$\begin{array}{cc}\begin{array}{cc}\left[{\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)},\hat{\Omega }\right]=0,& \left[\hat{A},\hat{\Omega }\right]\end{array}=0.& \left(42\right)\end{array}$

If (40) is fulfilled, Û_actual becomes:

$\begin{array}{cc}{\hat{U}}_{\mathrm{actual}}=e^{i\chi {\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}+i\lambda \hat{A}\hat{\Omega }}=\cos \left(\hat{\phi }\right)+i\left(\frac{\sin \left(\hat{\phi }\right)}{\hat{\phi }}\right)\left[{\mathrm{\chi \sigma }}_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}+i\lambda \hat{A}\hat{\Omega }\right],& \left(43\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}\hat{\phi }=\sqrt{{\chi }^2+{\lambda }^2{\hat{A}}^2{\hat{\Omega }}^2}.& \left(44\right)\end{array}$

Then, up to second order in the error strength λ, the state fidelity F_S becomes:

$\begin{array}{cc}{F}_S={❘1+i\left(\frac{\lambda }{\chi }\right)\cos \left(\chi \right)\sin \left(\chi \right)\langle {\psi }_0❘\hat{A}\hat{\Omega }❘{\psi }_0\rangle +\left(\frac{\lambda }{\chi }\right){\sin }^2\left(\chi \right)\langle {\psi }_0❘{\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}\hat{A}\hat{\Omega }❘{\psi }_0\rangle -\left(\frac{{\lambda }^2}{2{\chi }^2}\right){\sin }^2\left(\chi \right)\langle {\psi }_0❘{\hat{A}}^2{\hat{\Omega }}^2❘{\psi }_0\rangle ❘}^2.& \left(45\right)\end{array}$

Now, (40) and (42) lead to:

$\begin{array}{cc}\langle {\psi }_0❘{\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}\hat{A}\hat{\Omega }❘{\psi }_0\rangle =-\langle {\psi }_0❘\hat{A}{\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}\hat{\Omega }❘{\psi }_0\rangle =-{\langle {\psi }_0❘\hat{\Omega }{\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}\hat{A}❘{\psi }_0\rangle }^{*}=-{\langle {\psi }_0❘{\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}\hat{A}\hat{\Omega }❘{\psi }_0\rangle }^{*}.& \left(46\right)\end{array}$

This means that ψ₀|σ_x⁽¹⁾σ_x⁽²⁾Â{circumflex over (Ω)}|ψ₀ is purely imaginary, and thus:

$\begin{array}{cc}\langle {\psi }_0❘{\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}\hat{A}\hat{\Omega }❘{\psi }_0\rangle =i\left\{\langle {\psi }_0❘{\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}\hat{A}\hat{\Omega }❘{\psi }_0\rangle \right\}.& \left(47\right)\end{array}$

With this result, the state fidelity F_S can be written, up to second order in the error strength λ, as:

$\begin{array}{cc}{F}_S=1-\left(\frac{{\lambda }^2}{{\chi }^2}\right){\sin }^2\left(\chi \right)\langle {\psi }_0❘{\hat{A}}^2{\hat{\Omega }}^2❘{\psi }_0\rangle +\left(\frac{{\lambda }^2}{{\chi }^2}\right){{\sin }^2\left(\chi \right)\left[\cos \left(\chi \right)\langle {\psi }_0❘\hat{A}\Omega ❘{\psi }_0\rangle -i\sin \left(\chi \right)\langle {\psi }_0❘{\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}\hat{A}\Omega ❘{\psi }_0\rangle \right]}^2=1-{\lambda }^2{\sigma }_{\hat{C}}^2,& \left(48\right)\end{array}$ $\begin{array}{cc}C=\left(\frac{\sin \left(\chi \right)}{\chi }\right)\left[\cos \left(\chi \right)A\Omega -i\sin \left(\chi \right){\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}\hat{A}\hat{\Omega }\right]=\left(\frac{\sin \left(\chi \right)}{\chi }\right)e^{-i\chi {\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}}\hat{A}\hat{\Omega }.& \left(49\right)\end{array}$

It should be noted that, because of (42), and since the error operator  in the computational space is assumed to fulfill (40), iσ_x⁽¹⁾σ_x⁽²⁾Â{circumflex over (Ω)} is a Hermitian operator. So, Ĉ is Hermitian, as it should be, and Ĉ²=Ĉ^†Ĉ. With this result, the state fidelity {circumflex over (F)}_S can be explicitly written as:

$\begin{array}{cc}{F_S={\left(\frac{\mathrm{\lambda sin}\left(\chi \right)}{\chi }\right)}^2\{\langle {\psi }_0❘{\hat{A}}^2{\hat{\Omega }}^2❘{\psi }_0\rangle -\langle {\psi }_0❘e^{-i\chi {\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}}\hat{A}\hat{\Omega }❘{\psi }_0\rangle }^2\}.& \left(50\right)\end{array}$

The explicit form (48) of the state fidelity F_S confirms the general result (39) in cases where the condition (40) is met. Many of the most important error operators occurring satisfy (40) and thus (50) is applicable.

In the case of a two-qubit gate, the error operator Ê takes the form Ê=Σ_j=1,2 Â^(j){circumflex over (Ω)}^(j). In this case (48), (49), and (50) may immediately be generalized using the substitution:

$\begin{array}{cc}\hat{A}\hat{\Omega }\to {\sum }_{j=1,2}{\hat{A}}^{\left(j\right)}{\hat{\Omega }}^{\left(j\right)}.& \left(51\right)\end{array}$

The output of the gate simulator code is the complete state |ψ_out^M(τ), which includes computational states and phonon states. This way, the complete state information is available, which can be used to compute the state fidelity defined in (34) with |ψ_actual=|ψ_out^M(τ). In Section VII the state fidelity F_S is used to assess the quality of the various model Hamiltonians Ĥ_M^Nc,Ns).

The state fidelity F_S, as defined in (34), depends on the initial state |ψ₀. A more global measure of the fidelity of the quantum process that implements the two-qubit XX gate is the process fidelity

$\begin{array}{cc}{F}_P=\frac{1}{16}\mathrm{Tr}\left[E_{\mathrm{exact}}^{†}{E}_{\mathrm{actual}}\right],& \left(52\right)\end{array}$

where E_exact is the exact XX gate process and E_actual is the actual two-qubit XX-gate process as computed with the gate simulator as described in Section V.

A related fidelity measure, also used in Section VII, is the average gate fidelity F_G, which is defined as:

$\begin{array}{cc}{F}_G=\frac{1}{80}{\sum }_{j=1}^{16}Tr\left[{\widehat{U}}_{exact}{\widehat{U}}_j^{†}{\widehat{U}}_{exact}^{†}E\left({\hat{U}}_j\right)\right],& \left(53\right)\end{array}$

where Û_j, j=1, . . . ,16, is an operator basis in the computational space.

In addition to the state infidelity F_S=1−F_S, defined in (37), the average gate infidelity F_G and the process gate infidelity F_P are defined as:

$\begin{array}{cc}\begin{array}{c}{\overline{F}}_G=1-F_G,\\ {\overline{F}}_P=1-F_P.\end{array}& \left(54\right)\end{array}$

For characterizing the quality of the various Hamiltonians, it is useful to define the gate angle error:

$\begin{array}{cc}\mathrm{\Delta \chi }=\chi -\frac{\pi }{4},& \left(55\right)\end{array}$

where χ is the actual gate angle computed by running the XX gate simulator (see Sections V and VII). A positive gate angle error Δχ corresponds to an over-rotated gate angle χ, while a negative Δχ corresponds to an under-rotated gate angle χ. The error operator associated with the gate angle error Δχ is Ê_χ=Δχ{circumflex over (σ)}⁽¹⁾{circumflex over (σ)}⁽²⁾. Then, according to (35), for |ψ₀=|00|ph, for instance, where |00 is the computational state and |ph is any normalized phonon state, the state infidelity F_S^(χ) caused by Ê_χ is

$\begin{array}{cc}{{\stackrel{\_}{F}}_S}^{\left(\chi \right)}={\left(\Delta \chi \right)}^2.& \left(56\right)\end{array}$

VII. Numerical Results

As described in Section IV, the inventors have computed AMFM pulse functions g(t) for gate times ranging from τ=100 μs to τ=600 μs in steps of 100 μs, and used these pulse functions in the XX gate simulator (see Section V). As described in Section V, using the full Hamiltonian Ĥ_MS defined in (5), the full state function |ψ(τ) that includes both the computational space and the phonon space were computed. From |ψ(τ), as described in Section VI, the state infidelity F_S, the gate infidelity F_G, the process infidelity F_P, and the gate angle error Δχ were computed. The results are shown in Table III.

TABLE III
Phonon schemes and infidelities for XX gate between ion pair (2,5) as
a function of gate time τ in μs. All infidelities are in units of 10−4.
τ	100	200	300	400	500	600
phonon	2255111	4522111	2621111	2621111	2621111	2621111
scheme
F S	3.1	2.6	4.0	1.9	2.5	1.7
F G	2.4	2.1	3.1	1.5	2.0	1.3
F P	3.0	2.6	3.8	1.8	2.5	1.6
Δχ	−0.012	−0.011	−0.011	−0.011	−0.011	−0.011
F S c	1.6	1.5	2.9	0.78	1.4	0.61
F G c	1.3	1.2	2.2	0.56	1.1	0.47
F P c	1.6	1.5	2.7	0.70	1.4	0.58
F S Φ, c	0.88	0.14	0.62	0.41	0.34	0.45
F G Φ, c	0.67	0.10	0.36	0.28	0.26	0.34
F P Φ, c	0.84	0.12	0.45	0.34	0.32	0.42

Accepting Ĥ_MS as a good approximation of Ĥ_R, the most important result from Table III is that a gate infidelity F_G<10⁻⁴ cannot be achieved if the control pulses g(t) are constructed on the basis of the standard Hamiltonian Ĥ_S. Nevertheless, while not quite meeting the goal of ≤10⁻⁴, the infidelities obtained are very close to this goal. It can also be seen that the three different infidelity measures, i.e., the state infidelity F_S, the average gate infidelity F_G, and the process infidelity F_P yield similar results and any one of them may be used as a proxy for assessing the infidelity of the two-qubit XX gate.

Next, the control pulse g(t) will be fixed using the 300 μs pulse from Table III and the quality of the various approximations Ĥ_M^(Nc,Ns) is determined with respect to the full Hamiltonian Ĥ_MS by computing the process infidelity {circumflex over (F)}_P for some of the Hamiltonians Ĥ_M^(Nc,Ns). The result is shown in Table IV. It can be seen that, expectedly, F_P<<10⁻⁴ for Ĥ_M^−2,1)=Ĥ_S, since in this case the pulse is generated on the basis of the standard Hamiltonian Ĥ_S and the gate simulator is controlled by the standard Hamiltonian Ĥ_S as well. So, perfect fidelity is ideally expected. The difference from zero in Table IV in this case is not due to the accuracy of the numerical integrator, which, as stated in Section V is of the order of 10⁻⁷, but is due to the phonon scheme. Including more phonon states in the basis increases the accuracy of the simulations and drives the infidelity in the case of Ĥ_M^{Nc=−2,Ns=1)}=Ĥ_S closer to zero. The table entry for (N_c=−2, N_s=1) also provides us with an estimate of the accuracy of the infidelity entries in Tables III and IV. As indicated by the (N_c=−2, N_s=1) entry in Table IV, the chosen phonon schemes guarantee an accuracy of the computed infidelities of approximately 3×10⁻⁵.

TABLE IV
Process infidelities FP in units of 10−4 for phonon
scheme 2621111 and ion pair (2, 5) produced by model
Hamiltonians HM(Nc, Ns). A dash in the table means that
the corresponding quantity was not computed.
	Ns = 1	Ns = 3	Ns = 5	Ns = ∞
	Nc = −2	0.31	1.6	1.6	1.6
	Nc = 0	2.5	3.8	3.8	—
	Nc = 2	2.5	3.8	3.8	—
	Nc = ∞	—	—	—	3.8

Looking at the infidelity results in Table IV for the Hamiltonians expanded to zeroth and second orders in the cos-term in (5), it can be seen that neglecting the zeroth-order term in (5) (the carrier term) is not justified. But it can also be seen that, apparently, expansion to second order of the cos-term of (5) is not necessary. This is shown analytically in more detail in Sections IX. B. and IX. C.

Turning now to the expansions of the sin-term of (5), Table IV shows that truncating this expansion at the first order of the sin function, i.e., linearizing the sin-term in the Lamb-Dicke parameters, is not accurate enough. As a consequence, the expansion of the sin term in (5) has to be carried to the third order in the Lamb-Dicke parameters, but expansion to the 5th order is not necessary.

As an overall result of the performance tests for different (N_c, N_s) model Hamiltonians, Ĥ_M^(0,3) is a good enough approximation of Ĥ_MS on the ≤10⁻⁴ fidelity level. Conversely, the important result is that pulse construction on the basis of the standard Hamiltonian Ĥ_S=Ĥ_M^(−2,1) does not yield infidelities <10⁻⁴. Thus, to improve control-pulse construction, at a minimum, the carrier term [zeroth-order cos term in (5)] and the third-order sin term in (5) needed to be included. In Sections IX. B. and IX. C., the effects of these two additional Hamiltonian terms are further explored. It can also be seen analytically that inclusion of the second- and fourth-order terms of the Hamiltonian (5) are not needed.

VIII. Calibration

Table III shows that the gate angle error Δχ is quite large and may make a significant contribution to the infidelity. However, by slightly adjusting the amplitude of the control pulse g(t), the gate angle error Δχ can be completely eliminated and thus any contribution to the infidelity that otherwise would be due to the gate angle error Δχ can be eliminated. Adjusting the pulse function g(t) to result in Δχ=0 is called calibration. This procedure, extensively used in the laboratory, is an attractive way of reducing the infidelity, since the phase-space closure condition (14) linear in the amplitude of the control pulse g(t), is invariant under a change in the control-pulse amplitude. So, despite calibration of the pulse, phase-space closure is always guaranteed exactly, independent of the pulse amplitude.

The error operator for the gate angle error Δχ is Ê=Δχ{circumflex over (σ)}_x⁽¹⁾{circumflex over (σ)}_x⁽²⁾, which trivially commutes with {circumflex over (σ)}_x⁽¹⁾{circumflex over (σ)}_x⁽²⁾. Therefore, for all entries in Table III, the contribution of the gate angle error Δχ to the state infidelity F_S may be estimated according to (35) and (36) as:

$\begin{array}{cc}{\overline{F_S}}^{\left(\Delta \chi \right)}={\left(\Delta \chi \right)}^2\left[1-{\langle {\psi }_0❘{\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}❘{\psi }_0\rangle }^2\right].& \left(57\right)\end{array}$

Based on this result, and for all pulse lengths listed in Table III, F_S^(Δχ)˜1.2×10⁻⁴ is obtained, which makes a significant contribution to the infidelities listed in Table III. However, as mentioned above, this infidelity may be eliminated by calibrating the control pulse. Denoting by χ_target the desired gate angle and by χ_g=χ_target+Δχ the gate angle actually obtained with the control pulse g(t), the calibration factor c, i.e., the (real) factor by which the amplitude of g(t) has to be multiplied to eliminate the gate angle error Δχ is obtained explicitly, with (15), as:

$\begin{array}{cc}c={\left({\chi }_{\mathrm{target}}/{\chi }_g\right)}^{1/2}={\left(\frac{{\chi }_{target}}{{\chi }_{tar\mathrm{get}}+\mathrm{\Delta \chi }}\right)}^{1/2}.& \left(58\right)\end{array}$

The entries for the state infidelity F_S^c, the gate infidelity F_G^c, and the process infidelity F_P^c in Table III represent the results for the respective infidelities computed with the same control pulses used for the corresponding entries F_S, F_G, and F_P, but multiplied (calibrated) with the calibration factor c, computed according to (58) with the corresponding gate angle error Δχ values from Table III and χ_target=π/4. As expected, the result is a significant reduction of the infidelity. Since the starting state |ψ₀ for the results in Table III is |ψ₀=|00|0_ph, a reduction by Δχ² can be expected. This agrees well with the results in Table III. Since they can no longer be due to the gate angle error Δχ, the rest of the infidelities in Table III have to come from other sources. To look for these sources, the XX-gate propagator can be computed via a Magnus expansion as outlined in the following section.

IX. Magnus Expansion

Given the time-dependent Schrödinger equation:

$\begin{array}{cc}i\frac{\partial }{\partial t}❘\psi \left(t\right)\rangle =\widehat{H}\left(t\right)❘\psi \left(t\right)\rangle ,& \left(59\right)\end{array}$

the time evolution operator Û(τ) of (59) over the time interval τ may be constructed systematically and analytically, using a Magnus expansion, i.e., up to 3rd order in the Hamiltonian:

$\begin{array}{cc}\widehat{U}\left(\tau \right)=\exp \left[i{\widehat{W}}_1\left(\tau \right)+i{\widehat{W}}_2\left(\tau \right)+i{\widehat{W}}_3\left(\tau \right)+\dots \right],& \left(60\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}{\widehat{W}}_1\left(\tau \right)=\frac{1}{i}{\int }_0^{\tau }\left(-\right)\widehat{H}\left(t_1\right)d{t}_1,& \left(61\right)\end{array}$ ${\widehat{W}}_2\left(\tau \right)=\frac{1}{2i}{\int }_0^{\tau }d{t}_1{\int }_0^{t_1}d{t}_2\left[\left(-\right)\widehat{H}\left(t_1\right),\left(-\right)\widehat{H}\left(t_2\right)\right],$ ${\widehat{W}}_3\left(\tau \right)=\frac{1}{6i}{\int }_0^{\tau }d{t}_1{\int }_0^{t_1}d{t}_2{\int }_0^{t_2}d{t}_3\left\{\left[\left(-\right)\widehat{H}\left(t_1\right),\left[\left(-\right)\widehat{H}\left(t_2\right),\left(-\right)\widehat{H}\left(t_3\right)\right]\right]+\left[\left(-\right)\widehat{H}\left(t_3\right),\left[\left(-\right)\widehat{H}\left(t_2\right),\left(-\right)\widehat{H}\left(t_1\right)\right]\right]\right\}$

are Hermitian operators.

In this section, a consistent Magnus expansion up to 4th order in the Lamb-Dicke parameters η_p^j is shown. Therefore, expanding the cosine-term in (5) up to fourth order in {circumflex over (V)}² and the sine-term in (5) up to third order in {circumflex over (V)}_j, the Hamiltonian:

$\begin{array}{cc}\widehat{H}\left(t\right)=g\left(t\right){\sum }_{\alpha =0}^4{\hat{h}}_{\alpha }\left(t\right)& \left(62\right)\end{array}$

is used in the Magnus expansion (61), where

$\begin{array}{cc}\begin{array}{c}{\hat{h}}_0\left(t\right)={\sum }_j{\sigma }_y^{\left(j\right)},\\ {\hat{h}}_1\left(t\right)={\sum }_j\widehat{V_j}\left(t\right){\hat{\sigma }}_x^{\left(j\right)},\\ {\hat{h}}_2\left(t\right)=-\frac{1}{2}{\sum }_j{\hat{V}}_j^2\left(t\right){\hat{\sigma }}_y^{\left(j\right)},\\ {\hat{h}}_3\left(t\right)=-\frac{1}{6}{\sum }_j{\hat{V}}_j^3\left(t\right){\hat{\sigma }}_x^{\left(j\right)},\\ {\hat{h}}_4\left(t\right)=\frac{1}{24}{\sum }_j{\hat{V}}_j^4\left(t\right){\hat{\sigma }}_y^{\left(j\right)}.\end{array}& \left(63\right)\end{array}$

It should be noted that ĥ₀(t) is actually time independent, but the time argument is formally kept for notational convenience. Based on the results listed in Table IV, it is known that the zeroth order of cos[{circumflex over (V)}_j(t)] contributes substantially to the infidelity, while the 2nd and higher orders of cos[{circumflex over (V)}_j(t)] do not, and that the first and third orders of sin[{circumflex over (V)}_j(t)] contribute substantially, while the 5th and higher orders do not. Thus, the Hamiltonian (62) covers all these important cases. According to Table IV, it may not be necessary to include the second- and fourth-order expansion terms of the cosine function in (5). However, for a consistent expansion up to fourth order in the Lamb-Dicke parameters, and also to show that these terms are negligible, they are included in the expansion (62) of (5).

In Section IX. A., Ŵ₁(τ) is computed on the basis of (62). In Sections IX.B. and IX.C., Ŵ₂(τ) and Ŵ₃(τ), respectively are constructed. Many of the resulting error terms are found to be negligible. Further, the most significant terms that make significant contributions to the infidelity are identified.

With the operators defined in (63), the following Hermitian operators are defined:

$\begin{array}{cc}{\hat{T}}_{\alpha }=-{\int }_0^{\tau }{\hat{h}}_{\alpha }\left(t\right)g\left(t\right)dt,\alpha =\mathrm{0,1,2,3,4,}& \left(64\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{\mathrm{\alpha \beta }}=-\frac{1}{2i}{\int }_0^{\tau }d{t}_1{\int }_0^{t_1}d{t}_{2}g\left(t_1\right)g\left(t_2\right)\left[{\hat{h}}_{\alpha }\left(t_1\right),{\hat{h}}_{\beta }\left(t_2\right)\right],\alpha ,\beta =\mathrm{0,1,2,3,4,}& \left(65\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}{\hat{T}}_{a\beta \gamma }=\left(\frac{1}{6}\right){\int }_0^{\tau }d{t}_1{\int }_0^{t_1}d{t}_2{\int }_0^{t_2}d{t}_{3}g\left(t_1\right)g\left(t_2\right)g\left(t_3\right)\left\{\left[{\hat{h}}_{\alpha }\left(t_1\right),\left[{\hat{h}}_{\beta }\left(t_2\right)\right],{\hat{h}}_{\gamma }\left(t_3\right)\right]+\left[{\hat{h}}_a\left(t_3\right),\left[{\hat{h}}_{\beta }\left(t_2\right)\right],{\hat{h}}_{\gamma }\left(t_1\right)\right]\right\},\alpha ,\beta ,\gamma =\mathrm{0,1,2,3,4}.& \left(66\right)\end{array}$

In the following sections, these operators are used (i) to compute Ŵ₁, Ŵ₂, and Ŵ₃ and (ii) to compute infidelity contributions to the gate evolution operator Û(τ) defined in (60).

IX. A. First Order

With the definition in (64) Ŵ₁(τ) becomes:

$\begin{array}{cc}{\widehat{W}}_1\left(\tau \right)={\sum }_{\alpha =0}^4{\hat{T}}_{\alpha }\left(\tau \right),& \left(67\right)\end{array}$

where the carrier condition (13) together with

$\begin{array}{cc}{\int }_0^{\tau }g\left(t\right)e^{i{\omega }_{p}t}dt=0\Rightarrow {\int }_0^{\tau }g\left(t\right){\hat{V}}_j\left(t\right)dt=0,& \left(68\right)\end{array}$

is used to arrive at

$\begin{array}{cc}{\hat{T}}_0\left(\tau \right)=0,& \left(69\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_1\left(\tau \right)=0,& \left(70\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_2\left(\tau \right)=\left(\frac{1}{2}\right){\sum }_{pqj}{\eta }_p^j{\eta }_q^j{\hat{\sigma }}_y^{\left(j\right)}\left[Q\left({\omega }_p+{\omega }_q\right){\hat{a}}_p^{†}{\hat{a}}_q^{†}+h.c.\right]+{\sum }_{pqj,p\ne q}{\eta }_p^j{\eta }_q^j{\hat{\sigma }}_y^{\left(j\right)}Q\left({\omega }_p-{\omega }_q\right){\hat{a}}_p^{†}{\hat{a}}_q,& \left(71\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_3\left(\tau \right)=\frac{1}{6}{\sum }_{j=1,2}{\hat{\sigma }}_x^{\left(j\right)}{\int }_0^{\tau }g\left(t\right){\widehat{V_j}}^3\left(t\right)\mathrm{dt}=\frac{1}{6}{\sum }_{p,q,r,j}{\eta }_p^j{\eta }_q^j{\eta }_r^j{\hat{\sigma }}_x^{\left(j\right)}\left\{Q\left({\omega }_p+{\omega }_q+{\omega }_r\right){\hat{a}}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_r^{†}+Q\left({\omega }_p+{\omega }_q-{\omega }_r\right){\hat{a}}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_r+Q\left({\omega }_p-{\omega }_q+{\omega }_r\right){\hat{a}}_p^{†}{\hat{a}}_q{\hat{a}}_r^{†}+Q\left(-{\omega }_p+{\omega }_q+{\omega }_r\right){\hat{a}}_p{\hat{a}}_q^{†}{\hat{a}}_r^{†}+h.c.\right\},& \left(72\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}{\hat{T}}_4\left(\tau \right)=\left(-\frac{1}{24}\right){\sum }_j{\hat{\sigma }}_y^{\left(j\right)}{\int }_0^{\tau }g\left(t\right){\widehat{V_j}}^4\left(t\right)\mathrm{dt}=\left(-\frac{1}{24}\right){\sum }_j{\hat{\sigma }}_y^{\left(j\right)}{\int }_0^{\tau }g\left(t\right){\sum }_{pqrs}{\eta }_p^j{\eta }_q^j{\eta }_r^j{\eta }_s^j\left\{{\hat{a}}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_r^{†}{\hat{a}}_s^{†}{e}^{i\left({\omega }_p+{\omega }_q+{\omega }_r+{\omega }_s\right)t}+{\hat{a}}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_r^{†}{\hat{a}}_s{e}^{i\left({\omega }_p+{\omega }_q+{\omega }_r-{\omega }_s\right)t}\dots \right\},& \left(73\right)\end{array}$

where Q(w), defined in (22), is used, in particular, with the carrier condition (13), Q(0)=0.

While {circumflex over (T)}₀(τ) and {circumflex over (T)}₁(τ) vanish, {circumflex over (T)}₂(τ), {circumflex over (T)}₃(τ), and {circumflex over (T)}₄(τ) need to be investigated further since they may produce undesirable contributions to the infidelity of the XX gate.

According to (71), the size of {circumflex over (T)}₂(τ) is controlled by

$\begin{array}{cc}{\gamma }_2=\underset{\mathrm{pqj},\sigma =\pm 1}{\max }❘{\eta }_p^j{\eta }_q^{j}Q\left({\omega }_p+\sigma {\omega }_q\right)❘.& \left(74\right)\end{array}$

Numerically, for the 300 μs test pulse,

$\begin{array}{cc}{\gamma }_2=3.5\times 10^{-6},& \left(75\right)\end{array}$

was obtained. Since γ₂ is significantly smaller than 10⁻⁴, {circumflex over (T)}₂(τ) may be neglected. According to (72), the size of {circumflex over (T)}₃(τ) is controlled by

$\begin{array}{cc}{\gamma }_3^{(+)}=\left(\frac{1}{6}\right)\underset{\mathrm{pqrj}}{\max }❘{\eta }_p^j{\eta }_q^j{\eta }_r^{j}Q\left({\omega }_p+{\omega }_q+{\omega }_r\right)❘& \left(76\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}\begin{array}{c}{\gamma }_3^{(-)}=\left(\frac{1}{6}\right)\underset{\mathrm{pqrj}}{\max }❘{\eta }_p^j{\eta }_q^j{\eta }_r^{j}Q\left({\omega }_p+{\omega }_q-{\omega }_r\right)❘\\ =\left(\frac{1}{6}\right)\underset{\mathrm{pqrj}}{\max }❘{\eta }_p^j{\eta }_q^j{\eta }_r^{j}Q\left({\omega }_p-{\omega }_q+{\omega }_r\right)❘\\ =\left(\frac{1}{6}\right)\underset{\mathrm{pqrj}}{\max }❘{\eta }_p^j{\eta }_q^j{\eta }_r^{j}Q\left(-{\omega }_p+{\omega }_q+{\omega }_r\right)❘.\end{array}& \left(77\right)\end{array}$

For the 300 μs test pulse,

$\begin{array}{cc}{\gamma }_3^{(+)}=5.8\times 10^{-8},& \left(78\right)\end{array}$ ${\gamma }_3^{(-)}=3.1\times 10^{-4},$

was obtained. Since {circumflex over (σ)}_x^(j) commutes with {circumflex over (σ)}x⁽¹⁾{circumflex over (σ)}_x⁽²⁾, (35) may be used to estimate the infidelity F_S^(T3)) caused by {circumflex over (T)}₃(τ). Since, according to (78), both γ₃⁽⁺⁾ and γ₃⁽⁻⁾ are smaller than 10⁻³, and since F_S^(T3), according to (35), involves the square of {circumflex over (T)}₃(τ), the infidelity caused by {circumflex over (T)}₃(τ) is negligible on the level of 10⁻⁴.According to (73), the size of {circumflex over (T)}₄(τ) is controlled by

$\begin{array}{cc}{\gamma }_4=\underset{{\sigma }_p{\sigma }_q{\sigma }_r{\sigma }_s=\pm 1}{\underset{j,\mathrm{pqrs}}{\max }}❘{\eta }_p^j{\eta }_q^j{\eta }_r^j{\eta }_s^{j}Z\left({\sigma }_p{\omega }_p,{\sigma }_q{\omega }_q+{\sigma }_r{\omega }_r+{\sigma }_s{\omega }_s\right)❘.& \left(79\right)\end{array}$

For the 300 μs test pulse,

$\begin{array}{cc}{\gamma }_4=1.1\times 10^{-7},& \left(80\right)\end{array}$

was obtained. Therefore, {circumflex over (T)}₄(τ) can be neglected.

As a result of this section it has been found that, on the 10⁻⁴ level, the first-order terms in the Magnus expansion do not contribute to the infidelity. However, it should be noted that this is only true if the pulse function g(t) satisfies the carrier condition (13). Thus, the carrier condition (13) is an important condition that needs to be required for high-fidelity XX gates.

IX. B. Second Order

Now, the second-order terms (65) are evaluated in the Magnus expansion, i.e., Ŵ₂(τ) is computed analytically for Ĥ(t) defined in (62) up to fourth order in the Lamb-Dicke parameters η. With the definitions (65) of the operators {circumflex over (T)}_αβ(τ), the definition (15) of the gate angle, making use of the carrier condition (13), and the functions defined in (16)-(28),

$\begin{array}{cc}{\hat{W}}_2\left(\tau \right)={\sum }_{\alpha ,\beta =0}^4{\hat{T}}_{\mathrm{\alpha \beta }}\left(\tau \right),& \left(81\right)\end{array}$

was obtained, where the operators {circumflex over (T)}_αβ, listed only up to fourth order in η, are

$\begin{array}{cc}{\hat{T}}_{00}=0,& \left(82\right)\end{array}$ $\begin{array}{cc}\begin{array}{c}{\hat{T}}_{01}=-{\sum }_{j=1,2}{\hat{\sigma }}_z^{\left(j\right)}{\int }_0^{\tau }g\left(t\right)G\left(t\right)\widehat{V_j}\left(t\right)\mathrm{dt}\\ =-{\sum }_{\mathrm{jp}}{\eta }_p^j{\hat{\sigma }}_z^{\left(j\right)}\left[f\left({\omega }_p\right){\hat{a}}_p^{†}+f^{*}\left({\omega }_p\right){\hat{a}}_p\right],\end{array}& \left(83\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{02}=0,& \left(84\right)\end{array}$ $\begin{array}{cc}\begin{array}{c}{\hat{T}}_{03}=\left(\frac{1}{6}\right)\sum _{j=1,2}{\hat{\sigma }}_z^{\left(j\right)}{\int }_0^{\tau }g\left(t\right)G\left(t\right){\widehat{V_j}}^3\left(t\right)\mathrm{dt}\\ =\left(\frac{1}{6}\right)\sum _{\mathrm{pqrj}}{\hat{\sigma }}_z^{\left(j\right)}{\eta }_p^j{\eta }_q^j{\eta }_r^j\{f\left({\omega }_p+{\omega }_q+{\omega }_r\right){\hat{a}}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_r^{†}+\\ f\left({\omega }_p+{\omega }_q-{\omega }_r\right){\hat{a}}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_r^{†}\\ +f\left({\omega }_p-{\omega }_q+{\omega }_r\right){\hat{a}}_p^{†}{\hat{a}}_q{\hat{a}}_r^{†}+f\left(-{\omega }_p+{\omega }_q+{\omega }_r\right){\hat{a}}_p{\hat{a}}_q^{†}{\hat{a}}_r^{†}+h.c.\},\end{array}& \left(85\right)\end{array}$ ${\hat{T}}_{04}=0,$ $\begin{array}{cc}{\hat{T}}_{10}={\hat{T}}_{01},& \left(87\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{11}=\tilde{\chi }+\chi {\hat{\sigma }}_x^{\left(1\right)}{\hat{\sigma }}_x^{\left(2\right)},& \left(88\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{12}=\left(\frac{1}{4}\right){\sum }_j{\hat{\sigma }}_z^{\left(j\right)}{\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)\left\{{\hat{V}}_j\left(t_1\right),{\widehat{V_j}}^2\left(t_2\right)\right\}-{\sum }_{\mathrm{pq},j\ne k}{\eta }_p^j{\eta }_p^k{\eta }_q^k{\hat{\sigma }}_x^{\left(j\right)}{\hat{\sigma }}_y^{\left(k\right)}\left[S_p\left({\omega }_q\right){\hat{a}}_q^{†}+S_p^{*}\left({\omega }_q\right){\hat{a}}_q\right],& \left(89\right)\end{array}$ $\begin{array}{cc}\begin{array}{c}{\hat{T}}_{13}=\left(-\frac{1}{2}\right){\sum }_{\mathrm{jkp}}{\hat{\sigma }}_x^{\left(j\right)}{\hat{\sigma }}_x^{\left(k\right)}{\eta }_p^j{\eta }_p^k{\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}\\ {\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)\sin \left[{\omega }_p\left(t_1-t_2\right)\right]{\hat{V}}_j^2\left(t_1\right)\\ =\left(-\frac{1}{2}\right){\sum }_{\mathrm{jkpqr}}{\hat{\sigma }}_x^{0)}{\hat{\sigma }}_x^{\left(k\right)}{\eta }_p^j{\eta }_p^k{\eta }_q^j{\eta }_r^j[S_p\left({\omega }_q+{\omega }_r\right){\hat{a}}_q^{†}{\hat{a}}_r^{†}\\ +S_p\left({\omega }_q-{\omega }_r\right){\hat{a}}_q^{†}{\hat{a}}_r+S_p^{*}\left({\omega }_q-{\omega }_r\right){\hat{a}}_q{\hat{a}}_r^{†}+S_p^{*}\left({\omega }_q+{\omega }_r\right){\hat{a}}_q{\hat{a}}_r],\end{array}& \left(90\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{20}=0,& \left(91\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{21}={\hat{T}}_{12},& \left(92\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{22}=\left(\frac{1}{2}\right){\sum }_{\mathrm{jkp}}{\eta }_p^j{\eta }_p^k{\hat{\sigma }}_y^{\left(j\right)}{\hat{\sigma }}_y^{\left(k\right)}{\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)\sin \left[{\omega }_p\left(t_1-t_2\right)\right]\left\{{\hat{V}}_j\left(t_1\right),{\hat{V}}_k\left(t_2\right)\right\},& \left(93\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{30}={\hat{T}}_{03},& \left(94\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{31}={\hat{T}}_{13},& \left(95\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}{\hat{T}}_{40}=0.& \left(96\right)\end{array}$

The above explicit, analytical results are computed making explicit use of the carrier condition (13).

Since {circumflex over (χ)} in (88) is only a c-number, which causes only a global phase, {circumflex over (T)}₁₁ generates the desired XX gate. All the other operators {circumflex over (T)}_αβ in (82)-(96), if not identically zero, are unwanted error operators that generate infidelities. Because of symmetries and the operators that are identically zero, the sizes of only five remaining, non-zero error operators, i.e., {circumflex over (T)}₀₁, {circumflex over (T)}₀₃, {circumflex over (T)}₁₂, {circumflex over (T)}₁₃, and {circumflex over (T)}₂₂ need to be investigated.

The size of {circumflex over (T)}₀₁ is controlled by

$\begin{array}{cc}{\gamma }_{01}=\underset{\mathrm{pj}}{\max }❘{\eta }_p^{j}f\left({\omega }_p\right)❘.& \left(97\right)\end{array}$

For the 300 μs test pulse

$\begin{array}{cc}{\gamma }_{01}=1.8\times 10^{-6},& \left(98\right)\end{array}$

was obtained. Thus, on the 10⁻⁴ level, {circumflex over (T)}₀₁ can safely be neglected.

Next, {circumflex over (T)}₀₃ is evaluated, which proceeds in analogy to the evaluation of {circumflex over (T)}₃(τ) in Section IX. A. The size of {circumflex over (T)}₀₃ is controlled by

$\begin{array}{cc}{\gamma }_{03}^{(+)}=\left(\frac{1}{6}\right)\underset{\mathrm{pqrj}}{\max }❘{\eta }_p^j{\eta }_q^j{\eta }_r^{j}f\left({\omega }_p+{\omega }_q+{\omega }_r\right)❘& \left(99\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}\begin{array}{c}{\gamma }_{03}^{(-)}=\left(\frac{1}{6}\right)\underset{\mathrm{pqrj}}{\max }❘{\eta }_p^j{\eta }_q^j{\eta }_r^{j}f\left({\omega }_p+{\omega }_q-{\omega }_r\right)❘\\ =\left(\frac{1}{6}\right)\underset{\mathrm{pqrj}}{\max }❘{\eta }_p^j{\eta }_q^j{\eta }_r^{j}f\left({\omega }_p-{\omega }_q+{\omega }_r\right)❘\\ =\left(\frac{1}{6}\right)\underset{\mathrm{pqrj}}{\max }❘{\eta }_p^j{\eta }_q^j{\eta }_r^{j}f\left(-{\omega }_p+{\omega }_q+{\omega }_r\right)❘.\end{array}& \left(100\right)\end{array}$

For the 300 μs test pulse,

$\begin{array}{cc}{\gamma }_{03}^{(+)}=7.4\times 10^{-10},& \left(101\right)\end{array}$ ${\gamma }_{03}^{(-)}=9.3\times 10^{-10},$

was obtained. Thus, {circumflex over (T)}₀₃ is negligible.According to (89), the operator {circumflex over (T)}₁₂ consists of two parts, an anti-commutator part and a part that originated from a commutator. The size of the anti-commutator part of {circumflex over (T)}₁₂ is controlled by

$\begin{array}{cc}{\gamma }_{12}^{\left(a\right)}=\underset{{\sigma }_p{\sigma }_q{\sigma }_r=\pm 1}{\underset{j,\mathrm{pqr}}{\max }}❘{\eta }_p^j{\eta }_q^j{\eta }_r^{j}Z\left({\sigma }_p{\omega }_p,{\sigma }_q{\omega }_q+{\sigma }_r{\omega }_r\right)❘.& \left(102\right)\end{array}$

The commutator part is controlled by

$\begin{array}{cc}{\gamma }_{12}^{\left(c\right)}=\underset{\mathrm{pq},j\ne k}{\max }❘{\eta }_p^j{\eta }_p^k{\eta }_q^k{S}_p\left({\omega }_q\right)❘.& \left(103\right)\end{array}$

For the 300 μs test pulse both turned out to be very small. Thus, {circumflex over (T)}₁₂ can be neglected.

The operator {circumflex over (T)}₁₃ in (90) can be split into a diagonal part

$\begin{array}{cc}{\hat{T}}_{13}^{\left(d\right)}=\left(-\frac{1}{2}\right){\sum }_{\mathrm{jpqr}}{\left({\eta }_p^j\right)}^2{\eta }_q^j{\eta }_r^j\left[S_p\left({\omega }_q+{\omega }_r\right){\hat{a}}_q^{†}{\hat{a}}_r^{†}+S_p\left({\omega }_q-{\omega }_r\right){\hat{a}}_q^{†}{\hat{a}}_r+S_p^{*}\left({\omega }_q-{\omega }_r\right){\hat{a}}_q{\hat{a}}_r^{†}+S_p^{*}\left({\omega }_q+{\omega }_r\right){\hat{a}}_q{\hat{a}}_r\right]& \left(104\right)\end{array}$

and an off-diagonal part

$\begin{array}{cc}{\hat{T}}_{13}^{\left(o\right)}=\left(-\frac{1}{2}\right){\sum }_{\mathrm{jkpqr}}{\hat{\sigma }}_x^{\left(1\right)}{\hat{\sigma }}_x^{\left(2\right)}{\eta }_p^{\left(1\right)}{\eta }_p^{\left(2\right)}\left[{\eta }_q^{\left(1\right)}{\eta }_r^{\left(1\right)}+{\eta }_q^{\left(2\right)}{\eta }_r^{\left(2\right)}\right]\left[S_p\left({\omega }_q+{\omega }_r\right){\hat{a}}_q^{†}{\hat{a}}_r^{†}+S_p\left({\omega }_q-{\omega }_r\right){\hat{a}}_q^{†}{\hat{a}}_r+S_p^{*}\left({\omega }_q-{\omega }_r\right){\hat{a}}_q{\hat{a}}_r^{†}+S_p^{*}\left({\omega }_q+{\omega }_r\right){\hat{a}}_q{\hat{a}}_r\right].& \left(105\right)\end{array}$

The size of {circumflex over (T)}₃^(d) is controlled by

$\begin{array}{cc}{\gamma }_{13}^{\left(d\right)}=\underset{\begin{array}{c}j,\mathrm{pqr}\\ \sigma =\pm 1\end{array}}{\max }❘{\left({\eta }_p^j\right)}^2{\eta }_q^j{\eta }_r^j{S}_p\left({\omega }_q+{\mathrm{\sigma \omega }}_r\right)❘.& \left(106\right)\end{array}$

For the 300 μs test pulse,

$\begin{array}{cc}{\gamma }_{13}^{\left(d\right)}=1.17\times 10^{-3},& \left(107\right)\end{array}$

were obtained. Thus, compared to the target infidelity of <<10⁻⁴, {circumflex over (T)}₁₃^(d) cannot be dismissed outright. However, taking {circumflex over (T)}₁₃^(d) as the error operator and since {circumflex over (T)}₁₃^(d) commutes with σ_x⁽¹⁾α_x⁽²⁾, the infidelity estimate (35) can be used, which involves the squares of {circumflex over (T)}₁₃^(d), i.e., the contribution of {circumflex over (T)}₁₃^(d) to the infidelity of Û(τ) is expected to be of the order of (γ₁₃^(d))²˜10⁻⁶. This means that the contribution of {circumflex over (T)}₁₃^(d) to the infidelity of Û(τ) can be neglected.

The term {circumflex over (T)}₁₃^(o) contains σ_x⁽¹⁾σ_x⁽²⁾ and thus contributes to over/under rotation of the XX gate, i.e., it contributes to Δχ (see Table III). The operator {circumflex over (T)}₁₃^(o) acts in the computational space but also produces phonon excitations. However, according to (105), the two-phonon excitation terms are proportional to S_p(ω_q+ω_r), which are nonresonant and very small. Assuming that the initial state starts out in the phonon ground state, |0>_ph, the operator {circumflex over (T)}₁₃^(o) is, to an excellent approximation,

$\begin{array}{cc}{\hat{T}}_{13}^{\left(0\right)}=\left(-\frac{1}{2}\right){\sum }_{pq}{\eta }_p^{\left(1\right)}{\eta }_p^{\left(2\right)}\left[{\left({\eta }_q^{\left(1\right)}\right)}^2+{\left({\eta }_q^{\left(2\right)}\right)}^2\right]{\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}{S}_p\left(0\right).& \left(108\right)\end{array}$

Now, with (15),

$\begin{array}{cc}{\sum }_p{\eta }_p^{\left(1\right)}{\eta }_p^{\left(2\right)}{S}_p\left(0\right)={\sum }_p{\eta }_p^{\left(1\right)}{\eta }_p^{\left(2\right)}{\int }_0^{\tau }d{t}_1{\int }_0^{t_1}d{t}_{2}g\left(t_1\right)g\left(t_2\right)\sin \left[{\omega }_p\left(t_1-t_2\right)\right]=\frac{\chi }{2}.& \left(109\right)\end{array}$

With this result,

$\begin{array}{cc}{\hat{T}}_{13}^{\left(0\right)}=\left[\left(-\frac{\chi }{4}\right){\sum }_{jp}{\left({\eta }_p^j\right)}^2\right]{\sigma }_x^{\left(1\right)}{\sigma }_x^{\left(2\right)}.& \left(110\right)\end{array}$

The prefactor of the α_x⁽¹⁾σ_x⁽²⁾ operator is a c-number. Therefore, {circumflex over (T)}₁₃^(o) produces a contribution to the gate angle χ. Since {circumflex over (T)}₃₁={circumflex over (T)}₁₃, the total contribution to the gate angle error Δχ is

$\begin{array}{cc}\Delta \chi =\left(-\frac{\chi }{2}\right){\sum }_{jp}{\left({\eta }_p^j\right)}^2.& \left(111\right)\end{array}$

Numerically, for the N=7 case (see Table II),

$\begin{array}{cc}{\sum }_{jp}{\left({\eta }_p^j\right)}^2=2.45\times 10^{-2},& \left(112\right)\end{array}$

was obtained. For χ=π/4, this results in

$\begin{array}{cc}\Delta \chi =-9.62\times 10^{-3}.& \left(113\right)\end{array}$

According to the numerical simulations (see Table III), Δχ≈−0.011 is obtained. Thus, the analytical calculations predict the correct sign (under-rotation) of χ. In addition, the magnitude of the relative error of the analytical prediction is (0.011−9.62×10⁻³)/0.011≈0.13, i.e., the analytical prediction is only of the order of 10% off.

According to (111), Δχ depends only on the Lamb-Dicke parameters, and not on the gate duration τ. This is reflected in Table III and explains why Δχ in Table III is approximately constant, independent of τ.

Summarizing the results obtained in this Section, it is found that the second-order Magnus-expansion operators yield only two substantial contributions, i.e., the operator {circumflex over (T)}₁₁, which generates the desired XX gate and the operator 2{circumflex over (T)}₁₃ which explains the under-rotation of the gate angle, its approximate independence of the gate time τ, and, approximately, its size, as listed in Table Ill.

At the end of Section VIII, it is noted that Δχ explains a significant contribution to the infidelity of Û(τ), but cannot explain the entire infidelity contribution. Thus, the additional sources of infidelity in the various orders of the Magnus expansion were explored. In this section the Inventors have found that the second order explains only the origin of Δχ, but does not reveal any additional significant sources of infidelity. Since the second order of the Magnus expansion did not reveal these sources, the third order of the Magnus expansion is now investigated, which, indeed, reveals the remaining significant sources of infidelity. In addition, it will be found that this source of infidelity can be eliminated by the addition of a single linear equation to the control-pulse construction protocol.

IX. C Third Order

In this section, Ŵ₃(τ) is computed. With the definitions (66) stated at the end of the introduction to Section IX,

$\begin{array}{cc}{\widehat{W}}_3\left(\tau \right)={\sum }_{jki}{\hat{T}}_{\mathrm{jkl}}& \left(114\right)\end{array}$

was obtained. All operators {circumflex over (T)}_jkl up to fourth order in η are listed in Section XIII. There is only a single operator of zeroth order in η, i.e., {circumflex over (T)}₀₀₀˜η⁰. It is trivially zero and does not contribute to the infidelity.

There are three operators of first order in η, i.e., {circumflex over (T)}₀₀₁, {circumflex over (T)}₀₁₀, {circumflex over (T)}₁₀₀, where [see Section XIII]{circumflex over (T)}₀₁₀=2{circumflex over (T)}₀₀₁, and {circumflex over (T)}₁₀₀=0. This means that only {circumflex over (T)}₀₀₁'s contribution to the infidelity has to be investigated. With (25), the infidelity caused by {circumflex over (T)}₀₀₁ is controlled by

$\begin{array}{cc}{\gamma }_{001}=\underset{\mathrm{jp}}{\max }❘{\eta }_p^j{J}_p❘.& \left(115\right)\end{array}$

For the 300 μs control pulse,

$\begin{array}{cc}{\gamma }_{001}=1.2\times 10^{-6},& \left(116\right)\end{array}$

was obtained. Therefore, {circumflex over (T)}₀₀₁, and with it {circumflex over (T)}₀₁₀, are negligible. Since {circumflex over (T)}₁₀₀=0, all operators ˜η¹ can be neglected.

There are six operators of second order in η, i.e., {circumflex over (T)}₀₀₂, {circumflex over (T)}₀₁₁, {circumflex over (T)}₀₂₀, {circumflex over (T)}₁₀₁, {circumflex over (T)}₁₁₀, and {circumflex over (T)}₂₀₀. Of these, according to Section XIII, only {circumflex over (T)}₀₁₁, {circumflex over (T)}₁₀₁, and {circumflex over (T)}₁₁₀ are nonzero, and of those, only {circumflex over (T)}₀₁₁ and {circumflex over (T)}₁₁₀ are non-negligible. While {circumflex over (T)}₀₁₁ acts only on the computational space, {circumflex over (T)}₁₁₀ also has a part that produces phonon excitations. This part, however, is negligibly small. As an overall result of the operators ˜η², it has been found that the effective error operator corresponding to the sum of the leading parts of {circumflex over (T)}₀₁₁ and {circumflex over (T)}₁₁₀ is:

$\begin{array}{cc}{\hat{E}}^{{\sigma }_x{\sigma }_z}=4\Phi \left[\eta ,g\right]\left[{\hat{\sigma }}_x^{\left(1\right)}{\hat{\sigma }}_z^{\left(2\right)}+{\hat{\sigma }}_x^{\left(2\right)}{\hat{\sigma }}_z^{\left(1\right)}\right].& \left(117\right)\end{array}$

The inventors have checked that all operators of order three and four in η are negligible. Thus, (117) is the only significant error operator that results from the third-order Magnus expansion. It should be noted that this operator acts only in the computational space, not in the phonon space. With (50) and λ=4Φ[η, g], the infidelity contribution of (117) for |ψ₀=|00|0_ph.

$\begin{array}{cc}{\overline{F_S}}^{\left({\sigma }_x{\sigma }_z\right)}={\left(\frac{\lambda \sin \left(\chi \right)}{\chi }\right)}^2\langle {\psi }_0❘\left[{\hat{\sigma }}_x^{\left(1\right)}{\hat{\sigma }}_z^{\left(2\right)}+{\hat{\sigma }}_x^{\left(2\right)}{\hat{\sigma }}_z^{\left(1\right)}\right]❘{\psi }_0\rangle ={\left(\frac{4\lambda }{\pi }\right)}^2,& \left(118\right)\end{array}$

was obtained. For the 300 μs test pulse,

$\begin{array}{cc}\lambda =1.17\times 10^{-2},& \left(119\right)\end{array}$

was obtained. Thus, with (118),

$\begin{array}{cc}{\overline{F_S}}^{\left({\sigma }_x{\sigma }_z\right)}=2.2\times 10^{-4},& \left(120\right)\end{array}$

was obtained. Together with (113), an estimate for the total infidelity of the 300 μs test pulse is obtained according to:

$\begin{array}{cc}{\stackrel{\_}{F}}_S={\left(\mathrm{\Delta \chi }\right)}^2+{\stackrel{\_}{F}}_S^{\left({\sigma }_x{\sigma }_z\right)}={\left(-9.62\times {10}^{-3}\right)}^2+2.2\times {10}^{-4}=3.1\times {10}^{-4}.& \left(121\right)\end{array}$

According to Table III, this accounts for about 80% of the infidelity of the 300 μs test pulse. Since, via calibration, Δχ can always be set to zero, generating control pulses that zero out the Φ functional may go a long way to reduce the infidelity. How to generate such control pulses, and that this recipe actually works to suppress the infidelity below 10⁻⁴, is shown in the following section.

X. Improved Control-Pulse Construction

In Section VIII, it has been argued and proved numerically, that a significant portion of the infidelity can be removed by calibration of the control pulse. Here, it is shown that by eliminating the infidelity due to the Φ-functional (see Section IX. C.), the infidelity can be further suppressed below the level of 10⁻⁴. To eliminate Φ, the auxiliary condition:

$\begin{array}{cc}{\sum }_n\frac{B_n}{n}=0& \left(122\right)\end{array}$

is added to the AMFM pulse-solver code and obtain new pulses {tilde over (g)}(t) that zero out Φ[η, {tilde over (g)}], defined in (27). That the AMFM code with the condition (122) added produces control pulses {tilde over (g)}(t) with Φ[η, {tilde over (g)}]=0 was confirmed explicitly. As described in Section VIII, {tilde over (g)}(t) may be renormalized (calibrated) such that {tilde over (g)}(t) not only produces Φ[η, {tilde over (g)}]=0, but simultaneously produces Δχ=0. Running the gate-simulator code (see Section V) with the calibrated pulses {tilde over (g)}(t), the infidelities F_S^Φ,c, F_G^Φ,c, and F_Φ^Φ,c were found as shown in Table III. It can be sees that in all cases the calibrated pulses {tilde over (g)}(t) produce infidelities below 10⁻⁴. It should be noted that the calibrated pulses {tilde over (g)}(t) require only insignificantly larger power compared with the original pulses g(t). This is as expected, since only one additional condition, i.e., the auxiliary condition (122), was added to the original set of linear phase-space closure conditions (14). Thus, the computation of the pulse function g(t) involves a linear process of solving a set of linear equations with respect to the Fourier coefficients B_n. It should be noted that, in conjunction with calibration, the construction of the improved AMFM pulses {tilde over (g)}(t) is still a linear process. Nonlinear optimizer codes are not required.

XI. Scaling

So far the 7-ion case has been discussed, for which a complete analysis machinery is in place consisting of pulse construction, analytical formulas for error estimates, and a gate simulator that includes all the relevant phonon states (see Section V). The analytical results do not depend on the number of ions N in the chain, i.e., the analytical results are valid for any number of qubits. In particular, as soon as a control pulse is generated (it does not matter whether this is an AM pulse, FM pulse, AMFM pulse, or any other type of pulse), this pulse can immediately be inserted into the analytical formulas, which then may be used to obtain infidelity estimates for this particular N-ion control pulse. This method is illustrated for N=36 by computing the infidelity (118), i.e., the leading source of infidelity, for N=36 uncalibrated control pulses for 2-qubit XX gates between all possible gate combinations (i₀, j₀), j₀=1, . . . , i₀−1, i₀=2, . . . ,36. This results in 630 gate combinations. The infidelities obtained are displayed in FIG. 8 in the form of a bar graph, where the height of a bar shows the frequency of occurrence of infidelities within the width of the bar. The infidelities in FIG. 8 are in units of 10⁻⁴=1 pptt, where the unit “pptt” denotes one part per ten thousand. It can be seen that the infidelities generated by N=36 control pulses constructed on the basis of the Standard Hamiltonian Ĥ_S are significant, and in many cases they are much larger than 10 pptt. This particular error source, in conjunction with calibration, can now be eliminated completely by using the linear construction technique outlined in Section X.

The histogram in FIG. 8 was made for 300 μs AMFM pulses. The question is: How does this scale with the gate time τ. To answer this question, 36-ion AMFM histograms at 700 μs were also computed. The result is that most of the Φ-infidelities for the 700 μs AMFM histogram are below 5 pptt. About half of the gates are good gates with infidelities less than 1 pptt and most of the other gates have Φ-infidelities less than 5 pptt. This indicates that the Φ-infidelity is a sensitive function of gate time τ.

XII: Gate Simulator Matrix Elements

In this section, the matrix elements of the full Hamiltonian Ĥ_MS are presented as well as the ones of the model Hamiltonians Ĥ^(Nc,Ns). The matrix elements of the full Hamiltonian, Ĥ_MS are calculated.

The matrix elements of the cosine part of Ĥ_MS are

$\begin{array}{cc}{C}_{n_1{n}_2\dots ;m_1{m}_2\dots }^{\left(j\right)}=\langle n_1{n}_2\dots ❘\cos \left[{\sum }_p{\eta }_p^j\left({\hat{a}}_p^{†}+{\hat{a}}_p\right)\right]❘m_1{m}_2\dots \rangle =\frac{1}{2}\langle n_1{n}_2\dots ❘e^{i{\sum }_p{\eta }_p^i\left({\hat{a}}_p^{†}+{\hat{a}}_p\right)}+e^{-i{\sum }_p{\eta }_p^i\left({\hat{a}}_p^{†}+{\hat{a}}_p\right)}❘m_1{m}_2\dots \rangle =\frac{1}{2}\left\{{\prod }_p\langle n_p❘e^{i{\eta }_p^j\left({\hat{a}}_p^{†}+{\hat{a}}_p\right)}❘m_p\rangle +{\prod }_p\langle n_p❘e^{i{\eta }_p^j\left({\hat{a}}_p^{†}+{\hat{a}}_p\right)}❘m_p\rangle \right\}.& \left(123\right)\end{array}$

It can be seen that C_{n1n2. . . ;m1m2. . .} ^(j) is real and symmetric in n_p↔m_p. With n_p≥m_p and

$\begin{array}{cc}\langle n❘e^{\lambda {\hat{a}}^{†}}❘m\rangle =\{\begin{array}{cc}0& \mathrm{if}m>n\\ {\frac{{\lambda }^{n-m}}{\left(n-m\right)!}\left[\frac{n!}{m!}\right]}^{1/2}& \mathrm{if}m\le n,\end{array}& \left(124\right)\end{array}$ $\begin{array}{cc}{C}_{n_1{n}_2\dots ;m_1{m}_2\dots }^{\left(j\right)}=\frac{1}{2}{\prod }_p{e^{-{\left({\eta }_p^j\right)}^2/2}\left(\frac{m_p!}{n_p!}\right)}^{1/2}{\left({\eta }_p^j\right)}^{n_p-m_p}{L}_{m_p}^{\left(n_p-m_p\right)}\left[{\left({\eta }_p^j\right)}^2\right]\left\{{\prod }_p{i}^{\left(n_p-m_p\right)}+{\prod }_p{\left(-i\right)}^{\left(n_p-m_p\right)}\right\},& \left(125\right)\end{array}$

are obtained. Now, define

$\begin{array}{cc}\sigma =\left[{\sum }_p\left(n_p-m_p\right)\right]\mathrm{mod}4.& \left(126\right)\end{array}$

Then:

$\begin{array}{cc}{C}_{n_1{n}_2\dots ;m_1{m}_2\dots }^{\left(j\right)}={\prod }_p{e^{-{\left({\eta }_p^j\right)}^2/2}\left(\frac{m_p!}{n_p!}\right)}^{1/2}{\left({\eta }_p^j\right)}^{n_p-m_p}{L}_{m_p}^{\left(n_p-m_p\right)}\left[{\left({\eta }_p^j\right)}^2\right]\{\begin{array}{cc}1,& \mathrm{if}\sigma =0,\\ 0,& \mathrm{if}\sigma =1,\\ -1,& \mathrm{if}\sigma =2,\\ 0,& \mathrm{if}\sigma =3.\end{array}& \left(127\right)\end{array}$

Because of

$\begin{array}{cc}{C}_{n_1{n}_2\dots ;m_1{m}_2\dots }^{\left(j\right)}~{\prod }_p{\left({\eta }_p^j\right)}^{n_p-m_p},& \left(128\right)\end{array}$

and |η_p^j|<<1, it can be seen that C_{n1n2. . . ;m1m2. . .} ^(j) is very close to diagonal, i.e., only the first few off-diagonals are significantly different from zero. This fact can be used to speed up the numerical integration of the system of linear equations significantly.

Similarly,

$\begin{array}{cc}{S}_{n_1{n}_2\dots ;m_1{m}_2\dots }^{\left(j\right)}=\langle n_1{n}_2\dots ❘\sin \left[{\sum }_p{\eta }_p^j\left({\hat{a}}_p^{†}+{\hat{a}}_p\right)\right]❘m_1{m}_2\dots \rangle ={\prod }_p{e^{-{\left({\eta }_p^j\right)}^2/2}\left(\frac{m_p!}{n_p!}\right)}^{1/2}{\left({\eta }_p^j\right)}^{n_p-m_p}{L}_{m_p}^{\left(n_p-m_p\right)}\left[{\left({\eta }_p^j\right)}^2\right]\{\begin{array}{cc}0,& \mathrm{if}\sigma =0,\\ 1,& \mathrm{if}\sigma =1,\\ 0,& \mathrm{if}\sigma =2,\\ -1,& \mathrm{if}\sigma =3,\end{array}& \left(129\right)\end{array}$

are obtained.

XIII: Third-Order Commutators

In this section, the results of some third-order commutators T_jkl as defined in (66) with η orders η^m, m≤4 are listed, grouped according to their η order m.

Commutators {circumflex over (T)}_jkl˜η⁰:

$\begin{array}{cc}{\hat{T}}_{000}=0,& \left(130\right)\end{array}$

Commutators {circumflex over (T)}_jkl˜η¹:

$\begin{array}{cc}{\hat{T}}_{001}=\left(\frac{2}{3}\right){\sum }_{j=1,2}\left[{\int }_0^{\tau }g\left(t\right)G^2\left(t\right){\hat{V}}_j\left(t\right)\mathrm{dt}\right]{\hat{\sigma }}_x^{\left(j\right)},& \left(131\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{010}=2{\hat{T}}_{001},& \left(132\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{100}=0,& \left(133\right)\end{array}$

Commutators {circumflex over (T)}_jkl˜η²:

$\begin{array}{cc}{\hat{T}}_{002}=0,& \left(134\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{011}=\left(\frac{8}{3}\right)\Phi \left[\eta ,g\right]\left[{\hat{\sigma }}_x^{\left(1\right)}{\hat{\sigma }}_z^{\left(2\right)}+{\hat{\sigma }}_x^{\left(2\right)}{\hat{\sigma }}_z^{\left(1\right)}\right],& \left(135\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{020}=0,& \left(136\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{101}=\left(-\frac{2}{3}\right){\sum }_{j=1,2}{\hat{\sigma }}_y^{\left(j\right)}{\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)\left[G\left(t_1\right)-G\left(t_2\right)\right]\left[{\hat{V}}_j\left(t_1\right){\hat{V}}_j\left(t_2\right)+{\hat{V}}_j\left(t_2\right){\hat{V}}_j\left(t_1\right)\right],& \left(137\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{110}=\left(-\frac{2}{3}\right){\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)G\left(t_1\right){\sum }_j\left[{\hat{V}}_j\left(t_1\right){\hat{V}}_j\left(t_2\right)+{\hat{V}}_j\left(t_2\right){\hat{V}}_j\left(t_1\right)\right]{\hat{\sigma }}_y^{\left(j\right)}+\left(\frac{4}{3}\right)\Phi \left[\eta ,g\right]\left[{\hat{\sigma }}_x^{\left(1\right)}{\hat{\sigma }}_z^{\left(2\right)}+{\hat{\sigma }}_x^{\left(2\right)}{\hat{\sigma }}_z^{\left(1\right)}\right],& \left(138\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{200}=0.& \left(139\right)\end{array}$

Commutators {circumflex over (T)}_jkl˜η³:

$\begin{array}{cc}{\hat{T}}_{003}=\left(-\frac{1}{9}\right){\sum }_j\left\{{\int }_0^{\tau }g\left(t\right)G^2\left(t\right){\hat{V}}_j^3\left(t\right)\mathrm{dt}\right\}{\hat{\sigma }}_x^{\left(j\right)},& \left(140\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{012}=\left(\frac{1}{6}\right){\sum }_k{\hat{\sigma }}_x^{\left(k\right)}{\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)\left\{G\left(t_2\right)\left[{\hat{V}}_k\left(t_2\right){\hat{V}}_k^2\left(t_1\right)+{\hat{V}}_k^2\left(t_1\right){\hat{V}}_k\left(t_2\right)\right]-G\left(t_1\right)\left[{\hat{V}}_k\left(t_1\right){\hat{V}}_k^2\left(t_2\right)+{\hat{V}}_k^2\left(t_2\right){\hat{V}}_k\left(t_1\right)\right]\right\}-\left(\frac{2}{3}\right){\sum }_p{\eta }_p^1{\eta }_p^2{\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)\sin \left[{\omega }_p\left(t_1-t_2\right)\right]\left\{\left[G\left(t_2\right){\hat{V}}_1\left(t_1\right)+G\left(t_1\right){\hat{V}}_1\left(t_2\right)\right]{\hat{\sigma }}_y^{\left(1\right)}{\hat{\sigma }}_z^{\left(2\right)}+\left[G\left(t_2\right){\hat{V}}_2\left(t_1\right)+G\left(t_1\right){\hat{V}}_2\left(t_2\right)\right]{\hat{\sigma }}_y^{\left(2\right)}{\hat{\sigma }}_z^{\left(1\right)}\right\},& \left(141\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{021}=\left(\frac{1}{3}\right){\sum }_k{\hat{\sigma }}_x^{\left(k\right)}{\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)G\left(t_1\right)\left[{\hat{V}}_k^2\left(t_1\right){\hat{V}}_k\left(t_2\right)+{\hat{V}}_k\left(t_2\right){\hat{V}}_k^2\left(t_1\right)\right]-\left(\frac{2}{3}\right){\sum }_p{\eta }_p^1{\eta }_p^2{\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)\sin \left[{\omega }_p\left(t_1-t_2\right)\right]\left\{\left[G\left(t_2\right){\hat{V}}_1\left(t_2\right)+G\left(t_1\right){\hat{V}}_1\left(t_1\right)\right]{\hat{\sigma }}_y^{\left(1\right)}{\hat{\sigma }}_z^{\left(2\right)}+\left[G\left(t_2\right){\hat{V}}_2\left(t_2\right)+G\left(t_1\right){\hat{V}}_2\left(t_1\right)\right]{\hat{\sigma }}_y^{\left(2\right)}{\hat{\sigma }}_z^{\left(1\right)}\right\},& \left(142\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{030}=\left(-\frac{2}{9}\right){\sum }_{j=1,2}{\hat{\sigma }}_x^{\left(j\right)}{\int }_0^{\tau }g\left(t\right)G^2\left(t\right){\hat{V}}_j^3\left(t\right)\mathrm{dt},& \left(143\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{102}=0,& \left(144\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{111}=0,& \left(145\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{120}=0,& \left(146\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{021}=\left(-\frac{1}{6}\right){\sum }_j{\hat{\sigma }}_x^{\left(j\right)}{\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)\left[G\left(t_1\right)-G\left(t_2\right)\right]\left\{{\hat{V}}_j^2\left(t_1\right){\hat{V}}_j\left(t_2\right)+{\hat{V}}_j\left(t_2\right){\hat{V}}_j^2\left(t_1\right)+{\hat{V}}_j^2\left(t_2\right){\hat{V}}_j\left(t_1\right)+{\hat{V}}_j\left(t_1\right){\hat{V}}_j^2\left(t_2\right)\right\}+\left(\frac{2}{3}\right){\sum }_{p,j\ne k}{\hat{\sigma }}_y^{\left(j\right)}{\hat{\sigma }}_z^{\left(k\right)}{\eta }_p^j{\eta }_p^k{\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)\left[G\left(t_1\right)-G\left(t_2\right)\right]\sin \left[{\omega }_p\left(t_1-t_2\right)\right]{\widehat{[V}}_j\left(t_1\right)-{\hat{V}}_j\left(t_2\right)],& \left(147\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{210}=\left(\frac{1}{6}\right){\sum }_j{\hat{\sigma }}_x^{\left(j\right)}{\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)\left\{G\left(t_2\right)\left[{\hat{V}}_j^2\left(t_1\right){\hat{V}}_j\left(t_2\right)+{\hat{V}}_j\left(t_2\right){\hat{V}}_j^2\left(t_1\right)\right]-G\left(t_1\right)\left[{\hat{V}}_j^2\left(t_2\right){\hat{V}}_j\left(t_1\right)+{\hat{V}}_j\left(t_1\right){\hat{V}}_j^2\left(t_2\right)\right]\right\}-\left(\frac{2}{3}\right){\sum }_{p,j\ne k}{\hat{\sigma }}_y^{\left(j\right)}{\hat{\sigma }}_z^{\left(k\right)}{\eta }_p^j{\eta }_p^k{\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)\sin \left[{\omega }_p\left(t_1-t_2\right)\right]\left[G\left(t_2\right){\hat{V}}_j\left(t_1\right)+G\left(t_1\right){\hat{V}}_j\left(t_2\right)\right],& \left(148\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{300}=0.& \left(149\right)\end{array}$

Commutators {circumflex over (T)}_jkl˜η⁴:

$\begin{array}{cc}{\hat{T}}_{004}=0,& \left(150\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{013}=\left(-\frac{1}{3}\right)\left[{\hat{\sigma }}_x^{\left(1\right)}{\hat{\sigma }}_z^{\left(2\right)}+{\hat{\sigma }}_x^{\left(2\right)}{\hat{\sigma }}_z^{\left(1\right)}\right]{\sum }_{p,j=1,2}{\eta }_p^1{\eta }_p^2{\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)\sin \left[{\omega }_p\left(t_1-t_2\right)\right]\left[G\left(t_1\right){\hat{V}}_j^2\left(t_2\right)+G\left(t_2\right){\hat{V}}_j^2\left(t_1\right)\right],& \left(151\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{022}=0& \left(152\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{031}=\left(-\frac{2}{3}\right)\left[{\hat{\sigma }}_x^{\left(1\right)}{\hat{\sigma }}_z^{\left(2\right)}+{\hat{\sigma }}_x^{\left(2\right)}{\hat{\sigma }}_z^{\left(1\right)}\right]{\sum }_{p,j=1,2}{\eta }_p^1{\eta }_p^2{\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)G\left(t_1\right)\sin \left[{\omega }_p\left(t_1-t_2\right)\right]{\hat{V}}_j^2\left(t_1\right),& \left(153\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{040}=0& \left(154\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{103}={\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)\left[G\left(t_1\right)-G\left(t_2\right)\right]\left\{\left(\frac{1}{18}\right){\sum }_{j=1,2}{\hat{\sigma }}_y^{\left(j\right)}\left[{\hat{V}}_j\left(t_1\right){\hat{V}}_j^3\left(t_2\right)+{\hat{V}}_j^3\left(t_2\right){\hat{V}}_j\left(t_1\right)+{\hat{V}}_j\left(t_2\right){\hat{V}}_j^3\left(t_1\right)+{\hat{V}}_j^3\left(t_1\right){\hat{V}}_j\left(t_2\right)\right]+\left(\frac{1}{3}\right){\sum }_{j\ne k}{\hat{\sigma }}_x^{\left(j\right)}{\hat{\sigma }}_z^{\left(k\right)}{\sum }_p{\eta }_p^j{\eta }_p^k\left[{\hat{V}}_k^2\left(t_2\right)-{\hat{V}}_k^2\left(t_1\right)\right]\sin \left[{\omega }_p\left(t_1-t_2\right)\right]\right\},& \left(155\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{112}=\left(\frac{1}{6}\right){\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_2{\int }_0^{t_2}{\mathrm{dt}}_{3}g\left(t_1\right)g\left(t_2\right)g\left(t_3\right)\{\left(-\frac{1}{2}\right){\sum }_j{\hat{\sigma }}_y^{\left(j\right)}\left\{{\hat{V}}_j\left(t_1\right)\left[{\hat{V}}_j\left(t_2\right){\hat{V}}_j^2\left(t_3\right)+{\hat{V}}_j^2\left(t_3\right){\hat{V}}_j\left(t_2\right)\right]+\left[{\hat{V}}_j\left(t_2\right){\hat{V}}_j^2\left(t_3\right)+{\hat{V}}_j^2\left(t_3\right){\hat{V}}_j\left(t_2\right)\right]{\hat{V}}_j\left(t_1\right)+{\hat{V}}_j\left(t_3\right)\left[{\hat{V}}_j\left(t_2\right){\hat{V}}_j^2\left(t_1\right)+{\hat{V}}_j^2\left(t_1\right){\hat{V}}_j\left(t_2\right)\right]+\left[{\hat{V}}_j\left(t_2\right){\hat{V}}_j^2\left(t_1\right)+{\hat{V}}_j^2\left(t_1\right){\hat{V}}_j\left(t_2\right)\right]{\hat{V}}_j\left(t_3\right)\right\}-2{\sum }_{p,j\ne k}{\hat{\sigma }}_x^{\left(j\right)}{\hat{\sigma }}_z^{\left(k\right)}{\eta }_p^j{\eta }_p^k\left\{{\hat{V}}_k^2\left(t_3\right)\sin \left[{\omega }_p\left(t_1-t_2\right)\right]-V_k^2\left(t_1\right)\sin \left[{\omega }_p\left(t_2-t_3\right)\right]+\left[{\hat{V}}_k\left(t_2\right){\hat{V}}_k\left(t_3\right)+{\hat{V}}_k\left(t_3\right){\hat{V}}_k\left(t_2\right)\right]\sin \left[{\omega }_p\left(t_1-t_3\right)\right]-\left[{\hat{V}}_k\left(t_2\right){\hat{V}}_k\left(t_1\right)+{\hat{V}}_k\left(t_1\right){\hat{V}}_k\left(t_1\right){\hat{V}}_k\left(t_2\right)\right]\sin \left[{\omega }_p\left(t_1-t_3\right)\right]\right\}+4{\sum }_{\mathrm{pq}}{\eta }_p^1{\eta }_p^2{\eta }_q^1{\eta }_q^2\left({\hat{\sigma }}_y^{\left(1\right)}+{\hat{\sigma }}_y^{\left(2\right)}\right)\sin \left[{\omega }_q\left(t_1-t_3\right)\right]\left\{\sin \left[{\omega }_p\left(t_2-t_3\right)\right]-\sin \left[{\omega }_p\left(t_2-t_1\right)\right]\right\}-2{\sum }_p{\eta }_p^1{\eta }_p^2{\hat{\sigma }}_x^{\left(1\right)}{\hat{\sigma }}_z^{\left(2\right)}\left[{\hat{V}}_2\left(t_1\right){\hat{V}}_2\left(t_3\right)+{\hat{V}}_2\left(t_3\right){\hat{V}}_2\left(t_1\right)\right]\left\{\sin \left[{\omega }_p\left(t_2-t_3\right)\right]+\left\{\sin \left[{\omega }_p\left(t_2-t_1\right)\right]\right\}\right\},& \left(156\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{121}=\left(-\frac{1}{12}\right){\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_2{\int }_0^{t_2}{\mathrm{dt}}_{3}g\left(t_1\right)g\left(t_2\right)g\left(t_3\right){\sum }_{\mathrm{jkl}}\left\{2\left[{\hat{V}}_j\left(t_1\right){\hat{V}}_k^2\left(t_2\right){\hat{V}}_l\left(t_3\right)+{\hat{V}}_j\left(t_3\right){\hat{V}}_k^2\left(t_2\right){\hat{V}}_l\left(t_1\right)\right]{\hat{\sigma }}_x^{\left(j\right)}{\hat{\sigma }}_y^{\left(k\right)}{\hat{\sigma }}_x^{\left(l\right)}-\left[{\hat{V}}_j\left(t_1\right){\hat{V}}_k\left(t_3\right){\hat{V}}_l^2\left(t_2\right)+{\hat{V}}_j\left(t_3\right){\hat{V}}_k\left(t_1\right){\hat{V}}_l^2\left(t_2\right)\right]{\hat{\sigma }}_x^{\left(j\right)}{\hat{\sigma }}_x^{\left(k\right)}{\hat{\sigma }}_y^{\left(l\right)}-\left[{\hat{V}}_j^2\left(t_2\right){\hat{V}}_k\left(t_3\right){\hat{V}}_l\left(t_1\right)+{\hat{V}}_j^2\left(t_2\right){\hat{V}}_k\left(t_1\right){\hat{V}}_l\left(t_3\right)\right]{\hat{\sigma }}_y^{\left(j\right)}{\hat{\sigma }}_x^{\left(k\right)}{\hat{\sigma }}_x^{\left(l\right)}\right\},& \left(157\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{130}=\left(\frac{1}{9}\right){\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right)G\left(t_1\right)\left\{{\sum }_{j=1,2}{\hat{\sigma }}_y^{\left(j\right)}\left[{\hat{V}}_j\left(t_2\right){\hat{V}}_j^3\left(t_1\right)+{\hat{V}}_j^3\left(t_1\right){\hat{V}}_j\left(t_2\right)\right]\right\}-\left(\frac{2}{3}\right){\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_{2}g\left(t_1\right)g\left(t_2\right){\sum }_{\begin{array}{c}j,k=1,2\\ j\ne k\end{array}}{\sum }_{p=1}^N{\eta }_p^j{\eta }_p^k{\sigma }_x^j{\sigma }_z^{k}G\left(t_1\right){\hat{V}}_k^2\left(t_1\right)\sin \left[{\omega }_p\left(t_1-t_2\right)\right],& \left(158\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{202}=0,& \left(159\right)\end{array}$ $\begin{array}{cc}{\hat{T}}_{211}=\left(\frac{4}{3}\right){\sum }_{p,k\ne j}{\eta }_p^1{\eta }_p^2{\hat{\sigma }}_z^{\left(j\right)}{\hat{\sigma }}_x^{\left(k\right)}{\int }_0^{\tau }{\mathrm{dt}}_1{\int }_0^{t_1}{\mathrm{dt}}_2{\int }_0^{t_2}{\mathrm{dt}}_{3}g\left(t_1\right)g\left(t_2\right)g\left(t_3\right){\hat{V}}_j^2\left(t_1\right)\sin \left[{\omega }_p\left(t_2-t_3\right)\right],& \left(160\right)\end{array}$

Claims

1. A method of performing a two-qubit gate operation, comprising: computing, by a classical computer, a control pulse to be applied to a pair of trapped ions in a plurality of trapped ions in a quantum processor, each of the plurality of trapped ions having two frequency-separated states defining a qubit, wherein computing the control pulse comprises: computing a pulse function of the control pulse based on a phase-space closure condition and an auxiliary condition; andcomputing the pulse function of the control pulse further based on a gate angle condition; andapplying, by a system controller, the control pulse, having the computed pulse function, to the pair of trapped ions.

2. The method of claim 1, wherein the computing of the pulse function based on the phase-space closure condition and the auxiliary condition comprises computing Fourier coefficients of the pulse function of the control pulse by solving a set of linear equations with respect to the Fourier coefficients.

3. The method of claim 2, wherein the phase-space closure condition requires the pair of trapped ions that are excited by a delivery of the control pulse return to their initial motional state.

4. The method of claim 3, wherein the auxiliary condition implies the vanishing of the Φ functional.

5. The method of claim 1, further comprising: prior to the application of the control pulse, calibrating the pulse function of the control pulse by multiplying the pulse function by a calibration factor.

6. The method of claim 5, wherein the calibration factor is based on a gate angle error and a target gate angle determined by the gate angle condition.

7. The method of claim 6, wherein the target gate angle is

8. A quantum computing system, comprising: a quantum processor comprising a plurality of trapped ions, wherein each of the trapped ions having two frequency-separated states defining a qubit;a classical computer configured to: compute a control pulse to be applied to a pair of trapped ions in the quantum processor, wherein computing the control pulse comprises: computing a pulse function of the control pulse based on a phase-space closure condition and an auxiliary condition; andcomputing the pulse function of the control pulse further based on a gate angle condition; anda system controller configured to: apply the control pulse, having the computed pulse function, to the pair of trapped ions.

9. The quantum computing system of claim 8, wherein the computing of the pulse function, based on the phase-space closure condition and the auxiliary condition, comprises computing Fourier coefficients of the pulse function of the control pulse by solving a set of linear equations with respect to the Fourier coefficients.

10. The quantum computing system of claim 9, wherein the phase-space closure condition requires the pair of trapped ions that are excited by a delivery of the control pulse return to their initial motional state.

11. The quantum computing system of claim 10, wherein the auxiliary condition causes the Φ functional to vanish.

12. The quantum computing system of claim 8, further comprising: prior to the applying of the control pulse, calibrating the pulse function of the control pulse by multiplying the pulse function by a calibration factor.

13. The quantum computing system of claim 12, wherein the calibration factor is based on a gate angle error and a target gate angle determined by the gate angle condition.

14. The quantum computing system of claim 13, wherein the target gate angle is

15. A quantum computing system comprising non-volatile memory having a number of instructions stored therein which, when executed by one or more processors, causes the quantum computing system to perform operations comprising: computing, by a classical computer, a control pulse to be applied to a pair of trapped ions in a plurality of trapped ions in a quantum processor, each of the plurality of trapped ions having two frequency-separated states defining a qubit, wherein computing the control pulse comprises: computing a pulse function of the control pulse based on a phase-space closure condition and an auxiliary condition; andcomputing the pulse function of the control pulse further based on the gate angle condition; andapplying, by a system controller, the control pulse, having the computed pulse function, to the pair of trapped ions.

16. The quantum computing system of claim 15, wherein the computing of the pulse function, based on the phase-space closure condition and the auxiliary condition, comprises computing Fourier coefficients of the pulse function of the control pulse by solving a set of linear equations with respect to the Fourier coefficients.

17. The quantum computing system of claim 16, wherein the phase-space closure condition requires the pair of trapped ions that are excited by a delivery of the control pulse return to their initial motional state.

18. The quantum computing system of claim 17, wherein the auxiliary condition implies the vanishing of the Φ functional.

19. The quantum computing system of claim 15, further comprising: prior to the applying of the control pulse, calibrating the pulse function of the control pulse by multiplying the pulse function by a calibration factor.

20. The quantum computing system of claim 19, wherein the calibration factor is based on a gate angle error and a target gate angle determined by the gate angle condition.