Fast separation of two trapped ions

Trapped cold ions provide a leading platform to implement quantum information processing. Separating ion chains is in the toolkit of basic operations required. (Merging chains is the corresponding reverse operation so we shall only refer to separation hereafter.) It has been used to implement two-qubit quantum gates [1]; also to purify entangled states [2, 3], or teleport material qubits [4]. Moreover, an architecture for processing information scalable to many ions could be developed based on shuttling, separating and merging ion crystals in multisegmented traps [5].

Ion-chain separation is known to be a difficult operation [6]. Experiments have progressed towards lower final excitations and shorter times but much room for improvement still remains [7–9]. Problems identified include anomalous heating, so devising short-time protocols via shortcuts to adiabaticity (STA) techniques was proposed as a way-out worth exploring [10]. STA methods intend to speed up different adiabatic operations [11, 12] without inducing final excitations. An example of an elementary (fast quasi-adiabatic) STA approach [11] was already applied for fast chain splitting in [7]. Here, we design, using a more general and efficient STA approach based on dynamical normal modes (NM) [13, 14], protocols to effectively separate two equal ions, initially in a common electrostatic linear harmonic trap, into a final configuration where each ion is in a different well. The motion is assumed to be effectively one-dimensional due to tight radial confinement. The external potential for an ion at q is approximated as

$$\begin{eqnarray}&&{V}_{\mathrm{ext}}=\alpha (t){q}^{2}+\beta (t){q}^{4},\end{eqnarray} \tag{ 1 }$$

which is experimentally realizable with state-of-art segmented Paul traps [6, 15].

Using dynamical NM [13, 14] a Hamiltonian will be set which is separable in a harmonic approximation around potential minima. By means of Lewis–Riesenfeld invariants [16] we shall design first the approximate dynamics of an unexcited splitting, taking into account anharmonicities in a perturbative manner, and from that inversely find the corresponding protocol, i.e. the α(t) and β(t) functions.

The Hamiltonian of the system of two ions of mass m and charge e is, in the laboratory frame

$$\begin{eqnarray}H & = & \displaystyle \frac{{p}_{1}^{2}}{2m}+\displaystyle \frac{{p}_{2}^{2}}{2m}+V,\\ V & = & \alpha (t)({q}_{1}^{2}+{q}_{2}^{2})+\beta (t)({q}_{1}^{4}+{q}_{2}^{4})+\displaystyle \frac{{C}_{c}}{{q}_{1}-{q}_{2}},\end{eqnarray} \tag{ 2 }$$

where p₁, p₂ are the momentum operators for both ions, q₁, q₂ their position operators, and ${C}_{c}=\frac{{e}^{2}}{4\pi {\epsilon }_{0}}$, ${\epsilon }_{0}$ being the vacuum permittivity. We use on purpose a c-number notation since we shall also consider classical simulations. The context will make clear if c-numbers or q-numbers are required. We suppose that, due to the strong Coulomb repulsion, ${q}_{1}\gt {q}_{2}$. By minimizing the potential part of the Hamiltonian V, we find for the equilibrium distance between the two ions, d(t), the quintic equation [6]

$$\begin{eqnarray}&&\beta (t){d}^{5}(t)+2\alpha (t){d}^{3}(t)-2{C}_{c}=0,\end{eqnarray} \tag{ 3 }$$

which will be quite useful for inverse engineering the ion-chain splitting, even without an explicit solution for d(t). At t = 0 a single external well is assumed, β(0) = 0 and $\alpha (0)\gt 0$, whereas in the final double-well configuration $\beta ({t}_{{\rm{f}}})\gt 0,\alpha ({t}_{{\rm{f}}})\lt 0$. At some intermediate time t_i the potential becomes purely quartic ($\alpha ({t}_{i})=0$). Our aim is to design the functions α(t) and β(t) so that each of the ions ends up in a different external well as shown in figure 1, in times as short as possible, and without any final excitation.

Zoom In Zoom Out Reset image size

To define dynamical NM coordinates, we calculate first at equilibrium (the point $\{{q}_{1}^{(0)},{q}_{2}^{(0)}\}$ in configuration space where the potential is a minimum, $\partial V/\partial {q}_{1}=\partial V/\partial {q}_{2}=0$) the matrix ${V}_{{ij}}=\frac{1}{m}\frac{{\partial }^{2}V}{\partial {q}_{i}\partial {q}_{j}}{| }_{{}_{\mathrm{eq}}}$.

The equilibrium positions are ${q}_{1}^{(0)}=\frac{d(t)}{2}$, ${q}_{2}^{(0)}=-\frac{{\rm{d}}(t)}{2}$, and the matrix takes the form

$$\begin{eqnarray}{V}_{{ij}}=\displaystyle \frac{1}{m}\left(\begin{array}{cc}2\alpha +12\beta \displaystyle \frac{{d}^{2}}{4}+\displaystyle \frac{2{C}_{c}}{{d}^{3}} & -\displaystyle \frac{2{C}_{c}}{{d}^{3}}\\ \displaystyle \frac{-2{C}_{c}}{{d}^{3}} & 2\alpha +12\beta \displaystyle \frac{{d}^{2}}{4}+\displaystyle \frac{2{C}_{c}}{{d}^{3}}\end{array}\right).\end{eqnarray} \tag{ 4 }$$

The eigenvalues are

$$\begin{eqnarray}{\lambda }_{-} & = & \displaystyle \frac{1}{m}(2\alpha +3\beta {d}^{2}),\\ {\lambda }_{+} & = & \displaystyle \frac{1}{m}\left(2\alpha +3\beta {d}^{2}+\displaystyle \frac{4{C}_{c}}{{d}^{3}}\right),\end{eqnarray} \tag{ 5 }$$

which define the NM frequencies as ${{\rm{\Omega }}}_{\pm }=\sqrt{{\lambda }_{\pm }}$ corresponding to center-of-mass (−) and relative (stretch) motions (+). These relations, with equation (3) written as

$$\begin{eqnarray}&&\beta (t)=\displaystyle \frac{2{C}_{c}}{{d}^{5}(t)}-\displaystyle \frac{2\alpha (t)}{{d}^{2}(t)},\end{eqnarray} \tag{ 6 }$$

allow us to write α(t) and d(t) as functions of the NM frequencies

$$\begin{eqnarray}&&\alpha (t)=\displaystyle \frac{1}{8}m\left(3{{\rm{\Omega }}}_{+}^{2}-5{{\rm{\Omega }}}_{-}^{2}\right),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\rm{d}}(t)=\sqrt[3]{\displaystyle \frac{4{C}_{c}}{m\left({{\rm{\Omega }}}_{+}^{2}-{{\rm{\Omega }}}_{-}^{2}\right)}}.\end{eqnarray} \tag{ 8 }$$

Substituting these expressions into equation (6), β(t) may also be written in terms of NM frequencies.

The normalized eigenvectors are

$$\begin{eqnarray}{v}_{-} & = & \displaystyle \frac{1}{\sqrt{2}}\displaystyle \left(\displaystyle \genfrac{}{}{0em}{}{1}{1}\right),\\ {v}_{+} & = & \displaystyle \frac{1}{\sqrt{2}}\displaystyle \left(\displaystyle \genfrac{}{}{0em}{}{1}{-1}\right),\end{eqnarray} \tag{ 9 }$$

which we denote as ${v}_{\pm }=\displaystyle \left(\displaystyle \genfrac{}{}{0em}{}{{a}_{\pm }}{{b}_{\pm }}\right)$. The (mass weighted) dynamical NM coordinates are defined in terms of the laboratory coordinates as

$$\begin{eqnarray}&&{{\mathsf{q}}}_{\pm }={a}_{\pm }\sqrt{m}({q}_{1}-{q}_{1}^{(0)})+{b}_{\pm }\sqrt{m}({q}_{2}-{q}_{2}^{(0)}).\end{eqnarray} \tag{ 10 }$$

The unitary transformation-of-coordinates matrix is

$$\begin{eqnarray}&&U=\langle {{\mathsf{q}}}_{+},{{\mathsf{q}}}_{-}| {q}_{1},{q}_{2}\rangle =\delta [{q}_{1}-{q}_{1}({{\mathsf{q}}}_{+},{{\mathsf{q}}}_{-})]\delta [{q}_{2}-{q}_{2}({{\mathsf{q}}}_{+},{{\mathsf{q}}}_{-})].\end{eqnarray} \tag{ 11 }$$

The Hamiltonian in the dynamical equation for $| \psi ^{\prime} \rangle =U| \psi \rangle$, where $| \psi \rangle$ is the lab-frame time-dependent wave function (in q₁, q₂ coordinates) evolving with H, is given by

$$\begin{eqnarray}H^{\prime} & = & {{UHU}}^{\dagger }-{\rm{i}}{\rm{\hslash }}U({\partial }_{t}{U}^{\dagger })\\ & = & \displaystyle \frac{{{\mathsf{p}}}_{+}^{2}}{2}+\displaystyle \frac{1}{2}{{\rm{\Omega }}}_{+}^{2}{{\mathsf{q}}}_{+}^{2}+\displaystyle \frac{\dot{d}}{\sqrt{2}}\sqrt{m}{{\mathsf{p}}}_{+}\\ & & +\displaystyle \frac{{{\mathsf{p}}}_{-}^{2}}{2}+\displaystyle \frac{1}{2}{{\rm{\Omega }}}_{-}^{2}{{\mathsf{q}}}_{-}^{2},\end{eqnarray} \tag{ 12 }$$

plus qubic and higher order terms in the potential that we neglect by now (they will be considered in section 4 below). Similarly to [13, 14], we apply a further unitary transformation ${\mathcal{U}}={{\rm{e}}}^{-{\rm{i}}\sqrt{m}\dot{d}{{\mathsf{q}}}_{+}/(\sqrt{2}{\rm{\hslash }})}$, to write down an effective Hamiltonian for $| {\psi }^{\prime\prime }\rangle ={\mathcal{U}}| \psi ^{\prime} \rangle$ with the form of two independent harmonic oscillators in NM space, $H^{\prime \prime} ={\mathcal{U}}H^{\prime} {{\mathcal{U}}}^{\dagger }-{\rm{i}}{\rm{\hslash }}{\mathcal{U}}({\partial }_{t}{{\mathcal{U}}}^{\dagger })$,

$$\begin{eqnarray}H^{\prime \prime} & = & \displaystyle \frac{{{\mathsf{p}}}_{+}^{2}}{2}+\displaystyle \frac{1}{2}{{\rm{\Omega }}}_{+}^{2}{\left({{\mathsf{q}}}_{+}+\displaystyle \frac{\sqrt{m}\ddot{d}}{\sqrt{2}{{\rm{\Omega }}}_{+}^{2}}\right)}^{2}\\ & & +\displaystyle \frac{{{\mathsf{p}}}_{-}^{2}}{2}+\displaystyle \frac{1}{2}{{\rm{\Omega }}}_{-}^{2}{{\mathsf{q}}}_{-}^{2}={H}_{+}^{\prime\prime }+{H}_{-}^{\prime\prime }.\end{eqnarray} \tag{ 13 }$$

These oscillators have dynamical invariants of the form [16]

$$\begin{eqnarray}{I}_{\pm } & = & \displaystyle \frac{1}{2}{[{\rho }_{\pm }({{\mathsf{p}}}_{\pm }-{\dot{x}}_{\pm })-{\dot{\rho }}_{\pm }({{\mathsf{q}}}_{\pm }-{x}_{\pm })]}^{2}\\ & & +\displaystyle \frac{1}{2}{{\rm{\Omega }}}_{0\pm }^{2}{\left(\displaystyle \frac{{{\mathsf{q}}}_{\pm }-{x}_{\pm }}{{\rho }_{\pm }}\right)}^{2},\end{eqnarray} \tag{ 14 }$$

where the auxiliary functions ${\rho }_{\pm }$ and ${x}_{+}$ satisfy

$$\begin{eqnarray}&&{\ddot{\rho }}_{\pm }+{{\rm{\Omega }}}_{\pm }^{2}{\rho }_{\pm }=\displaystyle \frac{{{\rm{\Omega }}}_{0\pm }^{2}}{{\rho }_{\pm }^{3}},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\ddot{x}}_{+}+{{\rm{\Omega }}}_{+}^{2}{x}_{+}=-\sqrt{\displaystyle \frac{m}{2}}\ddot{d},\end{eqnarray} \tag{ 16 }$$

with ${{\rm{\Omega }}}_{0\pm }={{\rm{\Omega }}}_{\pm }(0)$, and, due to symmetry, ${x}_{-}=0$.

The physical meaning of the auxiliary functions may be grasped from the solutions of the time-dependent Schrödinger equations (for each NM Hamiltonian ${H}_{\pm }^{\prime\prime }$ in equation (13)) proportional to the invariant eigenvectors [22]. They form a complete basis for the space spanned by each Hamiltonian ${H}_{\pm }^{\prime\prime }$ and take the form

$$\begin{eqnarray}&&| {\psi }_{n\pm }^{\prime\prime }(t)\rangle ={{\rm{e}}}^{\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}\left[[,\displaystyle \frac{{\dot{\rho }}_{\pm }{{\mathsf{q}}}_{\pm }^{2}}{2{\rho }_{\pm }}+({\dot{x}}_{\pm }{\rho }_{\pm }-{x}_{\pm }{\dot{\rho }}_{\pm })\displaystyle \frac{{{\mathsf{q}}}_{\pm }}{{\rho }_{\pm }},]\right]}\displaystyle \frac{{{\rm{\Phi }}}_{n}({\sigma }_{\pm })}{{\rho }_{\pm }^{1/2}},\end{eqnarray} \tag{ 17 }$$

where ${\sigma }_{\pm }=\frac{{{\mathsf{q}}}_{\pm }-{x}_{\pm }}{{\rho }_{\pm }}$, and ${{\rm{\Phi }}}_{n}({\sigma }_{\pm })$ are the eigenfunctions of the static harmonic oscillator at time t = 0. Thus ${\rho }_{\pm }$ are scaling factors proportional to the state 'width' in NM coordinates, whereas the ${x}_{\pm }$ are the dynamical-mode centers in the space of NM coordinates. Within the harmonic approximation there are dynamical states of the factorized form $| {\psi }^{\prime\prime }(t)\rangle =| {\psi }_{+}^{\prime\prime }(t)\rangle | {\psi }_{-}^{\prime\prime }(t)\rangle$ for the ion chain dynamics, where the NM wave functions $| {\psi }_{\pm }^{\prime\prime }(t)\rangle$ evolve independently with ${H}_{\pm }^{\prime\prime }$. They may be written as combinations of the form $| {\psi }_{\pm }^{\prime\prime }(t)\rangle ={\displaystyle \sum }_{n}{C}_{n\pm }| {\psi }_{n\pm }^{\prime\prime }(t)\rangle$, with constant amplitudes ${C}_{n\pm }$. The average energies of the nth basis states for the two NM are ${E}_{n\pm }^{\prime\prime }=\langle {\psi }_{n\pm }^{\prime\prime }| {H}_{\pm }^{\prime\prime }| {\psi }_{n\pm }^{\prime\prime }\rangle$,

$$\begin{eqnarray}{E}_{n-}^{\prime\prime } & = & \displaystyle \frac{(2n+1){\rm{\hslash }}}{4{{\rm{\Omega }}}_{0-}}\left({\dot{\rho }}_{-}^{2}+{{\rm{\Omega }}}_{-}^{2}{\rho }_{-}^{2}+\displaystyle \frac{{{\rm{\Omega }}}_{0-}^{2}}{{\rho }_{-}^{2}}\right),\\ {E}_{n+}^{\prime\prime } & = & \displaystyle \frac{(2n+1){\rm{\hslash }}}{4{{\rm{\Omega }}}_{0+}}\left({\dot{\rho }}_{+}^{2}+{{\rm{\Omega }}}_{+}^{2}{\rho }_{+}^{2}+\displaystyle \frac{{{\rm{\Omega }}}_{0+}^{2}}{{\rho }_{+}^{2}}\right)\\ & & +\displaystyle \frac{1}{2}{\dot{x}}_{+}^{2}+\displaystyle \frac{1}{2}{{\rm{\Omega }}}_{+}^{2}{\left({x}_{+}-\displaystyle \frac{\sqrt{m}\ddot{d}}{\sqrt{2}{{\rm{\Omega }}}_{+}^{2}}\right)}^{2}.\end{eqnarray} \tag{ 18 }$$

Once the Hamiltonian and Lewis–Riesenfeld invariants are defined, we proceed to apply the invariant based inverse engineering technique and design STA. The results for the simple harmonic oscillator in [11] serve as a reference but have to be extended here since the two modes are not really independent from the perspective of the inverse problem. This is because a unique protocol, i.e., a single set of α(t) and β(t) functions has to be designed.

We first set the initial and target values for the control parameters α(t) and β(t). At time t = 0, the external trap—for a single ion—is purely harmonic, with (angular) frequency ${\omega }_{0}$. From equation (5), we find that ${{\rm{\Omega }}}_{-}(0)={\omega }_{0}$ and ${{\rm{\Omega }}}_{+}(0)=\sqrt{3}{\omega }_{0}$. The equilibrium distance is $d(0)=\sqrt[3]{\frac{2{C}_{c}}{m{\omega }_{0}^{2}}}$. For the final time, we set a tenfold expansion of the equilibrium distance, ${\rm{d}}({\text{}}{t}_{{\rm{f}}})=10d(0)$, and ${{\rm{\Omega }}}_{-}({t}_{{\rm{f}}})={\omega }_{0}$. This also implies ${{\rm{\Omega }}}_{+}({t}_{{\rm{f}}})=\sqrt{1.002}\;{\omega }_{0}\approx {{\rm{\Omega }}}_{-}({t}_{{\rm{f}}})$, i.e., the final frequencies of both NM are essentially equal, the Coulomb interaction is negligible, and the ions can be considered to oscillate in independent traps.

The inverse problem is somewhat similar to the expansion of a trap with two equal ions in [14], but complicated by the richer structure of the external potential. The Hamiltonian (13) and the invariant (14) must commute at both boundary times $[H({t}_{{\rm{b}}}),I({t}_{{\rm{b}}})]=0$,

$$\begin{eqnarray}&&{t}_{{\rm{b}}}=0,{t}_{{\rm{f}}},\end{eqnarray} \tag{ 19 }$$

to drive initial levels into final levels via a one-to-one mapping. This is achieved by applying appropriate boundary conditions (BCs) to the auxiliary functions ${\rho }_{\pm }$, ${x}_{\pm }$ and their derivatives

$$\begin{eqnarray}&&{\rho }_{\pm }(0)=1,\ \ {\rho }_{\pm }({t}_{{\rm{f}}})={\gamma }_{\pm },\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{\dot{\rho }}_{\pm }({t}_{{\rm{b}}})={\ddot{\rho }}_{\pm }({t}_{{\rm{b}}})=0,\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{x}_{+}({t}_{{\rm{b}}})={\dot{x}}_{+}({t}_{{\rm{b}}})={\ddot{x}}_{+}({t}_{{\rm{b}}})=0,\end{eqnarray} \tag{ 22 }$$

where ${\gamma }_{\pm }=\sqrt{\frac{{{\rm{\Omega }}}_{\pm }(0)}{{{\rm{\Omega }}}_{\pm }({t}_{{\rm{f}}})}}$. Let us recall that ${x}_{-}=0$ for all times so this parameter does not have to be considered further.

Inserting the BC for ${x}_{+}({t}_{{\rm{b}}})$ and ${\ddot{x}}_{+}({t}_{{\rm{b}}})$ in equation (16) we find that $\ddot{d}({t}_{{\rm{b}}})=0$. Additionally, $\dot{d}({t}_{{\rm{b}}})=0$ is to be imposed so that ${\mathcal{U}}({t}_{{\rm{b}}})=1$. According to equation (8), $\dot{d}({t}_{{\rm{b}}})=\ddot{d}({t}_{{\rm{b}}})=0$ by imposing ${\dot{{\rm{\Omega }}}}_{\pm }({t}_{{\rm{b}}})={\ddot{{\rm{\Omega }}}}_{\pm }({t}_{{\rm{b}}})=0$. With $\dot{d}({t}_{{\rm{b}}})=\ddot{d}({t}_{{\rm{b}}})=0$ the Hamiltonians and wave functions coincide at the boundaries, $H^{\prime} ({t}_{{\rm{b}}})=H^{\prime \prime} ({t}_{{\rm{b}}})$, $| {\psi }^{\prime }({t}_{{\rm{b}}})\rangle =| {\psi }^{\prime\prime }({t}_{{\rm{b}}})\rangle$, which simplifies the calculation of the excitation energy.

From equation (15), the NM frequencies may be written as

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\pm }=\sqrt{\displaystyle \frac{{{\rm{\Omega }}}_{0\pm }^{2}}{{\rho }_{\pm }^{4}}-\displaystyle \frac{{\ddot{\rho }}_{\pm }}{{\rho }_{\pm }}}.\end{eqnarray} \tag{ 23 }$$

Thus the BC ${\dot{{\rm{\Omega }}}}_{\pm }({t}_{{\rm{b}}})={\ddot{{\rm{\Omega }}}}_{\pm }({t}_{{\rm{b}}})=0$ are satisfied by imposing on the auxiliary functions the additional BC

$$\begin{eqnarray}&&{\dddot{\rho }}_{\pm }({t}_{{\rm{b}}})=\;{\ddddot{\rho }}_{\pm }({t}_{{\rm{b}}})=0.\end{eqnarray} \tag{ 24 }$$

We may now design ansatzes for the auxiliary functions ${\rho }_{\pm }$ that satisfy the ten BC in equations (20), (21) and (24), plus the BC for ${x}_{+}({t}_{{\rm{b}}})$ and ${\dot{x}}_{+}({t}_{{\rm{b}}})$ in equation (22) (since $\ddot{d}({t}_{{\rm{b}}})=0$, ${\ddot{x}}_{+}({t}_{{\rm{b}}})=0$ is then automatically satisfied, see equation (16)). Finally, from the NM frequencies given by equation (23) we can inverse engineer the control parameters $\alpha (t)$ and $\beta (t)$ from equations (6), (7) and (8) .

A simple choice for ${\rho }_{-}(t)$ is a polynomial ansatz of 9th order ${\rho }_{-}={\displaystyle \sum }_{i=0}^{9}{b}_{i}{s}^{i}$, where $s=t/{t}_{{\rm{f}}}$. Substituting this form in the ten BC in equations (20), (21) and (24), we finally get

$$\begin{eqnarray}{\rho }_{-} & = & 126({\gamma }_{-}-1){s}^{5}-420({\gamma }_{-}-1){s}^{6}+540({\gamma }_{-}-1){s}^{7}\\ & & -315({\gamma }_{-}-1){s}^{8}+70({\gamma }_{-}-1){s}^{9}+1.\end{eqnarray} \tag{ 25 }$$

For ${\rho }_{+}$ we will use an 11th order polynomial ${\rho }_{+}={\displaystyle \sum }_{n=0}^{11}{a}_{n}{s}^{n}$ to satisfy as well ${x}_{+}({t}_{{\rm{b}}})={\dot{x}}_{+}({t}_{{\rm{b}}})=0$. The parameters ${a}_{0-9}$ are fixed so that the ten BC for ${\rho }_{+}$ are fulfilled (see the appendix), whereas ${a}_{10}$, ${a}_{11}$ are left free, and will be numerically determined by a shooting program [17] ('fminsearch' in MATLAB which uses the Nelder–Mead simplex method for optimization), so that the remaining BC for ${x}_{+}({t}_{{\rm{b}}})$ and ${\dot{x}}_{+}({t}_{{\rm{b}}})$ are also satisfied. Specifically, for each pair $\{{a}_{10},{a}_{11}\}$, ${{\rm{\Omega }}}_{\pm }(t)$ and d(t) are determined from equations (8) and (15), to solve equation (16) for ${x}_{+}(t)$ with initial conditions ${x}_{+}(0)={\dot{x}}_{+}(0)=0$. The free constants are changed until ${x}_{+}({t}_{{\rm{f}}})=0$ and ${\dot{x}}_{+}({t}_{{\rm{f}}})=0\;$ are satisfied. Numerically a convenient way to find the solution is to minimize the energy ${E}_{n+}^{\prime\prime }({t}_{{\rm{f}}})$ in equation (18).

Figures 2(a) and (b) depict the control parameters $\alpha (t)$ and $\beta (t)$ found with this method, using equations (6) and (7), for some value of t_f and ${\omega }_{0}$, see the caption, while figure 2(c) represents the equilibrium distance between ions as a function of time (8), and figure 2(d) the NM frequencies. In figure 3(a) the excitation energy is shown, versus final time, for optimized parameters given in figure 3(b). The initial state is the ground state of the two ions. It is calculated by propagating an initial guess of the wave function in imaginary time until it relaxes to the lowest eigenfunction [18]. The excitation energy is ${E}_{\mathrm{ex}}=E({t}_{{\rm{f}}})-{E}_{0}({t}_{{\rm{f}}})$, where $E({t}_{{\rm{f}}})$ is the final energy, calculated in the lab frame, and ${E}_{0}({t}_{{\rm{f}}})$ is the final ground-state energy. The wave function evolution is calculated using the 'split-operator method' with the full Hamiltonian (2). If the harmonic approximation were exact, there would not be any excitation with this STA method, $E({t}_{{\rm{f}}})={E}_{0+}^{\prime\prime }({t}_{{\rm{f}}})+{E}_{0-}^{\prime\prime }({t}_{{\rm{f}}})={E}_{0}({t}_{{\rm{f}}})$, see equation (18). The actual result is perturbed by the anharmonicities and NM couplings. The final ground state is also calculated with an 'imaginary-time evolution'. The corresponding final ground state energy is essentially two times a harmonic oscillator ground state energy plus the (negligible) Coulomb repulsion at distance ${\rm{d}}({t}_{{\rm{f}}})$. For the final times of all the examples, as it was noted in previous works [10, 13, 14, 19], classical simulations (solving Hamilton's equations from the equilibrium configuration instead of Schrödinger's equation) give indistinguishable results in the scale of figure 3(a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The excitation energy in figure 3(a) (solid line) increases at short times since the harmonic NM approximation fails [13, 14]. However, it goes down rapidly below one excitation quantum at times which are still rather small compared to experimental values used so far [7, 8]. In the following section we shall apply a perturbative technique to minimize the excitation further.

An improvement of the protocol is to consider the perturbation of the higher order terms neglected in the Hamiltonian (13). These 'anharmonicities' [20] are cubic and higher order terms in the Taylor expansion of the Coulomb term ${C}_{c}/({q}_{1}-{q}_{2})$,

$$\begin{eqnarray}\delta V & = & {\displaystyle \sum }_{j=3}^{\infty }\delta {V}^{(j)}\\ & = & {\displaystyle \sum }_{j=3}^{\infty }\displaystyle \frac{{(-1)}^{j}{C}_{c}}{{d}^{j+1}}{\left[\left({q}_{2}-{q}_{2}^{(0)}\right)-\left({q}_{1}-{q}_{1}^{(0)}\right)\right]}^{j}.\end{eqnarray} \tag{ 26 }$$

In NM coordinates the terms take the simple form

$$\begin{eqnarray}&&\delta {V}^{(j)}={(-1)}^{j+1}\displaystyle \frac{{C}_{c}}{{d}^{j+1}}{\left(\sqrt{\displaystyle \frac{2}{m}}{{\mathsf{q}}}_{+}\right)}^{j},\end{eqnarray} \tag{ 27 }$$

which may be regarded as a perturbation to be added to ${H}_{+}^{\prime\prime }$ in equation (13). (The perturbation does not couple the center-of-mass and relative subspaces.) To first order, the excess energy due to these perturbative terms at final time is given by

$$\begin{eqnarray}&&\delta {E}_{n+}^{(j)}=\langle {\psi }_{n+}^{\prime\prime }({t}_{{\rm{f}}})| \delta {V}^{(j)}| {\psi }_{n+}^{\prime\prime }({t}_{{\rm{f}}})\rangle ,\end{eqnarray} \tag{ 28 }$$

where the $| {\psi }_{n+}^{\prime\prime }\rangle$ are the unperturbed states in equation (17). Inverse engineering the splitting process may now be carried out by considering a 12th order polynomial for ${\rho }_{+}$ (see (A2)), with three free parameters so as to fix the BC for ${x}_{+}$ and also minimize the excitation energy. In practice we use MATLAB's 'fminsearch' function for the shooting to minimize ${E}_{0+}({t}_{{\rm{f}}})+\delta {E}_{0+}^{(3)}$ as no significant improvement occurs by including higher order terms. Figure 3(a) compares the performance of such a protocol with the simpler one with the 11th-order polynomial (A1). Figure 3(c) gives the values of optimized parameters at different final times.

A large quartic potential is desirable to control the excitations produced at the point where the harmonic term changes its sign [10]. At this point, the harmonic potential switches from confining to repulsive, which reduces the control of the system and potentially increases diabaticities and heating. In the inverse approach proposed here there is no special design of the protocol at this point, but the method naturally seeks high quartic confinements there. In figure 2(b) β reaches its maximum value right at the time where α changes sign (see figure 2(a)). However, the maximum value that β can reach will typically be limited in a Paul trap [6].

In table 1 we summarize the different maximal values of β and critical times (final times at which excitations below 0.1 quanta are reached) for different values of ${\omega }_{0}$ using the 11th order polynomial (A1) for ${\rho }_{+}$.

Table 1.  Maximum values of β, and critical times (final times at which excitations below 0.1 quanta are reached) for different values of ${\omega }_{0}$. The calculations were performed with the 11th order polynomial for ${\rho }_{+}$.

${\omega }_{0}$ (MHz)	${\beta }_{\mathrm{max}}$ (10−3N m−3)	tcrit (μs)
3	44.2	2.9
2	11.4	4.4
1.2	2.082	7.4
0.8	0.539	11.2

The maximum β decreases with t_f, such that the shortest possible t_f at a given maximum tolerable excitation energy is limited by the achievable β. The trap used in [8] yields a maximum β of about 10⁻⁴ N m⁻³, at ±10 V steering range. In a recent experiment reported in [23], the value used was β ≈ 5× 10⁻³ N m⁻³. The numbers reported in the table are thus within reach, as the β coefficients scale with the inverse 4th power of the overall trap dimension, and technological improvements on arbitrary waveform generators may allow for operation at an increased voltage range.

Another potential limitation the method could encounter in the laboratory is due to biases (a linear slope) in the trapping potential, ${V}_{\mathrm{ext}}=\alpha {q}^{2}+\beta {q}^{4}+\lambda q$, with λ constant and unknown [10]. Figure 4 represents the excitation energy versus the energy difference between the two final minima of the external potential, ΔE (also versus λ). To calculate the results, α(t) and β(t) are designed as if λ = 0, but the dynamics is carried out with a non-zero λ, in particular the initial state is the actual ground state, including the perturbation. Note that ΔE should be more than a thousand vibrational quanta to excite the final energy by one quantum. In [8] an energy increase of ten phonons at about 150 zN and 80 μs separation time was reported, so the STA ramps definitely improve the robustness against bias.

Zoom In Zoom Out Reset image size

Further experimental limitations may be due to random fluctuations in the potential parameters, or higher order terms in the external potential. We leave these important issues for a separate study but note that the structure of the STA techniques used here is well adapted to deal with noise or perturbations [24–26]. Combining STA with optimal-control-theory methods is also feasible, see e.g. [12].

Finally, we compare in figure 5 the performance of the protocols based on the polynomials (25) and (A1) with a simple non-optimized protocol based on those experimentally used in [8]. There, the equilibrium distance d is first designed as ${\rm{d}}(t)={\rm{d}}(0)+[{\rm{d}}({t}_{{\rm{f}}})-{\rm{d}}(0)]{s}^{2}\mathrm{sin}(s\pi /2)$, where $s=t/{t}_{{\rm{f}}}$. From the family of possible potential ramps consistent with this function, we chose a polynomial that drives α from $\alpha (0)={\alpha }_{0}$ to $\alpha ({t}_{{\rm{f}}})=-{\alpha }_{0}/2$ (as in figure 2) and whose first derivatives are 0 at both boundary times. β is given by equation (3). For the times analysed in figure 5, the method based on equations (25) and (A1) clearly outperforms the non-optimized ramp. To get excitations below the single motional excitation quantum with the non-optimized protocol, final times as long as ${t}_{{\rm{f}}}\sim 80\;\mu {\rm{s}}$ would be needed, which is in line with current experiments.

Zoom In Zoom Out Reset image size

We conclude that the method presented here, could bring a clear improvement with respect to the best results experimentally reported so far [7, 8]. The parameters required are realistic in current trapped ions laboratories. The simulations show that, under ideal conditions, the separation of two ions could be performed in a few oscillation periods, at times similar to those required for other operations as transport [13] or expansions [14], also studied with STA.

This work has been possible thanks to the close collaboration with three ion-trap groups at Boulder, Mainz and Zurich. We acknowledge in particular discussions with Didi Leibfried, Ludwig de Clercq, Joseba Alonso, and Jonathan Home that provided essential elements to develop the approach. It was supported by the Basque Country Government (Grant No. IT472-10), Ministerio de Economí a y Competitividad (Grant No. FIS2012-36673-C03-01), and the program UFI 11/55 of UPV/EHU. MP and SM-G acknowledge fellowships by UPV/EHU. This research was funded by the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), through the Army Research Office grant W911NF-10-1-0284. All statements of fact, opinion or conclusions contained herein are those of the authors and should not be construed as representing the official views or policies of IARPA, the ODNI, or the US Government.

The ansatz for ${\rho }_{+}$ that satisfies the BC ${\rho }_{+}(0)=1$, ${\rho }_{+}({t}_{{\rm{f}}})={\gamma }_{+}$, ${\dot{\rho }}_{+}({t}_{{\rm{b}}})={\ddot{\rho }}_{+}({t}_{{\rm{b}}})={\dddot{\rho }}_{+}({t}_{{\rm{b}}})={\ddddot{\rho }}_{+}({t}_{{\rm{b}}})=0$ with two free parameters takes the form

$$\begin{eqnarray}{\rho }_{+} & = & 1-(126-126{\gamma }_{+}+{a}_{10}+5{a}_{11}){s}^{5}\\ & & +(420-420{\gamma }_{+}+5{a}_{10}+24{a}_{11}){s}^{6}\\ & & -(540-540{\gamma }_{+}+10{a}_{10}+45{a}_{11}){s}^{7}\\ & & +(315-315{\gamma }_{+}+10{a}_{10}+40{a}_{11}){s}^{8}\\ & & -(70-70{\gamma }_{+}+5{a}_{10}+15{a}_{11}){s}^{9}\\ & & +{a}_{10}{s}^{10}+{a}_{11}{s}^{11}.\end{eqnarray} \tag{ A1 }$$

To minimize the perturbation energy in equation (28), three free parameters are introduced,

$$\begin{eqnarray}{\rho }_{+} & = & 1-(126-126{\gamma }_{+}+{c}_{10}+5{c}_{11}+15{c}_{12}){s}^{5}\\ & & +(420-420{\gamma }_{+}+5{c}_{10}+24{c}_{11}+70{c}_{12}){s}^{6}\\ & & -(540-540{\gamma }_{+}+10{c}_{10}+45{c}_{11}+126{c}_{12}){s}^{7}\\ & & +(315-315{\gamma }_{+}+10{c}_{10}+40{c}_{11}+105{c}_{12}){s}^{8}\\ & & -(70-70{\gamma }_{+}+5{c}_{10}+15{c}_{11}+35{c}_{12}){s}^{9}\\ & & +{c}_{10}{s}^{10}+{c}_{11}{s}^{11}+{c}_{12}{s}^{12}.\end{eqnarray} \tag{ A2 }$$