Shortcuts to adiabatic rotation of a two-ion chain

Trapped ions stand out as a flexible architecture to control internal and/or motional states and dynamics for fundamental research of quantum phenomena and technological applications. Pure motional control without internal state transitions is in particular crucial in proposals of two-qubit gates, see e.g. [1], or interferometry [2–4], as well as to scale up the number of ions for quantum information processing [5–12]. The toolbox of basic operations induced by controlling the voltage of electrodes in different Paul trap configurations or detuned laser fields includes transport, expansions and compressions, separation and merging of ion chains, and rotations, the latter being the central topic of this work.

Specific motivations to implement rotations are: reordering an ion chain (to scale up quantum information processing or to locate cooling ions at appropriate positions) [13, 14]; rotation sensing [2]; different simulations (e.g. of black holes [15] or diatomic molecules [16]); probing the exchange phase of quantum statistics [17]; or sorting ions according to charge and mass [18].

Trap rotations, to impart some angular momentum to an ion or ion chain, or to reorient the longitudinal axis of the trap, have been implemented in experiments with improving accuracy [13, 14, 16, 19], and investigated theoretically [20, 21].

Motional control operations, and rotations in particular, need in most applications to be fast, relative to adiabatic dynamics, but also gentle, avoiding final excitations, two requirements met with shortcut to adiabaticity (STA) driving protocols [22]. There are different STA techniques but, for trapped ion driving, invariant-based inverse engineering STA has proven useful [1, 3, 4, 21, 23–31], also to design trap rotations for a single ion [20].

In this paper we extend to a two-ion chain the design of STA one-dimensional (1D) trap rotations done in reference [20]. Our aim is to inverse engineer the time dependence of the rotation angle to implement a fast process, free from final excitations. The work in reference [20] was indeed presented as a preliminary step toward the more complex scenario of the chain rotation, which allows for different, and surely more relevant applications, in particular reordering. Engineering the two-ion rotation also entails non-trivial technical complications due to the increase in the number of equations to be solved, and also because, for some configurations, in particular for two different species in a harmonic trap, there is not in general a point-transformation that provides independent dynamical normal modes [21]. ⁴ Inverse engineering is much easier -to describe the motion and with respect to computational time-for independent modes than for a system which is not separable by point-transformations ⁵ .

We introduce now the basic model. We opt for a cavalier, idealised modeling in the main text where the trap is assumed for simplicity to be a 1D rotating line, e.g. a tightly confined trap in the radial and vertical directions, as depicted in figure 1. We leave aside in the main text peculiarities of the experimental settings, such as micromotion effects and detailed electrode configurations, that may vary significantly among different traps. Our solutions will therefore be guiding starting points for a realistic implementation [19, 33]. Steps toward a more realistic setting, specifically the effect of the transversal dimension and of a circularly symmetric ponderomotive pseudopotential, are explored in two appendices.

Zoom In Zoom Out Reset image size

The trapping line rotates in a horizontal plane in a time t_f up to a predetermined final angle, θ_f = π in all examples. We first find the classical Hamiltonian from the corresponding Lagrangian and then quantize the result. Let s_i , i = 1, 2, denote the points on the line where each ion lays. s_i may take positive and negative values. The Cartesian (lab frame) components of a trajectory s_i (t) are x_i = x_i (s, t), y_i = y_i (s, t),

$$\begin{equation}{x}_{i}={s}_{i}\enspace \mathrm{cos}(\theta ),\quad {y}_{i}={s}_{i}\enspace \mathrm{sin}(\theta ),\end{equation} \tag{ 1 }$$

where θ = θ(t) is the rotation angle. For two different ions in trap potentials f_i (s_i ) the Lagrangian is (we have considered that the magnetic interaction between the two moving charges can be safely neglected, see appendix A)

$$\begin{equation}L=\sum\limits _{i=1,2}\left[\frac{{m}_{i}}{2}{\dot {s}}_{i}^{2}-{f}_{i}({s}_{i})+\frac{{m}_{i}{\dot {\theta }}^{2}}{2}{s}_{i}^{2}\right]-\frac{{C}_{\text{c}}}{{s}_{2}-{s}_{1}},\end{equation} \tag{ 2 }$$

(ion 1 is to the left of ion 2 for θ = 0) with corresponding Hamiltonian

$$\begin{equation}H=\frac{{p}_{1}^{2}}{2{m}_{1}}+\frac{{p}_{2}^{2}}{2{m}_{2}}+V,\end{equation} \tag{ 3 }$$

$$\begin{equation}V=\sum\limits _{i=1,2}\left[{f}_{i}({s}_{i})-\frac{{m}_{i}}{2}{\dot {\theta }}^{2}{s}_{i}^{2}\right]+\frac{{C}_{\text{c}}}{{s}_{2}-{s}_{1}}.\end{equation} \tag{ 4 }$$

In the Coulomb repulsion term, C_c = e²/(4π₀), where ₀ is the vacuum permittivity, and e the electric charge of the electron. The fact that the potentials f_i may be different for different ions makes the model quite flexible, in the sense that the trap could have mass-dependent ponderomotive terms. This is made explicit in appendix C motivated by actual traps [16, 34].

The equilibrium positions $\left\{{s}_{i}^{(0)}\right\}$ of the ions are found by solving the set of equations {∂V/∂s_i = 0}. Since different external traps may be considered, the following equations are for a generic V, the results for harmonic traps are given later in section 2.1.

We define the equilibrium distance between ions as

$$\begin{equation}d={s}_{2}^{(0)}-{s}_{1}^{(0)}\end{equation} \tag{ 5 }$$

and expand V around the equilibrium positions, keeping terms up to second order. Using mass-weighted coordinates ${~{s}}_{i}=\sqrt{{m}_{i}}{s}_{i}$ and momenta ${~{p}}_{i}={p}_{i}/\sqrt{{m}_{i}}$, H is simplified to the quadratic form

$$\begin{equation}H=\frac{{~{p}}_{1}^{2}}{2}+\frac{{~{p}}_{2}^{2}}{2} + ({~{s}}_{1}-{s}_{1}^{(0)},{~{s}}_{2}-{s}_{2}^{(0)}){\mathsf {v}}\left(\begin{matrix}\hfill {~{s}}_{1}-{s}_{1}^{(0)}\hfill \\ \hfill {~{s}}_{2}-{s}_{2}^{(0)}\hfill \end{matrix}\right),\end{equation} \tag{ 6 }$$

where the matrix ${\mathsf {v}}$ has elements ${\mathsf{v}}_{ij}=\frac{1}{\sqrt{{m}_{i}{m}_{j}}}{\left.\frac{{\partial }^{2}V}{\partial {s}_{i}\partial {s}_{j}}\right\vert }_{\left\{{s}_{i},{s}_{j}\right\}=\left\{{s}_{i}^{(0)},{s}_{j}^{(0)}\right\}}$.

We may try to decouple the dynamics by diagonalising ${\mathsf {v}}$. As explained in reference [21], moving to a frame defined by the eigenvectors of ${\mathsf {v}}$ leads, after a classical canonical transformation or, equivalently, quantum unitary transformation, to the following effective Hamiltonian [21]

$$\begin{equation}{H}^{\prime }=\sum\limits _{\nu ={\pm}}\left[\frac{{p}_{\nu }^{2}}{2}+\frac{{{\Omega}}_{\nu }^{2}}{2}{\left({s}_{\nu }+\frac{{\dot {p}}_{0\nu }}{{{\Omega}}_{\nu }^{2}}\right)}^{ 2}\right] - \dot {\mu }({s}_{-}{p}_{+}-{s}_{+}{p}_{-}),\end{equation} \tag{ 7 }$$

where μ is the tilting angle of the potential in configuration space defined by the relation

$$\begin{equation}\mathrm{tan}\enspace 2\mu =\frac{2{v}_{12}}{{v}_{11}-{v}_{22}},\end{equation} \tag{ 8 }$$

and the momentum shifts

$$\begin{equation}{p}_{0+}={\dot {s}}_{1}^{(0)}\sqrt{{m}_{1}}\enspace \mathrm{cos}\enspace \mu +{\dot {s}}_{2}^{(0)}\sqrt{{m}_{2}}\enspace \mathrm{sin}\enspace \mu ,\end{equation} \tag{ 9 }$$

$$\begin{equation}{p}_{0-}=-{\dot {s}}_{1}^{(0)}\sqrt{{m}_{1}}\enspace \mathrm{sin}\enspace \mu +{\dot {s}}_{2}^{(0)}\sqrt{{m}_{2}}\enspace \mathrm{cos}\enspace \mu ,\end{equation} \tag{ 10 }$$

have been defined. The coordinates that diagonalise ${\mathsf {v}}$ are

$$\begin{equation} {s}_{+}=\sqrt{{m}_{1}}({s}_{1} - {s}_{1}^{(0)})\mathrm{cos}\enspace \mu + \sqrt{{m}_{2}}({s}_{2} - {s}_{2}^{(0)})\mathrm{sin}\enspace \mu ,\end{equation} \tag{ 11 }$$

$$\begin{equation} {s}_{-}=-\sqrt{{m}_{1}}({s}_{1} - {s}_{1}^{(0)})\mathrm{sin}\enspace \mu + \sqrt{{m}_{2}}({s}_{2} - {s}_{2}^{(0)})\mathrm{cos}\enspace \mu ,\end{equation} \tag{ 12 }$$

with conjugate momenta

$$\begin{equation}{p}_{+}=\frac{\mathrm{cos}\enspace \mu }{\sqrt{{m}_{1}}}{p}_{1}+\frac{\mathrm{sin}\enspace \mu }{\sqrt{{m}_{2}}}{p}_{2},\end{equation} \tag{ 13 }$$

$$\begin{equation}{p}_{-}=-\frac{\mathrm{sin}\enspace \mu }{\sqrt{{m}_{1}}}{p}_{1}+\frac{\mathrm{cos}\enspace \mu }{\sqrt{{m}_{2}}}{p}_{2}.\end{equation} \tag{ 14 }$$

The squares of the frequencies are

$$\begin{align}{{\Omega}}_{+}^{2}& ={v}_{11}\enspace {\mathrm{cos}}^{2}\enspace \mu +{v}_{22}\enspace {\mathrm{sin}}^{2}\enspace \mu +{v}_{12}\enspace \mathrm{sin}\enspace 2\mu ,\\ {{\Omega}}_{-}^{2}& ={v}_{11}\enspace {\mathrm{sin}}^{2}\enspace \mu +{v}_{22}\enspace {\mathrm{cos}}^{2}\enspace \mu -{v}_{12}\enspace \mathrm{sin}\enspace 2\mu .\end{align} \tag{ 15 }$$

s_± describe independent, dynamical normal modes whenever μ is time independent, see equation (7). In a quantum scenario this means that any wave-function dynamics can be decomposed in terms of the dynamics of two independent harmonic oscillators with time-dependent parameters. Different scenarios to achieve this decoupling will be discussed, for equal ions see section 3 and, for different ions, the first paragraph of section 4 and appendix B.

2.1. Results for harmonic traps

We consider now harmonic potentials, ${f}_{i}({s}_{i})={k}_{i}{s}_{i}^{2}/2$, for simplicity and because of the feasible experimental realization [16, 34]. Equation (4) takes the form

$$\begin{equation}V=\frac{1}{2}{u}_{1}{s}_{1}^{2}+\frac{1}{2}{u}_{2}{s}_{2}^{2}+\frac{{C}_{\text{c}}}{{s}_{2}-{s}_{1}},\end{equation} \tag{ 16 }$$

where

$$\begin{equation}{u}_{i}={m}_{i}({\omega }_{i}^{2}-{\dot {\theta }}^{2}),\hspace{25.0pt}{m}_{i}{\omega }_{i}^{2}={k}_{i}.\end{equation} \tag{ 17 }$$

The u_i are effective spring constants affected by the rotation speed. Unless m₁ = m₂, they are different for both ions. With this V we find the explicit relations

$$\begin{align}\begin{aligned}{s}_{i}^{(0)}& =-{\left[\frac{{C}_{\text{c}}{u}_{j}^{2}}{{u}_{i}{({u}_{i}+{u}_{j})}^{2}}\right]}^{ 1/3} ,\quad i\ne j,\\ d& ={\left[\frac{{C}_{\text{c}}({u}_{1}+{u}_{2})}{{u}_{1}{u}_{2}}\right]}^{1/3},\\ {\mathsf {v}}& =\left(\begin{matrix}\hfill \frac{\frac{2{C}_{\text{c}}}{{d}^{3}}+{u}_{1}}{{m}_{1}}\hfill & \hfill -\frac{2{C}_{\text{c}}}{{d}^{3}\sqrt{{m}_{1}{m}_{2}}}\hfill \\ \hfill -\frac{2{C}_{\text{c}}}{{d}^{3}\sqrt{{m}_{1}{m}_{2}}}\hfill & \hfill \frac{\frac{2{C}_{\text{c}}}{{d}^{3}}+{u}_{2}}{{m}_{2}}\hfill \end{matrix}\right),\\ \mathrm{tan}\enspace 2\mu & =\frac{4{C}_{\text{c}}\sqrt{{m}_{1}{m}_{2}}}{({m}_{1}-{m}_{2})2{C}_{\text{c}}+{d}^{3}({m}_{1}{k}_{2}-{m}_{2}{k}_{1})}.\end{aligned}\end{align} \tag{ 18 }$$

In the main text we shall only consider the simple case of a common harmonic trap, namely, k₁ = k₂ = k. In physical terms this configuration occurs for equal masses, irrespective of the origin (electrostatic, ponderomotive or mixed) of the longitudinal trapping potential, as well as for different masses if the longitudinal trapping is purely electrostatic. A configuration with k₁ ≠ k₂ due to ponderomotive pseudopotentials for different ions is discussed in appendix C.

If the ions are equal, m₁ = m₂ = m, the tilting angle takes the constant value μ = −π/4. The decoupling condition is therefore identically satisfied at all times. Also, u₁ = u₂ = u = mω², with

$$\begin{equation}{\omega }^{2}={\omega }_{0}^{2}-{\dot {\theta }}^{2}\end{equation} \tag{ 19 }$$

and ω₀ constant. The angular velocity of the rotation $\dot {\theta }(t)$ could be negative at some intervals, whereas ω² may also be positive or negative. The equilibrium positions are simplified to

$$\begin{equation}{s}_{2}^{(0)}=-{s}_{1}^{(0)}=\frac{{x}_{0}}{2}\quad \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\quad {x}_{0}={\left(\frac{2{C}_{\text{c}}}{m{\omega }^{2}}\right)}^{1/3},\end{equation} \tag{ 20 }$$

which are symmetrical with respect to the trap center s = 0. The decoupled, effective Hamiltonian is therefore

$$\begin{equation}{H}^{{\prime\prime}}=\sum\limits _{\nu ={\pm}}\left[\frac{{p}_{\nu }^{2}}{2}+\frac{{{\Omega}}_{\nu }^{2}}{2}{\left({s}_{\nu }+\frac{{\dot {p}}_{0\nu }}{{{\Omega}}_{\nu }^{2}}\right)}^{2}\right],\end{equation} \tag{ 21 }$$

with

$$\begin{align}\begin{aligned}{s}_{{\pm}}& =\sqrt{\frac{m}{2}}\left[\left({s}_{1}+\frac{{x}_{0}}{2}\right)\mp \left({s}_{2}-\frac{{x}_{0}}{2}\right)\right],\\ {p}_{0+}& =-\sqrt{\frac{m}{2}}{\dot {x}}_{0},\quad {p}_{0-}=0,\\ {{\Omega}}_{+}^{2}& =3{\omega }^{2},\hspace{25.0pt}{{\Omega}}_{-}^{2}={\omega }^{2}.\end{aligned}\end{align} \tag{ 22 }$$

We consider rotation protocols with a smooth behavior of θ at the boundary times t_b = 0, t_f,

$$\begin{equation}\theta (0)=0,\hspace{25.0pt}\theta ({t}_{\text{f}})={\theta }_{\text{f}},\end{equation} \tag{ 23 }$$

$$\begin{equation}\dot {\theta }({t}_{\text{b}})=\ddot {\theta }({t}_{\text{b}})=0.\end{equation} \tag{ 24 }$$

These conditions imply that

$$\begin{equation}\omega ({t}_{\text{b}})={\omega }_{0},\end{equation} \tag{ 25 }$$

$$\begin{equation}\dot {\omega }({t}_{\text{b}})=\ddot {\omega }({t}_{\text{b}})={\dot {p}}_{0\nu }({t}_{\text{b}})=0.\end{equation} \tag{ 26 }$$

The two independent harmonic oscillators expand or compress through the time dependence of Ω_ν and experiment a 'transport', in s_ν space, along ${\dot {p}}_{0\nu }/{{\Omega}}_{\nu }^{2}$. The Hamiltonian (21) has a dynamical invariant [35]

$$\begin{equation}I=\sum\limits _{\nu ={\pm}}\frac{1}{2}{[{b}_{\nu }({p}_{\nu }-{\dot {\alpha }}_{\nu })-{\dot {b}}_{\nu }({s}_{\nu }-{\alpha }_{\nu })]}^{2}+\frac{1}{2}{{\Omega}}_{0\nu }^{2}{\left(\frac{{s}_{\nu }-{\alpha }_{\nu }}{{b}_{\nu }}\right)}^{2},\end{equation} \tag{ 27 }$$

where Ω_0± = Ω_±(0), and b_± (scaling factors of the normal mode wavefunctions) and α_± (reference classical trajectories for each oscillator) are auxiliary functions that have to satisfy, respectively, the Ermakov and Newton equations,

$$\begin{equation}{\ddot {b}}_{{\pm}}+{{\Omega}}_{{\pm}}^{2}{b}_{{\pm}}=\frac{{{\Omega}}_{0{\pm}}^{2}}{{b}_{{\pm}}^{3}},\end{equation} \tag{ 28 }$$

$$\begin{equation}{\ddot {\alpha }}_{{\pm}}+{{\Omega}}_{{\pm}}^{2}{\alpha }_{{\pm}}={\dot {p}}_{0{\pm}}.\end{equation} \tag{ 29 }$$

The time-dependent Schrödinger equation can be solved by superposing, with constant coefficients, elementary solutions which are also eigenstates of the invariant, with the ('Lewis–Riesenfeld') phase adjusted to be also solutions of the Schrödinger equation [23],

$$\begin{equation}\vert {\psi }_{n{\pm}}^{{\prime\prime}}\rangle ={e}^{\frac{i}{\hslash }\left[\frac{{\dot {b}}_{{\pm}}{s}_{{\pm}}^{2}}{2{b}_{{\pm}}}+({\dot {\alpha }}_{{\pm}}{b}_{{\pm}}-{\alpha }_{{\pm}}{\dot {b}}_{{\pm}})\frac{{s}_{{\pm}}}{{b}_{{\pm}}}\right]}\frac{1}{\sqrt{{b}_{{\pm}}}}{{\Phi}}_{n}({\sigma }_{{\pm}}),\end{equation} \tag{ 30 }$$

where ${\sigma }_{{\pm}}=\frac{{s}_{{\pm}}-{\alpha }_{{\pm}}}{{b}_{{\pm}}}$ and Φ_n are the eigenfunctions for the static harmonic oscillators with frequencies Ω_0,±. The average energies for the nth elementary solution of each mode can be calculated [29, 30],

$$\begin{align}{E}_{n{\pm}}^{{\prime\prime}}& =\langle {\psi }_{n{\pm}}^{{\prime\prime}}\vert {H}^{{\prime\prime}}\vert {\psi }_{n{\pm}}^{{\prime\prime}}\rangle \\ & =\frac{(2n+1)\hslash }{4{{\Omega}}_{0{\pm}}}\left({\dot {b}}_{{\pm}}^{2}+{{\Omega}}_{{\pm}}^{2}{b}_{{\pm}}^{2}+\frac{{{\Omega}}_{0{\pm}}^{2}}{{b}_{{\pm}}^{2}}\right)+\frac{1}{2}{\dot {\alpha }}_{{\pm}}^{2}+\frac{1}{2}{{\Omega}}_{{\pm}}^{2}{\left({\alpha }_{{\pm}}-\frac{{\dot {p}}_{0{\pm}}}{{{\Omega}}_{{\pm}}^{2}}\right)}^{2}.\end{align} \tag{ 31 }$$

As ${\dot {p}}_{0\nu }({t}_{\text{f}})=0$, the final values are minimised when the only contribution is due to the eigenenergies for the oscillators, with

$$\begin{equation}{b}_{{\pm}}({t}_{\text{f}})=1,\alpha ({t}_{\text{f}})=\dot {\alpha }({t}_{\text{f}})={\dot {b}}_{{\pm}}({t}_{\text{f}})=0.\end{equation} \tag{ 32 }$$

3.1. Inverse engineering

Imposing commutativity between Hamiltonian and invariant at initial t = 0 and final times t = t_f, the invariant drives the initial eigenstates of H to corresponding final eigenstates along the elementary solutions (30), although there could be diabatic excitations at intermediate times, when the commutation between Hamiltonian and invariant is not guaranteed. By inspection of equation (30), commutativity at the boundary times is achieved if the conditions in equation (32) are satisfied, which occur automatically when the final energies (31) are minimised. To inverse engineer the rotation we proceed similarly to reference [20], with an ansatz for θ(t) that satisfies boundary conditions (23) and (24) with some free parameters. We use up to 4 free parameters,

$$\begin{align}\theta (t)& =\frac{1}{16}(32{c}_{3}+80{c}_{4}+144{c}_{5}+224{c}_{6}-9{\theta }_{\text{f}})\mathrm{cos}\left(\frac{\pi t}{{t}_{\text{f}}}\right)\\ & \quad -\frac{1}{16}(48{c}_{3}+96{c}_{4}+160{c}_{5}+240{c}_{6}-{\theta }_{\text{f}})\mathrm{cos}\left(\frac{3\pi t}{{t}_{\text{f}}}\right)\\ & \quad +{c}_{3}\enspace \mathrm{cos}\left(\frac{5\pi t}{{t}_{\text{f}}}\right)+{c}_{4}\enspace \mathrm{cos}\left(\frac{7\pi t}{{t}_{\text{f}}}\right)+{c}_{5}\enspace \mathrm{cos}\left(\frac{9\pi t}{{t}_{\text{f}}}\right)+{c}_{6}\enspace \mathrm{cos}\left(\frac{11\pi t}{{t}_{\text{f}}}\right)+\frac{{\theta }_{\text{f}}}{2}.\end{align} \tag{ 33 }$$

This gives an expression of $\dot {\theta }$, from which we find ω in equation (19). We introduce ω in (22) to get the normal mode angular frequencies Ω_± needed in the Ermakov equation (28). For a given set of values of these parameters we solve the 'direct problem' (Ermakov and Newton equations) with initial conditions

$$\begin{equation}\begin{gathered}{b}_{{\pm}}(0)=1,{\dot {b}}_{{\pm}}(0)=0,\\ {\alpha }_{{\pm}}(0)={\dot {\alpha }}_{{\pm}}(0)=0,\end{gathered}\end{equation} \tag{ 34 }$$

and compute easily the final energies with equation (31). The values of the parameters are varied with a subroutine that minimises the sum of the final mode energies (31) (we use the MATLAB 'fminsearch' and n = 0 but note that the optimal final values of b(t_f), α(t_f) and their derivatives would minimise the energies for any n). The excess energy found with the optimal parameters for the normal modes is negligible in the range of final times depicted in figure 2.

Zoom In Zoom Out Reset image size

Once the free parameters are defined such that the design of θ(t) minimises the excitation energy of the normal modes, we perform the quantum evolution driven by the full Hamiltonian with (16) to check the performance of the designed protocol. We use the 'split-operator method', and the initial ground state is found performing an evolution in imaginary time. Figure 2 shows the final excitation, i.e. the excess energy with respect to the initial energy after performing the evolution with the full Hamiltonian (3) using the potential (16). In figure 2(a) this excitation is depicted in a linear scale, and in figure 2(b) in a logarithmic scale. The results improve significantly by using more optimisation parameters. Even when using a single optimising parameter, the results are clearly better than the protocol without free parameters. Figure 3 shows some examples of the rotation protocols with 4 parameters for different rotation times.

Zoom In Zoom Out Reset image size

Let us first explore some possible manipulations to make the modes separable when the ions are different. For k₁ = k₂ = k, (18) gives

$$\begin{equation}\mathrm{tan}\enspace 2\mu =\frac{4{C}_{\text{c}}\sqrt{{m}_{1}{m}_{2}}}{({m}_{1}-{m}_{2})(2{C}_{\text{c}}+{d}^{3}k)},\end{equation} \tag{ 35 }$$

so d³ k should be constant to keep μ constant and have independent modes. If the only parameter that depends on time is d, this condition cannot be satisfied. (This is also true for k₁ ≠ k₂.) But if k is allowed to be a time-dependent controllable parameter, it would be in principle possible to keep d³ k constant. Using the expressions for d and the u_i , this condition may be satisfied for two values of ${\dot {\theta }}^{2}={a}_{l}k$ for each k, l = 1, 2, where the a_l are two constants. The proportionality between ${\dot {\theta }}^{2}$ and k, however, is problematic. If we wish to approach $\dot {\theta }=0$ smoothly at the time boundaries, then k → 0 there, which implies a vanishing trapping potential and d → ∞. A way out is explored in appendix B making use of a more complex external trap potential with linear and quartic terms added, as in reference [31]. In the main text we stay within the harmonic trap configuration with constant k and renounce to separate the modes. Thus a different, pragmatic strategy is adopted, minimising the excitation energy directly to find the rotation protocol. We use the same ansatz for the control parameter θ as in equation (33) and solve the full (quantum) dynamics for the potential (16) to find the final excess energy for specific values of the free parameters c_3–6. Then, as in section 3.1, we minimise the excess energy letting the MATLAB subroutine 'fminsearch' find the optimal parameters.

In figure 4 we depict this final excitation, optimising the result using from 1 to 4 free parameters for the θ, and compare it with the results for no free parameters. This direct minimisation provides even better results than the indirect one based on the normal mode energy in section 3.1. The best protocol (4 optimising parameters) gives an excitation below 0.1 quanta at a final time t_f = 0.56 μs. The price to pay though, is that the computational time required increases dramatically, as we have to solve the full dynamics of the system at every iteration of the shooting method, whereas in the method based on normal modes we only needed to solve four ordinary differential equations at each iteration. Figure 5 shows the equilibrium and dynamical positions of both ions during the evolution for t_f = 1 μs. The trajectories are not symmetric since the two ions experience different effective spring constants, see equation (17).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In configurations where mass-dependent ponderomotive pseudopotentials [36] affect the longitudinal trapping, k₁ and k₂ are no longer equal, but there are no substantial changes in the formal treatment, as discussed in appendix C, and the same methodology proposed here can be applied.

We have designed protocols to rotate a linear, transversally tight trap containing two ions, without final excitation. For two equal ions in a rotating, rigid harmonic trap, there are uncoupled dynamical normal modes. The separation facilitates inverse engineering since it is only necessary to solve ordinary differential equations for independent variables to minimise the final energy. These Ermakov and Newton equations are for the auxiliary functions in the invariants associated with the uncoupled Hamiltonians. Following this method and for a given ansatz for the rotation angle and for some allowed final excitation threshold, process-time lower limits are met due to the eventual failure of the small oscillation regime for very rapid rotations. Faster processes can be achieved by increasing the number of parameters in the ansatz. For two different ions in a harmonic trap, this method is not possible as the modes are coupled for a rigid trap, or can be uncoupled for a non-rigid trap but only for impractical boundary conditions for the trap. Instead we used direct optimisation of the rotation ansatz parameters with the full Hamiltonian. This direct approach is efficient to design fast processes but the computational effort is much more demanding.

The idealised conditions assumed-1D motion for each ion along a rotating line and a common trap even for different ions-are useful to set limits. Of course specific trap settings need a dedicated analysis. Steps toward experimental implementations are given in appendices C and D that discuss a 2D configuration with mass-dependent trap potentials due to ponderomotive pseudopotentials. In particular the harmonic model is very well adapted to the setting in references [16, 34]. As well, modest transverse to longitudinal frequency ratios are shown to be sufficient to render the 1D results a useful guidance. In this regard, progress on applying STA to systems of coupled oscillators in non-separable configurations [37, 38], may eventually lead to STA protocols that could take the actual dimensionality of the rotating trap for the ion chain into account from the start. At this point though, such a program has not yet been fulfilled, so that the reference of 1D approximation is very valuable.

Several natural extensions of this work are possible: for example, to consider different boundary conditions, such as a final rotating trap with $\dot {\theta }({t}_{\text{f}})\ne 0$, as in reference [16], to transfer an angular momentum to the chain. Another possible extension would be adding noises and perturbations in the model to make the protocols robust with respect to them [39, 40]. Finally, specific protocols could be designed to simultaneously rotate longer chains of ions, although it is possible to sequentially rotate them in groups of 2 using the protocols designed here.

We thank Uli Poschinger for discussions on the early stages of this paper. This work was supported by the Basque Country Government (Grant No. IT986-16), and by the Spanish Ministry of Science and Innovation through projects PGC2018-101355-B-I00 and PGC2018-095113-B-I00 (MCIU/AEI/FEDER, UE).

The data that support the findings of this study are available upon reasonable request from the authors.

Two charged particles moving in a direction perpendicular to the direction in which they are aligned experience a magnetic force, with magnitude

$$\begin{equation}{F}_{\text{mag}}=\frac{{\mu }_{0}}{4\pi }\frac{{e}^{2}}{{r}^{2}}\vert { \overrightarrow {v}}_{1}{\times}({ \overrightarrow {v}}_{2}{\times}\hat{r})\vert ,\end{equation} \tag{ A.1 }$$

where μ₀ is the permeability constant, ${ \overrightarrow {v}}_{i}$ the velocity vectors of each ion, $\overrightarrow {r}={ \overrightarrow {s}}_{2}-{ \overrightarrow {s}}_{1}$ the position vector of ion 2 with reference to ion 1, and $\hat{r}= \overrightarrow {r}/r$. The Coulomb interaction, which is the only one considered so far, gives a force of magnitude

$$\begin{equation}{F}_{\text{el}}=\frac{{C}_{\text{c}}}{{r}^{2}}.\end{equation} \tag{ A.2 }$$

The ratio of these two forces is, using μ_{0 0} = c⁻², where c is the speed of light,

$$\begin{equation}R=\frac{{F}_{\text{mag}}}{{F}_{\text{el}}}=\frac{\vert { \overrightarrow {v}}_{1}{\times}({ \overrightarrow {v}}_{2}{\times}\hat{r})\vert }{{c}^{2}}.\end{equation} \tag{ A.3 }$$

With $\vert { \overrightarrow {v}}_{1}\vert ,\vert { \overrightarrow {v}}_{2}\vert \approx \frac{r}{2}\dot {\theta }$ we get

$$\begin{equation}R\approx \frac{{r}^{2}{\dot {\theta }}^{2}}{4{c}^{2}}.\end{equation} \tag{ A.4 }$$

For the protocols designed in the main text, the maximum values during the simulations at the represented times are ${\dot {\theta }}_{\mathrm{max}}=5{\times}1{0}^{6}\enspace {\mathrm{s}}^{-1}$ and r_max = 5.5 × 10⁻⁶ m so the magnetic interaction is negligible with respect to the electric force.

If the matrix ${\mathsf {v}}$ is time dependent, the normal modes get decoupled if ${\mathsf{v}}_{11}={\mathsf{v}}_{22}$, see equation (8), i.e.

$$\begin{equation}{\left.{\frac{1}{m}}_{1}\frac{{\partial }^{2}V}{\partial {s}_{1}^{2}}\right\vert }_{{s}_{1}^{(0)}}={\left.{\frac{1}{m}}_{2}\frac{{\partial }^{2}V}{\partial {s}_{2}^{2}}\right\vert }_{{s}_{2}^{(0)}}.\end{equation} \tag{ B.1 }$$

As explained in the main text, the rotation of different ions trapped by harmonic potentials cannot be described in general in terms of dynamical normal modes. For a non-rigid trap there is a formal solution which does not lead to practically useful boundary conditions. Here we consider different confining potentials that obey equation (B.1), and thus allow us to inverse engineer the rotation using the Lewis–Riesenfeld family of invariants. We use for the equilibrium positions the parametrisation ${s}_{1}^{(0)}={s}_{0}-d/2$ and ${s}_{2}^{(0)}={s}_{0}+d/2$, where s₀ is the middle point between them.

Specifically we consider a tilted double well potential, which combines a repulsive harmonic potential with the confinement provided by the quartic term and a linear term [31],

$$\begin{equation}V=\gamma (t)({s}_{1}+{s}_{2})+\frac{1}{2}{u}_{1}(t){s}_{1}^{2}+\frac{1}{2}{u}_{2}(t){s}_{2}^{2}+{\beta }_{1}{s}_{1}^{4}+{\beta }_{2}{s}_{2}^{4}+\frac{{C}_{\text{c}}}{{s}_{2}-{s}_{1}}.\end{equation} \tag{ B.2 }$$

As for its feasibility, 'Mexican hat' ponderomotive pseudopotentials to provide the basic double well shape have been already realized [41], whereas the tilt and rotation could be implemented by a biased rotating quadrupole field, see also appendix C. This form of V gives the potential matrix

$$\begin{equation}{\mathsf {v}}=\left(\begin{matrix}\hfill \frac{\frac{2{C}_{\text{c}}}{{d}^{3}}+{u}_{1}+12{\left(-\frac{d}{2}+{s}_{0}\right)}^{ 2}{\beta }_{1}}{{m}_{1}}\hfill & \hfill -\frac{2{C}_{\text{c}}}{{d}^{3}\sqrt{{m}_{1}{m}_{2}}}\hfill \\ \hfill -\frac{2{C}_{\text{c}}}{{d}^{3}\sqrt{{m}_{1}{m}_{2}}}\hfill & \hfill \frac{\frac{2{C}_{\text{c}}}{{d}^{3}}+{u}_{2}+12{\left(\frac{d}{2}+{s}_{0}\right)}^{ 2}{\beta }_{2}}{{m}_{2}}\hfill \end{matrix}\right).\end{equation} \tag{ B.3 }$$

Equating the diagonal terms leads to an equation from which the needed γ(t) may be found. The general case is however very cumbersome to treat. To illustrate the details of the procedure with the simplest possible example, we shall set β = β₁ = β₂, and also take k = k₁ = k₂ for the quadratic term in (17). This idealised case would correspond to a purely electrostatic origin of the tilted double well potential.

The main-text equations from equations (7) to (15) are still valid here. We assume that the controllable parameters are the linear potential and the rotation speed. Equating diagonal terms in (B.3), d is found to obey

$$\begin{align}& \left\{A({m}_{1}-{m}_{2})({u}_{1}-{u}_{2})+{d}^{2}\left[6A\beta ({m}_{1}+{m}_{2})+({m}_{1}-{m}_{2}){({u}_{1}-{u}_{2})}^{2}\right]\right.\\ & \hspace{25.0pt}\left.+24\beta {C}_{\text{c}}d({m}_{1}-{m}_{2})+12{\beta }^{2}{d}^{6}({m}_{1}-{m}_{2})\right\}=0,\end{align} \tag{ B.4 }$$

where we have defined

$$\begin{equation}A=\sqrt{{d}^{3} \left[24\beta {C}_{\text{c}} - 12{\beta }^{2}{d}^{5} - 12\beta {d}^{3}({u}_{1} + {u}_{2}) + d{({u}_{1} - {u}_{2})}^{2}\right]}.\end{equation} \tag{ B.5 }$$

The force γ(t) that would produce the desired evolution for d is, using (B.1) and (B.2),

$$\begin{align}\gamma & =\frac{1}{108{\beta }^{2}{d}^{6}}\left\{18\beta {C}_{\text{c}}{d}^{2}({u}_{2}-{u}_{1})-24{\beta }^{2}{d}^{5}A-d{({u}_{1}-{u}_{2})}^{2}A\right.\\ & \quad -\left.6\beta {C}_{\text{c}}A+36{\beta }^{2}{d}^{7}({u}_{1}-{u}_{2})+{d}^{3}\left[-6\beta A({u}_{1}+{u}_{2})-{({u}_{1}-{u}_{2})}^{3}\right]\right\},\end{align} \tag{ B.6 }$$

and the corresponding evolution for the middle point between the ions is

$$\begin{equation}{s}_{0}=\frac{A+{d}^{2}({u}_{1}-{u}_{2})}{12\beta {d}^{3}}.\end{equation} \tag{ B.7 }$$

The frequencies of the normal modes Ω_± can be analytically expressed in terms of d, the parameters that define the potential (u₁, u₂ and β) and the masses m₁ and m₂, but they are too lengthy to be reproduced here. Provided that equation (B.4) is satisfied, the rotation of the potential in equation (B.2) is governed by an uncoupled Hamiltonian of the form (21), with the corresponding frequencies Ω_± and momentum shifts that read

$$\begin{equation}{p}_{0{\pm}}=\frac{1}{\sqrt{2}}\left[(\sqrt{{m}_{1}}{\pm}\sqrt{{m}_{2}})\dot {s}+(\sqrt{{m}_{2}}\mp \sqrt{{m}_{1}})\frac{\dot {d}}{2}\right].\end{equation} \tag{ B.8 }$$

From here on the procedure to design the protocol is similar to the one explained in section 3.1. We start from the same ansatz for θ(t), see equation (33), which satisfies the boundary conditions (23) and (24) by design, and search for the values of the free parameters that minimise the final excitation. Decoupling the dynamics of the system into independent dynamical normal modes, however, is more demanding here than for equal ions. We compute the necessary force, see equation (B.6), and equilibrium positions, see equations (B.4) and (B.7), for each test value of the free parameters in θ(t).

Figure B1(a) shows that, for a rotation of a ⁴⁰Ca⁺ and a ⁹Be⁺ ion chain, any of the protocols produce no excitations in the normal modes for processes as fast as 0.4 μs. It also illustrates the improvement of the results by increasing the number of free parameters for θ(t). Normal mode excitation is an approximation of the exact excitation, nevertheless, our results suggest that performing the rotation with the double well may provide excitationless protocols at short time scales.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures B1(b) and (c) depict, respectively, the initial potential and the required force γ(t) for a specific rotation protocol using the tilted double well potential in equation (B.2). Notice that even the lowest value of the force, at boundary times, produces a considerable bias with little to none barrier potential between the two wells. Despite this, each equilibrium position, whose evolution is depicted in figure B1(c), initially lays in its own well. This unusual potential shape would be the price to pay for mode separability. We note that a potential bias may be imposed or canceled using STA methods as well [42].

In the main text we use k₁ = k₂ for a centered, harmonic trap and different ions, which occurs if the longitudinal trapping is purely electrostatic. This is admittedly challenging, but, as we shall see in this appendix, nothing fundamental changes in the structure of the equations for a more feasible setting in which a ponderomotive potential with cylindrical symmetry [16] combines with an electrostatic one in the longitudinal direction, so that k₁ ≠ k₂ because of the mass dependences.

For a single ion, we assume an 'electrostatic' quadrupole potential with principal axes s, s', z,

$$\begin{equation}{V}_{\text{DC}}={\alpha }_{s}{s}^{2}+{\alpha }_{{s}^{\prime }}{{s}^{\prime }}^{2}+{\alpha }_{z}{z}^{2}\end{equation} \tag{ C.1 }$$

where the s and s' axes rotate about z with respect to the laboratory frame axes by an angle θ. Rotating frame and lab frame coordinates on the rotation plane are related by (figure C1)

$$\begin{equation}s=x\enspace \mathrm{cos}\enspace \theta +y\enspace \mathrm{sin}\enspace \theta ,\end{equation} \tag{ C.2 }$$

$$\begin{equation}{s}^{\prime }=-x\enspace \mathrm{sin}\enspace \theta +y\enspace \mathrm{cos}\enspace \theta .\end{equation} \tag{ C.3 }$$

Laplace's equation implies

$$\begin{equation}{\alpha }_{s}+{\alpha }_{{s}^{\prime }}+{\alpha }_{z}=0,\end{equation} \tag{ C.4 }$$

and we assume the simple configuration

$$\begin{equation}{\alpha }_{{s}^{\prime }}=-{\alpha }_{s}=\alpha ;{\alpha }_{z}=0,\end{equation} \tag{ C.5 }$$

In addition let us suppose that there is a radio-frequency (rf) potential with cylindrical symmetry -the simplest symmetry to implement rotations-of the form Φ(z, r)cos(Ωt), where r is the radial coordinate $r={({x}^{2}+{y}^{2})}^{1/2}$. The rf potential is treated approximately by a static ponderomotive pseudopotential [36],

$$\begin{equation}{\Psi}(z,r)=\frac{{e}^{2}}{4m{{\Omega}}^{2}}\vert \nabla {\Phi}(z,r){\vert }^{2}.\end{equation} \tag{ C.6 }$$

A set of planar, annular, concentric electrodes allows for many different 3D trap designs [43]. We shall assume assume a strongest confinement in the z direction so that the motion is effectively 2D, at some distance from the electrodes, as in [16]. The dependence of Ψ(z, r) on r may be adjusted to different forms. For the sake of concreteness and for comparison with the harmonic model in the main text we assume here a 'point' harmonic trap (i.e. with harmonic terms in the three directions) with the radial part characterized as

$$\begin{equation}{V}_{\text{radial}}=\frac{{\lambda }^{2}}{m{{\Omega}}^{2}}({s}^{2}+{{s}^{\prime }}^{2}).\end{equation} \tag{ C.7 }$$

This form has been implemented by Urban et al experimentally [16, 34] by three rings with the rf voltage applied to the second one and the others grounded. See [16, 34] for technical details on details on ring sizes and applied voltages.

To write down the effective Lagrangian for the two ions, a subindex i = 1, 2 is set for each ion's coordinates and mass. The resulting (classical) Hamiltonian, written for simplicity in terms of s and $\dot {s}$ takes the form

$$\begin{equation}H=\frac{1}{2}\sum\limits _{i}\;{m}_{i}({{\dot {s}}_{i}}^{2}+{\dot {s}}_{i}^{\prime 2})+{V}_{\text{trap}}-\frac{{\dot {\theta }}^{2}}{2}\sum\limits _{i}\;{m}_{i}({s}_{i}^{2}+{{s}_{i}^{\prime }}^{2})+{V}_{\text{Coul}},\end{equation} \tag{ C.8 }$$

where

$$\begin{equation}{V}_{\text{trap}}=\frac{1}{2}\sum\limits _{i}\;{m}_{i}({\omega }_{{s}_{i}}^{2}{s}_{i}^{2}+{\omega }_{{s}_{i}^{\prime }}^{2}{{s}_{i}^{\prime }}^{2}),\end{equation} \tag{ C.9 }$$

$$\begin{equation}{\omega }_{{s}_{i}}^{2}=\frac{2{\lambda }^{2}}{{m}_{i}^{2}{{\Omega}}^{2}}-\frac{2\alpha }{{m}_{i}},\end{equation} \tag{ C.10 }$$

$$\begin{equation}{\omega }_{{s}_{i}^{\prime }}^{2}=\frac{2{\lambda }^{2}}{{m}_{i}^{2}{{\Omega}}^{2}}+\frac{2\alpha }{{m}_{i}},\end{equation} \tag{ C.11 }$$

$$\begin{equation}{V}_{\text{Coul}}=\frac{{C}_{\text{c}}}{{[{({s}_{2}-{s}_{1})}^{2}+{({s}_{2}^{\prime }-{s}_{1}^{\prime })}^{2}]}^{1/2}}.\end{equation} \tag{ C.12 }$$

In terms of the conjugate momenta (we skip ion labels)

$$\begin{equation}{P}_{s}=m\dot {s}-m{s}^{\prime }\dot {\theta },\end{equation} \tag{ C.13 }$$

$$\begin{equation}{P}_{{s}^{\prime }}=m{\dot {s}}^{\prime }+ms\dot {\theta },\end{equation} \tag{ C.14 }$$

the Hamiltonian takes a form that depends on 'angular momentum' couplings between longitudinal and transversal motions,

$$\begin{equation}H=\sum\limits _{i}\frac{{P}_{{s}_{i}}^{2}}{2{m}_{i}}+\frac{{P}_{{s}_{i}^{\prime }}^{2}}{2{m}_{i}}-\dot {\theta }\sum\limits _{i}({s}_{i}{P}_{{s}_{i}^{\prime }}-{s}_{i}^{\prime }{P}_{{s}_{i}})+{V}_{\text{trap}}+{V}_{\text{Coul}}\end{equation} \tag{ C.15 }$$

A 1D version is found in the limit of tight transversal confinement in direction s'. Setting ${s}^{\prime }={\dot {s}}^{\prime }=0$, the formal results and treatment in sections 1 and 2.1 are still valid, with k_i = 2λ²/(m_i Ω²) − 2α.

In this appendix we test the range of validity of the 1D trap model and the optimised rotation protocols found for that model when transversal motion is allowed. We use the Hamiltonian in appendix C and equation (C.15), for equal masses, so that the effective frequencies do not depend on the specific ion. As a simple test of the effect of the transversal dimension we solve Hamilton's equations of motion starting in the equilibrium configuration and compute the final excitation energy.

Figure D1 shows the excess energy for the protocols with 4 optimising parameters used in figure 2 for different final times. The excitations decay quite rapidly for tighter confinements. For modest, feasible frequency ratios (∼2), and for relatively short final times (∼2 μs), we can already get excitations below 1 quantum.

Zoom In Zoom Out Reset image size