Robust optical clock transitions in trapped ions

We present a novel method for engineering an optical clock transition that is robust against external field fluctuations and is able to overcome limits resulting from field inhomogeneities. The technique is based on the application of continuous driving fields to form a pair of dressed states essentially free of all relevant shifts. Specifically, the clock transition is robust to magnetic shifts, quadrupole and other tensor shifts, and amplitude fluctuations of the driving fields. The scheme is applicable to either a single ion or an ensemble of ions, and is relevant for several types of ions, such as $^{40}\mathrm{Ca}^{+}$, $^{88}\mathrm{Sr}^{+}$, $^{138}\mathrm{Ba}^{+}$ and $^{176}\mathrm{Lu}^{+}$. Taking a spherically symmetric Coulomb crystal formed by 400 $^{40}\mathrm{Ca}^{+}$ ions as an example, we show through numerical simulations that the inhomogeneous linewidth of tens of Hertz in such a crystal together with linear Zeeman shifts of order 10~MHz are reduced to form a linewidth of around 1~Hz. We estimate a two-order-of-magnitude reduction in averaging time compared to state-of-the art single ion frequency references, assuming a probe laser fractional instability of $10^{-15}$. Furthermore, a statistical uncertainty reaching $2.9\times 10^{-16}$ in 1~s is estimated for a cascaded clock scheme in which the dynamically decoupled Coulomb crystal clock stabilizes the interrogation laser for an $^{27}\mathrm{Al}^{+}$ clock.


I. INTRODUCTION
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 −18 [1][2][3][4]. Taking advantage of these record uncertainties for applications ranging from relativistic geodesy [5][6][7][8] over fundamental physics [9][10][11] to frequency metrology [12][13][14][15][16] 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 × 10 −16 / √ τ [17] and by implementing an effectively dead-time-free clock consisting of two independent clocks probed in an interleaved fashion [18], reaching a statistical uncertainty of 0.6 × 10 −16 / √ τ. 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 [19] to levels of a few parts in 10 −15 / √ τ [3,20,21]. The statistical uncertainty can be improved by probing for longer times, ultimately limited by the excited clock state lifetime or the laser coherence time [22,23]. Alternatively, the number of probed ions can be increased.
Recently, multi-ion clock schemes have been proposed to address this issue [24][25][26][27]. 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 [28][29][30]. 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 ma-chined linear Paul trap with negligible excess micromotion from trap imperfections [26,31]. 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 In + or Al + , or by employing ring traps in which the QPS is the same for all ions [27]. A high-accuracy multiion 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 −19 level [24].
An alternative approach based on large 3d Coulomb crystals of ions has been theoretically investigated [25] for ion species for which the micromotion-induced Doppler shift and the ac Stark shift, both driven by the rf trapping field, can be made to cancel at a "magic" rf drive frequency [28]. This cancellation has been employed for single Sr + [32], single Ca + [33], and is currently being investigated for Lu + [25,34,35]. However, electronic states with J > 1/2 are subject to rank 2 tensor shifts, such as rf electric field-induced tensor ac Stark (TASS) and QPS [28,30]. For clouds of 100s 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 [25], or by operating at a judiciously chosen magnetic-field-insensitive point [35] for ions with hyperfine structure.
CDD or dressed-state engineering in the context of clocks has been proposed for linear Zeeman shift cancellation in neu-  (c) Illustration of the robust clock transition. While the unperturbed clock transition is subject to magnetic shifts, quadrupole shifts, and tensor shifts, the dressing fields substantially mitigate these shifts and result in a robust clock transition.
tral atom clocks through rf dressing of Zeeman substates in a regime where the drive Rabi frequency is larger than the drive frequency [43]. However, this scheme is unable to cancel tensorial shifts, such as the QPS or tensor Stark shifts. In a similar approach, magnetic field noise suppression up to second order in radio-frequency clocks via weak rf dressing has been proposed [44] and experimentally implemented to engineer synthetic clock states for rf spectroscopy in ultra-cold rubidium [63,64].
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 ion crystals containing on the order of 100 ions or more. We demonstrate the basic principle, which is illustrated in Fig. 1, on the 2 S 1/2 ↔ 2 D 5/2 clock transition in 40 Ca + , but the scheme is directly applicable to other systems as well. It therefore allows the operation of a multi-ion clock [24] 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 [22,23].
The paper is organized as follows. In Section II we show how robustness to QPS and TASS can be achieved by the application of a continuous detuned driving field. In Section III we present the general CDD scheme for the construction of a robust optical clock transition. Then, in Section IV we consider the implementation of the scheme in the case of a multiion crystal clock and analyze the performance of the scheme in terms of the expected statistical uncertainties of the robust optical clock transition. We end with the conclusions in Section V.

II. ROBUSTNESS TO TENSOR SHIFTS
In the absence of hyperfine structure, tensor shifts including the quadrupole shift are proportional to Q J,m j = J(J +1)− 3m 2 j [30], 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 ∑ J m J =−J Q J,m j = 0 [65]. 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 these type of shifts is achieved by the application of a detuned driving field, mixing all m J states to form a 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 D 5/2 states of e.g. Ca + ions, where g d µ B B is the Zeeman splitting due to the static magnetic field B, g d is the gyromagnetic ratio of the D 5/2 states, Ω 1 is the Rabi frequency of the driving field, and δ 1 is the detuning. Moving to the interaction picture (IP) with respect to H 0 = (g d µ B B − δ 1 )S z and taking the rotating-waveapproximation (RWA) ((g d µ B B − δ 1 ) Ω 1 ), results in Since the bare D 5/2 states have tensor shifts which are proportional to Q 5/2,±5/2 = −10, Q 5/2,±3/2 = +2, and Q 5/2,±1/2 = +8, the tensor shifts of the dressed states (the eigenstates of H I ) are proportional to with Hence, by choosing δ 1 = ± 1 8 (g d Ω 1 ), all of the dressed states have a zero (first order) tensor shift, Q = 0 (see Fig. 2). This can also be understood in the lab frame, in which the coupling Ω 1 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.
The cancellation of the tensor shift to first order can also be understood as follows. The tensor and quadrupole shift operatorQ = S 2 − 3S 2 z can also be written asQ = S 2 x + S 2 y − 2S 2 z . Moving to a rotated basis (dressed states basis) defined by z → cos(θ )z + sin(θ )x, x → cos(θ )x − sin(θ )z, and y → y, and neglecting purely off-diagonal terms we obtain thatQ ≈ S 2 z 1 − 3 cos 2 (θ ) + S 2 x 1 − 3 sin 2 (θ ) + S 2 y . To first order (diagonal terms) we have thatQ ≈ (S + S − + S − S + ) 2 − 3 sin 2 (θ ) + S 2 z 1 − 3 cos 2 (θ ) . Hence to first order, the tensor shift vanishes for cos (θ ) = 1 √ 3 , which is analogous to magic-angle spinning in solid-state NMR spectroscopy [66]. The detuned driving field results in dressed states, which are the eigenstates of Eq. 2. The rotation to the dressed states basis is given by U = e iθ S y , where cos (θ ) =

III. THE SCHEME
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 S 1/2 and D 5/2 states, the usual clock states in Ca + optical clocks. Four driving fields are employed, where two driving fields operate on the S 1/2 states, and two driving fields operate on the D 5/2 states. The first driving field of the D 5/2 states mitigates tensor shifts, including the QPS, as shown in the previous section. However, the D 5/2 states, as well as the S 1/2 states, are still sensitive to the linear Zeeman shift. Moreover, the 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 S 1/2 ↔ 2 D 5/2 transition is achieved. Hence, the driving scheme results in doublydressed S 1/2 and D 5/2 states where one S 1/2 ↔ D 5/2 transition (between the doubly-dressed states) is robust to magnetic shifts, quadrupole shifts, tensor shifts, and driving amplitude shifts (see Fig. 3). For the D 5/2 (S 1/2 ) states, the Rabi frequencies, Ω k , and the detunings, δ k , of the first and second driving fields are denoted by {Ω 1 , δ 1 } and {Ω 2 , δ 2 } ({Ω 3 , δ 3 } and {Ω 4 , δ 4 }), respectively. Similar to the derivation in Section II, 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 where S i and s i are the spin matrices of the D 5/2 and S 1/2 states and g d = 6/5 and g s = 2 are the gyromagnetic ratios of the D 5/2 and 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 A detailed derivation is given in Appendix I. We consider the doubly-dressed S 1/2 and D 5/2 states with the smallest positive eigenvalue as the robust optical clock states. By setting δ 1 = 1 8 (g d Ω 1 ) robustness to quadrupole and tensor shifts is attained. We denote by δ b and δ a magnetic shift and as the relative driving amplitude shift respectively. In order to achieve robustness to shifts in the magnetic and driving fields, we first add to the Hamiltonian the magnetic shift terms δ bS z + δ bs z , and driving shifts terms by Ω k → (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 (Eq. 6) now includes both δ b and δ , which modify the eigenvalues (the energies) of the double-dressed states, and hence the optical transition frequency. We then calculate the power series expansion of the eigenvalues to orders of δ b i and δ i . We denote the series expansion terms of the magnetic shift of the S 1/2 and D 5/2 states by Z S i δ b i and Z D i δ b i respectively. Similarly, we denote the series expansion terms of the amplitude driving shift of the S 1/2 and D 5/2 states by O S i δ i and O D i δ i respectively. The expansion terms of the correlated shifts are denoted by ZO S i δ b i δ and ZO D i δ b i δ for the S 1/2 and D 5/2 states respectively. We calculate the magnetic energy shifts, Z S i δ b i and Z D i δ b i , up to fourth order (i = 1, ..., 4), driving energy shifts, O S i δ i and O D i δ i , up to second order (i = 1, 2), and the correlated shifts, ZO S i δ b i δ and ZO D i δ b i δ (i = 1, 2) as function of the driving parameters, Ω i and δ i . A detailed derivation is given in Appendix II. We continue by defining a goal function, Given distributions of the shifts, δ b and δ , of a specific experimental set-up, we then define an averaged goal function, where δ b m and δ 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.
In the derivation of the optimal driving parameters we assumed the RWA and neglected cross-driving effects (for example, the effect of the S 1/2 driving field on the D 5/2 states). In the numerical analysis of the resulting shift distribution we took the counter-rotating terms of the driving fields and the cross-driving effect into account. In addition, the simulations where performed using the full driving Hamiltonian without making any approximations. For more details see Appendix sections C, D, and E.

IV. ROBUST MULTI-ION CRYSTAL CLOCK
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 Ca + as a widely used clock species [67][68][69] 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 [26,31]. 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 Ca + was measured to be Θ(3d, 5/2) = 1.83(1)ea 2 0 [70]. 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 [25,30,70] with the Euler angles {α, β , γ} defined as in [30]. The state dependence can be described by with a total QPS of ∆ f QPS = Θ h × f (α, β , γ) × 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 cancelled by the dynamical decoupling scheme, see below.
The shift originating from the space charge will in general lead to inhomogeneous line broadening. The space chargeinduced 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 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 [25]. For this reason, we consider a spherical trap configuration in which it is straightforward to obtain near isotropic trap frequencies ω r,x ≈ ω r,y ≈ ω a ≈ 2π×1 MHz. A slight deviation from complete symmetry is desirable to pin the orientation of the crystal and allow efficient Doppler cooling. 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 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 crystal ions by minimizing the pseudo-potential energy in the linear Paul trap. We next calculate the effective electric field tensor ∂ E i ∂ x j and the resulting quadrupole shift. In Fig. 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.
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 tensor ac Stark shift. While for 40 Ca + the scalar shift can be canceled by operating the trap at the "magic" drive frequency [33], the tensor ac Stark shift 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 ... over one cycle with the dc tensor polarizability α dc from [71]. Fig. 5 shows the result for our configuration, exhibiting an approximately uniform distribution with about 25 Hz width.
For the numerical optimization of the goal function Eq. 7, we consider a static magnetic field such that the Zeeman splitting of the 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 shifts of the driving fields δ to be normally distributed with a zero mean and a width of 4 × 10 −4 .
In order to test the performance of the scheme in the multiion spherical crystal configuration, we have realized the averaged goal function G A , which is given by Eq. 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): 6, Ω 2 = 2π × 13.6, δ 2 = 2π × 5, Ω 3 = 2π × 93.6, δ 3 = 2π × 27.2, and Ω 4 = 2π × 14.8, δ 4 = 2π × 25.6. Then, 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 took into account the effect of the S 1/2 D 5/2 drive on the D 5/2 S 1/2 states, as well as a correction to the counter-rotating terms of the driving fields (see Appendix Sections C, D, and E). The simulation results, shown in Fig. (6), indicate that the shift distribution of the robust transition has a narrow width of ∼ 1 Hz. 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.119 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 S 1/2 , m J = 0 ↔ 2 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 Fig. 7 were performed with the driving parameters used in the simulations of Fig. 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 D 5/2 states (assuming the RWA and no effect of the S 1/2 drive) we ex-pect a quadratic scaling. However, the S 1/2 drive and the counter-rotating terms of the D 5/2 drive result in an amplitude mixing between the ideal dressed 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. 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.
The coupling between the dressed 2 S 1/2 states and the dressed 2 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 S 1/2 , m J = + 1 2 ↔ 2 D 5/2 , m J = + 5 2 transition is reduced by a factor of 0.51 (0.3) for the transition between a single (double) dressed 2 S 1/2 state and a single (double) dressed 2 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.

V. CONCLUSIONS
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 containing on the order of 100 ions or more. We analyzed the proposed scheme in the case of a multi-ion crystal of 400 Ca + ions and showed that the shift of the robust transition is δ f 1 Hz with a width of ∼ 1 Hz, which is close to the observed linewidth when probing the transition for a few hundred milliseconds. Using the results of [23], we estimate the achievable statistical uncertainty of a 400 ion Ca + clock when probed with a flicker floor-limited probe laser at 10 −14 (10 −15 ) to be 5.8 × 10 −16 / τ/s (1.8 × 10 −16 / τ/s). This represents an order of magnitude improvement in instability over current single ion clocks [3,20,21], corresponding to a reduction in averaging time by a factor of 100. Through appropriate characterization of the residual line center shift away from an effective 2 S 1/2 , m J = 0 ↔ 2 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 our analysis we assumed that all driving fields suffer from driving fluctuations. Generating the second driving fields (Ω 2 and Ω 4 ) by a phase modulation, as proposed in [39], should result in stable driving fields with negligible amplitude fluctuations. In this case we expect a further improvement in the performance of the scheme.
One of the many potential applications of the proposed scheme is a cascaded clock [72][73][74] in which a clock laser σ l T Ca (s) σ (τ) Ca τ/s T Al (s) σ (τ) Al τ/s T CaAl (s) σ (τ) CaAl τ/s G  Table I. 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 Ca + ions using the described DD scheme (k = Ca), a single Al + ion (k = Al), and a cascaded scheme in which a single Al + ion is probed by a laser pre-stabilized through a cloud of 400 Ca + ions (k = CaAl). The reduction in averaging time to achieve a certain statistical measurement uncertainty is given by G. is first stabilized to an ensemble of Ca + ions to improve its phase coherence time and thus allow extended probing times [22,23] in a high-accuracy single-ion clock (e.g. Al + or Yb + ). To bridge the difference in clock transition frequencies, a transfer scheme using a frequency comb would be employed [74][75][76]. Table I shows the achievable instabilities of an Al + clock for different initial clock laser instabilities. We have used the results from Leroux et al. [23] 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 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 Ca + ions and the gain is reduced. Even higher gain can be obtained from larger crystals or reference atoms with narrower linewidth than Ca + , such as lutetium [25,34]. The cascaded clock scheme enables short averaging times with lasers that are commercially available. Furthermore, 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. We would like to note that during the preparation of this manuscript we became aware of a related independent work by Shaniv et al. [77].

VI. ACKNOWLEDGEMENTS
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.
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 rotating-wave-approximation (RWA) is valid, and that the driving fields of the S 1/2 D 5/2 do not operate on the D 5/2 S 1/2 states.

The D 5/2 states
The driving Hamiltonian of the D 5/2 states is given by where g d µ B B is the Zeeman splitting due to the static magnetic field B, g d = 6/5 is the gyromagnetic ratio of the D 5/2 states, S z and S x are the z and x spin-5/2 matrices, Ω 1 , Ω 2 , δ 1 = 1 8 g d Ω 1 and δ 2 are the Rabi frequencies and the detunings of the driving fields, respectively. Moving to the interaction picture (IP) with respect to the first drive (Ω 1 ) with H 01 = (g d µ B B − δ 1 ) S z and assuming the RWA (g d µ B B − δ 1 Ω 1 ) we get We continue by moving to the basis of the dressed states with U 1 = e iθ d S y , where θ d = arccos  , which leads to and then to the second IP with respect to H 02 = g d Ω 1 2 2 + δ 2 1 − δ 2 S z . Assuming the RWA, The eigenstates of H I 2 D are the double-dressed D 5/2 states. The eigenstate with the smallest positive eigenvalue of 1 2 is used for the robust optical transition.

The S 1/2 states
The driving Hamiltonian of the S 1/2 states is given by where g s µ B B is the Zeeman splitting due to the static magnetic field B, g s = 2 is the gyromagnetic ratio of the S 1/2 states, s z and s x are the z and x spin-1/2 matrices, Ω 3 , Ω 4 , δ 3 and δ 4 are the Rabi frequencies and the detunings of the driving fields, respectively. Moving to the interaction picture (IP) with respect to the first drive (Ω 3 ) with H 01 = (g s µ B B − δ 3 ) s z and assuming the RWA (g s µ B B − δ 3 Ω 3 ) we get We continue by moving to the basis of the dressed states with U 1 = e iθ s s y , where θ s = arccos   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 Sections C,D). The simulations were performed using the full driving Hamiltonian without making any approximations (see Appendix Section E). For simplicity we'll show the derivation for the S 1/2 states. The derivation for the D 5/2 follows the same calculations.
We start by adding to the driving Hamiltonian of the S 1/2 states, Eq. (A5), a magnetic noise term, which is given by g s δ bs z . The drive shift is introduced by replacing Ω 3 and Ω 4 by Ω 3 (1 + δ ) and Ω 4 (1 + δ ), where δ represents a relative error shift 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 (g s µ B B − δ 3 Ω 3 ) we now obtain We continue by moving to the basis of the dressed states, including the shifts, with U 1 = e iθ s s y , where θ s = arccos δ 3 +g s δ b √ (δ 3 +g s δ b) 2 +(Ω 3 (1+δ )) 2 , which leads to and then to the second IP with respect to H 02 = δ 2 3 + Ω 2 3 − δ 4 s z . Assuming the RWA, δ 2 3 + Ω 2 3 − δ 4 Ω 4 , we obtain The positive eigenvalue, which we consider for the robust clock transition, is given by Following the same calculations for the D 5/2 states we obtain the lowest positive eigenvalue of the double-dressed D 5/2 states, which is given by 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 δ b i and δ i .
Appendix C: Modified energy gaps from cross-driving A drive of the S 1/2 D 5/2 states also drives the D 5/2 S 1/2 states off-resonantly and results in a Stark shift of the initial sub-levels energy gap. For simplicity we'll show the derivation for the S 1/2 states. The derivation for the D 5/2 follows the same calculations. Consider the off-resonant drive of the S 1/2 states, where ω 0 = g s µ B B is the Zeeman splitting and δ is the detuning of the drive. We first move the the IP of the counter-rotating terms of the drive with respect to H 01 = − (ω 0 − δ ) s z . This results in We continue by moving to the diagonalized basis of the time-independent part of H I 1 s with U 1 = e iθ s s y , where   , and then to a second IP of the rotating terms of the drive with respect to H 02 = 2 (ω 0 − δ ) s z .
The time-independent part of H I 2 s is given by and hence, the eigenvalues are equal to ± , which gives the modified energy gap. Plugging in the parameters of the first D 5/2 drive, (Ω 1 and δ = ω 0 − ω 1 , where ω 1 = g d µ B B − δ 1 and g s = 2 we obtained the modified energy gap of the S 1/2 states, which is given by For the D 5/2 states the modified energy gap is given by where here ω 0 = g d µ B B. Expanding the modified energy gaps in a power series of Ω 1 and Ω 3 results in E S ≈ ω 0 + . 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. The correction can be understood as follows. 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. This results 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 Instead of moving to the interaction picture (IP) of the rotating frame we first move to the IP of the counter-rotating frame with respect to H 0 = − Ω 1 2 σ x and obtain We continue by moving to the diagonal basis of the time-independent part of H I , If we choose 2ω 2 = 1 2 4(Ω 1 + ω 2 ) 2 + Ω 2 2 the rotating terms are on-resonance with the energy gap of the modified bare states. The on-resonance condition is given by 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, Ω 2 →Ω 2 , wherẽ BecauseΩ 2 corresponds to an optimal driving parameter, we must take that into account and set the initial Rabi frequency, Ω 2 , accordingly.

Appendix E: Numerical analysis
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 IV of the main text.

Fixing the Driving parameters
When fixing the driving parameters we must take into account the effect of the S 1/2 (D 5/2 ) drive on the D 5/2 (S 1/2 ) states (Section C), as well as the effect of the counter-rotating terms of a drive (Section D) simultaneously. We start by parametrically solving for a general detuned drive the Bloch-Siegert correction equation, 2ω i = 1 2 (ω i + ω 0 ) 2 + Ω 2 i − x iΩi , where x i would correspond to the ratio between an optimal drive detuning and Rabi frequency, that is x i = δ i Ω i . We denote the parametric solution of ω i by ω sol . Then, in order to fix the driving parameters of the first driving fields, ω 1 ,Ω 1 ,ω 3 , and Ω 3 , we numerically solve the following four equations: For the case of µ B B = 5 MHz (so g s µ B B = 10 MHz), this results in ω 1 = 5.904881 MHz, ω 3 = 9.980794 MHz, Ω 1 = 225.3107 kHz, and Ω 3 = 93.6311 kHz. The ideal driving parameters, without the above corrections are given by ω 1 = 5.904411 MHz, ω 3 = 9.972789 MHz, Ω 1 = 225.3035 kHz, and Ω 3 = 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 = 160.5892 kHz, ω 4 = 72.0503 kHz, Ω 2 = 13.6373 kHz, and Ω 4 = 14.8157 kHz.