Robust optical clock transitions in trapped ions using dynamical decoupling

Optical clocks based on neutral atoms trapped in optical lattices and single trapped ions have reached estimated systematic uncertainties of a few parts in 10⁻¹⁸ [1–4] or even below [5]. Taking advantage of these record uncertainties for applications ranging from relativistic geodesy [6–9] over fundamental physics [10–12] to frequency metrology [13–17] requires achieving statistical measurement uncertainties of the same level within practical averaging times τ (given in seconds). This has been achieved with single-ensemble optical lattice clocks in self-comparison experiments up to a level of $1.6\times {10}^{-16}/\sqrt{\tau }$ [18] and by implementing an effectively dead-time-free clock consisting of two independent clocks probed in an interleaved fashion [19, 20], reaching a statistical uncertainty in the range of $5\times {10}^{-17}/\sqrt{\tau }$. In contrast to neutral atom lattice clocks, which are typically probed with hundreds to thousands of atoms, single ion clocks are currently limited in their statistical uncertainty by quantum projection noise [21] to levels of a few parts in ${10}^{-15}/\sqrt{\tau }$ [3, 22, 23]. The statistical uncertainty can be improved by probing for longer times, ultimately limited by the excited clock state lifetime or the laser coherence time [24, 25]. Alternatively, the number of probed ions can be increased.

Recently, multi-ion clock schemes have been proposed to address this issue [26–29]. However, several challenges have to be overcome to maintain and transfer the small and very well characterizable systematic shifts achievable with single trapped ions to larger ion crystals. The oscillating rf field in Paul traps results in ac Stark and second order Doppler shifts through micromotion [30–32]. Furthermore, electric field gradients from the trapping fields and the surrounding ions couple to atomic quadrupole moments, resulting in an electric quadrupole shift (QPS). The effects of micromotion can be avoided by trapping strings of ions in a precision-machined linear Paul trap with negligible excess micromotion from trap imperfections [28, 33]. The QPS in such chains can be avoided by choosing an ion species with negligible differential electric quadrupole moment between the clock states, such as ${\mathrm{In}}^{+}$ or ${\mathrm{Al}}^{+}$ , or by employing ring traps in which the QPS is the same for all ions [29]. A high-accuracy multi-ion clock based on ion chains containing on the order of tens of ions in a linear quadrupole trap has been proposed and is expected to achieve trap-induced fractional systematic uncertainties at the 10⁻¹⁹ level [26].

An alternative approach based on large 3d Coulomb crystals of ions has been theoretically investigated [27] for ion species for which the micromotion-induced Doppler shift and the scalar ac Stark shift, both driven by the rf trapping field, can be made to cancel at a 'magic' rf drive frequency [30]. This cancellation has been employed for single ${\mathrm{Sr}}^{+}$ [34], single ${\mathrm{Ca}}^{+}$ [35], and is currently being investigated for ${\mathrm{Lu}}^{+}$ [27, 36, 37]. However, electronic states with J > 1/2 are subject to rank 2 tensor shifts, such as rf electric field-induced tensor ac Stark shift (TASS) and QPS [30, 32]. For 3d Coulomb crystals with tens up to several hundreds of ions, this results in position-dependent shifts, since the electric field environment of the ions differ. It has been proposed to reduce this inhomogeneous broadening across the ion crystal by employing a spherical ion crystal to minimize the QPS and by adding a compensating laser field for the TASS [27], or by operating at a judiciously chosen magnetic-field-insensitive point [37] for ions with hyperfine structure.

Achieving insensitivity of atomic energy levels to external field fluctuations has been theoretically investigated using pulsed or continuous dynamical decoupling (CDD) [38–48]. Such schemes have been experimentally implemented in various systems ranging from NV centers and solid state spin systems [49–57] to neutral atoms [58] and trapped ions [59–64].

CDD or dressed-state engineering in the context of clocks has been proposed for linear Zeeman shift cancellation in neutral atom clocks through rf dressing of Zeeman substates in a regime where the drive Rabi frequency is larger than the drive frequency [45]. However, this scheme is unable to cancel tensorial shifts, such as the QPS or TASS. In a similar approach, magnetic field noise suppression up to second order in radio-frequency clocks via weak rf dressing has been proposed [46] and experimentally implemented to engineer synthetic clock states for rf spectroscopy in ultra-cold rubidium [65, 66].

Here, we show through numerical simulations that CDD using four rf frequencies significantly suppresses all relevant homogeneous (linear Zeeman) and inhomogeneous (micromotion-induced second order Doppler, QPS, scalar and tensor ac Stark) frequency shifts on an optical clock transition for an ion crystal containing 400 ions. We note that the scheme in principle works for any ion number and crystal size. We demonstrate the basic principle, which is illustrated in figure 1, on the ${}^{2}{{\rm{S}}}_{1/2}$ $\leftrightarrow$ ${}^{2}{{\rm{D}}}_{5/2}$ clock transition in ${}^{40}{\mathrm{Ca}}^{+}$ , but the scheme is directly applicable to other systems as well. It therefore allows the operation of a multi-ion clock [26] using ion species whose clock transitions have a non-vanishing differential electric quadrupole moment. One of the many possible applications of such a multi-ion frequency reference is the phase stabilization of a probe laser for a single ion clock to allow near-lifetime-limited probe times and correspondingly reduced statistical uncertainties [24, 25].

Zoom In Zoom Out Reset image size

The paper is organized as follows. In section 2 we first show how robustness to QPS and TASS can be achieved by the application of a single continuous detuned driving field. The discussion of only one driving field gives a simple and intuitive understanding of our approach. In section 3 we extend the discussion to the full CDD scheme employing four driving fields for the construction of a robust optical clock transition. For simplicity and clarity in the discussion we ignore effects of the counter-rotating terms of the driving fields by assuming the rotating-wave-approximation (RWA), as well as the cross-driving effects (for example, the effect of the ${}^{2}{{\rm{S}}}_{1/2}$ driving field on the ${}^{2}{{\rm{D}}}_{5/2}$ states). These two effects are, however, taken into account when numerically optimizing the driving parameters as well as in the simulations presented. In appendices C–E, we show in detail how this is achieved by incorporating these effects in the optimization. Then, in section 4 we consider the implementation of the scheme in the case of a multi-ion crystal clock and analyze the performance of the scheme in terms of the expected statistical uncertainties of the robust optical clock transition. After discussing possible applications in section 5, we end with the conclusions in section 6.

In the absence of hyperfine structure, tensor shifts including the QPS are proportional to ${Q}_{J,{m}_{j}}=J(J+1)-3{m}_{j}^{2}$ [32], where J is the total angular momentum and m_J the magnetic quantum number. These shifts reduce the precision of atomic clocks when not suppressed by suitable averaging schemes. Previously employed schemes to suppress the quadrupole shift include averaging the transition frequency over all m_J states, since ${\sum }_{{m}_{J}=-J}^{J}{Q}_{J,{m}_{j}}=0$ [67]. Such an average will also eliminate the linear Zeeman shift, assuming the field does not change between the frequency measurements contributing to the average. We propose a novel dynamical decoupling scheme in which robustness to this type of shifts is achieved by the application of a detuned driving field, mixing all m_J states to form dressed states with effective ${Q}_{J,{m}_{j}}=0$. While the cancellation scheme is general and applies to tensor shifts of arbitrary electronic states, we consider in the following the Hamiltonian of the ${}^{2}{{\rm{D}}}_{5/2}$ states of e.g. ${\mathrm{Ca}}^{+}$ ions

$$\begin{eqnarray}&&H={g}_{d}{\mu }_{B}{{BS}}_{z}+{g}_{d}{{\rm{\Omega }}}_{1}\cos \left[({g}_{d}{\mu }_{B}B-{\delta }_{1})t\right]{S}_{x},\end{eqnarray} \tag{ 1 }$$

where g_d μ_BB is the Zeeman splitting due to the static magnetic field B, g_d is the gyromagnetic ratio of the ${}^{2}{{\rm{D}}}_{5/2}$ states, S_z and S_x are the z and x spin-5/2 matrices, ${{\rm{\Omega }}}_{1}$ is the Rabi frequency of the driving field, and δ₁ is the detuning. Moving to the interaction picture (IP) with respect to H₀ = (g_dμ_BB − δ₁)S_z and taking the RWA ($({g}_{d}{\mu }_{B}B-{\delta }_{1})\gg {{\rm{\Omega }}}_{1}$), results in

$$\begin{eqnarray}&&{H}_{I}={\delta }_{1}{S}_{z}+\displaystyle \frac{{g}_{d}{{\rm{\Omega }}}_{1}}{2}{S}_{x},\end{eqnarray} \tag{ 2 }$$

where H_I is the Hamiltonian in the IP. Since the bare ${}^{2}{{\rm{D}}}_{5/2}$ states have tensor shifts which are proportional to ${Q}_{5/2,\pm 5/2}=-10$, ${Q}_{5/2,\pm 3/2}=+2$, and ${Q}_{5/2,\pm 1/2}=+8$, the tensor shifts of the dressed states (the eigenstates of H_I) are proportional to

$$\begin{eqnarray}&&Q=-5q({\delta }_{1},{{\rm{\Omega }}}_{1}),\quad +q({\delta }_{1},{{\rm{\Omega }}}_{1}),\quad +4q({\delta }_{1},{{\rm{\Omega }}}_{1}),\end{eqnarray} \tag{ 3 }$$

with

$$\begin{eqnarray}&&q({\delta }_{1},{{\rm{\Omega }}}_{1})=\displaystyle \frac{-8{\delta }_{1}^{2}+{\left({g}_{d}{{\rm{\Omega }}}_{1}\right)}^{2}}{4{\delta }_{1}^{2}+{\left({g}_{d}{{\rm{\Omega }}}_{1}\right)}^{2}}.\end{eqnarray} \tag{ 4 }$$

Hence, by choosing ${\delta }_{1}=\pm \sqrt{\tfrac{1}{8}}\left({g}_{d}{{\rm{\Omega }}}_{1}\right)$, all of the dressed states have a zero (first order) tensor shift, Q = 0 (see figure 2). Similarly, for a fixed detuning δ₁ choosing ${{\rm{\Omega }}}_{1}=\sqrt{8}{\delta }_{1}/{g}_{d}$ results in a zero tensor shift. This can also be understood in the lab frame, in which the coupling Ω₁ drives a rotation among the bare states, averaging their shifts in time just as in the m_J averaging scheme mentioned above. However, here the averaging takes place at a rate corresponding to the Rabi frequency rather than the experiment repetition rate, allowing it to suppress much faster field fluctuations.

Zoom In Zoom Out Reset image size

The cancellation of the tensor shift to first order can also be understood as follows. The tensor and quadrupole shift operator $\hat{Q}={S}^{2}-3{S}_{z}^{2}$ can also be written as $\hat{Q}={S}_{x}^{2}+{S}_{y}^{2}-2{S}_{z}^{2}$. Moving to a rotated basis (dressed states basis) defined by $z\to \cos (\theta )z+\sin (\theta )x$, $x\to \cos (\theta )x-\sin (\theta )z$, and $y\to y$, and neglecting purely off-diagonal terms we obtain that $\hat{Q}\approx {S}_{z}^{2}\left(1-3{\cos }^{2}\left(\theta \right)\right)+{S}_{x}^{2}\left(1-3{\sin }^{2}\left(\theta \right)\right)+{S}_{y}^{2}$. To first order (diagonal terms) we have that $\hat{Q}\approx \left({S}_{+}{S}_{-}+{S}_{-}{S}_{+}\right)\left(2-3{\sin }^{2}\left(\theta \right)\right)+{S}_{z}^{2}\left(1-3{\cos }^{2}\left(\theta \right)\right),$ where S_± = S_x ± iS_y are the spin ladder operators. Hence to first order, the tensor shift vanishes for $\cos \left(\theta \right)=\tfrac{1}{\sqrt{3}}$, which is analogous to magic-angle spinning in solid-state NMR spectroscopy [68]. The detuned driving field results in dressed states, which are the eigenstates of equation (2). The rotation to the dressed states basis is given by $U={{\rm{e}}}^{{\rm{i}}\theta {S}_{y}}$, where $\cos \left(\theta \right)=\tfrac{{\delta }_{1}^{2}}{\sqrt{{\delta }_{1}^{2}+\tfrac{{\left({g}_{d}{{\rm{\Omega }}}_{1}\right)}^{2}}{4}}}$. For ${\delta }_{1}=\sqrt{\tfrac{1}{8}}\left({g}_{d}{{\rm{\Omega }}}_{1}\right)$ we obtain $\cos \left(\theta \right)=\tfrac{1}{\sqrt{3}}$.

Our scheme is based on the application of continuous driving fields for the construction of a robust optical clock transition. We consider the sub-levels of the ${}^{2}{{\rm{S}}}_{1/2}$ and ${}^{2}{{\rm{D}}}_{5/2}$ states, the usual clock states in ${\mathrm{Ca}}^{+}$ optical clocks. Four driving fields are employed, where two driving fields operate on the ${}^{2}{{\rm{S}}}_{1/2}$ states, and two driving fields operate on the ${}^{2}{{\rm{D}}}_{5/2}$ states. The rf driving fields continuously drive the bare ${}^{2}{{\rm{S}}}_{1/2}$ and ${}^{2}{{\rm{D}}}_{5/2}$ states, that is, continuously couple between the Zeeman sub-levels, which results in the desired robust optical transition of the dressed states. The first driving field of the ${}^{2}{{\rm{D}}}_{5/2}$ states mitigates tensor shifts, including the QPS, as shown in the previous section. However, the ${}^{2}{{\rm{D}}}_{5/2}$ states, as well as the ${}^{2}{{\rm{S}}}_{1/2}$ states, are still sensitive to the linear Zeeman shift. Moreover, the ${}^{2}{{\rm{D}}}_{5/2}$ states are also sensitive to amplitude fluctuations of the driving field. The purpose of adding three more driving fields is to have enough control degrees of freedom that can be tuned such that the suppression of both linear Zeeman shift and amplitude fluctuations of a single (dressed) ${}^{2}{{\rm{S}}}_{1/2}$ $\leftrightarrow$ ${}^{2}{{\rm{D}}}_{5/2}$ transition is achieved. Hence, the driving scheme results in doubly-dressed ${}^{2}{{\rm{S}}}_{1/2}$ and ${}^{2}{{\rm{D}}}_{5/2}$ states where one ${}^{2}{{\rm{S}}}_{1/2}\leftrightarrow {}^{2}{{\rm{D}}}_{5/2}$ transition (between the doubly-dressed states) is robust to magnetic shifts, quadrupole shifts, tensor shifts, and driving amplitude shifts caused by amplitude fluctuations of the driving fields (see figure 3). For the ${}^{2}{{\rm{D}}}_{5/2}$ (${}^{2}{{\rm{S}}}_{1/2}$) states, the Rabi frequencies, ${{\rm{\Omega }}}_{k}$, and the detunings, δ_k, of the first and second driving fields are denoted by {Ω₁, δ₁} and {Ω₂, δ₂} ({Ω₃, δ₃} and {Ω₄, δ₄}), respectively. We assume that a static magnetic field B is applied, which results in the Zeeman splitting of the ${}^{2}{{\rm{S}}}_{1/2}$ and ${}^{2}{{\rm{D}}}_{5/2}$ states. The rf dressing is based on magnetic dipole coupling between the Zeeman sub-levels, where the Rabi frequency of a driving field is related to the magnetic dipole moment by ${\hslash }{\rm{\Omega }}=\vec{\mu }\cdot {\vec{B}}_{{\rm{rf}}}$, where $\vec{\mu }$ is the magnetic dipole moment and ${\vec{B}}_{{\rm{rf}}}$ is the amplitude vector of the oscillating rf field.

Zoom In Zoom Out Reset image size

For simplicity and clarity, in the following analysis we neglect the effect of the counter-rotating terms of the driving fields by assuming the RWA, as well as the cross-driving effects. Both the counter-rotating terms and the cross-driving fields result in small energy shifts, which are taken into account when optimizing the driving parameters. This, however, does not change the overall physical picture, and results only in small variations of the optimal driving parameters. In appendices C–E we show how this can be done by incorporating these effects in the following optimization.

Similar to the derivation in section 2, by moving to the first IP with respect to the frequencies of the first driving fields we obtain the dressed states, which are the eigenstates of

$$\begin{eqnarray}&&{H}_{{\rm{I}}}={\delta }_{1}{S}_{z}+\displaystyle \frac{{g}_{d}{{\rm{\Omega }}}_{1}}{2}{S}_{x}+{\delta }_{3}{s}_{z}+\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}{s}_{x},\end{eqnarray} \tag{ 5 }$$

where S_i and s_i are the spin matrices of the ${}^{2}{{\rm{D}}}_{5/2}$ and ${}^{2}{{\rm{S}}}_{1/2}$ states and g_d = 6/5 and g_s = 2 are the gyromagnetic ratios of the ${}^{2}{{\rm{D}}}_{5/2}$ and ${}^{2}{{\rm{S}}}_{1/2}$ states, respectively. We proceed by moving to the basis of the dressed states and then to the IP with respect to the frequencies of the second driving fields and obtain the doubly-dressed states, which are the eigenstates of

$$\begin{eqnarray}&&{H}_{{\rm{II}}}={\delta }_{2}{S}_{\tilde{z}}+\displaystyle \frac{{g}_{d}{{\rm{\Omega }}}_{2}}{4}{S}_{y}+{\delta }_{4}{s}_{\tilde{\tilde{z}}}+\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{4}}{4}{s}_{y},\end{eqnarray} \tag{ 6 }$$

where we denote by $\tilde{z}$ and $\tilde{\tilde{z}}$ the diagonalized basis of the dressed ${}^{2}{{\rm{D}}}_{5/2}$ and ${}^{2}{{\rm{S}}}_{1/2}$ states, respectively. A detailed derivation is given in appendix A.

We consider the doubly-dressed ${}^{2}{{\rm{S}}}_{1/2}$ and ${}^{2}{{\rm{D}}}_{5/2}$ states with the smallest positive eigenvalue as the robust optical clock states. By setting ${\delta }_{1}=\sqrt{\tfrac{1}{8}}\left({g}_{d}{{\rm{\Omega }}}_{1}\right)$ robustness to quadrupole and tensor shifts is attained. We denote by δ_B = μ_Bδ B and δ_Ω a magnetic field shift and the relative driving amplitude shift respectively (please recall that the detunings of the driving fields are denoted by δ_k, see above). In order to achieve robustness to shifts in the magnetic and driving fields, we first add to the Hamiltonian the magnetic shift terms g_dδ_BS_z + g_sδ_Bs_z, and the driving amplitude shift terms δ_ΩΩ_k (so the driving amplitude including the shift is given by (1 + δ_Ω)Ω_k). We use a single parameter δ_Ω to describe the field amplitude fluctuations, which we assume to be correlated across all four dressing fields. This describes the experimental situation where the dominant variations in dressing-field amplitude are due to spatial inhomogeneities or to a common rf amplifier through which all four signals pass. The Hamiltonian of the double-dressed states H_II (equation (6)) now includes both δ_B and δ_Ω, which modify the eigenvalues (the energies) of the double-dressed states, and hence the optical transition frequency. Our aim is to mitigate the sensitivity of the optical transition frequency to δ_B and δ_Ω. We achieve this by tuning the driving parameters such that the leading order contributions of δ_B and δ_Ω to the eigenvalues of the double-dressed ${}^{2}{{\rm{S}}}_{1/2}$ and ${}^{2}{{\rm{D}}}_{5/2}$ states will be as close as possible (and ideally identical). We therefore calculate the power series expansion of the eigenvalues to orders of ${\delta }_{{\rm{B}}}^{i}$ and ${\delta }_{{\rm{\Omega }}}^{i}$ (i = 1, 2, ...). We denote the series expansion terms of the magnetic shift of the ${}^{2}{{\rm{S}}}_{1/2}$ and ${}^{2}{{\rm{D}}}_{5/2}$ states by ${Z}_{{S}_{i}}{\delta }_{{\rm{B}}}^{i}$ and ${Z}_{{D}_{i}}{\delta }_{{\rm{B}}}^{i}$ respectively. Similarly, we denote the series expansion terms of the amplitude driving shift of the ${}^{2}{{\rm{S}}}_{1/2}$ and ${}^{2}{{\rm{D}}}_{5/2}$ states by ${O}_{{S}_{i}}{\delta }_{{\rm{\Omega }}}^{i}$ and ${O}_{{D}_{i}}{\delta }_{{\rm{\Omega }}}^{i}$ respectively. The expansion terms of the correlated shifts are denoted by ${{ZO}}_{{S}_{i}}{\delta }_{{\rm{B}}}^{i}{\delta }_{{\rm{\Omega }}}$ and ${{ZO}}_{{D}_{i}}{\delta }_{{\rm{B}}}^{i}{\delta }_{{\rm{\Omega }}}$ for the ${}^{2}{{\rm{S}}}_{1/2}$ and ${}^{2}{{\rm{D}}}_{5/2}$ states respectively. We calculate the magnetic energy shifts, ${Z}_{{S}_{i}}{\delta }_{{\rm{B}}}^{i}$ and ${Z}_{{D}_{i}}{\delta }_{{\rm{B}}}^{i}$, up to fourth order (i = 1, ..., 4), driving energy shifts, ${O}_{{S}_{i}}{\delta }_{{\rm{\Omega }}}^{i}$ and ${O}_{{D}_{i}}{\delta }_{{\rm{\Omega }}}^{i}$, up to second order (i = 1, 2), and the correlated shifts, ${{ZO}}_{{S}_{i}}{\delta }_{{\rm{B}}}^{i}{\delta }_{{\rm{\Omega }}}$ and ${{ZO}}_{{D}_{i}}{\delta }_{{\rm{B}}}^{i}{\delta }_{{\rm{\Omega }}}$ (i = 1, 2) as function of the driving parameters, Ω_k and δ_k. A detailed derivation is given in appendix B. We continue by defining a goal function

$$\begin{eqnarray}&&G=\displaystyle \sum _{i=1,j=1}^{i=4,j=2}| {Z}_{{S}_{i}}-{Z}_{{D}_{i}}| {\delta }_{{\rm{B}}}^{i}+| {O}_{{S}_{j}}-{O}_{{D}_{j}}| {\delta }_{{\rm{\Omega }}}^{j}+| {{ZO}}_{{S}_{j}}-{{ZO}}_{{D}_{j}}| {\delta }_{{\rm{B}}}^{j}{\delta }_{{\rm{\Omega }}}.\end{eqnarray} \tag{ 7 }$$

Given distributions of the shifts, δ_B and δ_Ω, of a specific experimental set-up, we then define an averaged goal function

$$\begin{eqnarray}&&{G}_{A}=\langle G\left({\delta }_{{\rm{B}}}^{m},{\delta }_{{\rm{\Omega }}}^{n}\right){\rangle }_{m,n},\end{eqnarray} \tag{ 8 }$$

where ${\delta }_{{\rm{B}}}^{m}$ and ${\delta }_{{\rm{\Omega }}}^{n}$ are chosen randomly according to the given distributions, and numerically minimize G_A over the driving parameters, Ω_k and δ_k. The numerical minimization results in optimal sets of values of the driving fields, for which robustness to shifts in the magnetic and driving fields is obtained. We note that because it is experimentally possible to generate very precise and stable frequencies we do not consider uncertainties in the frequencies of the driving fields.

As already mentioned, in this derivation of the optimal driving parameters we assumed the RWA and neglected cross-driving effects. However, these effects, which result in small energy shifts, can be integrated in the optimization as we show in appendices C–E. Indeed, in the numerical analysis of the resulting shift distribution in section 4 we took the counter-rotating terms of the driving fields and the cross-driving effect into account and integrated them in the optimization. In addition, the simulations where performed using the full driving Hamiltonian without making any approximations.

We now consider the implementation of a robust multi-ion clock employing the proposed scheme and realizable with current ion trap technology. We estimate the dominating field inhomogeneities (QPS and TASS) for ${}^{40}{\mathrm{Ca}}^{+}$ as a widely used clock species [69–71] with convenient properties, and discuss effects of micromotion.

As a trap platform, we consider a linear Paul trap, in which axial micromotion can be made sufficiently small [28, 33]. In large three-dimensional ion crystals, each ion will experience micromotion driven by the rf-field of the trap with an amplitude proportional to the distance from the trap's symmetry axis. Due to the construction of the trap, this micromotion is oriented along the radial axis (perpendicular to the trap axis). When probing along the radial degrees of freedom, the ions' coupling to the probe laser is diminished by the strong Doppler modulation of the laser in the ions' frame. We therefore choose to probe along the trap symmetry axis z with a magnetic field oriented along the same direction.

The QPS of an ion in a multi-ion crystal is caused by electrical field gradients originating from the space charge of all other ions, and the electrical field gradients of the trap itself. The overall scale of the QPS is determined by the quadrupole moment expressed as a reduced matrix element, which for ${}^{40}{\mathrm{Ca}}^{+}$ was measured to be ${\rm{\Theta }}(3d,5/2)=1.83(1){{ea}}_{0}^{2}$ [72]. The QPS depends on the angle between the quantization axis and the electric field gradient as well as the state of the ion. The angle dependence is given by [27, 32, 72]

$$\begin{eqnarray}\begin{array}{rcl}f\left(\alpha ,\beta ,\gamma \right) & = & \displaystyle \frac{\partial {E}_{z}}{\partial z}\displaystyle \frac{1}{4}\left(3\cos {\left(\beta \right)}^{2}-1\right)+\displaystyle \frac{1}{2}\sin (2\beta )\left(\displaystyle \frac{\partial {E}_{x}}{\partial z}\cos \alpha +\displaystyle \frac{\partial {E}_{y}}{\partial z}\sin \alpha \right)\\ & & +\displaystyle \frac{1}{4}{\sin }^{2}\beta \left[\left(\displaystyle \frac{\partial {E}_{x}}{\partial x}-\displaystyle \frac{\partial {E}_{y}}{\partial y}\right)\cos (2\alpha )+2\displaystyle \frac{\partial {E}_{x}}{\partial y}\sin (2\alpha )\right]\end{array}\end{eqnarray} \tag{ 9 }$$

with the Euler angles {α, β, γ} defined as in [32]. The state dependence can be described by

$$\begin{eqnarray}&&g(J,{m}_{J})=\displaystyle \frac{J(J+1)-3{m}_{J}^{2}}{J(2J-1)},\end{eqnarray} \tag{ 10 }$$

with a total QPS of ${\rm{\Delta }}{f}_{\mathrm{QPS}}=\tfrac{{\rm{\Theta }}}{h}\times f(\alpha ,\beta ,\gamma )\times g(J,{m}_{J})$.

The Paul trap features two types of electrical field gradients: the rf-field employed for radial confinement and the static field gradient for axial confinement. The rf-field is averaged out and therefore does not lead to a quadrupole shift of the transition. The static gradient for axial trapping is constant and the same for all ions, causing no line broadening but a constant shift for each individual ion of the crystal. Provided the axial trap voltages are well controlled, this constant shift does not pose a limit on clock stability. Furthermore, this constant shift is canceled by the dynamical decoupling scheme, see below.

The shift originating from the space charge will in general lead to inhomogeneous line broadening. The space charge-induced quadrupole shift falls off cubicly with distance and is therefore dominated by the ion's local environment. In a linear chain of ions, for instance, the quadrupole shift caused by the space charge results in a significant shift depending on the position of the ion, since contributions from the individual ions add up. For example, a chain of 30 ${}^{40}{\mathrm{Ca}}^{+}$ ions exhibits an inhomogeneous shift of about 80 Hz across the chain (radial trap frequency ω_r = 2π × 1 MHz and axial trap frequency ω_a = ω_r/12). In a spherically symmetric crystal configuration, in contrast, the symmetry suppresses the quadrupole shift due to space charge to a large extent. In the limit of a crystal of infinite size a body-centered-cubic lattice forms and the quadrupole shift vanishes [27].

For this reason, we consider a spherical trap configuration in which it is straightforward to obtain near isotropic trap frequencies ${\omega }_{r,x}\approx {\omega }_{r,y}\approx {\omega }_{a}\approx 2\pi \times 1$ MHz. A slight deviation from complete symmetry is desirable to pin the orientation of the crystal and allow efficient Doppler cooling. For our scheme, we assume ions cooled to the ground state, as for instance shown experimentally for large crystals in [73]. It is advantageous to work with as many ions as possible in order to improve the signal-to-noise ratio and hence the stability of the clock. On the other hand, cooling and trapping large crystals can be challenging and the size of the crystal needs to be reasonable so that the probe, cooling, and control lasers can address all the ions. Also, the probability of a background gas collision during clock interrogation grows with the number of ions. For our purpose, we consider a realistic implementation employing 400 ${}^{40}{\mathrm{Ca}}^{+}$ ions, resulting in an isotropic crystal with approximately 30 μm radius for our proposed trap parameters.

For this configuration, we estimate the QPS by first finding the equilibrium positions of the ions in the crystal by minimizing the pseudo-potential energy in the linear Paul trap. We next calculate the effective electric field tensor $\tfrac{\partial {E}_{i}}{\partial {x}_{j}}$ and the resulting quadrupole shift. In figure 4 we show the resulting distribution of the QPS. The distribution shows a mean shift of about 10 Hz due to the static axial field gradient of the trap and a standard deviation of 0.74 Hz.

Zoom In Zoom Out Reset image size

In a spherical ion crystal, most of the ions are trapped away from the symmetry axis and the resulting rf-field at the individual ions' position causes a scalar and TASS. While for ${}^{40}{\mathrm{Ca}}^{+}$ the scalar shift can be canceled by operating the trap at the 'magic' drive frequency [35], the TASS will result in shifted transition frequencies for each individual ion depending on its position. We estimate the TASS in our system by first calculating the amplitude of the rf-field due to the trap drive at each ion's position, employing a numerical simulation. The resulting TASS is then estimated by evaluating the time average $\langle ...\rangle$ over one cycle of the squared electrical field E [27]:

$$\begin{eqnarray}&&\displaystyle \frac{\delta f}{f}=\displaystyle \frac{g(J,{m}_{J})}{4}\displaystyle \frac{{\alpha }_{{\rm{dc}}}}{h\nu }\langle {E}^{2}\rangle ,\end{eqnarray} \tag{ 11 }$$

with the dc tensor polarizability α_dc from [74]. Figure 5 shows the result for our configuration, exhibiting an approximately uniform distribution with about 25 Hz width (see appendix F).

Zoom In Zoom Out Reset image size

With the distributions as described above, we now turn to the implementation and optimization of the CDD sequence.

For the numerical optimization of the goal function equation (7), we consider a static magnetic field such that the Zeeman splitting of the ${}^{2}{{\rm{S}}}_{1/2}$ states is equal to 10 MHz, where the uncertainty of the Zeeman splitting g_sδ_B is normally distributed with a zero mean and a width of 1 kHz. There is also a small contribution due to the second order Zeeman effect. For the assumed magnetic field noise, this shifts amounts to an additional broadening of about 0.5 mHz (averaged over all relevant levels) and can therefore be neglected. In a real experiment, the rf driving fields used for dressing the states will be imperfect and show some fluctuations. We estimate the corresponding relative amplitude shifts of the driving fields δ_Ω to be normally distributed with a zero mean and a width of 4 × 10⁻⁴.

In order to test the performance of the scheme in the multi-ion spherical crystal configuration, we have realized the averaged goal function G_A, which is given by equation (8), with an average over 100 pairs of {δ_B, δ_Ω} chosen randomly according to the above distributions, and then numerically minimized G_A to obtain optimal driving parameters (all in units of kHz): Ω₁ = 2π × 225.3, δ₁ = 2π × 95.6, Ω₂ = 2π × 13.6, δ₂ = 2π × 5, Ω₃ = 2π × 93.6, δ₃ = 2π × 27.2, and Ω₄ = 2π × 14.8, δ₄ = 2π × 25.6. We then integrated the effect of the counter-rotating terms of the driving fields and the cross driving fields (the effect of the ${}^{2}{{\rm{S}}}_{1/2}$ $\left({}^{2}{{\rm{D}}}_{5/2}\right)$ drive on the ${}^{2}{{\rm{D}}}_{5/2}$ $\left({}^{2}{{\rm{S}}}_{1/2}\right)$ states) and adjusted the optimal driving parameters accordingly. We simulated 4927 trials of the robust clock-transition with these optimal driving parameters and the above distributions of the Zeeman, quadrupole, tensor and driving amplitude shifts. In the simulations we used the full driving Hamiltonian without making any approximations. For more details see appendices C–E. The simulation results, shown in figure 6, indicate that the shift distribution of the robust transition has a narrow width of ∼1 Hz.

Zoom In Zoom Out Reset image size

The distribution shown in the figure is slightly asymmetric and shifted away from zero. The asymmetry is not a problem when using Ramsey interrogation with short broadband π/2 pulses which address all the ions at once, yielding a symmetric (cosine) central resonance. However, the overall shift and the asymmetry lead to a probe-time dependent fractional bias in the clock frequency of around 0.12 Hz for a probe time of 150 ms, which would have to be evaluated and corrected in absolute frequency measurements.

This shift from the non-perturbed virtual (dressed) ${}^{2}{{\rm{S}}}_{1/2},{m}_{J}=0\leftrightarrow {}^{2}{{\rm{D}}}_{5/2},{m}_{J}={Q}_{J,{m}_{J}}=0$ transition depends on the initial distribution of inhomogeneously shifted line centers and the shift suppression achievable with the selected dynamical decoupling parameters. A qualitative understanding of the influence of different parameters on this shift can be developed from simulation results for the scaling of the shift of the robust transition as a function of the magnitude of the individual shifts, namely, the Zeeman shift uncertainty, quadrupole and tensor shifts, and driving amplitude shift. The simulations shown in figure 7 were performed with the driving parameters used in the simulations of figure 6. For each individual shift, several simulations were conducted with different values of that individual shift while all other shift contributions were set to zero. When varying the magnetic field uncertainty, we see some higher-order (quadratic and cubic) contributions but no linear dependence, showing that the decoupling fields have greatly suppressed the first-order Zeeman shift. This can be understood as an avoided crossing of the dressed levels with respect to the magnetic shift uncertainty δ_B. The shift of the robust transition as function of the quadrupole and tensor shifts scales linearly. For a perfect drive of the ${}^{2}{{\rm{D}}}_{5/2}$ states (assuming the RWA and no effect of the ${}^{2}{{\rm{S}}}_{1/2}$ drive) we expect a quadratic scaling. However, the ${}^{2}{{\rm{S}}}_{1/2}$ drive and the counter-rotating terms of the ${}^{2}{{\rm{D}}}_{5/2}$ drive result in an amplitude mixing between the ideal dressed ${}^{2}{{\rm{D}}}_{5/2}$ states, which reintroduces a first order contribution to the quadrupole and tensor shifts. The shift of the robust transition as a function of the drive amplitude error scales linearly as the robust transition frequency depends linearly on the drive amplitude. In an experiment the described scaling behavior could be exploited to verify the correct parameter set for the dynamical decoupling scheme by measuring the expected scaling and identifying the correct operating point. Note that the simulations were performed with fixed driving parameters, optimized as described above for a particular distribution of field shifts and drive amplitudes. For different assumptions on the distribution of field and amplitude fluctuations, different driving parameters could be derived that improve the scheme's performance case-by-case. Moreover, as in other CDD schemes, the performance of the scheme could, in principle, be improved by adding more driving fields.

Zoom In Zoom Out Reset image size

In our analysis we assumed that all driving fields suffer from driving fluctuations. Generating the second driving fields (Ω₂ and Ω₄) by a phase modulation, as proposed in [41], should result in stable driving fields with negligible amplitude fluctuations. In this case we expect a further improvement in the performance of the scheme.

The coupling strength between the dressed ${}^{2}{{\rm{S}}}_{1/2}$ states and the dressed ${}^{2}{{\rm{D}}}_{5/2}$ states is achieved via the laser coupling between the bare states that have a non-vanishing amplitude in the desired dressed states. Hence, the effective laser coupling strength is modified by the overlap between the bare states and the single or double dressed states. In our case the laser coupling of the bare ${}^{2}{{\rm{S}}}_{1/2},{m}_{J}=+\tfrac{1}{2}\leftrightarrow {}^{2}{{\rm{D}}}_{5/2},{m}_{J}=+\tfrac{5}{2}$ transition is reduced by a factor of 0.51 (0.3) for the transition between a single (double) dressed ${}^{2}{{\rm{S}}}_{1/2}$ state and a single (double) dressed ${}^{2}{{\rm{D}}}_{5/2}$ state. Note that the effective laser coupling strength should in any case be smaller than the energy gap of the double dressed states, which in our case is ∼6 kHz.

Using the results of [25], we estimate the achievable statistical uncertainty of a 400 ion ${\mathrm{Ca}}^{+}$ clock when probed with a flicker floor-limited probe laser at 10⁻¹⁴ (10⁻¹⁵) to be $5.8\times {10}^{-16}/\sqrt{\tau /{\rm{s}}}$ ($1.8\times {10}^{-16}/\sqrt{\tau /{\rm{s}}}$). This represents an order of magnitude improvement in instability over current single ion clocks [3, 22, 23], corresponding to a reduction in averaging time by a factor of 100.

One of the many potential applications of the proposed scheme is a cascaded clock [75–77] in which a clock laser is first stabilized to an ensemble of ${\mathrm{Ca}}^{+}$ ions to improve its phase coherence time and thus allow extended probing times [24, 25] in a high-accuracy single-ion clock (e.g. ${\mathrm{Al}}^{+}$ or ${\mathrm{Yb}}^{+}$ ). To bridge the difference in clock transition frequencies, a transfer scheme using a frequency comb would be employed [77–79]. Table 1 shows the achievable instabilities of an ${\mathrm{Al}}^{+}$ clock for different initial clock laser instabilities. We have used the results from [25] to determine the optimum probe times assuming flicker-floor limited laser instability and neglecting spontaneous emission from the excited clock state. Note that after stabilization to the ${\mathrm{Ca}}^{+}$ ion crystal the laser exhibits a white frequency noise spectrum. As expected, the reduction in required averaging time is largest when the initial laser instability is large. As the laser improves, it approaches the quantum projection noise limit of the 400 ${\mathrm{Ca}}^{+}$ ions and the gain is reduced. Even higher gain can be obtained from larger crystals or reference atoms with narrower linewidth than ${\mathrm{Ca}}^{+}$ , such as ${}^{176}{\mathrm{Lu}}^{+}$ [27, 36], ${}^{88}{\mathrm{Sr}}^{+}$ [80, 81], or ${}^{138}{\mathrm{Ba}}^{+}$ [82, 83]. Our scheme works in the same way for other atomic species, however, the contribution of the different broadening mechanisms depend on the employed species and specific properties of the experimental setup. We note that in our case a narrower atomic linewidth (as for example provided by a different atomic species such as ${}^{88}{\mathrm{Sr}}^{+}$ ) would not necessarily lead to a narrower linewidth of the robust transition, due to the remnant effects of the magnetic noise, the amplitude drive noise and the QPS and TASS. Another reason for using ${\mathrm{Ca}}^{+}$ is the advantageous mass ratio relative to ${\mathrm{Al}}^{+}$ , facilitating efficient sympathetic cooling in an ${\mathrm{Al}}^{+}$ clock [84].

Table 1.  Estimated statistical uncertainties. The flicker-floor limited performance of the clock laser is denoted by σ_l, which is assumed to be independent of the averaging time. The Allan deviations and optimum probe times for different configurations k are denoted by σ_k and T_k, respectively. The investigated systems are 400 ${\mathrm{Ca}}^{+}$ ions using the described CDD scheme (k = Ca), a single ${\mathrm{Al}}^{+}$ ion (k = Al), and a cascaded scheme in which a single ${\mathrm{Al}}^{+}$ ion is probed by a laser pre-stabilized through a cloud of 400 ${\mathrm{Ca}}^{+}$ ions (k = CaAl). The reduction in averaging time to achieve a certain statistical measurement uncertainty is given by G.

σl	TCa (s)	$\sigma {\left(\tau \right)}_{\mathrm{Ca}}\sqrt{\tau /{\rm{s}}}$	TAl (s)	$\sigma {\left(\tau \right)}_{\mathrm{Al}}\sqrt{\tau /{\rm{s}}}$	TCaAl (s)	$\sigma {\left(\tau \right)}_{\mathrm{CaAl}}\sqrt{\tau /{\rm{s}}}$	G
10−14	0.015	5.8 × 10−16	0.003 5	3.6 × 10−15	0.25	5.1 × 10−16	50
10−15	0.15	1.8 × 10−16	0.035	1.1 × 10−15	0.78	2.9 × 10−16	14

The cascaded clock scheme enables short averaging times with lasers that are commercially available. For state-of-the-art multi-segment ion traps, it is conceivable to trap a large ${\mathrm{Ca}}^{+}$ crystal in one segment and a clock ion in a different segment, strongly reducing experimental overhead. Furthermore, for the case of additional experimental constraints such as a transportable setup or the lack of a cryogenic system, for instance, our scheme could offer an advantage over state-of-the art optical resonators.

Beyond the cascaded clock scheme, many applications in fundamental physics, navigation and industry do not require ultimate accuracy [1, 4], but rather high stability as provided by a dynamically-decoupled Coulomb crystal clock. Through appropriate characterization of the residual line center shift away from an effective ${}^{2}{{\rm{S}}}_{1/2},{m}_{J}=0\leftrightarrow {}^{2}{{\rm{D}}}_{5/2},{m}_{J}=0$ transition, it is conceivable to not only obtain a reference with small statistical uncertainty, but also a low systematic uncertainty. In this case, the dynamically-decoupled Coulomb crystal clock would provide a reasonable accurate reference with superior stability.

We have proposed a CDD scheme that significantly suppresses the Zeeman shift as well as the quadrupole and tensor ac Stark frequency shifts of an optical clock transition for ion crystals. We analyzed the proposed scheme in the case of a multi-ion crystal of 400 ${\mathrm{Ca}}^{+}$ ions and showed that the shift of the robust transition is $\delta f\ll 1\,\mathrm{Hz}$ with a width of ∼1 Hz, which is close to the observed linewidth when probing the transition for a few hundred milliseconds. Our approach allows to exploit the improved stability from the higher ion number without suffering from the line broadening mechanisms associated with large ion crystals. Our scheme is applicable to other atomic species and experimental setups, paving the way for dynamically-decoupled Coulomb crystal clocks as references with high stability. We would like to note that during the preparation of this manuscript we became aware of a related independent work by Shaniv et al [85].

We thank M Barrett for stimulating discussions. We acknowledge support from the DFG through CRC 1128 (geo-Q), project A03 and CRC 1227 (DQ-mat), project B03. This joint research project was financally supported by the State of Lower-Saxony, Hannover, Germany. A R acknowledges the support of ERC grant QRES, project No. 77092.

In this section we present the construction of the (ideal) dressed states. By ideal we mean that there are no noise, uncertainties, or systematic shifts, the RWA is valid, and that the driving fields of the ${}^{2}{{\rm{S}}}_{1/2}$ $\left({}^{2}{{\rm{D}}}_{5/2}\right)$ do not operate on the ${}^{2}{{\rm{D}}}_{5/2}$ $\left({}^{2}{{\rm{S}}}_{1/2}\right)$ states.

A.1. The ${}^{2}{{\rm{D}}}_{5/2}$ states

The driving Hamiltonian of the ${}^{2}{{\rm{D}}}_{5/2}$ states is given by

$$\begin{eqnarray}\begin{array}{rcl}{H}_{D} & = & {g}_{d}{\mu }_{B}{{BS}}_{z}+{g}_{d}{{\rm{\Omega }}}_{1}\cos \left[\left({g}_{d}{\mu }_{B}B-{\delta }_{1}\right)t\right]{S}_{x}+{g}_{d}{{\rm{\Omega }}}_{2}\cos \left[\left({g}_{d}{\mu }_{B}B-{\delta }_{1}\right)t+\displaystyle \frac{\pi }{2}\right]\\ & & \times \cos \left[\left(\sqrt{{\left(\displaystyle \frac{{g}_{d}{{\rm{\Omega }}}_{1}}{2}\right)}^{2}+{\delta }_{1}^{2}}-{\delta }_{2}\right)t\right]{S}_{x},\end{array}\end{eqnarray} \tag{ A1 }$$

where g_d μ_BB is the Zeeman splitting due to the static magnetic field B, g_d = 6/5 is the gyromagnetic ratio of the ${}^{2}{{\rm{D}}}_{5/2}$ states, S_z and S_x are the z and x spin-5/2 matrices, Ω₁, Ω₂, ${\delta }_{1}=\sqrt{\tfrac{1}{8}}{g}_{d}{{\rm{\Omega }}}_{1}$ and δ₂ are the Rabi frequencies and the detunings of the driving fields, respectively. Moving to the IP with respect to the first drive (Ω₁) with ${H}_{01}=\left({g}_{d}{\mu }_{B}B-{\delta }_{1}\right){S}_{z}$ and assuming the RWA $\left({g}_{d}{\mu }_{B}B-{\delta }_{1}\gg {{\rm{\Omega }}}_{1}\right)$ we get

$$\begin{eqnarray}&&{H}_{D}^{{I}_{1}}={\delta }_{1}{S}_{z}+\displaystyle \frac{{g}_{d}{{\rm{\Omega }}}_{1}}{2}{S}_{x}+\displaystyle \frac{{g}_{d}{{\rm{\Omega }}}_{2}}{2}\cos \left[\left(\sqrt{{\left(\displaystyle \frac{{g}_{d}{{\rm{\Omega }}}_{1}}{2}\right)}^{2}+{\delta }_{1}^{2}}-{\delta }_{2}\right)t\right]{S}_{y}.\end{eqnarray} \tag{ A2 }$$

We continue by moving to the basis of the dressed states with ${U}_{1}={{\rm{e}}}^{{\rm{i}}{\theta }_{d}{S}_{y}}$, where ${\theta }_{d}=\arccos \left(\tfrac{{\delta }_{1}}{\sqrt{{\delta }_{1}^{2}+{\left(\tfrac{{g}_{d}{{\rm{\Omega }}}_{1}}{2}\right)}^{2}}}\right)$, which leads to

$$\begin{eqnarray}&&{H}_{D}^{{I}_{1}}=\sqrt{{\delta }_{1}^{2}+{\left(\displaystyle \frac{{g}_{d}{{\rm{\Omega }}}_{1}}{2}\right)}^{2}}{S}_{z}+\displaystyle \frac{{g}_{d}{{\rm{\Omega }}}_{2}}{2}\cos \left[\left(\sqrt{{\left(\displaystyle \frac{{g}_{d}{{\rm{\Omega }}}_{1}}{2}\right)}^{2}+{\delta }_{1}^{2}}-{\delta }_{2}\right)t\right]{S}_{y},\end{eqnarray} \tag{ A3 }$$

and then to the second IP with respect to ${H}_{02}=\left(\sqrt{{\left(\tfrac{{g}_{d}{{\rm{\Omega }}}_{1}}{2}\right)}^{2}+{\delta }_{1}^{2}}-{\delta }_{2}\right){S}_{z}$. Assuming the RWA, $\left(\sqrt{{\left(\tfrac{{g}_{d}{{\rm{\Omega }}}_{1}}{2}\right)}^{2}+{\delta }_{1}^{2}}-{\delta }_{2}\right)\gg {{\rm{\Omega }}}_{2}$, we obtain

$$\begin{eqnarray}&&{H}_{D}^{{I}_{2}}={\delta }_{2}{S}_{z}+\displaystyle \frac{{g}_{d}{{\rm{\Omega }}}_{2}}{4}{S}_{y}.\end{eqnarray} \tag{ A4 }$$

The eigenstates of ${H}_{D}^{{I}_{2}}$ are the double-dressed ${}^{2}{{\rm{D}}}_{5/2}$ states. The eigenstate with the smallest positive eigenvalue of $\tfrac{1}{2}\sqrt{{\delta }_{2}^{2}+{\left(\tfrac{{g}_{d}{{\rm{\Omega }}}_{2}}{4}\right)}^{2}}$ is used for the robust optical transition.

A.2. The ${}^{2}{{\rm{S}}}_{1/2}$ states

The driving Hamiltonian of the ${}^{2}{{\rm{S}}}_{1/2}$ states is given by

$$\begin{eqnarray}\begin{array}{rcl}{H}_{S} & = & {g}_{s}{\mu }_{B}{{Bs}}_{z}+{g}_{s}{{\rm{\Omega }}}_{3}\cos \left[\left({g}_{s}{\mu }_{B}B-{\delta }_{3}\right)t\right]{s}_{x}+{g}_{s}{{\rm{\Omega }}}_{4}\cos \left[\left({g}_{s}{\mu }_{B}B-{\delta }_{3}\right)t+\displaystyle \frac{\pi }{2}\right]\\ & & \times \cos \left[\left(\sqrt{{\left(\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}\right)}^{2}+{\delta }_{3}^{2}}-{\delta }_{4}\right)t\right]{s}_{x},\end{array}\end{eqnarray} \tag{ A5 }$$

where g_sμ_BB is the Zeeman splitting due to the static magnetic field B, g_s = 2 is the gyromagnetic ratio of the ${}^{2}{{\rm{S}}}_{1/2}$ states, s_z and s_x are the z and x spin-1/2 matrices, Ω₃, Ω₄, δ₃ and δ₄ are the Rabi frequencies and the detunings of the driving fields, respectively. Moving to the IP with respect to the first drive (Ω₃) with ${H}_{01}=\left({g}_{s}{\mu }_{B}B-{\delta }_{3}\right){s}_{z}$ and assuming the RWA $\left({g}_{s}{\mu }_{B}B-{\delta }_{3}\gg {{\rm{\Omega }}}_{3}\right)$ we get

$$\begin{eqnarray}&&{H}_{S}^{{I}_{1}}={\delta }_{3}{s}_{z}+\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}{s}_{x}+\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{4}}{2}\cos \left[\left(\sqrt{{\left(\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}\right)}^{2}+{\delta }_{3}^{2}}-{\delta }_{4}\right)t\right]{s}_{y}.\end{eqnarray} \tag{ A6 }$$

We continue by moving to the basis of the dressed states with ${U}_{1}={{\rm{e}}}^{{\rm{i}}{\theta }_{s}{s}_{y}}$, where ${\theta }_{s}=\arccos \left(\tfrac{{\delta }_{3}}{\sqrt{{\delta }_{3}^{2}+{\left(\tfrac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}\right)}^{2}}}\right)$, which leads to

$$\begin{eqnarray}&&{H}_{S}^{{I}_{1}}=\sqrt{{\delta }_{3}^{2}+{\left(\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}\right)}^{2}}{s}_{z}+\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{4}}{2}\cos \left[\left(\sqrt{{\left(\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}\right)}^{2}+{\delta }_{3}^{2}}-{\delta }_{4}\right)t\right]{s}_{y},\end{eqnarray} \tag{ A7 }$$

and then to the second IP with respect to ${H}_{02}=\left(\sqrt{{\delta }_{3}^{2}+{\left(\tfrac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}\right)}^{2}}-{\delta }_{4}\right){s}_{z}$. Assuming the RWA, $\left(\sqrt{{\delta }_{3}^{2}+{\left(\tfrac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}\right)}^{2}}-{\delta }_{4}\right)\gg {{\rm{\Omega }}}_{4}$, we obtain

$$\begin{eqnarray}&&{H}_{S}^{{I}_{2}}={\delta }_{4}{s}_{z}+\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{4}}{4}{s}_{y}.\end{eqnarray} \tag{ A8 }$$

The eigenstates of ${H}_{S}^{{I}_{2}}$ are the double-dressed ${}^{2}{{\rm{S}}}_{1/2}$ states. The eigenstate with the positive eigenvalue of $\tfrac{1}{2}\sqrt{{\delta }_{4}^{2}+{\left(\tfrac{{g}_{s}{{\rm{\Omega }}}_{4}}{4}\right)}^{2}}$ is used for the robust optical transition.

In this section we show how the expansions of the magnetic shifts, ${Z}_{{S}_{i}}$ and ${Z}_{{D}_{i}}$, the drive shifts, ${O}_{{S}_{i}}$ and ${O}_{{D}_{i}}$, and the correlated shifts, ${{ZO}}_{{S}_{i}}$ and ${{ZO}}_{{D}_{i}}$ are derived. In the derivation we assumed the RWA and neglected cross-driving effects. In the numerical simulations we took the counter-rotating terms of the driving fields and the cross-driving effect into account (see appendix C, D). The simulations were performed using the full driving Hamiltonian without making any approximations (see appendix E). For simplicity we will show the derivation for the ${}^{2}{{\rm{S}}}_{1/2}$ states. The derivation for the ${}^{2}{{\rm{D}}}_{5/2}$ follows the same steps.

We start by adding to the driving Hamiltonian of the ${}^{2}{{\rm{S}}}_{1/2}$ states, equation (A5), a magnetic noise term, which is given by ${g}_{s}{\delta }_{{\rm{B}}}{s}_{z}$. The drive shift is introduced by replacing Ω₃ and Ω₄ by ${{\rm{\Omega }}}_{3}\left(1+{\delta }_{{\rm{\Omega }}}\right)$ and ${{\rm{\Omega }}}_{4}\left(1+{\delta }_{{\rm{\Omega }}}\right)$, where δ_Ω represents a relative shift error of the driving fields. We assume that the relative errors of the driving fields are correlated since we expect that these errors are mostly due to changes in the amplifier chain and antenna, which are common to all drives. Moving to the IP with respect to the first drive as before and assuming the RWA $\left({g}_{s}{\mu }_{B}B-{\delta }_{3}\gg {{\rm{\Omega }}}_{3}\right)$, we now obtain

$$\begin{eqnarray}&&{H}_{S}^{{I}_{1}}={g}_{s}{\delta }_{{\rm{B}}}{s}_{z}+{\delta }_{3}{s}_{z}+\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{3}\left(1+{\delta }_{{\rm{\Omega }}}\right)}{2}{s}_{x}+\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{4}\left(1+{\delta }_{{\rm{\Omega }}}\right)}{2}\cos \left[\left(\sqrt{{\left(\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}\right)}^{2}+{\delta }_{3}^{2}}-{\delta }_{4}\right)t\right]{s}_{y}.\end{eqnarray} \tag{ B1 }$$

We continue by moving to the basis of the dressed states, including the shifts, with ${U}_{1}={{\rm{e}}}^{{\rm{i}}{\theta }_{s}{s}_{y}}$, where ${\theta }_{s}\,=\arccos \left(\tfrac{{\delta }_{3}+{g}_{s}{\delta }_{{\rm{B}}}}{\sqrt{{\left({\delta }_{3}+{g}_{s}{\delta }_{{\rm{B}}}\right)}^{2}+{\left(\tfrac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}\left(1+{\delta }_{{\rm{\Omega }}}\right)\right)}^{2}}}\right)$, which leads to

$$\begin{eqnarray}&&{H}_{S}^{{I}_{1}}=\sqrt{{\left({\delta }_{3}+{g}_{s}{\delta }_{{\rm{B}}}\right)}^{2}+{\left(\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}\left(1+{\delta }_{{\rm{\Omega }}}\right)\right)}^{2}}{s}_{z}+\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{4}\left(1+{\delta }_{{\rm{\Omega }}}\right)}{2}\cos \left[\left(\sqrt{{\left(\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}\right)}^{2}+{\delta }_{3}^{2}}-{\delta }_{4}\right)t\right]{s}_{y},\end{eqnarray} \tag{ B2 }$$

and then to the second IP with respect to ${H}_{02}=\left(\sqrt{{\delta }_{3}^{2}+{\left(\tfrac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}\right)}^{2}}-{\delta }_{4}\right){s}_{z}$. Assuming the RWA, $\left(\sqrt{{\delta }_{3}^{2}+{\left(\tfrac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}\right)}^{2}}-{\delta }_{4}\right)\gg {{\rm{\Omega }}}_{4}$, we obtain

$$\begin{eqnarray}&&{H}_{S}^{{I}_{2}}=\left({\delta }_{4}+\sqrt{{\left({\delta }_{3}+{g}_{s}{\delta }_{{\rm{B}}}\right)}^{2}+{\left(\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}\left(1+{\delta }_{{\rm{\Omega }}}\right)\right)}^{2}}-\sqrt{{\delta }_{3}^{2}+{\left(\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{3}}{2}\right)}^{2}}\right){s}_{z}+\displaystyle \frac{{g}_{s}{{\rm{\Omega }}}_{4}\left(1+{\delta }_{{\rm{\Omega }}}\right)}{4}{s}_{y}.\end{eqnarray} \tag{ B3 }$$

Assuming now that g_s = 2, the positive eigenvalue that we consider for the robust clock transition is given by

$$\begin{eqnarray}\begin{array}{rcl}{e}_{s} & = & \displaystyle \frac{1}{4}\left[4{\delta }_{4}\left(2\sqrt{{\left({\delta }_{{\rm{\Omega }}}+1\right)}^{2}{{\rm{\Omega }}}_{3}^{2}+{\left({\delta }_{3}+2{\delta }_{{\rm{B}}}\right)}^{2}}-2\sqrt{{\delta }_{3}^{2}+{{\rm{\Omega }}}_{3}^{2}}+{\delta }_{4}\right)\right.\\ & & -8\sqrt{\left({\delta }_{3}^{2}+{{\rm{\Omega }}}_{3}^{2}\right)\left({\left({\delta }_{{\rm{\Omega }}}+1\right)}^{2}{{\rm{\Omega }}}_{3}^{2}+{\left({\delta }_{3}+2{\delta }_{{\rm{B}}}\right)}^{2}\right)}\\ & + & 4{\left.\left({\delta }_{{\rm{\Omega }}}\left({\delta }_{{\rm{\Omega }}}+2\right)+2\right){{\rm{\Omega }}}_{3}^{2}+{\left({\delta }_{{\rm{\Omega }}}+1\right)}^{2}{{\rm{\Omega }}}_{4}^{2}+8\left({\delta }_{3}^{2}+2{\delta }_{3}{\delta }_{{\rm{B}}}+2{\delta }_{{\rm{B}}}^{2}\right)\right]}^{\tfrac{1}{2}}.\end{array}\end{eqnarray} \tag{ B4 }$$

Following the same calculations for the ${}^{2}{{\rm{D}}}_{5/2}$ states and assuming that g_d = 6/5 we obtain the lowest positive eigenvalue of the double-dressed ${}^{2}{{\rm{D}}}_{5/2}$ states, which is given by

$$\begin{eqnarray}&&{e}_{d}=\displaystyle \frac{1}{2}\sqrt{\displaystyle \frac{1}{100}{\left(3\sqrt{2\left(2{\delta }_{{\rm{\Omega }}}^{2}+4{\delta }_{{\rm{\Omega }}}+3\right){{\rm{\Omega }}}_{1}^{2}+16{\delta }_{{\rm{B}}}^{2}+8\sqrt{2}{\delta }_{{\rm{B}}}{{\rm{\Omega }}}_{1}}+10{\delta }_{2}-3\sqrt{6}{{\rm{\Omega }}}_{1}\right)}^{2}+\displaystyle \frac{9}{100}{\left({\delta }_{{\rm{\Omega }}}+1\right)}^{2}{{\rm{\Omega }}}_{2}^{2}}.\end{eqnarray} \tag{ B5 }$$

The magnetic shifts, ${Z}_{{S}_{i}}$ and ${Z}_{{D}_{i}}$, the drive shifts, ${O}_{{S}_{i}}$ and ${O}_{{D}_{i}}$, and the correlated shifts, ${{ZO}}_{{S}_{i}}$ and ${{ZO}}_{{D}_{i}}$, are obtained by the power series expansion of e_s and e_d to orders of ${\delta }_{{\rm{B}}}^{i}$ and ${\delta }_{{\rm{\Omega }}}^{i}$.

A drive of the ${}^{2}{{\rm{S}}}_{1/2}$ $\left({}^{2}{{\rm{D}}}_{5/2}\right)$ states also drives the ${}^{2}{{\rm{D}}}_{5/2}$ $\left({}^{2}{{\rm{S}}}_{1/2}\right)$ states off-resonantly and results in a Stark shift of the initial sub-levels energy gap. For simplicity we will show the derivation for the ${}^{2}{{\rm{S}}}_{1/2}$ states. The derivation for the ${}^{2}{{\rm{D}}}_{5/2}$ follows the same steps. Consider the off-resonant drive of the ${}^{2}{{\rm{S}}}_{1/2}$ states

$$\begin{eqnarray}&&{H}_{s}={\omega }_{0}^{s}{s}_{z}+{g}_{s}{\rm{\Omega }}\cos \left[\left({\omega }_{0}^{s}-{\rm{\Delta }}\right)t\right]{s}_{x},\end{eqnarray} \tag{ C1 }$$

where ${\omega }_{0}^{s}={g}_{s}{\mu }_{B}B$ is the Zeeman splitting and Δ is the detuning of the drive. We first move to the IP of the counter-rotating terms of the drive with respect to ${H}_{01}=-\left({\omega }_{0}^{s}-{\rm{\Delta }}\right){s}_{z}$. This results in

$$\begin{eqnarray}&&{H}_{s}^{{I}_{1}}=\left(2{\omega }_{0}^{s}-{\rm{\Delta }}\right){s}_{z}+\displaystyle \frac{{g}_{s}{\rm{\Omega }}}{2}\left[{s}_{x}+\displaystyle \frac{1}{2}\left({{\rm{e}}}^{{\rm{i}}\left(2{\omega }_{0}^{s}-{\rm{\Delta }}\right)t}{\sigma }_{+}+{{\rm{e}}}^{-{\rm{i}}\left(2{\omega }_{0}^{s}-{\rm{\Delta }}\right)t}{\sigma }_{-}\right)\right].\end{eqnarray} \tag{ C2 }$$

We continue by moving to the diagonalized basis of the time-independent part of ${H}_{s}^{{I}_{1}}$ with ${U}_{1}={{\rm{e}}}^{{\rm{i}}{\theta }_{s}{s}_{y}}$, where ${\theta }_{s}=\arccos \left(\tfrac{\left(2{\omega }_{0}^{s}-{\rm{\Delta }}\right)}{\sqrt{{\left(2{\omega }_{0}^{s}-{\rm{\Delta }}\right)}^{2}+{\left(\tfrac{{g}_{s}{\rm{\Omega }}}{2}\right)}^{2}}}\right)$, and then to a second IP of the rotating terms of the drive with respect to ${H}_{02}=2\left({\omega }_{0}^{s}-{\rm{\Delta }}\right){s}_{z}$. The time-independent part of ${H}_{s}^{{I}_{2}}$ is given by

$$\begin{eqnarray}&&{H}_{s}^{{I}_{2}}\approx \left(2\left({\rm{\Delta }}-{\omega }_{0}^{s}\right)+\sqrt{{\left({\rm{\Delta }}-2{\omega }_{0}^{s}\right)}^{2}+{\left(\displaystyle \frac{{g}_{s}{\rm{\Omega }}}{2}\right)}^{2}}\right){s}_{z}+\displaystyle \frac{{g}_{s}{\rm{\Omega }}}{4}\left(1+\displaystyle \frac{2{\omega }_{0}^{s}-{\rm{\Delta }}}{\sqrt{{\left(2{\omega }_{0}^{s}-{\rm{\Delta }}\right)}^{2}+{\left(\tfrac{{g}_{s}{\rm{\Omega }}}{2}\right)}^{2}}}\right){s}_{x},\end{eqnarray} \tag{ C3 }$$

and hence, the eigenvalues are equal to $\pm {\left[{\left(\left({\rm{\Delta }}-{\omega }_{0}^{s}\right)+\sqrt{{\left({\rm{\Delta }}-2{\omega }_{0}^{s}\right)}^{2}+{\left(\tfrac{{g}_{s}{\rm{\Omega }}}{2}\right)}^{2}}\right)}^{2}+\tfrac{{g}_{s}{\rm{\Omega }}}{8}{\left(1+\tfrac{2{\omega }_{0}^{s}-{\rm{\Delta }}}{\sqrt{{\left(2{\omega }_{0}^{s}-{\rm{\Delta }}\right)}^{2}+{\left(\tfrac{{g}_{s}{\rm{\Omega }}}{2}\right)}^{2}}}\right)}^{2}\right]}^{\tfrac{1}{2}}$, which gives the modified energy gap. Plugging in the parameters of the first ${}^{2}{{\rm{D}}}_{5/2}$ drive, (${{\rm{\Omega }}}_{1}$ and Δ = ω₀^s − ω₁, where ${\omega }_{1}={g}_{d}{\mu }_{B}B-{\delta }_{1}$ and g_s = 2, we obtained the modified energy gap of the ${}^{2}{{\rm{S}}}_{1/2}$ states, which is given by

$$\begin{eqnarray}&&{E}_{S}={\omega }_{1}+\sqrt{\displaystyle \frac{1}{4}{{\rm{\Omega }}}_{1}^{2}{\left(\displaystyle \frac{{\omega }_{0}^{s}+{\omega }_{1}}{\sqrt{{({\omega }_{0}^{s}+{\omega }_{1})}^{2}+{{\rm{\Omega }}}_{1}^{2}}}+1\right)}^{2}+{\left(\sqrt{{({\omega }_{0}^{s}+{\omega }_{1})}^{2}+{{\rm{\Omega }}}_{1}^{2}}-2{\omega }_{1}\right)}^{2}}.\end{eqnarray} \tag{ C4 }$$

For the ${}^{2}{{\rm{D}}}_{5/2}$ states the modified energy gap is given by

$$\begin{eqnarray}&&{E}_{D}={\omega }_{3}-\displaystyle \frac{1}{5}\sqrt{9{{\rm{\Omega }}}_{3}^{2}{\left(\displaystyle \frac{5\left({\omega }_{0}^{d}+{\omega }_{3}\right)}{2\sqrt{25{({\omega }_{0}^{d}+{\omega }_{3})}^{2}+9{{\rm{\Omega }}}_{3}^{2}}}+\displaystyle \frac{1}{2}\right)}^{2}+{\left(\sqrt{25{\left({\omega }_{0}^{d}+{\omega }_{3}\right)}^{2}+9{{\rm{\Omega }}}_{3}^{2}}-10{\omega }_{3}\right)}^{2}},\end{eqnarray} \tag{ C5 }$$

where here ${\omega }_{0}^{d}={g}_{d}{\mu }_{B}B$ and g_d = 6/5. Expanding the modified energy gaps in a power series of Ω₁ and Ω₃ results in ${E}_{S}\approx {\omega }_{0}^{s}+\tfrac{{\omega }_{0}^{s}{{\rm{\Omega }}}_{1}^{2}}{{\omega }_{0}^{{s}^{2}}-{\omega }_{1}^{2}}$ and ${E}_{D}\approx {\omega }_{0}^{d}-\tfrac{9{\omega }_{0}^{d}{{\rm{\Omega }}}_{3}^{2}}{25\left({\omega }_{3}^{2}-{\omega }_{0}^{{d}^{2}}\right)}$. Note that the second order correction could also be calculated by an effective Hamiltonian approach.

In this section we give a detailed derivation of the correction of the Bloch–Siegert shift arising from the counter-rotating terms of the driving fields. Without the correction, we first consider the dressed states due to the rotating-terms of a driving field and then consider the effect of the off-resonance counter-rotating terms on the dressed states. These result in an energy shift of the dressed states and (a time-dependent) amplitude-mixing between the dressed states, which is detrimental to the scheme, because it modifies the shifts of the dressed states considered for the robust transition. To correct this effect, we first consider the effect of the counter-rotating terms on the bare states, and then fix the frequency of the drive accordingly, such that the rotating-terms will drive the modified bare states. Consider, for example, the on-resonance driving Hamiltonian

$$\begin{eqnarray}&&{H}_{d}=\displaystyle \frac{{{\rm{\Omega }}}_{1}}{2}{\sigma }_{z}+\displaystyle \frac{{{\rm{\Omega }}}_{2}}{2}\cos \left({\omega }_{2}t\right){\sigma }_{x}.\end{eqnarray} \tag{ D1 }$$

Instead of moving to the IP of the rotating frame, we first move to the IP of the counter-rotating frame with respect to ${H}_{0}=-\tfrac{{{\rm{\Omega }}}_{1}}{2}{\sigma }_{x}$ and obtain

$$\begin{eqnarray}&&{H}_{{\rm{I}}}=\displaystyle \frac{{{\rm{\Omega }}}_{1}+{\omega }_{2}}{2}{\sigma }_{z}+\displaystyle \frac{{{\rm{\Omega }}}_{2}}{4}{\sigma }_{z}+\displaystyle \frac{{{\rm{\Omega }}}_{2}}{4}\left({\sigma }_{+}{{\rm{e}}}^{-2{\rm{i}}{\omega }_{2}t}+{\sigma }_{e}^{+2{\rm{i}}{\omega }_{2}t}\right).\end{eqnarray} \tag{ D2 }$$

We continue by moving to the diagonal basis of the time-independent part of H_I,

$$\begin{eqnarray}&&{H}_{{\rm{I}}}\approx \displaystyle \frac{1}{4}\sqrt{4{\left({{\rm{\Omega }}}_{1}+{\omega }_{2}\right)}^{2}+{{\rm{\Omega }}}_{2}^{2}}{\sigma }_{z}+\displaystyle \frac{{\tilde{{\rm{\Omega }}}}_{2}}{2}\left({\sigma }_{+}{{\rm{e}}}^{-2{\rm{i}}{\omega }_{2}t}+{\sigma }_{e}^{+2{\rm{i}}{\omega }_{2}t}\right).\end{eqnarray} \tag{ D3 }$$

If we choose $2{\omega }_{2}=\tfrac{1}{2}\sqrt{4{\left({{\rm{\Omega }}}_{1}+{\omega }_{2}\right)}^{2}+{{\rm{\Omega }}}_{2}^{2}}$ the rotating terms are on-resonance with the energy gap of the modified bare states. The on-resonance condition is given by

$$\begin{eqnarray}&&{\omega }_{2}=\displaystyle \frac{1}{6}\left(2{{\rm{\Omega }}}_{1}+\sqrt{16{{\rm{\Omega }}}_{1}^{2}+3{{\rm{\Omega }}}_{2}^{2}}\right).\end{eqnarray} \tag{ D4 }$$

In addition, due to the basis change from the basis of the bare states to the basis of the modified bare states, the Rabi frequency of the drive is slightly modified, ${{\rm{\Omega }}}_{2}\to {\tilde{{\rm{\Omega }}}}_{2}$, where

$$\begin{eqnarray}&&{\tilde{{\rm{\Omega }}}}_{2}=\displaystyle \frac{1}{2}{{\rm{\Omega }}}_{2}\left(1+\displaystyle \frac{2({\omega }_{2}+{{\rm{\Omega }}}_{1})}{\sqrt{4{\left({\omega }_{2}+{{\rm{\Omega }}}_{1}\right)}^{2}+{{\rm{\Omega }}}_{2}^{2}}}\right)\approx {{\rm{\Omega }}}_{2}-\displaystyle \frac{{{\rm{\Omega }}}_{2}^{3}}{16{\left({\omega }_{2}+{{\rm{\Omega }}}_{1}\right)}^{2}}.\end{eqnarray} \tag{ D5 }$$

Because ${\tilde{{\rm{\Omega }}}}_{2}$ corresponds to an optimal driving parameter, we must take that into account and set the initial Rabi frequency, Ω₂, accordingly.

In this section we provide a more detailed description of the numerical analysis of the scheme in the case of a robust multi-ion crystal clock that is presented in section 4 of the main text.

E.1. Optimized driving parameters

For the optimization of the goal function (see section 3 in main text and appendix B)

$$\begin{eqnarray}&&G=\displaystyle \sum _{i=1,j=1}^{i=4,j=2}| {Z}_{{S}_{i}}-{Z}_{{D}_{i}}| {\delta }_{{\rm{B}}}^{i}+| {O}_{{S}_{j}}-{O}_{{D}_{j}}| {\delta }_{{\rm{\Omega }}}^{j}+| {{ZO}}_{{S}_{j}}-{{ZO}}_{{D}_{j}}| {\delta }_{{\rm{B}}}^{j}{\delta }_{{\rm{\Omega }}},\end{eqnarray} \tag{ E1 }$$

we assume that the magnetic shift uncertainty δ_B is normally distributed with a zero mean and a width of 0.5 kHz (so the width of the ${}^{2}{{\rm{S}}}_{1/2}$ Zeeman splitting is 1 kHz), ${\delta }_{{\rm{B}}}\sim { \mathcal N }\left(0,0.5\right)$ kHz, and that the relative drive shift δ_Ω is normally distributed with a zero mean and a width of 4 × 10⁻⁴, ${\delta }_{{\rm{\Omega }}}\sim { \mathcal N }\left(0,4\times {10}^{-4}\right)$. Given these distributions of δ_B and δ_Ω, we defined an averaged goal function, ${G}_{A}=\langle G\left({\delta }_{{\rm{B}}}^{m},{\delta }_{{\rm{\Omega }}}^{n}\right){\rangle }_{m=1,n=1}^{m=100,n=100}$ over 100 realizations of δ_B and δ_Ω, which are chosen randomly according to the above distributions. We then numerically minimized G_A over the driving parameters, Ω_k and δ_k. The numerical minimization resulted in the following driving parameters (all in units of kHz): ${{\rm{\Omega }}}_{1}^{\star }=2\pi \times 225.3$, ${\delta }_{1}^{\star }=2\pi \times 95.6$, ${{\rm{\Omega }}}_{2}^{\star }=2\pi \times 13.6$, ${\delta }_{2}^{\star }=2\pi \times 5$, ${{\rm{\Omega }}}_{3}^{\star }=2\pi \times 93.6$, ${\delta }_{3}^{\star }=2\pi \times 27.2$, and ${{\rm{\Omega }}}_{4}^{\star }=2\pi \times 14.8$, ${\delta }_{4}^{\star }=2\pi \times 25.6$.

E.2. Fixing the driving parameters

When fixing the driving parameters we must take into account the effect of the ${}^{2}{{\rm{S}}}_{1/2}$ (${}^{2}{{\rm{D}}}_{5/2}$) drive on the ${}^{2}{{\rm{D}}}_{5/2}$ (${}^{2}{{\rm{S}}}_{1/2}$) states (appendix C), as well as the effect of the counter-rotating terms of a drive (appendix D) simultaneously. We start by parametrically solving the Bloch–Siegert correction equation, $2{\omega }_{i}=\tfrac{1}{2}\sqrt{{\left({\omega }_{i}+{\omega }_{0}\right)}^{2}+{{\rm{\Omega }}}_{i}^{2}}-{x}_{i}{\tilde{{\rm{\Omega }}}}_{i}$, for a general detuned drive, where x_i would correspond to the ratio between an optimal drive detuning and Rabi frequency, that is ${x}_{i}=\tfrac{{\delta }_{i}^{\star }}{{{\rm{\Omega }}}_{i}^{\star }}$. We denote the parametric solution of ω_i by ω_sol. Then, in order to fix the driving parameters of the first driving fields, ω₁, Ω₁, ω₃, and Ω₃, we numerically solve the following four equations:

$$\begin{eqnarray}&&{\omega }_{1}={\omega }_{{\rm{sol}}}\left({E}_{D}\left({\omega }_{3},{{\rm{\Omega }}}_{3}\right),{{\rm{\Omega }}}_{1}^{\star },{\delta }_{1}^{\star },{{\rm{\Omega }}}_{1}\right),\end{eqnarray} \tag{ E2 }$$

$$\begin{eqnarray}&&{\omega }_{3}={\omega }_{{\rm{sol}}}\left({E}_{S}\left({\omega }_{1},{{\rm{\Omega }}}_{1}\right),{{\rm{\Omega }}}_{3}^{\star },{\delta }_{3}^{\star },{{\rm{\Omega }}}_{3}\right),\end{eqnarray} \tag{ E3 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{1}^{\star }={\tilde{{\rm{\Omega }}}}_{1}\left({\omega }_{1},{E}_{D}\left({\omega }_{3},{{\rm{\Omega }}}_{3}\right),{{\rm{\Omega }}}_{1}^{\star },{\delta }_{1}^{\star },{{\rm{\Omega }}}_{1}\right),\end{eqnarray} \tag{ E4 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{3}^{\star }={\tilde{{\rm{\Omega }}}}_{3}\left({\omega }_{3},{E}_{S}\left({\omega }_{1},{{\rm{\Omega }}}_{1}\right),{{\rm{\Omega }}}_{2}^{\star },{\delta }_{3}^{\star },{{\rm{\Omega }}}_{3}\right).\end{eqnarray} \tag{ E5 }$$

For the case of μ_BB = $2\pi \times$ 5 MHz (so g_sμ_BB = $2\pi \times$ 10 MHz), this results in ω₁ = $2\pi \times$ 5.904 881 MHz, ω₃ = $2\pi \times$ 9.980 794 MHz, Ω₁ = $2\pi \times$ 225.3107 kHz, and Ω₃ = $2\pi \times$ 93.6311 kHz. The ideal driving parameters, without the above corrections are given by ω₁ = $2\pi \times$ 5.904 411 MHz, ω₃ = $2\pi \times$ 9.972 789 MHz, Ω₁ = $2\pi \times$ 225.3035 kHz, and Ω₃ = $2\pi \times$ 93.6306 kHz. For the second driving fields we neglect the cross-driving effect (modification of the energy gap) and take into account only the Bloch–Siegert corrections. We obtain that ω₂ = $2\pi \times$ 160.5892 kHz, ω₄ = $2\pi \times$ 72.0503 kHz, Ω₂ = $2\pi \times$ 13.6373 kHz, and Ω₄ = $2\pi \times$ 14.8157 kHz.

E.3. Simulations

We simulated the full driving Hamiltonian of the system with

$$\begin{eqnarray}\begin{array}{rcl}H & = & {g}_{s}{\mu }_{B}{{Bs}}_{z}+{g}_{s}{\delta }_{{\rm{B}}}{s}_{z}+{g}_{s}{{\rm{\Omega }}}_{3}\left(1+{\delta }_{{\rm{\Omega }}}\right)\cos \left[{\omega }_{3}t\right]{s}_{x}+{g}_{s}{{\rm{\Omega }}}_{4}\left(1+{\delta }_{{\rm{\Omega }}}\right)\cos \left[\left({\omega }_{3}t\right)+\displaystyle \frac{\pi }{2}\right]\cos \left[{\omega }_{4}t\right]{s}_{x}\\ & & +{g}_{s}{{\rm{\Omega }}}_{1}\left(1+{\delta }_{{\rm{\Omega }}}\right)\cos \left[{\omega }_{1}t\right]{s}_{x}+{g}_{s}{{\rm{\Omega }}}_{2}\left(1+{\delta }_{{\rm{\Omega }}}\right)\cos \left[\left({\omega }_{1}t\right)+\displaystyle \frac{\pi }{2}\right]\cos \left[{\omega }_{2}t\right]{s}_{x}\\ & & +{g}_{d}{\mu }_{B}{{BS}}_{z}+{g}_{d}{\delta }_{{\rm{B}}}{S}_{z}+{g}_{d}{{\rm{\Omega }}}_{1}\left(1+{\delta }_{{\rm{\Omega }}}\right)\cos \left[{\omega }_{1}t\right]{S}_{x}+{g}_{d}{{\rm{\Omega }}}_{2}\left(1+{\delta }_{{\rm{\Omega }}}\right)\cos \left[\left({\omega }_{1}t\right)+\displaystyle \frac{\pi }{2}\right]\cos \left[{\omega }_{2}t\right]{S}_{x}\\ & & +{g}_{d}{{\rm{\Omega }}}_{3}\left(1+{\delta }_{{\rm{\Omega }}}\right)\cos \left[{\omega }_{3}t\right]{S}_{x}+{g}_{d}{{\rm{\Omega }}}_{4}\left(1+{\delta }_{{\rm{\Omega }}}\right)\cos \left[\left({\omega }_{3}t\right)+\displaystyle \frac{\pi }{2}\right]\cos \left[{\omega }_{4}t\right]{S}_{x}\\ & & +\left[3{S}_{z}^{2}-{\rm{I}}\displaystyle \frac{5}{2}\left(\displaystyle \frac{5}{2}+1\right)\right]\left(Q+{Q}_{T}\right),\end{array}\end{eqnarray} \tag{ E6 }$$

where Q and Q_T represent the quadrupole and tensor shifts respectively. First, we simulated the system with δ_B = 0, δ_Ω = 0, and Q = Q_T = 0 for a time duration of 0.5 s, where we initialized the system in the equal superposition of $\tfrac{1}{\sqrt{2}}\left(| \,{S}_{1}\,\rangle +| \,{D}_{1}\,\rangle \right)$, where $| \,{S}_{1}\,\rangle$ and $| \,{D}_{1}\,\rangle$ are the double dressed ${}^{2}{{\rm{S}}}_{1/2}$ and ${}^{2}{{\rm{D}}}_{5/2}$ states with the smallest positive eigenvalue, respectively. This gave us a reference for the non-shifted transition frequency, ${\nu }_{S-D}$ . We then simulated 4927 realizations with ${\delta }_{{\rm{B}}}\sim { \mathcal N }\left(0,0.5\right)$ kHz, ${\delta }_{{\rm{\Omega }}}\sim { \mathcal N }\left(0,4\times {10}^{-4}\right)$, $Q\sim { \mathcal N }\left(-1,1\right)$ Hz, and ${Q}_{T}\sim { \mathcal U }\left(-3,0\right)$ Hz. Each simulation resulted in a frequency shift Δν with respect to ${\nu }_{S-D}$. The histogram ofΔν is shown in figure 6 in the main text.

In this section we give some details for the distribution of QPS and TASS.

In the Paul trap we consider for our scheme, the mean rf electrical field at the position of each ion increases linearly in dependence of the radial distance from the Paul trap symmetry axis. The TASS depends quadratically on the mean electrical field amplitude $(\tfrac{\delta f}{f}\propto \langle {E}^{2}\rangle \propto {r}^{2}$, see equation (11) in the main text ), so that overall we get a quadratic scaling with distance from the symmetry axis r. We plot the TASS in dependence on the radial distance from the Paul trap symmetry axis in figure F1, showing the expected quadratic behavior. The resulting distribution is shown in figure 5 of the main text (see also figure 3 in [27]). Note that while the electric field is linear in the radius, the frequency shifts, being quadratic in the field strength, vary slowly near the center of the crystal and more rapidly as one moves away from the line of symmetry. Thus, the width of an annulus corresponding to a given frequency interval decreases with increasing radius, roughly compensating for the increasing radial density of ions and leading to an approximately uniform density in frequency space as shown in figure 5 in the main text. For the purpose of implementation of the numerical simulation, we approximately describe the data with a uniform distribution.

Zoom In Zoom Out Reset image size

In the numerical simulation of the CDD we neglect possible correlations between TASS and QPS by modeling our system with random samples from corresponding independent distributions. In figure F2 we plot the QPS versus the TASS for each individual ion of our ion crystal. It is evident that the shifts are only weakly correlated. Furthermore, we note that the continuous decoupling drive renders any such correlations insignificant because the remnant effect of the QPS and TASS is a second order effect, and particulary in our case, this effect is much smaller than the effect of the shifts in the amplitudes of the drive, which are the dominant factor that determined the resulting total shift distribution of the robust optical transition (figure 6 of the main text).

Zoom In Zoom Out Reset image size