Fast separation of two trapped ions

We design fast protocols to separate or recombine two ions in a segmented Paul trap. By inverse engineering the time evolution of the trapping potential composed of a harmonic and a quartic term, it is possible to perform these processes in a few microseconds without final excitation. These times are much shorter than the ones reported so far experimentally. The design is based on dynamical invariants and dynamical normal modes. Anharmonicities beyond the harmonic approximation at potential minima are taken into account perturbatively. The stability versus an unknown potential bias is also studied.


I. INTRODUCTION
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][8][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 [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 which is experimentally realizable with state-of-art segmented Paul Traps [6,15]. Using dynamical normal modes (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, where p 1 , p 2 are the momentum operators for both ions, q 1 , q 2 their position operators, and C c = e 2 4πǫ0 , ǫ 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 > 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] 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 α(0) > 0, whereas in the final double-well configuration β(t f ) > 0, α(t f ) < 0. At some intermediate time t i the potential becomes purely quartic (α(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 Fig. 1, in times as short as possible, and without any final excitation.

II. DYNAMICAL NORMAL MODES
To define dynamical NM coordinates, we calculate first at equilibrium (the point {q , and the matrix takes the form The eigenvalues are which define the NM frequencies as Ω ± = λ ± corresponding to center-of-mass (−) and relative (stretch) motions (+). These relations, with Eq. (3) written as allow us to write α(t) and d(t) as functions of the NM frequencies, Substituting these expressions into Eq. (6), β(t) may also be written in terms of NM frequencies.
The normalized eigenvectors are which we denote as v ± = a± b± . The (mass weighted) dynamical NM coordinates are defined in terms of the laboratory coordinates as The unitary transformation of coordinates is where q + , q − |q 1 , The Hamiltonian in the dynamical equation for |ψ ′ = U |ψ , where |ψ is the lab-frame time-dependent wave function evolving with H, is given by plus qubic and higher order terms in the potential that we neglect by now (they will be considered in Sec. IV below). Similarly to [13,14], we apply a further unitary transformation U = e −i √ mḋq+/( √ 2 ) , to write down an effective Hamiltonian for |ψ ′′ = U|ψ ′ with the form of two independent harmonic oscillators in NM space, These oscillators have dynamical invariants of the form [16] where the auxiliary functions ρ ± and x + satisfyρ with Ω 0± = Ω ± (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 ′′ ± in Eq. (13)) proportional to the invariant eigenvectors [21]. They form a complete basis for the space spanned by each Hamiltonian H ′′ ± and take the form where σ ± = q±−x± ρ± , and Φ n (σ ± ) are the eigenfunctions of the static harmonic oscillator at time t = 0. Thus ρ ± are scaling factors proportional to the state "width" in NM coordinates, whereas the x ± are the dynamical-mode centers in the space of NM coordinates. Within the harmonic approximation there are dynamical states of the factorized form |ψ ′′ (t) = |ψ ′′ + (t) |ψ ′′ − (t) for the ion chain dynamics, where the NM wave functions |ψ ′′ ± (t) evolve independently with H ′′ ± . They may be written as combinations of the form |ψ ′′ ± (t) = n C n± |ψ ′′ n± (t) , with constant amplitudes C n± . The average energies of the n-th basis states for the two NM are

III. DESIGN OF THE CONTROL PARAMETERS
Once the Hamiltonian and Lewis-Riesenfeld invariants are defined, we proceed to apply the invariant based inverse engineering technique and design shortcuts to adiabaticity. 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 ω 0 . From Eq. (5), we find that Ω − (0) = ω 0 and Ω + (0) = √ 3ω 0 . The equilibrium distance is d(0) = 3 2Cc . For the final time, we set a tenfold expansion of the equilibrium distance, d(t f ) = 10d(0), and Ω − (t f ) = ω 0 . This also implies Ω + (t f ) = √ 1.002 ω 0 ≈ Ω − (t 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 Ref. [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 b ), I(t b )] = 0, to drive initial levels into final levels via a one-to-one mapping. This is achieved by applying appropriate boundary conditions (BC) to the auxiliary functions ρ ± , x ± and their derivatives, where γ ± = Ω±(0) Ω±(t 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 b ) andẍ + (t b ) in Eq. (16) we find thatd(t b ) = 0. Additionally,ḋ(t b ) = 0 is to be imposed so that U(t b ) = 1. According to Eq. From Eq. (15), the NM frequencies may be written as Thus the BCΩ ± (t b ) =Ω ± (t b ) = 0 are satisfied by imposing on the auxiliary functions the additional BC ...
A simple choice for ρ − (t) is a polynomial ansatz of 9th order ρ − = 9 i=0 b i s i , where s = t/t f . Substituting this form in the ten BC in Eqs. (20,21,24), we finally get For ρ + we will use an 11-th order polynomial ρ + = 11 n=0 a n s n to satisfy as well x + (t b ) =ẋ + (t b ) = 0. The parameters a 0−9 are fixed so that the 10 BC for ρ + 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 b ) andẋ + (t b ) are also satisfied. Specifically, for each pair {a 10 , a 11 }, Ω ± (t) and d(t) are determined from Eqs. (15) and (8), to solve Eq. (16) for x + (t) with initial conditions x + (0) =ẋ + (0) = 0. The free constants are changed until x + (t f ) = 0 andẋ + (t f ) = 0 are satisfied. Numerically a convenient way to find the solution is to minimize the energy E ′′ n+ (t f ) in Eq. (18). Fig. 2 (a) and (b) depict the control parameters α(t) and β(t) found with this method, using Eqs. (7) and (6), for some value of t f and ω 0 , see the caption, while Fig. 2 (c) represents the equilibrium distance between ions as a function of time (8), and Fig. 2 (d) the NM frequencies. In Fig. 3 (a) the excitation energy is shown, versus final time, for optimized parameters given in Fig. 3 (b). The initial state is the ground state of the two ions. It is calculated with an "imaginary time evolution". This is a variation of any numerical time evolution method (here we used a variation of the split-operator method), where instead of the real time one defines imaginary time in the evolution operator. By letting an ansatz wave function evolve for the static initial potential, it will eventually converge to the ground state of the system. The excitation energy is E ex = E(t f ) − E 0 (t f ), where E(t f ) is the final energy, calculated in the lab frame, and E 0 (t 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, (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 d(t f ). For the final times of all the examples, as it was noted in previous works [10,13,14,18], 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).
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.

IV. BEYOND THE HARMONIC APPROXIMATION
An improvement of the protocol is to consider the perturbation of the higher order terms neglected in the Hamiltonian (13). These "anharmonicities" [19] are cubic and higher order terms in the Taylor expansion of the Coulomb term C c /(q 1 − q 2 ), In NM coordinates the terms take the simple form which may be regarded as a perturbation to be added to H ′′ + in Eq. (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 where the |ψ ′′ n+ are the unperturbed states in Eq. (17). Inverse engineering the splitting process may now be carried out by considering a 12th order polynomial for ρ + (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 f ) + δE (3) 0+ as no significant improvement occurs by including higher order terms. Fig. 3 (a) compares the performance of such a protocol with the simpler one with the 11th-order polynomial (A1). Fig. 3 (c) gives the values of optimized parameters at different final times.

V. DISCUSSION
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 Fig. 2 (b) β reaches its maximum value right at the time where α changes sign (see Fig. 2 (a)). However, the maximum value that β can reach will typically be limited in a Paul trap [6]. In Table I 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 ω 0 using the 11-th order polynomial (A1) for ρ + . 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 Ref. [8] yields a maximum β of about 10 −4 N/m 3 , at ±10 V steering range. In a recent experiment reported in [22], where although the purpose was not ion separation a double well potential was produced, the value used was β ≈ 5 × 10 −3 N/m 3 . 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 ext = αq 2 + βq 4 + λq, with λ constant and unknown [10]. Fig. 4 represents the excitation energy versus the energy difference between the two final minima of the external potential, ∆E (also vs. λ). 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 Ref. [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.
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 [23][24][25].
Finally, we compare in Fig. 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 d(t) = d(0) + [d(t f ) − d(0)]s 2 sin(sπ/2), where s = t/t f . From the family of possible potential ramps consistent with this function, we chose a polynomial that drives α from α(0) = α 0 to α(t f ) = −α 0 /2 (as in Fig. 2) and whose first derivatives are 0 at both boundary times. β is given by Eq. (3). For the times analysed in Fig.  5, the method based on Eqs. (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 f ∼ 80 µs would be needed, which is in line with current experiments.
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.