Patent Publication Number: US-11645569-B2

Title: Stabilization of entangling gates for trapped-ion quantum computers

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
     This application is a continuation of co-pending U.S. patent application Ser. No. 16/854,033, filed Apr. 21, 2020, which claims the benefit to U.S. Provisional Application No. 62/851,280, filed May 22, 2019, and the benefit to U.S. Provisional Application No. 62/851,394, filed May 22, 2019. Each of the aforementioned related patent applications are incorporated by reference herein. 
    
    
     BACKGROUND 
     Field 
     The present disclosure generally relates to a method of generating an entangling gate in an ion trap quantum computer, and more specifically, to a method of stabilizing the entangling gate against fluctuations in frequencies of collective motional modes of trapped ions. 
     Description of the Related Art 
     In quantum computing, quantum bits or qubits, which are analogous to bits representing a “0” and a “1” in a classical (digital) computer, are required to be prepared, manipulated, and measured (read-out) with near perfect control during a computation process. Imperfect control of the qubits leads to errors that can accumulate over the computation process, limiting the size of a quantum computer that can perform reliable computations. 
     Among physical systems upon which it is proposed to build large-scale quantum computers, is a chain of ions (e.g., 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 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 their Coulombic interaction between the ions. As the size of a quantum computer increases, implementation of two-qubit gate operations between a pair of ions increases complexity, and thus errors associated with the implementation and resources required for the implementation, such as frequencies of the collective motional modes and laser powers, increase. 
     To increase the size of a quantum computer that may be able to implement algorithms to solve problems otherwise intractable in classical computer, there is a need for a procedure to accurately control qubits to perform in the presence of errors associated with the implementation and resources. 
     SUMMARY 
     Embodiments of the present disclosure generally relate to a method of performing a computation using a quantum computer includes generating a first laser pulse and a second laser pulse to cause entanglement interaction between a first trapped ion and a second trapped ion of a plurality of trapped ions that are aligned in a first direction, each of the plurality of trapped ions having two frequency-separated states defining a qubit, and applying the generated first laser pulse to the first trapped ion and the generated second laser pulse to the second trapped ion. Generating the first laser pulse and the second laser pulse includes stabilizing the entanglement interaction between the first and second trapped ions against fluctuations in frequencies of collective motional modes of the plurality of trapped ions in a second direction that is perpendicular to the first direction. 
    
    
     
       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 partial view of an ion trap quantum computer according to one embodiment. 
         FIG.  2    depicts a schematic view of an ion trap for confining ions in a chain according to one embodiment. 
         FIGS.  3 A,  3 B, and  3 C  depict a few schematic collective transverse motional mode structures of a chain of five trapped ions. 
         FIG.  4    depicts a schematic energy diagram of each ion in a chain of trapped ions according to one embodiment. 
         FIG.  5    depicts a qubit state of an ion represented as a point on a surface of the Bloch sphere. 
         FIGS.  6 A and  6 B  depict schematic views of motional sideband spectrum of each ion and a motional mode according to one embodiment. 
         FIG.  7 A  depicts infidelity of an XX-gate operation for a chain of five trapped ions according to one embodiment. 
         FIG.  7 B  depicts width of infidelity in  FIG.  7 A  of an XX-gate operation. 
         FIG.  7 C  depicts power requirement of stabilized pulses according to one embodiment. 
         FIG.  7 D  depicts effect of stabilization according to one embodiment. 
     
    
    
     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 
     An overall system that is able to perform quantum computations using trapped ions will include a classical computer, a system controller, and a quantum register. The classical computer performs supporting and system control tasks including selecting a quantum algorithm to be run by use of a user interface, such as graphics processing unit (GPU), compiling the selected quantum algorithm into a series of universal logic gates, translating the series of universal logic gates into laser pulses to apply on the quantum register, and pre-calculating parameters that optimize the laser pulses by use of a central processing unit (CPU). A software program for performing the task of decomposing and executing the quantum algorithms is stored in a non-volatile memory within the classical computer. The quantum register includes trapped ions that are coupled with various hardware, including lasers to manipulate internal hyperfine states (qubit states) of the trapped ions and an acousto-optic modulator to read-out the internal hyperfine states (qubit states) of the trapped ions. The system controller receives from the classical computer the pre-calculated parameters for power-optimal pulses at the beginning of running the selected algorithm on the quantum register, controls various hardware associated with controlling any and all aspects used to run the selected algorithm on the quantum register, and returns a read-out of the quantum register and thus output of results of the quantum computation(s) at the end of running the algorithm to the classical computer. 
     Among several known sets of universal logic gates by which any quantum algorithm can be decomposed, a set of universal 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. For those of ordinary skill in the art, it should be clear the R gate can be implemented with near perfect fidelity, while the formation of the XX gate is complex and requires optimization for a given type of trapped ions, number of ions in a chain of trapped ions, and the hardware and environment in which the trapped ions are trapped, to name just a few factors, such that the fidelity of the XX gate is increased and computational errors within a quantum computer are avoided or decreased. 
     In the following discussion, methods of generating and optimizing a pulse used to perform computations based the formation of an XX gate that has an improved fidelity in the presence of errors associated with the implementation and resources required for the implementation will be described. 
     General Hardware Configurations 
       FIG.  1    is a partial view of an ion trap quantum computer, or system  100 , according to one embodiment. The system  100  includes a classical (digital) computer  101 , a system controller  118  and a quantum register that is a chain  102  of trapped ions (i.e., five shown) that extend along the Z-axis. The classical computer  101  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  104 , 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)  106  for measurement of individual ions. Non-copropagating Raman laser beams from a laser  108 , which are provided along the X-axis, perform operations on the ions. A diffractive beam splitter  110  creates an array of static Raman beams  112  that are individually switched using a multi-channel acousto-optic modulator (AOM)  114  and is configured to selectively act on individual ions. A global Raman laser beam  116  illuminates all ions at once. The system controller (also referred to as a “RF controller”)  118  controls the AOM  114 . The system controller  118  includes a central processing unit (CPU)  120 , a read-only memory (ROM)  122 , a random access memory (RAM)  124 , a storage unit  126 , and the like. The CPU  120  is a processor of the RF controller  118 . The ROM  122  stores various programs and the RAM  124  is the working memory for various programs and data. The storage unit  126  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  120 , the ROM  122 , the RAM  124 , and the storage unit  126  are interconnected via a bus  128 . The RF controller  118  executes a control program which is stored in the ROM  122  or the storage unit  126  and uses the RAM  124  as a working area. The control program will include one or more software applications that include program code (e.g., instructions) that may be executed by 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  102  according to one embodiment. The confining potential is exerted by both static (DC) voltage and radio frequency (RF) voltages. A static (DC) voltage V s  is applied to end-cap electrodes  210  and  212  to confine the ions along the Z-axis (also referred to as an “axial direction,” “longitudinal direction,” or “first direction”). The ions in the chain  102  are nearly evenly distributed 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 1  (with an amplitude V RF /2) is applied to an opposing pair of the electrodes  202 ,  204  and a sinusoidal voltage V 2  with a phase shift of 180° from the sinusoidal voltage V 1  (and the amplitude V RF /2) is applied to the other opposing pair of the 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 the 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,” “transverse direction,” or “second direction”) for each of the trapped ions, which is proportional to a 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, as is discussed in greater detail below. 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. 
     Trapped Ion Configuration and Quantum Bit Information 
       FIGS.  3 A,  3 B, and  3 C  depict a few schematic structures of collective transverse motional modes (also referred to simply as “motional mode structures”) of a chain  102  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  102  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 energy is hereinafter referred to as |n ph     p , 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  102 .  FIGS.  3 A- 3 C  schematically illustrates examples of different types of collective transverse motional modes that may be experienced by five trapped ions that are positioned in a chain  102 .  FIG.  3 A  is a schematic view of a common motional mode |n ph     P  having the highest energy, where P is the number of motional modes. In the common motional mode |n   P , all ions oscillate in phase in the transverse direction.  FIG.  3 B  is a schematic view of a tilt motional mode |n ph     P−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.  3 C  is a schematic view of a higher-order motional mode |n ph     P−3  which has a lower energy than that of the tilt motional mode |n ph     P−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. 
       FIG.  4    depicts a schematic energy diagram  400  of each ion in the chain  102  of trapped ions according to one embodiment. In one example, each ion may be a positive Ytterbium ion,  171 Yb + , which has the  2 S 1/2  hyperfine states (i.e., two electronic states) with an energy split corresponding to a frequency difference (referred to as a “carrier frequency”) of ω 01 /2π=12.642821 GHz. A qubit is formed with the two hyperfine states, denoted as |0  and |1 , where the hyperfine ground state (i.e., the lower energy state of the  2 S 1/2  hyperfine states) is chosen to represent |0 . Hereinafter, the terms “hyperfine states,” “internal hyperfine states,” and “qubits” may be interchangeably used to represent |0  and |1 . Each ion may be cooled (i.e., kinetic energy of the ion may be reduced) to near the motional ground state |0   p  for any motional mode p 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 p denotes the motional ground state for a motional mode p of a chain  102  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  2 P 1/2  level (denoted as |e ). As shown in  FIG.  4   , a laser beam from the laser may be split into a pair of non-copropagating laser beams (a first laser beam with frequency ω 1  and a second laser beam with frequency ω 2 ) in the Raman configuration, and detuned by a one-photon transition detuning frequency Δ=ω 1 −ω 0e  with respect to the transition frequency ω 0e  between |0  and |e , as illustrated in  FIG.  4   . 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 a two-photon transition detuning frequency (also referred to simply as “detuning frequency”) δ=ω 1 −ω 2 −ω 01  (hereinafter denoted as ±μ,μ being a positive value), 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 a 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 configuration to manipulate internal hyperfine states of qubits (qubit states) may be referred to as a “composite pulse” or simply as a “pulse,” and the resulting time-dependent pattern of the two-photon Rabi frequency Ω(t) may be referred to as an “amplitude” of a pulse or simply as a “pulse,” which are illustrated and further described below. The detuning frequency δ=ω 1 −ω 2 −ω 01  may be referred to as detuning frequency of the composite pulse or 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. 
     It should be noted that the particular atomic species used in the discussion provided herein is just one example of atomic species which has 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.  5    is provided to help visualize a qubit state of an ion is represented as a point on a surface of the Bloch sphere  500  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 amplitudes of 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 amplitudes of 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 without loss of generality) 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 qubits 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  (e.g., |0 ±i|1  for ϕ=±π/2). 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. 
     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.  6 A and  6 B  schematically depict views of a motional sideband spectrum for an ion in the chain  102  in a motional mode |n ph     p  having frequency ω p  according to one embodiment. As illustrated in  FIG.  6 B , when the detuning frequency of the composite pulse is zero (i.e., a frequency difference between the first and second laser beams is tuned to the carrier frequency, δ=ω 1 −ω 2 −ω 01 =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, δ=ω 1 −ω 2 −ω 01 =μ&gt;0, referred to as a blue sideband), Rabi flopping between combined qubit-motional states |0 |n ph     p  and |1  |n ph +1   p occurs (i.e., a transition from the p-th motional mode with n-phonon excitations denoted by |n ph     p  to the p-th motional mode with (n ph +1)-phonon excitations denoted by |n ph +1   p  occurs when 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 the frequency ω p  of the motional mode |n ph     p , δ=ω 1 −ω 2 −ω 01 =−μ&lt;0, referred to as a red sideband), Rabi flopping between combined qubit-motional states |0 |n ph     p  and |1 |n ph −1   p  occurs (i.e., a transition from the motional mode |n ph     p  to the motional mode |n ph −1   p  with one less phonon excitations occurs when 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 |0 |n ph     p  into a superposition of |0 |n ph     p  and |1 |n ph −1   p . When the two-photon Rabi frequency Ω(t) is smaller as compared to the detuning frequency δ=ω 1 −ω 2 −ω 01 =±μ, 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. Entanglement between qubits 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 , and |1   i |1   j  as follows:
 
|0   i |0   j →|0   i |0   j   −i| 1   i |1   j  
 
|0   i |1   j →|0   i |1   j   −i| 1   i |0   j  
 
|1   i |0   j   →−i |0   i |1   j +|1   i |0   j  
 
|1   i |1   j   →−i |0   i |0   j +|1   i |1   j  
 
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 ph     p  is transformed into a superposition of |0   i |n ph     p  and |1   i |n ph +1   p , and thus the combined state of the two qubits and the motional mode is transformed into a superposition of |0   i |0   j |n ph     p  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 ph     p  is transformed to a superposition of |0   j |n ph     p  and |1   j |n ph +1   p , and the combined state |0   j |n ph +1   p  is transformed into a superposition of |0   j |n ph +1   p  and |1   j n ph     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 |0   j |n ph     p  into a superposition of |0   i |0   j |n ph     p  and |1   i |1   j |n ph     p , 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 motional mode having a different number of phonon excitations from the initial number of phonon excitations n ph  (i.e., |1   i |0   j |n ph +1   p  and |0   i |1   j |n ph −1   p ) 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 pulses on the sidebands for duration τ (referred to as a “gate duration”), having amplitude Ω (i) (t) and detuning frequency μ (i) (t), and amplitude Ω (i) (t) and detuning frequency μ (j) (t), respectively, can be described in terms of an entanglement interaction χ ij (τ) as follows:
 
|0   i |0   j →cos(χ ij (τ)/2)|0   i |0   j   −i sin(χ ij (τ)/2)|1   i |1   j  
 
|0   i |1   j →cos(χ ij (τ)/2)|0   i |1   j   −i sin(χ ij (τ)/2)|1   i |0   j  
 
|1   i |0   j   →−i sin(χ ij (τ)/2)|0   i |1   j +cos(χ ij (τ)/2)|1   i |0   j  
 
|1   i |1   j   →−i sin(χ ij (τ)/2)|0   i |0   j +cos(χ ij (τ)/2)|1   i |1   j  
 
where,
 
χ ij (τ)=4Σ p=1   P η p   (i) η p   (j) ∫ 0   τ   dt   2 ∫ 0   t     2     dt   1 Ω (i) ( t   2 )Ω (j) ( t   1 )sin(ψ (i) ( t   2 ))sin(ψ (j) ( t   1 ))sin[ω p ( t   2   −t   1 )]
 
and η p   (l) (l=i,j) is the Lamb-Dicke parameter that quantifies the coupling strength between the l-th ion and the p-th motional mode having the frequency ω p ,ψ (l) (t) is an accumulated phase ψ (l) (t)=ψ 0   (l) +∫ 0   t μ (l) (t′)dt′ of the pulse, ψ 0   (l)  is an initial phase which may be assumed to be zero (0) hereinafter for simplicity without loss of generality, and P is the number of the motional modes (equal to the number N of ions in the chain  102 ).
 
     Construction of Pulses for Entangling Gate Operations 
     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 operation (R gate) forms a set of universal gates {R, XX} that can be used to build a quantum computer that is configured to perform desired computational processes. To perform an XX-gate operation between m-th and n-th qubits, pulses that satisfy the condition χ ij  (τ)=θ ij (0&lt;θ ij ≤π/2) (i.e., the entanglement interaction χ ij (τ) has a desired value θ ij , referred to as condition for a non-zero entanglement interaction) are constructed and applied to the i-th and the j-th qubits. The transformations of the combined state of the i-th and the j-th qubits described above corresponds to the XX-gate operation with maximal entanglement when θ ij =π/2. The amplitudes and the detuning frequencies, (Ω (i) (t),μ (i) (t)) and (Ω (j) (t),μ (j) (t)), of the pulses to be applied to the i-th and the j-th qubits are control parameters that can be adjusted to ensure a non-zero tunable entanglement of the i-th and the j-th qubits to perform a desired XX gate operation on i-th and j-th qubits. 
     The control parameters, the amplitudes and detuning frequencies (Ω (l) (t),μ (l) (t))(l=i,j) of the pulses, must also satisfy conditions that all of the N trapped ions in the chain  102  that are displaced from their initial positions as the motional modes are excited by the delivery of the pulse return to the initial positions. The l-th qubit in a superposition state |0 ±|1  is displaced due to the excitation of the p-th motional mode during the gate duration τ and follows the trajectories ±α p   (l) (t′) in phase space (position and momentum) of the p-th motional mode. The phase space trajectories α p   (l) (t′)=η p   (l) ∫ 0   t′ g (l) (t)e iω     p     t dt are determined by the amplitude Ω (l) (t) and the accumulated phase ψ (l) (t)=ψ 0   (l) +∫ 0   t μ (l) (t′)dt′ of the pulse, where g (l) (t) is a pulse function defined as g (l) (t)=Ω (l) (t) sin(ψ (l) (t)). Thus, for the chain  102  of N trapped ions, the condition π p   (l) (τ)=0 (i.e., the phase space trajectories π p   (l) (τ) of l-th qubit must be closed, referred to as a condition for returning of trapped ions to their original positions and momentum values (or closure of phase space trajectories) must be imposed for all the P motional modes (p=1, 2, . . . , P) in addition to the condition for non-zero entanglement interaction, χ ij (τ)=θ ij (0&lt;θ ij ≤π/2). 
     The control parameters, the amplitudes and detuning frequencies (Ω (l) (t),μ (l) (t)) (l=i,j) of the pulses are also adjusted such that the resulting pulse is power-optimal, in which the required laser power is minimized (referred to as condition for minimized power). Since the required laser power is inversely proportional to the gate duration τ, a power-optimal pulse implements an XX gate operation with minimum power requirement if gate duration τ is fixed, or with shortest gate duration τ if a laser power budget is fixed. 
     In some embodiments, the amplitudes and the detuning frequencies (Ω (l) (t),μ (l) (t)) (l=i,j) are chosen to be symmetric or anti-symmetric in time with respect to a middle point of the gate duration, t=τ/2, i.e., 
                         Ω       (   l   )     ⁢     (   ±   )         (       τ   2     -   t     )     =     ±       Ω       (   l   )     ⁢     (   ±   )         (       τ   2     +   t     )         ,        
         μ       (   l   )     ⁢     (   ±   )         (       τ   2     -   t     )     =     ±         μ       (   l   )     ⁢     (   ±   )         (       τ   2     +   t     )     .                             
In the example described below, the amplitudes and the detuning frequencies (Ω (l) (t), μ (l) (t)) are chosen to be symmetric (Ω (l)(+) (t) and μ (l)(+) (t)) for simplicity and may be referred to as Ω (l) (t) and μ (l) (t) without the subscript (+). With the symmetric detuning frequency μ (l) (t), the accumulated phase ψ (l) (t) is anti-symmetric, i.e.,
 
                       sin   ⁡   (       ψ     (   l   )       (       τ   2     -   t     )     )     =     -     sin   ⁡   (       ψ     (   l   )       (       τ   2     +   t     )     )         .                       
The condition for returning of trapped ions to their original positions and momentum values can be rewritten in terms of the anti-symmetric component g (l)(−) (t) of the pulse function g (l) (t) (simply referred to as “pulse function” and denoted by g (l) (t) without the script (−) hereinafter), as
 
     
       
         
           
             
               
                 
                   
                     
                       α 
                       
                         l 
                         , 
                         p 
                       
                     
                     ( 
                     τ 
                     ) 
                   
                   = 
                   
                     
                       
                         η 
                         p 
                         
                           ( 
                           l 
                           ) 
                         
                       
                       ⁢ 
                       
                         
                           ∫ 
                           0 
                           τ 
                         
                         
                           
                             
                               g 
                               
                                 ( 
                                 l 
                                 ) 
                               
                             
                             ( 
                             t 
                             ) 
                           
                           ⁢ 
                           
                             sin 
                             [ 
                             
                               
                                 ω 
                                 p 
                               
                               ( 
                               
                                 
                                   τ 
                                   2 
                                 
                                 - 
                                 t 
                               
                               ) 
                             
                             ] 
                           
                           ⁢ 
                           d 
                           ⁢ 
                           t 
                         
                       
                     
                     = 
                     
                       0 
                       . 
                     
                   
                 
               
               
                   
               
             
           
         
       
     
     The pulse functions g (l) (t) that are determined by the amplitudes and the detuning frequencies (Ω (l) (t),μ (l) (t)) of the pulse (l=i,j) are derived such that these conditions are satisfied, by dividing the gate duration into N A  equi-spaced segments (n=1, 2, . . . , N A ), and varying the amplitudes and the detuning frequencies (Ω (l) (t), μ (l) (t)) of the pulse, thus varying the pulse functions g (l) (t), from one segment to another. Here, for example, the pulse function g (l) (t) is expanded in Fourier-sine basis, 
                       g     (   l   )       (   t   )     =       ∑     n   =   1       N   A           g   n     (   l   )       ⁢     sin   ⁡   (     2   ⁢   π   ⁢   n   ⁢     t   τ       )                             
over the gate duration τ, using basis function
 
             sin   ⁡   (     2   ⁢   π   ⁢   n   ⁢     t   τ       )         
and pulse function coefficients g n   (l) . The condition for returning of trapped ions to their original positions and momentum values can be rewritten as
 
                       α   p     (   l   )       (   τ   )     =       0   ↦       ∑     n   =   1       N   A           M     p   ⁢   n       ⁢     g   n     (   l   )             =     0   ⁢     (       p   =   1     ,   2   ,   …         ,   P     )                             
where M pn is defined as
 
     
       
         
           
             
               M 
               
                 p 
                 ⁢ 
                 n 
               
             
             = 
             
               
                 ∫ 
                 0 
                 τ 
               
               
                 
                   sin 
                   ⁡ 
                   ( 
                   
                     2 
                     ⁢ 
                     π 
                     ⁢ 
                     n 
                     ⁢ 
                     
                       t 
                       τ 
                     
                   
                   ) 
                 
                 ⁢ 
                 sin 
                 ⁢ 
                 
                   ( 
                   
                     
                       ω 
                       p 
                     
                     ( 
                     
                       
                         τ 
                         2 
                       
                       - 
                       t 
                     
                     ) 
                   
                   ) 
                 
                 ⁢ 
                 d 
                 ⁢ 
                 
                   
                     t 
                     ⁡ 
                     ( 
                     
                       
                         p 
                         = 
                         1 
                       
                       , 
                       2 
                       , 
                       … 
                            
                       , 
                       P 
                       , 
                         
                       
                         n 
                         = 
                         1 
                       
                       , 
                       2 
                       , 
                       … 
                           
                       , 
                       
                         N 
                         A 
                       
                     
                     ) 
                   
                   . 
                 
               
             
           
         
       
     
     Equivalently, the condition for returning of trapped ions to their original positions and momentum values (e.g., closure of phase space trajectories) can be written as M   (l) =0 in a matrix form, where M is a P×N A  coefficient matrix of M pn  and    (l)  is a N A  pulse function coefficient vector of g n   (l) . The number of the segments N A  is chosen to be larger than the number of motional modes P. Thus, there are N 0 (=N A −P) non-trivial (i.e., at least one of the pulse function coefficients A n   (l)  is non-zero) pulse function vectors    (l)[α] (α=1, 2, . . . , N 0 )that satisfy the condition for returning of trapped ions to their original positions and momentum values. The non-trivial pulse function vectors    (l)[α] (α=1, 2, . . . , N 0 ) are null-space vectors of the matrix M (i.e., M   (l)[α] =0). 
     The conditions for non-zero entanglement interaction and decoupling can be rewritten 
                       χ     i   ⁢   j       (   τ   )     =         θ     i   ⁢   j       ↦       ∑     n   =   1       N   A           ∑       n   ′     =   1       N   A           g   n     (   i   )       ⁢     D     nn   ′       (     i   ,   j     )       ⁢     g     n   ′       (   j   )               =     θ     i   ⁢   j                             
where D nn′   (i,j)  is defined as
 
                       D     nn   ′       (     i   ,   j     )       =     4   ⁢       ∑     p   =   1     P         η   p     (   i   )       ⁢     η   p     (   j   )       ⁢       ∫   0   τ       d   ⁢     t   2     ⁢       ∫   0     t   2         d   ⁢     t   1     ⁢     sin   ⁡   (     2   ⁢   π   ⁢   n   ⁢       t   2     τ       )     ⁢     sin   ⁡   (     2   ⁢   π   ⁢     n   ′     ⁢       t   1     τ       )     ⁢     sin   [       ω   p     (       t   1     -     t   2       )     ]                     ,                       
or equivalently, (   (i) ) T D (i,j)     (j) =θ ij  in a matrix form, where D (i,j)  is a N A ×N A  coefficient matrix of D nn′   (i,j)  and (   (l) ) T  is a transposed vector of    (l) . Since (   (i) ) T D (i,j)     (j)  is a scalar, the condition for non-zero entanglement interaction can be further rewritten in terms of a symmetric matrix S (i,j) =(D (i,j) +D (i,j)     T   ) /2, as (   (i) ) T S (i,j)     (j) =θ ij .
 
     The condition for minimized power corresponds to minimizing a power function,
 
 P ( t )=∥ g   (l) ( t )∥ 2 =2/τ∫ 0   96    [g   (l) ( t )] 2   dt =Σ   n=1   N     A     [g   n   (l) ] 2  
 
that is the absolute square value of the pulse function g (l) (t) averaged over the gate duration τ. Accordingly, a power-optimal pulse can be constructed by computing a linear combination ({right arrow over (g (l) )}=Σ α=1   N     0   Λ α   (l)  {right arrow over (g)} α   (i) ) of the non-trivial pulse function vectors {right arrow over (g α   (l) )} (the null-space vectors of the matrix M) that satisfy the condition for returning of trapped ions to their original positions and momentum values, in which the coefficients Λ α   (l)  are determined such that the power function P(t) is minimized, while the condition for non-zero entanglement interaction is also satisfied. The condition for non-zero entanglement interaction leads to additional constrains on the coefficients Λ α   (l) , as {right arrow over (Λ)} (i)     T   R (i,j) {right arrow over (Λ)} (j) =θ i,j , where {right arrow over (Λ)} (l) (l=i,j) is a N 0  coefficient vector of Λ α   (l) (α=1, 2, . . . , N 0 ), and R (i,j)  is a real and symmetric matrix having matrix elements R αβ   (i,j) ={right arrow over (g α   (l) )} T S (i,j) {right arrow over (g β   (J) )}(α, β=1, 2, . . . , N 0 ). Thus, the amplitudes and the detuning frequencies (Ω (l) (t), μ (l) (t)) of the pulses (l=i,j) can be computed based on the pulse functions g (l) (t) having the pulse function coefficients g n   (l)  (n=1, 2, . . . , N A ) or equivalently the pulse function coefficient vector    (l) , that satisfy the condition for non-zero entanglement interaction and the condition for minimized power. It should be noted the condition for returning of trapped ions to their original positions and momentum values and the condition for non-zero entanglement interaction are in linear algebraic forms in terms of the pulse function coefficient vector    (l) . Thus, the pulse function coefficients g n   (l)  that satisfy these conditions, along with the condition for minimized power, can be computed by known linear algebraic computational methods without approximation or iterations.
 
     The expansion of the pulse function g (l) (t) in terms of pulse function coefficients g n   (l)  corresponds to construction of a pulse in a frequency domain (with a frequency 2πn/τ), and thus a pulse constructed by the pulse function g (l) (t) may be directly implemented by multi-tone lasers (i.e., laser beams having multiple tones, each tone having a distinct amplitude and a corresponding frequency). That is, N A -tone laser beams each having a frequency 2πn/τ and amplitude g  n   (l) (n=1, 2, . . . , N A ) , with the phases of the laser beams fixed, may directly perform an XX gate operation. The pulse function may be expanded using any functions that form a complete set or an incomplete set over the gate duration. However, when the pulse function is expanded in an incomplete set, there is no guarantee that the pulse function g (l) (t) computed by the method described above is power-optimal. 
     It should be noted that the particular example embodiments described above are just some possible examples of a method of construction of pulse functions according to the present disclosure and do not limit the possible configuration, specifications, or the like of methods of construction of pulse functions. For example, the symmetry of the amplitudes and the detuning frequencies (Ω l (t), μ (l) (t)) (l=i,j) may be selected to be anti-symmetric (having a negative parity), or have mixed symmetry (having a mixed parity), based on convenience related to configurations, specifications, or the like of the system  100 . However, imposing a symmetry in the amplitudes and the detuning frequencies (Ω l (t), μ (l) (t)) may lead to eliminating errors that have certain symmetries. 
     Stabilization Against Mode Frequency Fluctuations and Calibration of Xx Gate 
     In the construction of a power-optimal pulse, additional conditions can be added without substantially increasing the complexity in computing the Fourier coefficient vector    (l) (l=i,j) due to the linear algebraic forms of the conditions imposed on the coefficient vector    (l) (l=i,j). For example, a condition for stabilizing the pulse against external errors, such as fluctuations Δω p  in the frequencies ω p  of the motional modes and the intensities of laser beams, can be additionally imposed while maintaining the linear algebraic forms of all the conditions. In the ion trap quantum computer, or system  100 , there can be fluctuations Δω p  in the frequencies ω p  of the motional modes due to stray electric fields, build-up charges in the ion trap  200  caused by photoionization or temperature fluctuations. Typically over a time span of minutes, the frequencies ω p  of the motional modes drift with excursion of Δω p /(2π)≈1 kHz. The conditions for non-zero entanglement interaction, returning of trapped ions to their original positions and momentum values, and minimized power based on the frequencies ω p  of the motional modes are no longer satisfied when the frequencies ω p  of the motional modes have drifted to ω p +Δω p , resulting in a reduction of the fidelity of the XX gate operation. It is known the infidelity 1−F of an XX gate operation between i-th and j-th qubits at zero temperature of the motional-mode phonons, is given by 
               1   -   F     =       4   5     ⁢       ∑   p         (           ❘   &#34;\[LeftBracketingBar]&#34;       α     i   ,   p         ❘   &#34;\[RightBracketingBar]&#34;       2     +         ❘   &#34;\[LeftBracketingBar]&#34;       α     j   ,   p         ❘   &#34;\[RightBracketingBar]&#34;       2       )     .               
This suggests the XX-gate operation can be stabilized against fluctuations Δω p  in the frequencies ω p  of the motional modes by requiring that k-th order variations of the phase space trajectories α l,p  of the l-th qubit (l=i,j) with respect to the fluctuations Δω p  in ω p ,
 
                           ∂   k         α     l   ,   p       (   τ   )         ∂     ω   p   k         =         ∫   0   τ           (     i   ⁢   t     )     k     ⁢       g     (   l   )       (   t   )     ⁢     e     i   ⁢     ω   p     ⁢   t           =   0       ⁢         (       l   =   i     ,   j   ,         p   =   1     ,   2       ,   …         ,   P   ,         k   =   1     ,   2   ,   …         ,     K   α       )                         
be zero (referred to as k-th order stabilization), where K α  is a maximal desired degree of stabilization for α. The pulse computed by requiring this condition for stabilization can perform an XX gate operation that is resilient against noise (i.e., fluctuations Δω p  in the frequencies ω p  of the motional modes).
 
     The fluctuations Δω p  in the frequencies ω p  of the motional modes may also affect the value of the entanglement interaction χ ij (τ), since the entanglement interaction χ ij (τ) is related to the frequencies ω p  of the motional modes. That is, the resulting entanglement interaction χ ij (τ) may have a different value from the desired value θ ij  that is set in the condition for non-zero entanglement interaction, as discussed above. Thus, in some embodiments, the condition for stabilizing the entanglement interaction against fluctuations Δω p  in the frequencies ω p  of the motional modes may also require that k-th order variations of the entanglement interaction χ ij (τ) with respect to the fluctuations Δω p  in the frequencies ω p·  of the motional modes, 
                 χ       i   ⁢   j     ,   p       (   k   )       =           ∂   k       χ     i   ⁢   j           ∂     ω   p   k         =     0   ⁢     (       k   =   1     ,   2   ,     …   ⁢           K   χ         )           ,         
be zero (referred to as k-th order stabilization), where K χ  is a maximal desired degree of stabilization.
 
     As noted above, the condition for non-zero entanglement interaction for a power-optimal pulse is written as {right arrow over (Λ)} (i)     T   R (i,j) {right arrow over (Λ)} (j) =θ ij , where the coefficients Λ α   (l)  (the elements of the coefficient vector {right arrow over (Λ)} (l) ) in the null-space expansion ({right arrow over (g (l) )}=Σ α=1   N     0   Λ α   (l) {right arrow over (g)} α   (l) ) of the pulse function    (l)  of the power-optimal pulse (l=i,j). Thus, the condition for the stabilization of the entanglement interaction against fluctuations in the frequencies ω p  of the motional modes up to k-th order with respect to fluctuations Δω p  in ω p· (k-th order stabilization) is written in a matrix form,
 
χ ij,p   (k) ={right arrow over (Λ)} (i)     T     R   p   (ij)(k) {right arrow over (Λ)} (j) =0( k =1, 2, . . .  K   χ )
 
where the matrix R p   (i,j)(k)  has elements R αβ,p   (i,j)(k) ={right arrow over (g)} α   (i)     T   S p   (i,j)(k) {right arrow over (g)} β   (j)  and the matrix S p   (i,j)(k)  is defined as
 
     
       
         
           
             
               S 
               p 
               
                 
                   ( 
                   
                     i 
                     , 
                     j 
                   
                   ) 
                 
                 ⁢ 
                 
                   ( 
                   k 
                   ) 
                 
               
             
             = 
             
               
                 
                   
                     ∂ 
                     k 
                   
                   
                     S 
                     
                       ( 
                       
                         i 
                         , 
                         j 
                       
                       ) 
                     
                   
                 
                 
                   ∂ 
                   
                     ω 
                     p 
                     k 
                   
                 
               
               . 
             
           
         
       
     
     Single-Pulse Stabilization 
     In some embodiments, an identical pulse is applied to both the i-th qubit and the j-th qubit (i.e., the coefficient vectors {right arrow over (Λ)} (i)  and {right arrow over (Λ)} (j)  are the same) to perform an XX gate operation between the i-th qubit and the j-th qubit. Then, the matrix R p   (i,j)(k)  and the coefficient vector {right arrow over (Λ)} (i) (={right arrow over (Λ)} (j) ) may be spectrally decomposed using the v-th eigenvalue λ v,p   (i,j)(k)  of R p   (i,j)(k)  and the corresponding eigenvector |λ v,p   (i,j)(k)   ,as R p   (i,j)(k) =Σ v=1   N     0   λ v,p   (i,j)(k)     λ v,p   (i,j)(k) |,and {right arrow over (Λ)} (i) ={right arrow over (Λ)} (j) =Σ v=1   N     0   c v,p   (i,j)(k) |λ v,p    (k)   , respectively. In this spectral decomposition, the condition for the stabilization of a power-optimal pulse may be written as χ ij,p   (k) =Σ v=1   N     0   |c v,p   (i,j)(k) | 2 λ v,p   (i,j)(k) =0. Thus, components in the null-space expansion {right arrow over (g (l) )}=Σ α=1   N     0   Λ α   (l) {right arrow over (g)} α   (l)  of a power-optimal pulse, which has a large eigenvalue λ v,p   (i,j)(k)  un-stabilize against fluctuations in the frequencies ω p  of the motional modes. 
     To stabilize the entanglement interaction against fluctuations Δω p  in the frequencies ω p  of the motional modes in the first order to a certain extent, the components in the power-optimal pulse {right arrow over (g (l) )}=Σ α=1   N     0   Λ α   (l) {right arrow over (g)} α   (l)  that have large eigenvalues λ v,p   (i,j)(k)  are removed, and the resulting pulse is renormalized to obtain the entanglement condition χ ij =π/2. 
     Two-Pulse Stabilization 
     In some embodiments, a pulse applied to the i-th ion is chosen that is different from a pulse applied to the j-th ion such that the coefficient vectors {right arrow over (Λ)} (i)  and {right arrow over (Λ)} (j)  are orthogonal. Specifically, a pulse having the coefficient vector {right arrow over (Λ)} (i)  to be applied to the i-th qubit to perform is chosen. Using the chosen coefficient vector {right arrow over (Λ)} (i) , the coefficient vector {right arrow over (Λ)} (j)  of a pulse to be applied to the j-th qubit is determined such that coefficient vector {right arrow over (Λ)} (j)  satisfies the condition for the stabilization of the entanglement interaction χ ij,p   (k) ={right arrow over (Λ)} (i)     T   R p   (i,j)(k) {right arrow over (Λ)} (j) =0(k=1, 2, . . . K χ ) as described above and the pulse is power- optimal. In some embodiments, the coefficient vectors {right arrow over (Λ)} (i)  and {right arrow over (Λ)} (j)  are chosen such that each of the coefficient vectors {right arrow over (Λ)} (i)  and {right arrow over (Λ)} (j)  satisfies the condition for the stabilization of the entanglement interaction in the single-pulse stabilization described above. 
     Fluctuations in intensities of the laser beams and Lamb-Dicke parameters η i,p  may also affect the value of the entanglement interaction χ ij (τ), since the entanglement interaction χ ij (τ) is related to the amplitude A n (n=1, 2, . . . , N A ) of each tone. That is, the resulting entanglement interaction χ ij (τ) may have a different value from the desired value θ ij  that is set in the condition for non-zero entanglement interaction. Thus, in some embodiments, known broadband laser pulse sequences, typically applicable for single-qubit gate operations, such as Solovay-Kitaev (SK) sequence and Suzuki-Trotter sequences, can be used to mitigate errors in the entanglement interaction χ ij (τ) with respect to, for instance, offsets in Lamb-Dicke parameters η i,p . The same technique may be used to stabilize the entanglement interaction χ ij (τ) against any error sources that perturbs the entanglement interaction χ ij (τ) value. 
     Alternatively or additionally to the stabilization of the entanglement interaction χ ij (τ), the resulting entanglement interaction χ ij (τ) can be calibrated to the desired value θ ij  by modifying the amplitude Ω(t) of the pulse. 
     Demodulation of Pulses 
     To apply the power-optimal and error-resilient pulse on the l-th qubit (l=i,j), the amplitude and the detuning frequency (Ω (l) (t), μ (l) (t)) of the power-optimal pulse need to be demodulated (i.e., the amplitude and the detuning frequency (Ω (l) (t), μ (l) (t)) are extracted and the pulse function g (l) (t) is converted into a pulse having a series of time-dependent pulse segments of a single laser beam) from the determined pulse function g (l) =Ω (l) (t)sin(ψ (l) (t))(l=i,j) where ψ (l) (t)=∫ 0   t μ (l) (t′)dt′ is the phase accumulated due to the detuning frequency μ (l) (t). If this demodulation process is performed with a fixed detuning frequency, i.e., μ (l) (t)=μ 0 , the resulting pulse is amplitude-modulated (AM) pulse, in which the amplitude Ω (l) (t) is modulated. If the demodulation process is performed with a fixed amplitude, i.e., Ω (l) (t)=Ω 0 , the resulting pulse is a phase-modulated (PM) pulse, in which the phase ψ (l) (t) is modulated. If the phase ψ (l) (t) is implemented by modulating detuning frequency μ (l) (t), the resulting pulse is a frequency-modulated (FM) pulse. The demodulation process can be performed in any combined modulation of amplitude Ω (l) (t), phase ψ (l) (t) (thereby the detuning frequency μ (l) (t)), and frequency to construct a power-optimal pulse by conventional demodulation methods known in the art of signal processing. 
     The first step of an exemplary demodulation process is to find zeros of the pulse function g (t)=Ω (l) (t)sin(ψ (l) (t)) at t=ζ j (j=0,1, . . . , N Z −1) (i.e. ,g(ζ j )=0). Here, N Z  is a total number of zeros of the pulse function g (l) (t). The amplitude Ω (l) (t) can be chosen such that the amplitude Ω (l) (t) does not have zeros. Thus, the pulse function g (l) (t) zero when sin(ψ (l) (t)) is zero (i.e., sin(ψ (l) (ζ j ))=0). Due to the nature of the sine function, sin (ψ (l) (ζ j ))=0 when ψ (l) (ζ j )=jπ(j=0,1, . . . , N Z −1) , including the zeros at the beginning and the end of the gate duration T of the pulse (i.e. t=ζ 0 =0 and t=ζ N     Z     −1 =τ). 
     The second step of the demodulation process is to compute the detuning frequency μ (l) (t) based on the zeros of the pulse function g (l) (t). In some embodiments, the detuning frequency μ (l) (t) is approximated as a constant value between adjacent zeros of the pulse function g (l) (t) (i.e., μ (l) (t)≈μ j for ζ j−1 &lt;t&lt;ζ j, j= 1, 2, . . . , N Z −1). Since the phase ψ (l) (t) is accumulated due to the detuning frequency μ (l) (t) as in ψ (l) (t)=∫ 0   t μ (l) (t′)dt′, the difference between the phase at t=ζ j  and t=ζ j−1  is ψ (l) (ζ j )−ψ (l) (ζ j−1 )=∫ ζ     j−1     ζ     j   μ (l) (t′)dt=μ j (ζ j −ζ j−1 )=π. Consequently, the detuning frequency μ j between t=ζ j−1  and t=ζ j  is determined as μ j =π/(ζ j −ζ j−1 ). The third step of the demodulation process is compute the amplitude Ω (l) (t). A time derivative of the pulse function g (l) (t)=Ω (l) (t)sin(ψ (l) (t)) at t=ζ j  is
 
( g   (l) )′(ζ j )=(Ω (l) )′(ζ j )sin (ψ (l) (ζ j )cos(ψ(ζ j ))(ψ (l) )′(ζ j )=(−1) j Ω (l) (ζ j )μ (l) (ζ j ),
 
where ψ (l) (ζ j )=jπand ψ (l) (t)=∫ 0   t μ (l) (t′)dt′ are used. Thus, the amplitude Ω (l) (t) at t=ζ j  is computed as Ω (l) (ζ j )=(−1) j (g (l) )′(ζ j )/μ (l) (ζ j ) using the time derivative of the computed pulse function g (l) (t)=Σ n=1   N     A   g k   (l) sin(2πnd/τ) (i.e.,
 
     
       
         
           
             
               
                 
                   
                     
                       
                         
                           
                             ( 
                             
                               g 
                               
                                 ( 
                                 l 
                                 ) 
                               
                             
                             ) 
                           
                           ′ 
                         
                         ⁢ 
                         
                           ( 
                           t 
                           ) 
                         
                       
                       = 
                       
                         
                           ∑ 
                           
                             n 
                             = 
                             1 
                           
                           
                             N 
                             A 
                           
                         
                         
                           
                             g 
                             n 
                             
                               ( 
                               l 
                               ) 
                             
                           
                           ⁢ 
                           
                             
                               2 
                               ⁢ 
                               π 
                               ⁢ 
                               n 
                             
                             τ 
                           
                           ⁢ 
                           cos 
                           ⁢ 
                           
                             ( 
                             
                               2 
                               ⁢ 
                               π 
                               ⁢ 
                                  
                               n 
                               ⁢ 
                               
                                 t 
                                 / 
                                 τ 
                               
                             
                             ) 
                           
                         
                       
                     
                     ) 
                   
                   . 
                 
               
               
                   
               
             
           
         
       
     
     In some embodiments, a set of the computed detuning frequencies μ j  (j=1, 2, . . . , N Z −1) is interpolated with splines (e.g., functions defined piecewise by one or more polynomials or other algebraic expressions) and the interpolated values of the detuning frequency μ (l) (t) is used for μ (l) (ζ j ) for computing the amplitude Ω (l) (ζ j ). In some embodiments, μ (l) (ζ j ) is (i) μ j , (ii) μ j+1 , or (iii) (μ j +μ j+1 )/2 is used as μ (l) (ζ j ) for computing the amplitude Ω (l) (ζ j ). 
     In some embodiments, a set of the computed amplitude Ω (l) (ζ j ) is also interpolated with splines to compute the time-dependent amplitude Ω (l) (t). 
     If the demodulation process for a phase-modulated (PM) pulse, a set of the computed phase ψ (l) (ζ j ) may be interpolated with splines to compute the time-dependent phase ψ (l) (t). 
     EXAMPLES 
     In the following, examples of stabilization against fluctuations Δω in the frequencies ω p  of the motional modes according to the embodiments described above are shown. 
     In  FIG.  7 A , infidelity 1−F of an XX gate with k-th order stabilization (k=0, 1, . . . , 8) is shown as a function of fluctuation Δω in the frequencies ω p  of the motional modes, where the frequencies ω p  of the p-th motional mode (p=1, 2, . . . , 5) are equally shifted to ω p +Δω. In this example, the gate duration τ is chosen as 300 μs over a chain  102  of five trapped ions. 
     The infidelity 1−F of the XX gate operation without the stabilization  702  increases substantially as a fluctuation Δω in the frequencies ω p  of the motional modes increases, as the phase space trajectories α p   (l) (τ) are related to the frequencies ω p  of the p-th motional mode and thus sensitive to a fluctuation Δω in the frequencies ω p . Thus, the pulse determined using inaccurate frequencies ω p  (i.e., the frequencies that are different from the actual frequencies) causes an XX gate operation to be inaccurate (i.e., different from the desired XX gate operation). The infidelity 1−F of an XX gate operation in which the first-order stabilization  704  is applied will stay low for small fluctuation Δω of the frequencies ω p  (up to about 0.5 kHz). That is, the XX gate operation can be robust against fluctuations in the frequency of the motional mode with the first-order stabilization. The widths of the infidelity curves in  FIG.  7 A  extracted at infidelity of 0.001 (i.e., the range of fluctuations Δω in which the infidelity 1−F can be made smaller than 0.001) increase from about 0.1 kHz to about 13 kHz as the order of stabilization k increases from 0 to 8 as shown in  FIG.  7 B . 
     However, with the stabilization, the required laser power increases, for example, up to 40% in the example described herein. Hence, there is a trade-off between the degree of stabilization and the optimization of the required laser power. The power requirement of the stabilized pulses for each order of stabilization k the requirement scales linearly in the order of stabilization k as shown in  FIG.  7 C . The effect of the stabilization increases inversely proportional to the gate duration τ as shown in  FIG.  7 D . 
     As described above, methods of generating power-optimum laser pulses to perform an entangling gate operation between two qubits such that the entangling gate operation can be stabilized against external errors, leading to an improved fidelity of the entangling gate operation. 
     Additionally, determining control parameters in generating power optimum and error-resilient laser pulses to perform an entangling gate operation includes solving a set of linear equations. Thus, determining the control parameters and subsequently constructing a power-optimal and error-resilient laser pulse can be performed in an efficient manner. At the end of running the desired quantum algorithm, the population of the qubit states (trapped ions) within the quantum register is measured (read-out), so that the results of the quantum computation(s) with the desired quantum algorithm can be determined and provided to the classical computer for use to obtain solutions to the problems that may be intractable by the classical computer. 
     While the foregoing is directed to specific embodiments, other and further embodiments may be devised without departing from the basic scope thereof, and the scope thereof is determined by the claims that follow.