Robust dynamical exchange cooling with trapped ions

We investigate theoretically the possibility for robust and fast cooling of a trapped atomic ion by transient interaction with a pre-cooled ion. The transient coupling is achieved through dynamical control of the ions’ equilibrium positions. To achieve short cooling times we make use of shortcuts to adiabaticity by applying invariant-based engineering. We design these to take account of imperfections such as stray fields, and trap frequency offsets. For settings appropriate to a currently operational trap in our laboratory, we find that robust performance could be achieved down to 6.3 motional cycles, comprising 14.2 μs for ions with a 0.44 MHz trap frequency. This is considerably faster than can be achieved using laser cooling in the weak coupling regime, which makes this an attractive scheme in the context of quantum computing.


Introduction
One of the major challenges in quantum computing is to realise fast operations, since these affect both the clock speed as well as the ability to preserve coherence in the presence of decoherence mechanisms. For trapped-ion approaches, direct operations on the qubits include single and multi-qubit gates and state detection. However in order to implement these processes in the flexible manner and with the high reliability required for quantum error-correction, transport of ions and re-cooling are expected to play an important role [1][2][3]. Thus to increase the speed of a trapped-ion quantum information processor all of these processes must be improved. While recent work has demonstrated impressive progress in the speed of one and two-qubit gates [4], detection [33] and transport [5,6], in many recent demonstrations in multi-zone chips the primary speed limitation was due to laser-recooling of ions, either after imperfect transport or following detection events, which heat the ions via photon recoil [3,7]. Laser cooling close to the ground state is performed by resolved-sideband cooling, while methods such as Sisyphus cooling [8] and electromagnetically-induced transparency cooling [9] offer higher rates. These methods are limited in rate by working in the weak coupling regime and by fundamental features of the atom, including finite decay rates for spontaneous emission of around 10 8 s −1 and the imperfect transfer of momentum between the atom and the light field, leading to cooling cycles of several hundred trap periods in current experiments [10]. While it may be possible to perform laser cooling on timescales faster than the ion oscillation period [11], the recoil rate presents a hard limit.
The premise for our work is that exchange of energy of ions via the Coulomb interaction can be used to extract excess energy from a hot ion by bringing it into resonance with a pre-cooled ion in a nearby potential well. Two ions held close to each other in an external potential experience a mutual repulsive force. This modifies the ion equilibrium positions relative to the minima of the external potential, but also couples the vibrations of the two different ions. For two ions of mass m 1 and m 2 separated by a distance d, which are held in harmonic potential wells in which they oscillate with frequencies ω 1 and ω 2 respectively, energy exchange between the oscillations of each ion occurs at a frequency This exchange means that energy can be removed from an initially hot ion by placing it close to a pre-cooled 'coolant' ion for the exchange time t e = π/(2Ω). This becomes useful when the coupling can be turned on and off without inducing excitation. One way to do this is to tune the two potential wells such that the ions come into resonance for a fixed time period, and subsequently detune them from each other. While this has been performed previously in the adiabatic regime [12], for the purposes of re-cooling ions in quantum computers it is desirable to increase the speed with which such operations are implemented. Thus we consider instead dynamic schemes. Taking advantage of the strong 1/d 3 scaling of the exchange frequency, we design an explicitly time-dependent Hamiltonian to transport two ions closer together from an outer distance d 0 . The closest point is briefly reached at an inner distance d in , from where the ions are brought back to their starting position after a total run-time t f . As the ion distance changes throughout the scheme, the exchange frequency is also constantly varying. The aim is to achieve a situation in which the aggregated exchange leads to full energy transfer. This scheme is illustrated in figure 1. We explore the limits of the speed with which this can be carried out using trajectories designed with shortcuts to adiabaticity (STA) [13]. Following a number of previous works [14][15][16][17][18][19][20][21][22], we investigate the robustness of these methods to imperfections, and find solutions which are tolerant to these. Contrary to similar protocols in earlier work where strong approximations were made [23], we consider a realistic double-well trapping potential including terms beyond the harmonic approximation, and assess the performance based on the full Hamiltonian of the system. The resulting schemes should allow robust cooling on timescales of 6 trap cycles, which is competitive with the operation speeds of high-fidelity multi-qubit gates [4].
The challenge of such a method is to design appropriate trajectories which do not add excitations during execution, such that the target ion may end up in the motional ground state. This can be realised by optimising the scheme with ground-state ions, disregarding the motional exchange at first. Thereafter, cooling can be achieved simply by finding the correct timing that leads to a complete exchange of energy.
This work is organised as follows. In section 2, we detail the physical constraints that the scheme needs to adhere to in order to be implementable in a real ion trap. Next, we obtain trajectories which we optimise to minimise the excitation of a pair of ground-state ions starting in separate potential wells. The theoretical techniques towards this end are presented in section 3, whereas necessary numerical optimisation is carried out in section 4. In section 5, we add robustness against unwanted homogeneous electric fields. Having gained the ability to perform this basic protocol robustly and without additional excitations, we apply it in section 6 to find timings where two ions swap their energy. We then analyse the timing of these trajectories when varying the inner distance or the maximal quartic confinement in a given trap. In section 7, we consider the effect of a homogeneous field on the exchange of energy, and calculate the level of control of these fields that would be needed to successfully implement our scheme. Finally in section 8, the energy swapping protocol is generalised and optimised for two ions with unequal mass.

Design constraints
The Hamiltonian for a system of two ions of the same mass m is of the form where {x i , p i } are the canonical position and momentum coordinates of the ions, C C = e 2 /(4π 0 ) and we assume that x 2 > x 1 . In this equation the positions and momenta of the ions are time-dependent, although this is suppressed. We initially choose a symmetric quartic form for a double well potential [24] V el ( The assumption of symmetry will be relaxed in section 8, but provides a useful starting point to understand the methods that we use. The parameters {α(t), β(t)} describe the harmonic and quartic parts of the potential. As long as α(t) < 0 and β(t) > 0, V el is a double-well potential. This potential leads to an equilibrium ion separation d which can be obtained from the equation and normal mode frequencies for the excursions about the equilibrium positions of where Ω − is associated with the in-phase motion (centre-of-mass mode) and Ω + with the out-of-phase motion (stretch mode). Due to the symmetry of the potential, the potential curvatures at the ion positions are equal and given by where i = {1, 2}. We define the initial curvature as ω 0 = ω i (t = 0). The protocol should be designed to transport two ions from a separation of d 0 to a separation of d in and back out again. For simplicity, we assume that the ions start and end the protocol at time t f at the same separation The distance d 0 should be chosen such that if no dynamical changes were made to the potential and the ions were simply kept at that distance during the protocol, the resulting exchange time t e would be slow, which we take to mean that it is on the order of current cooling times and thus above 100 motional periods.
The ions are considered to achieve the distance of closest approach at the halfway point of the protocol t f /2, which we define as d in = d(t f /2). We also choose to constrain our protocol by limiting the value of β β max , which requires the most demanding voltage settings for traps of the scale in use today [24]. To maintain ions in separate wells, which is the aim of our protocol, we can then see that d in must be greater than which is obtained from (5) with α = 0. Attempting to aim for d in < d c would involve combining the two ions into the same potential well, which involves taking the potential wells to conditions producing the lowest trap frequency for a given β. We make the assumption that this would be more experimentally challenging than keeping the ions in separated potentials. From the considerations above, we can now set the parameters of the initial and final potential α 0 ≡ α(0) = α(t f ) and β 0 ≡ β(0) = β(t f ) as well as the intermediate potential given by α in ≡ α(t f /2) and β in ≡ β(t f /2). To maintain two separate wells when the ions are closest, we take β in to its maximum and thus β in = β max . From (5) it then follows that So far, the initial potential parameters α 0 , β 0 are only constrained by the selected outer distance d 0 and can otherwise be chosen freely. To fix them entirely, we choose the centre-of-mass frequency to be the same initially as at the half-way point of the protocol: Ω − (0) = Ω − (t f /2) = Ω − (t f ). This then allows the determination of α 0 , β 0 as well as the boundary values of the normal mode frequencies Ω 0± ≡ Ω ± (0), Ω in± ≡ Ω ± (t f /2) by use of equations (5)- (7). Note that the initial curvature ω 0 is then given by In this way, the beginning, midpoint and final states of the protocol are dictated solely by these physical constraints, which consist of the desired initial and intermediate ion distances d 0 and d in , the maximal quartic confinement β max that can be achieved in a given trap and the ion mass m.
The remaining task of designing suitable protocols is then to provide a transition between α 0 , β 0 and α in , β in , and then back again. We make a working assumption that smooth transitions will provide the lowest excitations and use the symmetry of the problem to confine our search for protocols which are also symmetric in time around the point t f /2. To design protocols which achieve short execution times, we make use of so-called shortcut-to-adiabaticity techniques, which are described in the next section.

Inverse engineering of a shortcut to adiabaticity
The goal of STA techniques is to take a Hamiltonian H(t) and design it in such a way that the populations in the instantaneous bases at t = 0 and t = t f are the same. This thus allows tasks to be executed in arbitrarily short times while yielding the same final result as an adiabatic evolution. STA methods have been proposed for many applications in trapped-ion QIP [14,15,19,[25][26][27][28]. Of the various STA methods, we chose to use invariant-based inverse engineering as described in [29]. This is an approach that involves first designing dynamical invariants I(t) which commute with a general form of the system Hamiltonian H(t), and then deducing the explicit form of the latter from the resulting conditions. Dynamical invariants I(t) are operators with constant expectation values where |Ψ are solutions to the Schrödinger equation for H(t). For a Hamiltonian H HO (t) with pure harmonic oscillator form, such invariants are known explicitly. This allows for the use of a result due to Lewis and Riesenfeld [30], which states that if the Hamiltonian H(t) and the corresponding invariants I(t) are known, the individual solutions |Ψ to the Schrödinger equation can be given as a superposition of eigenvectors |n; t of I(t): where the coefficients c n are constant and the phases α n (t) are fully determined. This result is used in the invariant-engineering approach according to the following reasoning. If the invariant commutes with the Hamiltonian ([H HO (t b ), I(t b )] = 0) at boundary times t b = {0, t f }, they share an eigenbasis then. This means that the initial populations in the eigenbasis of H HO (0) are also the populations of the invariant eigenvectors |n; 0 . Since Lewis-Riesenfeld theory tells us that the population numbers c n are constant, the system is still in the superposition n c n e iα n (t f ) |n; t f at the final time, while the eigenvectors of the invariant |n; t f have evolved. Since the Hamiltonian and the invariant commute again at t f , they again share an eigenbasis, meaning that the populations in the initial instantaneous basis of H HO (0) have been carried over to the new basis of H HO (t f ). This yields an evolution that recovers the initial populations at the final time. Note that the system may generally not follow the adiabatic evolution at intermediate times, but reaches the same final state nonetheless. No invariant is known for the Hamiltonian (2). Thus in order to make use of STA methods, we first make an approximate transformation to the set of dynamical normal modes with frequencies Ω ± , for which the Hamiltonian takes the form of two harmonic oscillators H 2HO = H (+) HO + H (−) HO at all times. The procedure is given in detail in appendix A.
We find the Lewis-Riesenfeld invariants I ± for each dynamical normal mode to be where X ± and P ± are the transformed position and momentum coordinates and the auxiliary functions ρ ± and q ± are defined byρ The physical meaning of the auxiliary variables q ± can be understood as the normal mode centres, while ρ ± correspond to the effective curvatures of the oscillator potential in the transformed co-ordinates. As an example, in the case of transport of one ion in a constant trapping well the function ρ can be set to 1, as the ion experiences the same potential curvature at all times. In the case we are currently considering, q − is zero at all times due to the spatial symmetry of the potential.
In order to make use of the invariants, we must first find a parametrisation which satisfies the commutation relation H (±) HO , I (±) at the boundary times t b = {0, t f }. The commutation can be ensured by setting the boundary conditions on the auxiliary functions. The conditions on the zeroth and first derivatives arise from the Lewis-Riesenfeld theory (see appendix A for further details). In addition, we constrain the first two derivatives of the distanceḋ(t b ) andd(t b ) to zero, so that the scheme starts and ends with stationary ions and to minimise the energy. This together with the auxiliary equations (15) and (16) leads to the remaining conditions. To fulfil the physical constraints at the midpoint of the protocol, the normal mode frequencies need to reach Ω in± at t f /2, which is achieved approximately by setting This can be verified by inserting the condition into (15) and neglecting the term inρ ± , which is justified as long as the protocol takes several motional cycles. Any choice of {ρ ± , q ± } satisfying the boundary conditions above leads to a shortcut to adiabaticity with respect to the Hamiltonian H 2HO , leaving flexibility to optimise the scheme for various purposes, such as cancelation of residual excitations or robustness to experimental imperfections. However choosing {ρ ± , q ± } to simultaneously fulfil the ODEs in (15)- (20) is hard. Therefore we follow [14] and design a general form for ρ ± which satisfies only (15), (18) and (19), but contains additional free parameters. We then perform a numerical search in the free parameter space in order to obtain solutions which yield minimal excitations.
We expect that a smoothly varying function for ρ ± (t) would be the most satisfactory experimentally. For this reason, we use a polynomial interpolation function for ρ ± , which contains only even orders due to the chosen symmetry of the protocol. The simplest polynomial satisfying the boundary conditions for the centre-of-mass mode is ρ − = 1 (22) due to having chosen Ω 0− = Ω in− earlier.
Setting ρ + (t) to be a polynomial of order 14 leaves two free parameters for numerical optimisation after fulfiling all the constraints. The resulting parametrisation for ρ + is with the normalised time s = t/t f and the two free parameters A and B. We choose the former to constrain the curvature of ρ + at s = 1/2, while the latter is the co-efficient of the 14th order term. The protocol is now fully defined. The auxiliaries ρ + and ρ − are inserted into (15), yielding an expression for the mode frequencies Ω ± . From these, the potential parameters α(t) and β(t) can in turn be obtained from the relations which are found by rewriting (5)- (7). In the following section, we utilise the free parameter A to find optimised trajectories that do not create residual excitations when being executed on two ground-state ions. This is a prerequisite for achieving motional energy swapping. Furthermore, by using the parameter B as well, the protocol can be made robust to common experimental imperfections, for which the necessary steps are worked out in section 5.

Numerical optimisation of the shortcut for low residual excitations
Two obstacles stand in the way of analytically designing a perfect shortcut to adiabaticity for the given double-well potential. First of all, the ansatz for the auxiliary function ρ + does not a priori guarantee that the remaining boundary conditions (20) are fulfilled. But more importantly, the shortcut was designed for the harmonically approximated Hamiltonian H 2HO and can therefore not result in an excitation-free scheme for the full Hamiltonian.
This motivates that an optimal ρ + should minimise the additional excitation of the full Hamiltonian, instead of adhering to the approximate picture. In this section, we thus attempt to optimise the ansatz for ρ + in this sense and numerically determine the parameter A. In this algorithm, we only consider ions that start in their motional ground state, expecting that the protocol also adds close to minimal excitations for ions with small initial motional energy.
In appendix B, we additionally give the results of optimising the ansatz such that the boundary conditions are fulfilled, more in line with previous work [14], but find the performance inferior to the aforementioned method. While the method presented here is numerically more costly, we found that it was not a significant problem for the trajectories which we considered.
Note again that, once an optimised protocol is found, we will aim in section 6 to find run-times t f that yield a complete energy swap. This means that we now need to optimise the shortcut ansatz separately for a range of run-times, for each of which different values for the parameter A may be found.

Numerical prerequisites
To study the effectiveness of the optimisation approach, we need to compute the final energy of the ions after a given protocol. Numerical integration of the classical equations of motion given by the full Hamiltonian in (2) is used to calculate this. However, the total energy increase cannot be divided between the individual ions, as they are coupled by the Coulomb interaction. We therefore expand the Coulomb potential in (2) to first order in the displacement of the ions from the equilibrium positions and calculate the energy excitation of one of the ions at time t f with respect to the equilibrium energy as where i = 1, 2 denotes the ion. The positions x i (t) and momenta p i (t) of the ions are found by integrating the equations of motion. Note that when using the symmetric electrical potential (3) and symmetric starting conditions, both ions have the same energy increase E ex,1 = E ex,2 . This will not hold in later sections when we consider asymmetric potentials. In this section, the ions are assumed to be in their ground state at t = 0, which is implemented by setting their initial momentum to zero and placing them at the positions of the trapping potential minima. The physical constraints are shown in table 1 and reflect a realistic experiment. The achievable β max (and with it, the critical distance d c ) is that of the Sandia HOA2 surface trap [31] recently used in several research groups [32][33][34], including our own. It is obtained from simulations of the trap together with the assumption that potentials with a voltage range of ±10 V are supplied to the electrodes.
The inner distance d in was chosen slightly above the critical distance d c , such that the ions never get merged into a single well. The outer distance d 0 was chosen such that the exchange interaction is slow compared to the targeted protocol run-time. At the outer distance given in table 1, the exchange time t e is about 442 μs. As the initial curvature ω 0 is about 0.45 MHz, this corresponds to 200 trap periods.  The excitation E ex,1 of one of the ions is normalised by the energy ω 0 of one phonon, where ω 0 is the initial potential curvature. The dotted line in a) corresponds to final excitations of 0.1 quanta and marks the energy level below which the scheme can in principle yield cooling of one of the ions to better than 0.1 quanta, which we consider the target. Note that E ex,1 = E ex,2 holds here due to the symmetry of the problem.

Minimisation of residual excitations
The results are shown in figure 2. For each run-time, the free parameter A has been optimised and the final excitation E ex,1 (t f ) is plotted. Numerically, the optimisation is performed with the Nelder-Mead algorithm provided by the NumPy package in Python with E ex,1 (t f ) being the cost function. We are able to reduce the excitation to the level of numerical tolerance and have thus demonstrated a way to numerically construct a shortcut to adiabaticity and perform the proposed protocol such that no additional excitations are added.

Specific example trajectory
For illustration, the parameter functions of a protocol optimised in this way is shown in figure 3. We pick the trajectory with a run-time of t f = 30 μs. Shown is the time evolution of the auxiliary functions ρ ± (t), the normal mode frequencies Ω ± (t), the potential constituents α(t) and β(t) as well as the resulting distance trajectory d(t). Note that the centre-of-mass mode frequency Ω − is constant as designed and that the physical constraints specified in table 1 are met. Despite the high-order ansatz, the optimised trajectories are smooth.
We could now go on to explore how to achieve cooling. However, the protocols found in this section are susceptible to experimental imperfections. We therefore defer the analysis of cooling solutions to section 6 and first address experimental robustness in the following section.

Robustness optimisation of the shortcut
A commonly encountered experimental imperfection in ion traps is an additional homogeneous electric field [35,36], which could arise due to a number of sources, such as stray charges in the trap or inaccuracies in the applied trap voltages. It can also be regarded as the first order expansion of a general perturbation potential. In this section, we modify the previously obtained schemes to achieve optimised robustness towards this type of experimental imperfection.

Effect of a homogeneous field perturbation
Adding an additional homogeneous field γ to the double-well potential (4) yields Variables of the perturbed system are indicated with a tilde in this section. The equilibrium positionsx (0) i of the ions are not symmetric anymore and we introduce the centre-of-mass shifts to denotex (0) 1 =s −d/2 andx (0) 2 =s +d/2. The ion distance is given byd, which may differ from the unperturbed d. In what follows we assume that the additional field is weak, such that the shifts is much smaller than the equilibrium distance, allowing us to include it only to first order. Under this assumption, we find thatd ≈ d ands To describe the strength of the imperfection in a way that allows a dimensionless small parameter expansion, we define the perturbation parameter which turns the linear approximation condition |s| d into η = |s|/d in 1. This assumption can be validated through reference to literature values. In recent work on ion separation [36], the residual field after compensation was found to be 0.1 Vm −1 . Together with a critical distance of about 25 μm and a mode frequency of 0.174 MHz, a value of η = 0.008 is obtained, justifying the use as a small perturbation parameter.

Frequency mismatch
Here we should point out that an additional homogeneous field does not only cause excess excitations, but also introduces a mismatch between the potential curvatures at the two ions due to the quartic term in the potential. From the perturbed normal mode calculation in appendix B, one can calculate the shifted curvaturesω i due to a perturbation η as making the mismatch δω =ω 2 −ω 1 ≈ −12βdd in η/ (mω i ) to first order. This is problematic since the motional exchange (1) is a resonant effect with respect to the potential curvatures, meaning that such a homogeneous field can suppress the desired energy exchange. Note that no protocol optimisation can mitigate this off-resonance effect. We deduce the required level of experimental accuracy that this imposes on the protocol in section 7.

Numerical optimisation for optimal robustness and low residual excitations
Since a real experiment will never be able to implement the protocols found in section 4 perfectly, we aim to modify the numerical optimisation to find protocols that incorporate robustness. The cost function for which the algorithm finds optimal protocols should now also contain E ex,i (t f , η), the final excitation of ion i when executing a protocol perturbed with strength η. We thus define  (32). The final excitation E ex,1 (t f , η) of one of the ions is calculated on a grid of run-times and perturbation parameters, assuming the ions to be in their ground state initially. The physical constraints for all shown schemes are those from table 1. The perturbation is also shown in terms of the homogeneous field γ, which is obtained from η using (30). Note that plotting E ex,2 results in the same data, but mirrored along the line η = 0. The inset above (b) shows a cut through the data (indicated by the white dashed lines), corresponding to the unperturbed final excitations E ex,1 (t f , η = 0) together with a fit to that data and an envelope function. The dotted line in the insets corresponds to final excitations of 0.1 quanta and the hatched area marks the range of run-times where the envelope indicates final excitations above this level.
a combination of the residual excitations caused by both the perturbed and the unperturbed protocol. We chose a value of η = 0.015 which, at the design parameters in table 1, corresponds to an electrical field value of about 0.77 Vm −1 . While a field compensation with a resolution of 0.1 Vm −1 was performed in [36] by biasing nearby electrodes, this strongly depends on the parameters of the trap and electronics equipment. It is challenging to state an absolute limit of control, for example electrodes further away from the ion could be used for compensation, in principle yielding even finer resolution.
To calculate the final excitation E ex,i (t f , η) of ion i, the equations of motion of the full Hamiltonian are integrated as before. The perturbed potentialṼ el is inserted into the full Hamiltonian (2) and the final excitations calculated as in section 4.1. The initial conditions are defined by placing the ions at rest at the equilibrium positions of the perturbed potential. Note that E ex,1 = E ex,2 does not hold anymore for non-zero η due to the asymmetric potential. As before, we use the Nelder-Mead algorithm to find the minimum, this time by varying both free parameters A and B. As in section 4, we re-perform this optimisation for each considered run-time t f .
The resulting performance is displayed in figure 4(b). For each run-time t f , the optimisation was run to obtain the optimal parameters A and B. Then, using these values, the final excitations E ex,1 (t f , η) are calculated for a range of perturbation parameters η at each run-time. This yields the plotted excitations on a grid of run-times and perturbation parameters. For comparison, figure 4(a)) shows the same plot based on the results of the non-robust optimisation presented earlier in figure 2.
The comparison between the non-robust and robust methods shows that the latter offer a clear improvement, yielding a much larger area where final excitations are at an acceptable level for our scheme. However, the robust method produces exponentially growing excitations as run-times are shortened, giving rise to a minimal time beyond which the scheme cannot be sufficiently optimised anymore.
This breakdown can be seen in the inset of figure 4, which shows the final excitations E ex,1 (t f , η = 0) for unperturbed protocols (cut indicated by a white dashed line). While the envelope of the excitations increases exponentially, periodic minima are are also present. To quantify the run-time beyond which the robust protocols do not yield optimal solutions, the excitations at η = 0 are fitted with the function f (t) = a exp(−bt)sin 2 (ct + d) . The non-periodic part a exp(−bt) can then be interpreted as an envelope function and its intersection time with an energy level ofn = 0.1 (indicated by a grey dotted line) is then used to estimate the run-time below which the protocol is not well optimised (indicated by a hatched area). We refer to this time as the critical time T crit and it is found at T crit = 14.2 μs, equivalent to ∼ 6.4 motional cycles. Note that the definition of T crit is a conservative one, as timings resulting in negligible excitations can be found below T crit due to the periodic nature of the excitations. Another interesting feature of both robust methods is the existence of stripes with low excitations at constant, periodic run-times, marking configurations where the scheme is ultra-robust even against strong perturbations.
The fact that the robust method works well despite the simple choice of the cost function is encouraging, as it is plausible that it can be improved in similarly simple ways. One could for example use more free parameters and choose more values of η in (32), such that the robustness range is increased. Additionally in appendix B, we have constructed a method where robustness is added by taking the perturbation into account in the dynamical normal modes, but this was found to perform worse.
We conclude that the presented robustness optimisation method is a useful tool to make this STA protocol able to withstand experimental imperfections, at the cost of introducing a limit to the achievable run-times. Having gained sufficient control to perform it in around 6 trap cycles without excess excitations, we apply this result to find timings where the ion motional excitations are swapped.

Cooling characteristics
After gaining the ability to optimise the protocol for low residual excitations and optimal robustness over a range of run-times when executed on cold ions, the existence of cooling solutions can now be demonstrated. To this end, we deploy the protocols that were optimised for robustness in section 5 and now apply them to the case where one ion is initially hot. In this scenario, motional energy is exchanged between the ions and we therefore expect to find a run-time where a complete energy swap takes place. This timing can be found by calculating the final energy of the initially hot ion after the protocol for a range of run-times. The desired cooling solution is then given by the trajectory at the run-time where this final energy reaches a minimum, indicating an energy swap.
In this section, we only consider the ideal case where schemes are executed with no additional electric field (η = 0). The effect of additional homogeneous fields on the energy exchange is then explored in the subsequent section.
We also explore the effect that varying the inner distance d in has on the timing of the energy swap. As the motional exchange scales strongly with the ion distance, we expect to find cooling solutions at shorter run-times when bringing the ions closer together. As we have chosen to consider only protocols that leave the ions in separate potential wells, the available quartic confinement in the trap under consideration limits the available range of d in > d c . We thus go on to consider how to decrease cooling times by varying β max .

Numerical prerequisites
The energy of each ion is calculated in the same way as in (27) and section 5.3. However, now we also need to take into account that the ions do not necessarily start in their ground state. The initial conditions are given by setting p i (0) and x i (0). This is done by choosing the initial energy E in,i and distributing it onto the kinetic and potential parts of the energy by choosing the initial motional phase φ. A phase φ = 0 is understood to mean that all initial energy is kinetic and thus x i (0) is equal to the equilibrium position, while φ = π/2 implies that initially all energy is stored in the potential and thus p i (0) = 0. To subsume the effect of φ, we define the average energyĒ ex,i (t f , η), which is obtained by averaging the resulting E ex,i (t f , η) for φ ∈ [0, 2π]. For numerical examples in this work, 25 uniformly distributed values of φ are used.
One ion being initially hot is defined here to mean that it is initialised with a motional energy corresponding ton = E in,i / ω 0 = 10, which is on the order of the Doppler limit (around 20.5 quanta for 40 Ca + at a mode frequency of ∼0.48 MHz [10]).

Cooling solutions
We now go on to demonstrate the energy exchange, while varying the inner distance d in from 1.0d c to 1.25d c . Otherwise, the physical constraints are kept identical to those chosen in section 4 (see table table 1), namely β max = 0.85 × 10 −3 Nm −3 , d 0 = 5d c = 70.1 μm and the mass of a 40 Ca + ion. To obtain optimised trajectories, the ansatz (23) is constructed for each resulting set of physical constraints, leaving the free parameters A and B. Optimal values for these, which depend on the run-time t f , are then found by applying the robustness cost function to initially cold ions as described in the previous section.
Finally, for each run-time, the obtained optimised trajectories are applied to two ions, one of them initialised ton = 10 and the other one being in the ground state. The final excitationsĒ ex,1 of the initially excited ion are then calculated and the results are shown in figure 5. We observe that the final energy has a minimum far belown = 0.1 for each shown value of d in except for d in = 1.0d c . We refer to the run-time of these minima as the cooling time T c . This run-time of complete energy exchange tends to lower values as the ions are brought closer together, as predicted from (1). The final excitations at T c differ from zero due to the schemes causing a finite energy increase even for ground-state ions, as well as due to the simulations only being performed for discrete t f .
By decreasing the minimal distance d in , the cooling time is eventually lowered into a range of run-times where the adiabatic shortcut cannot be perfectly optimised anymore. The critical time T crit was introduced in section 5 to quantify the onset of this regime. After optimising the trajectories applied in figure 5 using the robustness cost function (32), T crit is calculated for each of the choices of d in and indicated in figure 5 by the black markers. For the lines where d in 1.1d c , the cooling time comes to lie above T crit . When bringing the ions closer together and thus reducing T c below the critical time, one might expect that the energy swap is not completed to a satisfactory level anymore, due to the adiabatic shortcut potentially adding excitations above 0.1 quanta per ion. This does however not occur in figure 5 until d in = 1.0d c , which is explained by the conservative definition of T crit . As one recalls from figure 4, the excitations of two ground-state ions do not increase monotonically with decreasing run-time, but show periodic minima indicating near-perfect shortcuts. However, T crit is defined using the envelope of said excitations, disregarding the minima. It is then indeed the case for d in = 1.05d c that the timing of such a trajectory coincides with the timing required for an energy swap, allowing cooling to below 0.1 quanta. Only for d in = 1.0d c is this not true anymore. To visualise this effect, we have also plotted the excitations E ex,1 (t f , η = 0) added to an ion initially in the ground state for these two choices of d in (dashed lines).
As the final excitationĒ ex,1 is an average over many initial motional phases, we check if this aggregation is justified. We plot the dependence of the non-averaged excitations E ex,1 on φ in figure 6(a), the example being the protocol with d in = 1.1d c at the run-time t f = 16.6 μs, which corresponds to the timing of the cooling solution. The dependence has a sinusoidal shape around the average, making it a reasonable replacement.
Furthermore, the final excitation should not depend strongly on choosing an initial energy that differs from the exemplaryn = 10. This is confirmed in figure 6(b), whereĒ ex,1 is shown to be well belown = 0.1 for initial excitations ranging up to the Doppler limit.
The existence of trajectories that provide almost complete energy swapping is thus demonstrated, with cooling times T c on the order of 10 motional cycles. While varying d in is a viable option to fine-tune the cooling time, the maximal speed is fundamentally limited by two factors. First, as seen in section 5, the integration of robustness into the protocol optimisation introduced the critical time T crit , below which no satisfying adiabatic shortcut is found. This limit holds despite the existence of viable trajectories at shorter run-times under certain conditions. Secondly, the maximal quartic confinement β max dictates d c and thus the available range of d in .
Note that in this discussion the intrinsic heating rate in the trap has not been taken into account. This would excite both ions throughout the execution of the protocol, preventing the cooled ion from reaching an arbitrarily low motional state. Depending on the parameters of a given apparatus, heating rates from a few tens to tens of thousands of quanta per second have been observed [37,38]. At the lower end of this range, an ion would only be heated on the order of 10 −4 quanta during the protocols discussed here.
We will therefore go on to analyse the behaviour of the protocol when scaling β max in addition to d in /d c , expecting that cooling can be achieved faster in a trap that provides a larger quartic confinement.

Scaling behaviour with the maximal quartic confinement
To study the effects of scaling β max in an easily comparable way, we go over to express all quantities that describe the protocol in a dimensionless form. All energies are expressed on a scale of ω 0 as before, while values denoting time can be made dimensionless by dividing by a motional cycle (2π)/ω 0 . The distance trajectory d(t) can be made dimensionless through division by the critical distance d c , which is solely dependent on β max . To set the physical constraints, one then picks d 0 /d c (= 5 in our case) and d in /d c > 1. Rewritten as such, protocols that share the same d 0 /d c and d in /d c can be compared across varying β max .
To determine how the motional exchange rate Ω throughout the protocol depends on β max and with it the timing of the cooling solutions, we make Ω dimensionless by division with ω 0 and obtain We show in appendix C that this expression scales as ω 0 √ ω 1 ω 2 d 3 ∝ (β max ) 0 . We then expect that for a given choice of {d 0 /d c , d in /d c }, varying β max leaves Ω/ω 0 invariant and thus the cooling times can always be found at the same number of motional cycles.
To demonstrate this result, we pick four values of β max spanning three orders of magnitude. For each of these maximal confinements, the inner distance d in /d c is varied from 1.05 to 1.25. The outer distance is chosen to be d 0 /d c = 5 as in all examples so far. For each combination of {β max , d in /d c }, the protocol is optimised for a range of run-times using the robustness function as described in section 5.3 and the critical time T crit is determined. We then determine the run-time T c leading to a cooling solution in the same way as in the previous subsection. Note that for β max = 1 · β HOA , these computations are equal to the ones performed for figure 5.
The cooling times and critical times are plotted in figure 7(c) depending on the dimensionless distance. The red symbols show the cooling times T c obtained for all the combinations of {β max , d in /d c }. No dependence on β max is observable as expected. Choosing a smaller inner distance leads to a faster overall exchange, as already discussed in figure 5. The critical times T crit are depicted by the blue symbols and only a weak dependence on both the quartic confinement and the inner distance is observed, allowing us to state that the robustness cost function produces well-optimised trajectories down to 6-7 motional cycles.
As already observed in figure 5, cooling solutions are also found which have shorter duration than the critical time of the robustness cost function T crit . This is shown in detail again in figures 7(a) and (b) through the example of the trajectories corresponding to β max = 1 · β HOA and d in /d c = 1.05. Figure 7(a) depicts the final excitations of an ion that was initially excited by 10 quanta, depending on the run-time. Figure 7(b) on the other hand shows the excitations caused by the protocol when both ions are initially in their ground state, together with a fit to the simulated data and the excitation envelope as described in section 5. It becomes apparent that this particular cooling solution is enabled by the fact that the excitations in figure 7(b) show a minimum at 6.3 motional cycles, despite the envelope already indicating far higher excitations.
The results presented in figure 7 provide an easy way to select the parameters of a desired cooling protocol. The energy of two ions can be swapped at best within 6-7 motional cycles. The ideal choice of the inner distance is then also immediately clear as being ∼ 1.05d c , providing a cooling solution at ∼ 6.3 motional cycles. As calculated in appendix C, the trap period 2π/ω 0 scales with β −3/10 max and the same is then true for the cooling solutions T c in absolute time. However, also consider that we have chosen d 0 /d c = 5 in all presented examples. If a larger value is selected in a given implementation, we would expect the cooling minima and the critical times to tend to longer run-times, as the ions spend more time far away from each other and are more strongly accelerated.

Suppression of motional exchange by an additional homogeneous electrical field
As noted in the introduction, the motional exchange is a resonant effect with respect to the ion frequencies.
Since an additional homogeneous electrical field introduces a mismatch of the potential curvatures of each ion as calculated in subsection 5.2, this resonance condition is not perfectly observed anymore. We therefore proceed to estimate the range of tolerable perturbations η such that cooling is still achieved.
Due to being a resonance effect, we choose an ansatz in which the final excitation E ex,1 of the initially hot ion has a Lorentzian shape with respect to the curvature mismatch δω: where E in,1 is the initial ion energy, Ω is the exchange frequency (1) and δω the frequency mismatch (31), while the factor k is a constant that is to be determined. From this, the term k δω Ω is calculated as where we have replaced the time-dependent values d and β by their values at t f /2, under the assumption that most of the energy exchange happens at that point in time. We understand η 1/2 to be the Lorentzian half-width (HWHM) with respect to the perturbation η. Using (9), we obtain the simple expression which depends on k and on how close one brings the ions to the critical point, but not on β max itself. This means that the required experimental accuracy in terms of the perturbation parameter η cannot be reduced by using a trap with larger quartic confinement. However, it can be maximised by choosing a minimal value for d in /d c . Thus the goal of swapping the motional energy as fast as possible is compatible with finding a cooling solution that is maximally insensitive to additional homogeneous fields.
To illustrate these results, figure 8 recreates the robustness plots from section 5, but with one ion initially excited by 10 quanta. This is done for the two inner distances yielding the fastest cooling, d in /d c = 1.05 and d in /d c = 1.10. For each run-time shown, we apply the trajectory that was optimised using the robustness cost function for figure 5, but with a range of perturbations parametrised by η. For each η, one ion is initialised with an excitation of 10 quanta and the equations of motion are calculated while keeping the perturbation constant throughout the scheme. The resulting final excitationĒ ex,1 (t f , η) is then calculated in the same way as described in section 6.1 and plotted on a grid of run-times and perturbations. Note that E ex,1 is again an average over the initial motional phase. The red contour line shows the area where the final excitation of the initially hot ion decreases belown = 0.1. For comparison, we overlay another contour line in white, marking the area where the trajectories, applied to ground-state ions, cause a final excitation of the ion below 0.1 quanta. To show the resonance effect, we take a cut though the data at the cooling time T c = 14.2 μs and T c = 16.6 μs, respectively, which are familiar from figure 5. The cuts are indicated by the white dashed lines. This yields the insets on the left, showing the excitationĒ ex,1 depending on the perturbation.
To the data of these cuts, we fit the Lorentzian (34). Specifically, we insert (35) into the Lorentzian and find the optimal value for the remaining parameter k by applying the Python method curve_fit on the central part of the data that fulfilsĒ ex,1 / ω 0 < 2. The resulting fit functions are also displayed in the left insets in figure 8, showing a good match to the data in the central peak. A value of k = 0.679 was obtained in (a), corresponding to η 1/2 = 0.048 and k = 0.678 in (b), corresponding to η 1/2 = 0.038. In order to assure cooling to better thann = 0.1, we therefore find that η must not exceed values of ±0.0048 and ±0.0038 respectively. The scheme that brings the ions closer together thus yields better tolerance to stray fields in addition to faster cooling, confirming (36). We repeat this analysis for d in /d c = {1.15, 1.2, 1.25}, obtaining values of k close to 2/3 as well. Thus we can regard the Lorentzian ansatz and the form of η 1/2 in (36) as a reasonable estimate for the resonance behaviour close to η = 0, despite the crude approximation that led to this result.
Contrary to η 1/2 , the corresponding homogeneous field strength γ 1/2 does scale with β max . From the conversion given in (30) and the results of appendix C, we obtain that The tolerable field strength thus increases with the quartic confinement, whereas η 1/2 is invariant. This allows us to compare this result to experimental demands in other experiments. β max scales with the overall trap size scale a as β max ∝ a −4 [24]. The tolerable value of η also strongly depends on the desired initial and final excitation level. This can be seen by rewriting (34) and (35) to For example, if one intends to cool from 1 to 0.1 quanta instead of from 10 to 0.1 quanta as shown in figure 8, the tolerable range of η increases by a factor of 3.3. The nature of the motional exchange imposes much stricter conditions on the maximally tolerable perturbation η than the in-and-out transport of two ground-state ions, as can be seen from the contour lines in figure 8. To achieve the desired motional exchange in our example, η cannot exceed values of ±0.0048, corresponding to an electric field of ±0.22 Vm −1 . Previously, a stray field calibration in steps of 0.1 Vm −1 was reported [36], although in a trap where the available β max was about 17 times lower than in the one considered here. As noted before, the achievable accuracy in compensating such perturbations depends on the details of the experimental apparatus and the chosen approach and a definitive limit is challenging to state. Our simulations indicate that the resonance condition imposes a need for control of the trap voltages on a level of tenths of mV −1 and with sub-μs resolution. Shifts of the higher-order potential parameters α and β are less constraining, allowing fractional deviations on the order of 1%.
Note finally that the available range of parameters also depends on the choice of the optimisation cost function (32). In a previous iteration of this work, we had optimised the schemes with ground-state ions using η = 0.03 in the cost function (32). While this led to additional robustness for two ground-state ions, the critical time T crit was found at a larger value of 10 motional periods. Consequently, the current best solution where the energy is completely exchanged after 6-7 motional periods was out of range. Slower protocols with larger inner distances had to be chosen. Due to the relation (36), this also led to reduced robustness when one ion is initially excited. It also seems of little use to optimise the robustness of the protocol for unexcited ions to η = ±0.03 when the resonance effect discussed in this section then limits the acceptable perturbation to around |η| < 0.005.
In general, we believe that there is an optimal choice of η, depending on the targeted result. Choosing a smaller η in the cost function will allow accessing shorter run-times and inner distances, thus increasing η 1/2 . However, the range of acceptable perturbations due to the resonance condition will be reduced by the now diminished robustness. On the other hand, choosing η too large will reduce the range of available run-times where the protocol behaves robustly. Since we chose to not combine the ions into a single well, we believe the choice of η = 0.015 in (32) to be reasonable, as it allows us to access protocols where the ions are almost brought to the critical distance.

Cooling ions of unequal mass
Thus far, we have only considered exchange protocols where both ions are of equal mass. However a number of recent works have used ion chains containing ions of different mass [3,[39][40][41][42][43], justifying the need to implement cooling in such configurations. We therefore aim to generalise the cooling protocol to ions of unequal masses m 1 = m 2 in this section.
When trying to extend the cooling scheme to ions of unequal masses m 1 = m 2 , a fundamental problem arises: if we keep using a symmetric potential such as the harmonic-quartic double-well potential (4), then the ions will always sit at symmetric equilibrium positions. The frequencies resulting from the potential curvatures given by (8) are then not equal, as the second derivative of the potential at the ion position is the same for both ions, but the mass is not. The resonance condition ω 1 = ω 2 for motional exchange to take place is thus violated. As noted in [44], we are furthermore unable to find decoupled dynamical normal modes for the unequal mass case when using a harmonic-quartic potential. This is mitigated by introducing an asymmetric term to the potential V el , which can be used to force the motional frequencies to be equal throughout the scheme. The simplest choice is to add a linear term leading to the trapping potential The linear term is reminiscent of section 5.1, where such a linear term appeared as an undesirable perturbation. Here in contrast, the additional field is intentional. The dynamical normal modes are derived in the same way as in the equal mass case and a detailed account is given in appendix D. The full Hamiltonian is again separated into two harmonic oscillator Hamiltonians and we find that the auxiliary functions ρ ± need to fulfil (15) and the commutation is ensured by the same boundary conditions (18)- (20) as before. If the COM-mode frequency is chosen to be constant again, this means that the ansatz for ρ + given in (23) can be reused, as well as ρ − = 1. The protocol is then defined by ρ ± and (15) as In the same way as in the equal mass case before, we now optimise the total final excitations E ex calculated with the full Hamiltonian by varying the free parameters in ρ + . Figure 9(a) shows the final excitations of two ground-state ions after optimising only A vs optimising both A and B for each run-time. The physical constraints are chosen to be as in table 1, but with the second ion (coolant) being only half the mass of the first, which is chosen to be that of 40 Ca + .
When only using one parameter, the excitations increase exponentially with decreasing run-time despite using the full Hamiltonian, while negligible excitation levels are reached when optimising both A and B. This can be compared to the results for non-robust protocols with ions of equal mass in section 2, where only one free parameter had to be optimised to achieve comparably optimal trajectories. Since q − = 0 does not hold anymore in this version of the scheme, the parameter search needs to satisfy the boundary conditions for both q + and q − and two free parameters are required.
The achievable motional exchange in the unequal-mass case is demonstrated in figure 9(b). The first ion is again chosen to have the mass of 40 Ca + , while the mass m 2 of the second ion is lowered. The two-parameter optimisation is run for a wide range of mass ratios of m 1 /m 2 = {1.25, 2, 5, 10}. Shown is the final energyĒ ex,1 of the first ion, initially excited by 10 motional quanta, exhibiting minima where the two ions exchange their motional energy almost completely. The cooling solutions are shifted to shorter run-times with increasing mass ratio. This is explained by the mass dependence of the exchange frequency Ω.
We have therefore demonstrated the existence of trajectories that provide almost complete energy swapping for two ions with mass ratios up to 10. Note that while no robustness towards experimental imperfections has been built into the scheme in this work, this could be generalised in a straightforward way from section 5.

Conclusion and outlook
Our work demonstrates the possibility to perform fast cooling of a trapped ion by transient interaction with a pre-cooled ancillary ion. For 40 Ca + ions, we obtain trajectories which achieve this in 6.3 trap periods, corresponding to 14.2 μs for a trap frequency of 0.44 MHz in a currently operated ion trap. The required electric field control is at levels similar to those reported in recent experiments. Similar methods to those reported here would be applicable to cooling of radial modes of motion, although it appears challenging to simultaneously cool both axial and radial modes.
In addition to direct application to cooling, dynamic resonant swapping of motional states would also be a key ingredient for quantum information using oscillator codes rather than internal state qubits [45]. In addition to the requirements above, this would require consideration of the phase control of the oscillator through the transport.
, (A.10) of the potential V tot , from which we obtain upon inserting the equilibrium positions ±d/2. The eigenvalues of K are given by the squares of the dynamical normal mode frequencies Ω ± , with the eigenvectors The coordinate transformation yielding decoupled Hamiltonians is given by and includes the dynamic componentḋ. In terms of these coordinates the Hamiltonian can be written as which also includes the dynamic componentd. We see that both of these have the same form as the harmonic oscillator Hamiltonian H HO obtained in (A.2), once with F(t) = − m/2d and once with F(t) = 0 and the mass M set to 1.
Applying the results of the previous subsection, we obtain the Lewis-Riesenfeld invariants with the auxiliary functions ρ ± and q ± , defined bÿ The mode energies are given by Note that we have defined here the part E (+) q of the stretch-mode energy that pertains to the auxiliary function q + , such that we may minimise its influence later on.

Appendix B. Optimising the shortcut of the approximate Hamiltonian
In addition to the method presented in section 4, where protocols are optimised such that the residual excitations of the full Hamiltonian are minimal, we have also attempted to complete the shortcut to adiabaticity by fulfiling all boundary conditions (18)- (20). This method has previously been applied to obtain optimal trajectories for trapped-ion schemes [14,46].
In the case where no robustness is considered, one way of doing so is to minimise the part E (+) q (t f ) of the mode energy in (A.20) pertaining to the auxiliary q + , thus finding parameters that come close to fulfiling (20).
The results are shown in figure B1, overlaid over the results of the full-Hamiltonian optimisation presented in figure 2 for comparison. In this case, E ex,1 (t f ) is not zero as it should be if a shortcut has indeed been found, restricting the accessible run-times. Instead the excitations increase exponentially for shorter run-times, consistent with the behaviour reported in [14]. This is due to a break down of the harmonic approximation used in constructing the dynamical normal modes.
In figure B1(b), the optimised parameter A depending on the run-time is depicted. The two optimisation approaches yield similar values of A, yet the difference in final excitations grows to more than ten orders of magnitude at a run-time of 15 μs. While the optimised parameters differ by several percent, the parameter A can only vary by a fraction of around 10 −3 if the final excitations are not to exceed 0.1 quanta. The periodic jumps in the optimal A are a consequence of the numerical optimisation. There are periodic 'valleys' of minima in the A − t f -plane that seem to be spaced roughly by a trap period. However, allowing the algorithm to follow one of these 'valleys' for the entire range of t f eventually leads to negative values of A, which has undesirable physical consequences, such as there being multiple turning points in the ion trajectories. We thus choose to always start the optimisation from the same initial value of A, which makes the algorithm find the nearest minimum and leading to the periodic jumps in A.  figure B1 and plot the final excitations E ex,1 (t f , η). To obtain plots (b) and (c), the free parameters A and B are optimised anew for each run-time, using the two cost functions (B.9) and (32). The final excitation E ex,1 (t f , η) of one of the ions is calculated on a grid of run-times and perturbation parameters, assuming the ions to be in their ground state initially. Note that plotting E ex,2 results in the same data, but mirrored along the line η = 0. The physical constraints for all shown schemes are those from table 1. The perturbation is also shown in terms of the homogeneous field γ, which is obtained from η using (30). The insets above (b) and (c) show a cut through the data (indicated by the white dashed lines), corresponding to the unperturbed final excitations E ex,1 (t f , η = 0) together with a fit to that data and an envelope function. The dotted line in the insets corresponds to final excitations of 0.1 quanta and the hatched area marks the range of run-times where the envelope indicates final excitations above this level. modified for the perturbed system. By applying the results of appendix A, the definition of q − changes instead to beq By adding these constraints to the optimisation of the free parameters A and B, similar to the approach presented in section 5, we are able to achieve optimal robustness with respect to a linear potential perturbation. Analogous to the procedure done for figure B1, this is most easily achieved by minimising the partẼ (−) q of the final energyẼ (−) n in addition to minimising E (+) q as before. Equation (B.5) also contains a mode-coupling term H c , which could only be decoupled in a further transformation ifΔ was time-independent [47]. Thus we ignore it and accept that any protocol optimised in this way will display some degree of mode coupling in the presence of such a homogeneous field. We see numerically that this does not produce a great problem for the situations considered here.
To make sure that all the boundary conditions, including the ones onq − and q + are fulfilled, we minimise the parts of the mode energies related to these auxiliary functions by choosing the cost function This ensures that a shortcut is indeed constructed (in the harmonic approximation that yielded (A.15)) and the effects of a perturbative field are minimised. SinceẼ (−) q depends on the strength of the perturbation η, we choose a value of η = 0.015 as in the full-Hamiltonian cost function (32).
The resulting performance is displayed in figure B2(c), alongside the results already presented in figure 4. For each run-time t f , the optimisation method is run to obtain the optimal parameters A and B. Then, using these values, the final excitations E ex,1 (t f , η) are calculated for a range of perturbation strengths η and for the run-times for which the trajectories have been optimised. This yields a grid of excitations depending both on the run-time and the perturbation. Overlaid over the plots are orange contour lines at a level ofn = 0.1, outlining the areas where we consider excitations to become negligible.
A critical time T crit,app can also be determined for these results in the same way as T crit in section 5 by fitting an envelope function to the excitations at η = 0 (cut indicated by a white dashed line). Whereas T crit = 14.2 μs in section 5.1, we now find T crit,app = 27.5 μs, allowing access to a reduced range of run-times. Overall, we see that compared to the results of the full-Hamiltonian cost function, this approach yields a significantly smaller region of robustness.

Appendix C. Scaling of experimental parameters with the maximal quartic confinement
We want to determine the behaviour of the ion distance d(t) and the potential curvatures ω i (t) when scaling the value of the maximal quartic confinement β max . For this we rewrite the distance as where D(t) is the dimensionless distance trajectory that we assume to be the same for all values of β max . This turns out to be true after comparing the results of the numerical optimisations.
In the examples shown in figure 7, D(t) varied from D(0) = 5 to D(t f /2) = {1.0, 1.15, 1.20, 1.25} and back to D(t f ) = 5. We thus conclude that the distance scales as d(t) ∝ (β max ) −1/5 . The quartic potential can be easily rewritten in the same way to where we defined the dimensionless quartic potential term B(t).
Using these results, the scaling of the potential curvatures ω i can be found by rewriting (8)  where we have used (5) for the second equality. The curvatures therefore scale as ω i ∝ (β max ) 3/10 . One can easily see from (6) and (7) that the same is true for Ω − and Ω + . The scaling of the motional exchange time in (33) can now easily be proven by collecting the scaling of the constituting terms. We find that

Appendix D. Dynamical normal modes for ions of unequal mass
We again follow [44] in constructing the dynamical normal modes for ions of unequal mass. The equilibrium positions fulfil