Reducing the instability of an optical lattice clock using multiple atomic ensembles

The stability of an optical atomic clock is a critical figure of merit for almost all clock applications. To this end, much optical atomic clock research has focused on reducing clock instability by increasing the atom number, lengthening the coherent interrogation times, and introducing entanglement to push beyond the standard quantum limit. In this work, we experimentally demonstrate an alternative approach to reducing clock instability using a phase estimation approach based on individually controlled atomic ensembles in a strontium (Sr) optical lattice clock. We first demonstrate joint Ramsey interrogation of two spatially-resolved atom ensembles that are out of phase with respect to each other, which we call"quadrature Ramsey spectroscopy,"resulting in a factor of 1.36(5) reduction in absolute clock instability as measured with interleaved self-comparisons. We then leverage the rich hyperfine structure of ${}^{87}$Sr to realize independent coherent control over multiple ensembles with only global laser addressing. Finally, we utilize this independent control over 4 atom ensembles to implement a form of phase estimation, achieving a factor of greater than 3 enhancement in coherent interrogation time and a factor of 2.08(6) reduction in instability over an otherwise identical single ensemble clock with the same local oscillator and the same number of atoms. We expect that multi-ensemble protocols similar to those demonstrated here will result in reduction in the instability of any optical lattice clock with an interrogation time limited by the local oscillator.


I. INTRODUCTION
Thanks to their remarkable precision and accuracy, optical atomic clocks are rapidly advancing the frontiers of timekeeping, quantum science, and fundamental physics [1].State-of-the-art optical clocks have now reached the level of 10 −18 in both stability and accuracy [2][3][4][5][6], enabling novel emerging applications such as relativistic geodesy, searches for ultralight dark matter, and gravitational wave detection [7][8][9][10][11][12][13][14][15][16].Instability is a critical figure of merit for optical atomic clocks, determining the sensitivity of clock comparisons to target signals in fundamental physics applications and the achievable precision in applications to timekeeping, metrology, and navigation [1,13].The fundamental limit to the instability of an optical clock with unentangled atoms is given by the quantum projection noise (QPN), where ν is the clock transition frequency, C is the contrast, T is the coherent interrogation time, T d is the dead time, τ is the averaging time, and N is the number of atoms being interrogated.Equation 1 implies that the clock instability can be reduced through longer coherent interrogation times.However, in optical clocks both conventional Ramsey and Rabi spectroscopy are limited by frequency noise of the local oscillator (LO), or clock laser, used to interrogate the atoms.In the case of Ramsey spectroscopy, the interrogation time is limited by the requirement that the * kolkowitz@berkeley.eduphase accumulated by the atoms in the rotating frame must remain within [−π/2, π/2] in order to avoid phase slips.In addition, for atomic clock operations with duty cycles less than 1, noise on the LO during the dead time leads to Dick noise, an aliased noise at harmonic frequencies of 1/T c (T c = T + T d corresponds to the cycle time) that compromises the clock instability [17,18].In particular, the Dick-limited clock instability for Ramsey spectroscopy can be approximated as (see Appendix F for details) where σ LO is the flicker frequency noise floor limited instability of the LO.From Eq. 2, we can see that for an atomic clock limited by Dick noise, the instability no longer scales with atom number, but can still be reduced by increasing the interrogation time.These factors have motivated considerable efforts into the development of optical LOs with longer coherence times [19][20][21][22][23], the use of synchronous differential comparisons to push the interrogation time beyond the coherence time of the LO [24][25][26][27], the use of multiple separate optical atomic clocks to actively perform feedback or feedforward on the LO [30], as well as the use of entangled atomic states to push beyond the QPN limit for the same LO coherence time [31][32][33], although the latter case will only be useful in scenarios where the clock is QPN-limited.
A proposed alternative and complementary path to reducing the instability of optical atomic clocks is to more efficiently allocate the available atomic resources.For example, in Ref. [34], Rosenband and Leibrandt argued that if the atoms in a clock are split up into multiple ensembles and probed with different interrogation times in parallel, the clock instability can be reduced exponentially with atom number rather than with 1/ √ N as in Eq. 1 for uniform interrogation.Related work on efficient atomic clock operations with several atomic ensembles was also proposed by Borregaard and Sørensen [35].Similarly, in a recent theoretical work [36], Li et al. argued that through joint interrogation of two atom ensembles that are 90 degrees out of phase, it should be feasible to extend the interpretable phases accumulated by the atoms from [−π/2, π/2] to [−π, π], thereby doubling the achievable Ramsey interrogation time for a given LO and reducing the corresponding clock instability.
There have been several prior and recent works on combining separate optical atomic clocks to realize a single compound clock with reduced instability and longer coherent interrogation times [25,[28][29][30].For example, in Dörscher et al. [28], dynamical decoupling sequences were used along with near synchronous interrogation of two separate optical atomic clocks in order to overcome the laser coherence limit, extending the coherent interrogation time by a factor of approximately 6.5.Bowden et al. [29] performed QND measurements of one ensemble using an optical cavity to extend the interrogation time of a second ensemble in a separate apparatus by a factor of 7. A simplified three-ensemble version of the Rosenband and Leibrandt scheme was also recently demonstrated by Kim et al. [30] using two separate interleaved Yb optical lattice clocks to realize zero-dead-time operation and feed-forward on the interrogation of an additional Al + single ion clock, achieving lifetime limited interrogation times.Related techniques have also been implemented in atom interferometers to increase their dynamic range [37].However, the use of multiple independent op-tical clock apparatus, including vacuum chambers, atom sources, and associated optics, is extremely resource intensive and is not easily scaled to larger numbers of ensembles, and any compound clock would still benefit from reductions in the instability of the constituent clocks, motivating work into realizing similar gains within individual optical clocks through the use of multiple atomic ensembles.
In this work, we demonstrate a multi-ensemble strontium optical lattice clock that takes advantage of up to 4 spatially-resolved atom ensembles within a single vacuum chamber and the same optical lattice to increase the coherent interrogation time and reduce the absolute clock instability.We first experimentally demonstrate joint interrogation on 2 out-of-phase atom ensembles [36], which we call "quadrature Ramsey spectroscopy", which enables the measurement of phase shifts of up to [−π, π], and reduces the clock instability by a factor of 1.36 (5), as measured in a self-comparison.We then demonstrate a novel approach which leverages the rich hyperfine structure of 87 Sr (I = 9/2) to enable independent coherent control of four atom ensembles in order to implement a simple form of phase estimation algorithm [38][39][40].We prepare four atom ensembles in pairs of two distinct nuclear spin states, ( 1 S 0 , m F = 5/2 and 3 P 0 , m F = 3/2 where m F is the hyperfine sub-level), and demonstrate independent coherent control of each ensemble pair on two different clock transitions through frequency multiplexing.Combining the four-ensemble scheme and the quadrature Ramsey technique, we further improve the interrogation time by a factor of roughly 1.7, achieving a factor of 2 reduction in instability over a single ensemble clock with the same atom number, duty cycle, and LO.
We note that our multi-ensemble protocols will introduce additional systematic frequency shifts that will require additional characterization.However, none of these shifts are new or unique to optical lattice clocks, and we do not expect any of them to fundamentally limit accuracy above the 10 −18 level (see Appendix J for details).These additional systematics will likely preclude the near term adoption of multi-ensemble protocols in the most accurate state-of-the-art laboratory-based clocks.We instead anticipate that multi-ensemble protocols such as the ones we demonstrate here will be first be adopted in applications involving transportable clocks, where stability is a more critical figure of merit and where record setting accuracy is already impractical.

II. EXPERIMENTAL APPARATUS AND QUADRATURE RAMSEY SPECTROSCOPY
The experimental apparatus used here is similar to the multiplexed strontium optical lattice clock presented in prior works [27,41].A movable, one-dimensional optical lattice is loaded with multiple atom ensembles, for example two atom ensembles at a 1 cm separation along the lattice axis (ẑ) as shown in Fig. 1(a), followed by nuclear spin-polarization through optical pumping and in-lattice-cooling.A bias magnetic field applied along x defines the quantization axis, with a typical magnitude of 5.5 G.A tunable magnetic field gradient along ẑ (∆B(z)), ranging from 0 to 45 mG/cm, is applied to introduce a differential frequency shift, dominated by the first-order Zeeman shift, and thus a differential phase shift for a Ramsey sequence of duration T of ∆φ ≈ γ∆B(z)T , where γ is the differential gyromagnetic ratio for the addressed clock transition, between the two ensembles.To take full advantage of the magnetic field gradient, we probe on the magnetically sensitive 3 P 0 , m F = 7/2 → 1 S 0 , m F = 5/2 (noted as |e, 7/2 → |g, 5/2 ) transition, which has a magnetic field sensitivity of roughly 564 Hz/G.Synchronous differential Ramsey spectroscopy is performed to characterize the differential phase ∆φ [27], and the magnetic field gradient is fine-tuned so that ∆φ = π/2, where the Ramsey fringes for the two ensembles are 90 degrees out-ofphase as shown in Fig. 1(b) in order to realize quadrature Ramsey spectroscopy (Appendix B).A typical parametric plot taken at a short free evolution time (T = 7.5 ms) is shown in Fig. 1(c), where synchronous Ramsey spectroscopy is performed by randomly sampling the phase of the final global π/2-pulse and each labeled quadrature corresponds to the matching region in Fig. 1(b) (see Appendix C for details).The differential phase has a typical long-term differential instability below 10 −19 as demonstrated in prior work [27].
A standard Ramsey spectroscopy sequence begins with a π/2-pulse that prepares the atoms in a 50:50 superposition of |g and |e .The atoms then evolve freely for time T and accumulate a phase θ, |g + e iθ |e , in which θ = 2πT ∆f and ∆f is the detuning of the clock laser from atomic resonance.A final π/2-pulse converts θ into a population difference and generates the error signal for clock feedback (Fig. 2(a)).However, phases beyond [−π/2, π/2] result in population differences that are no longer unique, producing an erroneous error signal.The interrogation times for conventional Ramsey spectroscopy are therefore typically constrained to well below the LO coherence time in order to avoid a failure in the feedback servo due to phase slips [42], an example of which is shown in the inset to Fig. 2(b).Clock operation with Ramsey spectroscopy is performed on the 1 S 0 → 3 P 0 clock transition.Because we lack a second independent optical lattice clock to perform direct comparisons against, we instead perform a self-comparison, where two independent atomic servos using the same experimental sequences are interleaved and compared in order to determine the clock instability (prior examples can be found in Refs.[43][44][45]).Separate digital servos correct the clock laser frequency independently for each interleaved sequence and generate a difference frequency signal which is then used to extract the Allan deviation and corresponding instability.Ramsey spectroscopy on the |e, 7/2 → |g, 5/2 transition is performed simultaneously on two spatially resolved atom ensembles, each containing roughly 4000 atoms.The interrogation time in our system is primarily limited by frequency noise of the LO, which limits our measured atom-laser coherence time to roughly 100 ms [27].The typical dead time per cycle in our system is about 1.5 s, which is primarily limited by atom preparation (roughly 1 s), and camera readout (roughly 250 ms), so there is a total of roughly 3 s of dead time for each clock in the self-comparison.The process of loading multiple ensembles contributes only roughly 100 ms of extra dead time (see Appendix D for details) and is not a limiting factor for the duty cycle in our apparatus.
In order to perform a fair comparison between the quadrature Ramsey scheme and standard Ramsey spectroscopy using the same total atom number and experimental sequence, the magnetic field gradient ∆B(z) is tuned such that the differential phase between the two ensembles is either ∆φ = 0 (standard Ramsey) or ∆φ = π/2 (quadrature Ramsey) (Fig. 2(a)).In the case of ∆φ = 0, the measurements from the two ensembles are simply averaged together and standard clock operation using Ramsey spectroscopy is performed, with the measured phase given by where P = (P 1 + P 2 )/2 is the average of the excitation fractions from both clocks, and P 0 ≈ 0.5 and C ≈ 0.95 are the measured half-maximum and contrast of the Ramsey fringe, respectively.In this case the two ensembles are equivalent to a single ensemble clock with 2N total atoms.In the case of ∆φ = π/2, quadrature Ramsey spectroscopy is employed, with the excitation fractions (P 1 , P 2 ) used as inputs for the phase decoder R(P 1 , P 2 ), which determines the quadrature of the acquired phase according to [36], in which θ 1 = sin −1 (2(P 1 − P 0 )/C), θ 2 = cos −1 (2(P 2 − P 0 )/C).We perform a series of self-comparisons for both ∆φ = 0 and π/2 cases by varying the interrogation time T , with each measurement averaging down for roughly an hour.Each self-comparison yields Allan deviations which are linearly fit assuming only white frequency noise (that averages down as τ −1/2 ), which is on the order of 1 × 10 −14 / √ τ level and is dominated here by the Dick noise from the LO.The extracted Allan deviations are plotted as a function of T in Fig. 2(b).In the case of ∆φ = 0, we find that phase slips begin to occur somewhere between the measurements at T = 20 and T = 30 ms, indicated by the cutoff of the shaded region in Fig. 2(b), which we conservatively place near T = 30 ms.This is consistent with Monte Carlo simulations of the phase slip probabilities based on the measured LO noise spectrum (Appendix H).Phase slips lead to the sudden accumulation of frequency differences between the two servos at multiples of the Ramsey fringe linewidth in Fig. 2(b).As a result, outside the shaded region in Fig. 2(b) the self-comparison no longer averages down and no meaningful Allan deviation can be defined.In the case of quadrature Ramsey spectroscopy with ∆φ = π/2, we are able to extend the interrogation time to T = 50 ms before we begin to observe phase slips in the self-comparison.The optimal clock performance for quadrature Ramsey spectroscopy in our apparatus is achieved at T = 40 ms with an instability of 9.5(3) × 10 −15 / √ τ , a factor of 1.36(5) smaller than the instability measured at T = 20 ms with conventional Ramsey interrogation (Fig. 2(c)).

III. HYPERFINE STATE INITIALIZATION AND INDIVIDUAL COHERENT CONTROL OVER ENSEMBLE PAIRS THROUGH GLOBAL ADDRESSING
While quadrature Ramsey spectroscopy with two ensembles enables a doubling of the coherent interrogation time, it cannot be naively extended beyond this by adding more ensembles because Ramsey fringes are inherently 2π periodic.In order to take advantage of additional ensembles and further extend the coherent interrogation time, different interrogation times are required for different samples, as proposed in Refs.[34,35], which constitute a form of a phase estimation algorithm [38][39][40].The lattice is loaded with 4 evenly-spaced atomic ensembles with 0.33 cm nearest-neighbor separation, and a linear magnetic field.Right: Relevant energy levels and the pulse sequence used to prepare the ensembles in different hyperfine states.All four ensembles are initially prepared in the stretched hyperfine ground state 1 S0, mF = 9/2 , or |g, 9/2 .A first π-pulse (i) transfers the population into 3 P0, mF = 7/2 (or |e, 7/2 ), followed by a π/2-pulse (ii) on the |e, 7/2 → |g, 5/2 transition that prepares a 50:50 superposition.A second π-pulse (iii) prepares a superposition of two nuclear spins (|e, 3/2 and |e, 7/2 ).During the free evolution time Tp, an applied field gradient (about 45 G/cm) introduces a phase shift of 0, π, 2π, and 3π for each atom ensemble, respectively.Final π-pulse (iv) and π/2-pulse (v) prepare the ensembles into opposite states, with ensembles (1, 3) on |g, 5/2 and ensembles (2, 4) on |e, 3/2 .(b) Global Rabi spectroscopy of the four atom ensembles after state preparation at a bias magnetic field of about 5.5 G. (c) Independent coherent Rabi oscillations for ensembles (1, 3) (orange, red) on the |g, 5/2 → |e, 7/2 transition and for ensembles (2, 4) (blue, light blue) on the |e, 3/2 → |g, 1/2 transition, respectively.Each plot represents the average of 4 independent scans, with each data point indicating the mean.Error bars correspond to 1σ standard deviation.The solid curves correspond to fits to the expected sinusoidal Rabi oscillations, and are used to bound the cross-talk between the two hyperfine clock transitions.
Potential approaches would be to address different spatially resolved ensembles in multiple optical lattices with different focused clock beams, or even to have different ensembles in different vacuum chambers.However, both approaches are challenging and would suffer from added differential phase [46,47] and intensity noise on the clock beams, as well as additional frequency shifts due to the differing environments of each ensemble.It is therefore desirable to have all of the atomic ensembles share the same lattice, vacuum chamber, and clock beam.Fortunately, the magnetic field gradient that enables quadrature Ramsey spectroscopy can be combined with the rich hyperfine level structure of 87 Sr (I = 9/2) to achieve independent clock interrogation with multiple ensembles using only global addressing of all the ensembles with a single clock beam.
The basic principle and relevant pulse sequence is shown in Fig. 3(a).Similar to the two-ensemble case, we start with a sample of 4 evenly spaced atom ensembles, each containing roughly 2000 atoms, in the stretched state |g, 9/2 .We first transfer the population into |e, 7/2 with a clock π-pulse (i), then create a 50:50 superposition between |e, 7/2 and |g, 5/2 using a clock π/2-pulse (ii).We then transfer half of the population into |e, 3/2 (iii), and let the sample evolve for T p , such that each ensemble accumulates phase shifts of 0, π, 2π, and 3π, respectively, due to the applied linear magnetic field gradient (at about 45 mG/cm).We then transfer the populations back to |g, 5/2 (iv), and convert the accumulated phases into population in either the |g, 5/2 and |e, 3/2 states with a final clock π/2-pulse (v).Ensembles 1 and 3 are left in |g, 5/2 , while ensembles 2 and 4 are left in |e, 3/2 (Fig. 3(a)).This state preparation sequence typically takes about 250 ms yielding an overall state preparation fidelity of about 90% limited by the clock π-pulse fidelity (> 95% per pulse).Cleanup pulses are then applied to remove unwanted residual populations in each of ensemble pairs.

IV. PHASE ESTIMATION WITH FOUR ATOMIC ENSEMBLES
Having achieved independent control with global addressing through a combination of spatial and frequency multiplexing, we proceed to demonstrate a combination of quadrature Ramsey spectroscopy [36] with a simple phase estimation algorithm [34,35] using four atomic ensembles.This approach allows us to extend the coher-ent interrogation time by a factor of 1.7 over quadrature Ramsey with only two ensembles, and thereby further reduce the clock instability (see Appendix G 1).
We interrogate a first ensemble pair to keep track of the laser phase evolution within the range of [−π, π] during T A using quadrature Ramsey, while simultaneously interrogating a second pair of ensembles on a different hyperfine clock transition with quadrature Ramsey for T B (T B > T A ).At the end of the measurement we perform global readout on all four ensembles, and then use a simple phase estimation algorithm to unwrap the phase accumulated by the second pair during T B .As a result, we are able to unambiguously measure the phase accumulated for longer coherent interrogation times that would otherwise lead to phase slips due to the modulo 2π ambiguity of Ramsey fringes.
The basic principle is shown in Fig. 4(a), in which we first apply a global π/2-pulse on the 3/2 − 1/2 transition for ensemble pair (2,4), and subsequently apply another global π/2-pulse on the 7/2 − 5/2 transition for (1,3).The delay between the asynchronous pulses for (1, 3) and (2, 4) is a few ms, limited by the rise and fall times of the clock pulses.As the two transitions are spectrally resolved, each pulse only impacts the targeted pair, aside from a negligibly small probe shift (see Appendix J 2).
We let ensembles (1, 3) freely evolve for time T A , with the magnetic field gradient ∆B(z) tuned such that the differential phase shift between ensembles (1, 3) is π/2, and a second global π/2-pulse is applied to the 7/2 − 5/2 transition and the phase accumulated by the ensembles (1, 3) is stored in the populations of the clock states and is shelved there temporarily.After a total free precession time of T B we apply the second π/2-pulse on the 3/2−1/2 transition for ensembles (2,4).The free precession time T B is chosen such that the differential phase of ensembles (2, 4) is again π/2, with T B = 1.7 × T A due to the difference in magnetic moment for the two clock transitions, resulting in two quadrature Ramsey measurements with interrogation times T A and T B = 1.7 × T A .
The populations of all four ensembles are then simultaneously read out using global imaging pulses, where the excitation fractions are used as inputs to the phase estimator.In particular, we obtain two phases from the two quadrature Ramsey measurements: θ A = R(P 1 , P 3 ) and θ B = R(P 2 , P 4 ).We constrain θ A within [−π, π], and we have θ A = 2π∆f T A , in which ∆f is the LO detuning.The actual accumulated phase shift θ act,B during T B can fall outside the range of [−π, π].Consequently, it cannot be handled by the two-ensemble quadrature Ramsey decoder alone (Eq.4) due to the 2π ambiguity.To unambiguously determine θ act,B , we employ a phase estimation algorithm that uses the measurement of θ A to remove the 2π ambiguity in θ B .This estimator is valid in the limit that the LO noise is dominated by low frequency noise and for a low clock duty cycle (see Appendix E 2 for details).We then feedback on the LO according to ∆f = θ act,B /(2πT B ).
The use of two quadrature Ramsey sequences with dif- 4. Phase estimation using four atom ensembles with independent control through frequency multiplexing.(a) Basic principle and pulse sequence for clock operation using phase estimation with four ensembles [34][35][36].We interrogate on the |e, 7/2 → |g, 5/2 transition for clock pair (1, 3) for a free evolution time of TA, during which the magnetic field gradient is chosen such that the phase difference between ensembles 1 and 3 is π/2.Asynchronously, we interrogate on the |e, 3/2 → |g, 1/2 transition for TB = 1.7 × TA, which is chosen such that the phase difference between ensembles 2 and 4 is also π/2 due to the difference in magnetic sensitivities for the two transitions.An estimator extracts the LO detuning using information from all 4 clocks, which is then used to feed back on the LO.(b) Comparison of inferred phase shifts using various decoders for a long interrogation time.The data is taken at (TA, TB) = (50, 85) ms in a single four-ensemble experiment for 1000 measurement cycles.The inferred phase shifts are post-processed and extracted using various decoders.Top: inferred with standard Ramsey decoder using only ensemble 2; Middle: inferred with quadrature Ramsey decoder using ensembles ferent interrogation times allows us to perform phase estimation and discriminate accumulated phases during T B for ensembles (2,4) beyond [−π, π].A representative plot from 1,000 measurement cycles of the inferred phase shifts from the longest interrogation times (T A = 50 ms, T B = 85 ms) without measured phase slips in the selfcomparison using the four-ensemble approach is shown in Fig. 4(b).We note that the inferred phase shifts are postprocessed using different decoders applied to the same 1,000 measurements.The distributions of inferred phases acquired by ensemble 2 using information from only a single ensemble (ensemble 2) with a standard Ramsey decoder (top, red), two ensembles (ensembles 2 and 4) using a quadrature Ramsey decoder (middle, yellow), and all four ensembles (ensembles 1, 2, 3, and 4) using the full phase estimation decoder (bottom, blue) are compared.Instances of 199 inferred phases beyond [−π/2, π/2] in the case of two ensembles, and 8 inferred phases beyond [−π, π] in the case of four ensembles can be observed respectively, which would otherwise result in corresponding phase slips for the single ensemble or two-ensemble decoders for this interrogation time.
As shown in Fig. 4(c), the clock self-comparison Allan deviation at 1 s from all three approaches (standard Ramsey, quadrature Ramsey with two ensembles, and quadrature Ramsey combined with phase estimation with four ensembles) are compared using identical experimental sequences and atom numbers, with the only differences being the applied field gradients and how the accumulated phases are extracted by the estimator.We are able to extend the coherent interrogation time T B out to 85 ms (with T A = 50 ms) without observing phase slips, and achieve a factor of 1.2(1) reduction in instability when comparing the four-ensemble phase estimation scheme to the two-ensemble quadrature Ramsey scheme, and a factor of 2.08 (6) reduction with respect to standard Ramsey.We note that the instabilities are slightly degraded over those shown in Fig. 2(b) for the same interrogation times, which we attribute to the added time required for the hyperfine state preparation (250 ms), which contributes a total additional dead time of 500 ms in the self-comparisons.

V. DISCUSSION AND OUTLOOK
We note that all the instabilities and Allan deviations quoted in this work are the instabilities extracted from the self-comparisons, and that the actual single clock instabilities are expected to be at least a factor of 2 lower due to the product of the factor of √ 2 from comparing two clocks and the factor of √ 2 from the doubling in cycle time for each clock due to the interleaved nature of the self-comparison.While the instabilities demonstrated here are not competitive with the state-of-the-art optical lattice clocks to date [3,4,6,11,14], this is primarily due to the relatively short coherence time of our LO (100 ms), which is compounded by the need for a self-comparison as we do not have another independent clock to compare with, resulting in very low duty-cycles (< 2%).These factors are not fundamental, and we anticipate that our multiple-ensemble techniques can be extended to other existing state-of-the-art Sr optical lattice clocks [6,11,14] in a straightforward fashion and that similar reductions in clock instability would be expected.
For example, with state-of-the-art optical cavities that enabled atom-laser coherence time on the order of 600 ms and duty-cycle over 50% [19][20][21][22][23], a factor of 2 improvement in the coherence time using quadrature Ramsey spectroscopy with two ensembles should be fairly straightforward to implement (Appendix I), offering roughly a factor of √ 2 reduction in instability.However, given the resulting higher duty-cycle (> 50%), more of the frequency drift of the local oscillator will take place during the Ramsey interrogation time rather than during the dead time, so the four-ensemble phase estimation scheme employed here would likely be insufficient.We anticipate that 6 ensembles split into three quadrature Ramsey pairs, with one pair interrogated for 2T in parallel to two pairs each interrogated back-to-back for T , would be sufficient to achieve a comparable enhancement of roughly a factor of 4 in coherent interrogation time and a factor of 2 reduction in instability.However, this would require the use of an additional hyperfine clock transition or an alternative mechanism for achieving independent control over three ensemble pairs.
It is also natural to consider scaling up to many more atomic ensembles, which would eventually require an alternative approach to achieving independent control that does not require a unique clock transition for each pair of ensembles.This could be achieved in optical lattice clocks through the introduction of larger magnetic field gradients to spectroscopically resolve each ensemble on the same clock transition [48], or targeted light shifting beams to bring individual ensembles in and out of resonance with the global clock beam [49,50].The resulting scaling of clock instability with the number of ensembles M will depend on the noise spectrum of the LO, the number of atoms available, and the clock duty cycle, but both the coherent interrogation time and clock instability can be expected to continue to improve with the addition of more ensembles as long as the coherent interrogation time is limited by the LO.
While we leave the determination of the optimal algorithm and number of ensembles for a given optical atomic clock for future work, a naive extension of the 6 ensemble algorithm discussed in the previous paragraph can be expected to provide a scaling of the instability with 1/ √ M for large M , so long as the clock remains Dick noise limited.If combined with non-destructive mid-circuit measurement [51][52][53][54][55], the use of multiple independently controlled atom ensembles could ultimately offer exponential reductions in clock instability with atom and ensemble number using the scheme put forward by Rosenband and Leibrandt [34].Furthermore, we note that similar multi-ensemble phase estimation algorithms have been proposed to achieve Heisenber-limited clock performance using maximally entangled atomic Greenberger-Horne-Zeilinger (GHZ) states [38], hybrid coherent and squeezed states [39], and Gaussian spin-squeezed states [40].
Although the multiple-ensemble protocols demonstrated in this work improve the Ramsey interrogation time and reduce the clock instability, they do come with some trade-offs.First, in this work we incur an additional 100 ms of dead time due to multiple-ensemble preparation, although this could likely be further optimized and reduced.Second, our approach introduces additional systematic shifts that could impact the achievable clock accuracy.Examples of potential effects include Zeeman shifts arising from the magnetic field gradient applied to generate the π/2 differential phase shift for quadrature Ramsey spectroscopy, probe shifts from the off-resonant clock pulses in the four-ensemble phase estimation approach, and the servo error associated with the differential phase shift in the quadrature Ramsey decoder.Estimates of the magnitude of the associated systematic shifts and potential mitigation strategies are discussed in Appendix J.While a full systematic evaluation with multi-ensemble protocols is beyond the scope of this work, we estimate that the effects on achievable accuracy will be manageable at the 10 −18 level.The added systematic effects therefore likely preclude the near term adoption of these techniques in optical clocks that are pushing the ultimate limits of accuracy, such as for the redefinition of the SI second.Instead, we expect that multi-ensemble techniques will initially prove most useful in situations where improved stability is beneficial and ultra-high accuracy is already either not critical or not possible, for example applications involving portable clocks such as positioning, navigation, and communication [1,[56][57][58].
Finally, we note that while this work has focused exclusively on protocols for extending the coherent interrogation times for Ramsey spectroscopy, many state-ofthe-art optical clocks instead employ Rabi spectroscopy, and there are significant trade-offs between the two approaches.To the best of our knowledge, prior works have also focused on Ramsey-like protocols to extend the interrogation times and perform phase estimation both theoretically [34-36, 38-40, 42] and experimentally [25,[28][29][30], and we are not aware of similar proposals for protocols making use of Rabi spectroscopy.The closest example we are aware of is the combination of high-resolution imaging and Rabi spectroscopy, as demonstrated in Marti et al. [59], where imaging spectroscopy serves as a multiplexed measurement of laser frequency noise with a long Rabi pulse, but in this case many of the atoms in the spatially extended ensemble contribute zero information to the measurement.We therefore leave consideration of extensions of our approaches to Rabi spectroscopy for future work.

VI. CONCLUSION
In conclusion, in this work we experimentally demonstrated an optical lattice clock with enhanced Ramsey interrogation times and reduced absolute clock instability by harnessing multiple atomic ensembles in a single clock.We first demonstrated quadrature Ramsey spectroscopy, a simple protocol for jointly interrogating two atom ensembles that are 90 degrees out of phase in order to extend the phase interrogation window from [−π/2, π/2] to [−π, π] [36], extending the coherent interrogation time without phase slips by roughly a factor of 2× and measuring a corresponding reduction in clock instability by a factor of 1.36 (5).We then leveraged the hyperfine structure of 87 Sr and a linear magnetic field gradient to achieve independent state preparation of multiple ensembles with only global addressing.Finally, we demonstrated a combination of quadrature Ramsey spectroscopy and enhanced phase estimation [34,35] with four ensembles to further extend the coherent interrogation time by a factor of 1.7 compared to the two-ensemble scheme, resulting in a factor of 2.08(6) reduction in clock instability compared to a standard single ensemble clock with the same atom number, LO, and duty-cycle.Our multi-ensemble approach could be extended to other existing optical lattice clocks with state-of-the-art LOs [19][20][21][22][23], and similar reductions in instability would be anticipated.In addition, we anticipate that multi-ensemble phase estimation protocols like the ones demonstrated in this work will be required to take full advantage of the metrological gains offered by entanglement such as in the use of spin-squeezed states [31][32][33][38][39][40].
atoms to a temperature of 1 µK.We then load multiple ensembles into a movable one-dimensional (1D) optical lattice as described and demonstrated in prior works [27,41], followed by spin-polarization and inlattice-cooling both axially (sideband cooling) and radially (Doppler cooling) on the narrow-linewidth 1 S 0 → 3 P 1 transition at 689 nm.The typical lattice trap depth for loading and cooling is at 70 E rec , where E rec /h ≈ 3.5 kHz is the recoil energy of a lattice photon at the magic wavelength around 813.4 nm.The lattice is then adiabatically ramped down to the operational trap depth at 15 E rec for state preparation and clock interrogation, after which it is ramped back up to 70 E rec for read-out.

Appendix B: Differential phase shift generation and calibration
A magnetic field of roughly 5.5 G is applied along x to define the quantization axis using 3 pairs of orthogonal Helmholtz bias coils, (B x , B y , B z ).The coil set-points are optimized to zero out B y and B z at the waist of the lattice, z = 0. To introduce the π/2 differential phase shift between a pair of ensembles at z = ±0.5 cm, we employ a magnetic field gradient which results in a differential Zeeman shift between the two ensembles.The magnetic field gradient can be tuned by scanning an offset voltage applied to both B z bias coils, as shown in Fig. 5.To take full advantage of the magnetic field gradient, we interrogate on the |e, 7/2 → |g, 5/2 clock transition, which has a magnetic field sensitivity of about γ = 564 Hz/G.The differential frequency shift is dominated by the differential first-order Zeeman shift as, ∆f ZS, where ∆B is the magnetic field gradient.This corresponds to a typical differential frequency shift of 25 Hz for a magnetic field gradient of ∆B = 45 mG/cm.The magnetic field gradient is then fine-tuned such that the differential phase ∆φ = 2π∆f ZS,1 st T equals to π/2 at each interrogation time T .In particular, we scan the B z offset set-points and map out ∆φ using synchronous Ramsey spectroscopy [27], in which the two ensembles are prepared in a 50:50 superposition of 1 S 0 and 3 P 0 , freely evolve for T and project using a final global π/2-pulse with its phase (φ 2nd ) randomly sampled from [−π, π].An ellipse fitting is then fit to the data using the least-squares method to extract ∆φ, which has a typical uncertainty of ∼ 10 mrad, limited by the QPN.The resulting ∆φ as a function of B z offset voltage is fitted using a quadratic function, and the intercept of the fit with π/2 is used as the B z offset set-point for a specific Ramsey evolution time T , as shown in Fig. 6.
We note that the ∂B z /∂z gradient potentially rotates the quantization axis along the lattice direction, which would introduce sizable tensor Stark shift gradient.In our previous work [41], we have characterized the tensor Stark shift gradient to be −8(1) × 10 −20 /E rec /cm on the 1 S 0 , m F = ±5/2 ↔ 3 P 0 , m F = ±3/2 transition.This corresponds to a frequency gradient of −1.2(2) × 10 −18 /cm at the operational lattice depth of 15 E rec .In future work, it appears to be feasible to mitigate this effect by generating a more uniform linear B x gradient through the addition of Golay coils [61], requiring an upgrade to our experimental apparatus that is currently underway.In Fig. 1(c), we perform synchronous Ramsey spectroscopy by randomly sampling the phase (φ 2nd ) of the final π/2-pulse from a uniform distribution within [−π, π].This is equivalent to randomly sampling the LO detunings corresponding to ∆f i ∈ [−1/2T, 1/2T ] under a short Ramsey free evolution time (T = 7.5 ms), where each quadrature can be mapped to a specific region as shown in Fig. 1(b).To quantify the accuracy of the decoder, we perform a point-by-point comparison between the phase shifts extracted using the quadrature Ramsey decoder (Eq.4) and the known randomly sampled phase shifts applied to the final π/2-pulse, as shown in the red points in Fig. 7.The blue points in Fig. 7 correspond to Monte-Carlo simulations using QPN and the LO noise spectrum as input parameters, and are in good agreement with the experiment.

Appendix D: Self-comparison clock interrogation with quadrature Ramsey protocol
Before clock interrogation, the atomic population is transferred from 1 S 0 , m F = 9/2 → 3 P 0 , m F = 7/2 (note as |g, 9/2 → |e, 7/2 ) via a π-pulse, followed by a clean-up pulse on resonant with the 461 nm 1 S 0 → 1 P 1 transition to blow out residual populations in the ground state.Quadrature Ramsey spectroscopy is then per-formed by interrogating the |e, 7/2 → |g, 5/2 clock transition, which has a first-order Zeeman shift sensitivity of 564 Hz/G.To generate an error signal for frequency feedback, the phase of the final π/2-pulse is detuned by π/2 with respect to the phase of the first π/2-pulse.After spectroscopy, the lattice is adiabatically ramped back up to full depth for readout.The populations in the ground and excited clock states of both ensembles are read-out in parallel with imaging pulses along the lattice axis, with scattered photons collected on an EMCCD camera (Andor, iXon-888).The excitation fraction is extracted through P = (N e −N bg )/(N g +N e −2N bg ), where N g , N e and N bg are the ground state population, excited clock state population, and background counts without atoms, respectively.The excitation fractions from both ensembles are used as inputs to the phase estimator to extract the clock laser detuning, which is used to feed back to the acousto-optical modulator on the clock laser path that addresses the atomic resonance.
Due to the lack of an additional independent optical lattice clock to compare against, in order characterize the clock instability we perform a self-comparison.Two asynchronous independent servos sharing the same LO are interleaved on a single experiment apparatus.Separate digital servos correct the clock laser frequency independently for each servo and generate a difference frequency signal which is used to extract the Allan deviation and infer the clock instability.A typical experimental timing sequence for the self-comparison is shown in Fig. 8.The typical dead time is about 1.5 s within each servo, and thus 3 s per clock cycle when accounting for the other servo.The primary contribution comes from sample loading, which takes about 1 s to cool, trap, and optically pump the atoms, with an additional 100 ms required to prepare multiple ensembles rather than a single ensemble.The camera imaging and data processing takes about 250 ms, while an additional 250 ms is spent on the nuclear spin state preparation when performing phase estimation with four atom ensembles, which is not shown in Fig. 8 Appendix E: Phase estimation and quadrature Ramsey with four atom ensembles

Hyperfine state initialization
The experimental sequence used to initialize four atom ensembles is shown in Fig. 3(a) in the main text.We elaborate on some of the experimental details here.First, we start with a sample of spin-polarized atoms in the |g, 9/2 state with four ensembles spaced evenly over a 1 cm extent.A first π-pulse on the 9/2 − 7/2 transition transfers the population into |e, 7/2 , followed by a clean-up pulse on the 461 nm 1 S 0 → 1 P 1 transition.We then prepare the atoms in a 50:50 superposition of the |e, 7/2 and |g, 5/2 state with a clock π/2-pulse, and coherently transfer the population in |g, 5/2 into |e, 3/2 via a π-pulse on resonance with the 5/2 − 3/2 transition.This creates a superposition of two nuclear spin states in the excited clock state manifold, |e, 7/2 and |e, 3/2 , which freely evolve for roughly T p = 40 ms under a linear magnetic field gradient tuned to about 45 mG/cm to maximize the differential phase shifts between the adjacent atom ensembles.We then coherently transfer the populations and close the interferometer via a π-pulse on the 7/2−5/2 transition and a final π/2-pulse on 5/2−3/2 transition, respectively.Pulses on resonance with either the 1 S 0 → 1 P 1 transition or the 3 P 0,2 → 3 S 1 transitions are combined with clock π-pulses to clean up residual populations leftover in unwanted states in each pair of ensembles.T p and the laser phase of the last π/2-pulse are fine-tuned such that the four atom ensembles accumulate 0, π, 2π, and 3π phase shifts, respectively, enabling preparation of two ensemble pairs in distinct nuclear spin states, pair (1, 3) in |g, 5/2 and pair (2,4) in |e, 3/2 , using only global pulses without the need for individual ensemble addressing.
Relevant transitions, their differential magnetic moments, and the measured Rabi frequencies used for state initialization and clock interrogation are summarized in Table I.As described in Sec.IV of the main text, by probing asynchronously on two hyperfine clock transitions, we obtain two phases from the two quadrature Ramsey measurements: θ A = R(P 1 , P 3 ) and θ B = R(P 2 , P 4 ).We constrain θ A within [−π, π] by choosing T A ≤ 50 ms, and we have θ A = 2π∆f T A , in which ∆f is the LO detuning.Given that T B = 1.7 × T A , the actual phase shift θ act,B accumulated during T B may exceed the range of [−π, π], which cannot be handled by the quadrature Ramsey decoder (Eq.4) due to the 2π ambiguity.Consequently, this results in θ B = θ act,B − k × 2π(k = 0, ±1) that could lead to phase slips.
To ascertain θ act,B , we utilize the following phase estimator to determine the value of k that resolves the 2π ambiguity.In the limit that the LO noise is dominated by low frequency noise and for a low duty cycle, we can assume that the frequency of the LO did not significantly change during the second half of the interrogation period T B , and we approximate the estimated phase shift accumulated by ensembles (2,4) as θ est,B = 1.7 × θ A .In our estimator, we compare θ est,B with θ B + k × 2π(k = 0, ±1) and identify the particular value of k that yields the smallest difference.We then use θ act,B = θ B + k × 2π to feedback on the LO according to ∆f = θ act,B /(2πT B ).

Appendix F: Scaling of the Dick-noise-limited clock instability
The instability of a clock limited by Dick noise is given by [17] where T c is the cycle time, m is the m-th harmonic, S f LO (m/T c ) is the one-sided power spectral density of the relative frequency fluctuations of the LO at Fourier frequencies m/T c , and the parameters g 0 , g s m , and g c m are defined as: with g(ξ) being the sensitivity function of Ramsey spectroscopy sequence.
In the limit that T π/2 ≪ T (where T π/2 is the duration of π/2-pulse), we can approximate the sensitivity function of a Ramsey spectroscopy with free evolution time T as The LO noise is modeled as in which b 0 corresponds to white frequency noise, b −1 corresponds to flicker frequency noise, and b −2 corresponds to random walk frequency noise.
Here we derive the scaling of Dick-limited clock instability as function of interrogation time T for an arbitrary duty-cycle.As pointed out in Refs.[17,18], for an LO dominated by its flicker frequency noise floor, , where σ LO is the LO flicker frequency instability, the Dick noise limited clock instability can be approximated as: where T c is the cycle time.
We then have the parameters g 0,1 as and thus after some algebraic manipulation we have: Plugging this back into Eq.F4, we arrive at It's worth noting that taking only the lowest order term m = 1, we have Eq.F9 as which matches the approximation given in Eq. 12 of Ref. [18].We can thus fit our data (as taken in Fig. 2(b) and Fig. 4(c) in the main text) to the scaling in Eq.F9.However, we need to bound the maximum m max in the actual numerical fitting given limited computational resources, and due to the finite bandwidth of the Ramsey sensitivity function.The cutoff m max can be roughly bounded by the multiples of sine waves within the minimum Ramsey free evolution time T min as in which we conservatively choose n = 2.In our actual data taking, we have a typical π/2-pulse duration of T π/2 = 1.5 ms, and a minimum T min = 2 ms.To account for the finite pulse duration, we can take the average of sensitivity function g ′ (ξ) over the entire duration of (0, T + 2T π/2 ), where we have where ΩT π/2 = π/2, with Ω being the Rabi frequency.For T π/2 = 1.5 ms and T = 2 ms, this yields an effective duration of T ′ min = 3.9 ms, giving m max = 1700 in our numerical fitting.
In order to fit to the instability data shown in both Fig. 2(b) and Fig. 4(c), we use a simplified expression for the Dick noise that assumes only flicker frequency noise, and treat the instability of the LO as a free parameter.The extracted LO flicker frequency instability of σ FF = 3.2(2) × 10 −15 , which is roughly a factor of 2 − 3× larger than the value provided by the manufacturer of the LO (Menlo Systems).This discrepancy likely arises in part from not including the contributions of random walk and white frequency noise in the fit.In addition, the LO is now in a relatively noisy vibration and thermal environment, and we believe its noise spectrum is now degraded (see next section for details).
Appendix G: Reducing clock instability by increasing interrogation times with multiple ensembles

Theoretically expected clock instability
The expected clock instability is given by the quadrature sum of QPN and Dick noise, In particular, QPN scales with the Ramsey interrogation time (T ) and atom number (N , per ensemble) as (Eq. 1) The Dick noise limited clock instability scales with T roughly as (Eq.F9) in which m max ≈ 1700 in our case.Both Eq.G2 and Eq.G3 suggest that increasing the interrogation time T reduces the clock instability.
In our case, the clock instability is Dick noise limited due to the LO noise (σ LO ≈ 10 −15 ) and relatively low duty-cycle (T /T c < 2%), and thus has no atom number dependence given that N is sufficiently large (4000 atoms per ensemble in quadrature Ramsey with two ensembles, and 2000 atoms per ensemble in phase estimation protocol with four ensembles).We note that, however, even in the case of QPN limited, we are still taking advantage of the information from all the atoms and should not suffer from reduced atom numbers with the use of multiple ensembles.
For standard Ramsey spectroscopy with a single ensemble, the interrogation time is T 1 , limited by requiring that the phase of the LO remain within the range of [−π/2, π/2].For quadrature Ramsey protocol with two ensembles, the interrogation time is increased to T 2 = 2×T 1 by extending the interpretable LO phase shift to [−π, π].Experimentally, we achieve a reduction of 1.36 (5) in the measured clock instability that is roughly consistent with theoretically expected reduction in the Dick noise limit (Fig. 2(b) in the main text).
For the phase estimation protocol with four ensembles we demonstrate in Fig. 4, the interrogation time of the second ensemble pair can be further extended by a factor of 1.7 compared to the case of two ensembles as T 4 = 1.7 × T 2 , or a factor 3.4 compared to the case of a single ensemble as T 4 = 3.4 × T 1 .And we achieve a reduction of 1.2(1) in the measured clock instability when compared to the quadrature Ramsey protocol, and a reduction of 2.08(6) when compared to the standard Ramsey approach, both are again roughly consistent with the theoretical expectations (Fig. 4(c) in the main text).

Numerical calculations of the clock instability
The QPN limited instability (σ QPN ) for Ramsey spectroscopy is given in Eq. 1 in the main text.To calculate σ Dick (Eq.F9), the LO frequency noise spectrum is modeled with Eq.F3, and the coefficients b m (m = 0, −1, −2) are extracted from the measured cavity instabilities provided by the manufacturer by beating the LO against another independent cavity-stabilized laser, and are given as where σ WH , σ FF , and σ RW correspond to the estimated fractional white, flicker, and random walk frequency noise instabilities, respectively.These instabilities (at τ = 1 s) are then converted into b m as The calculated single clock instability is shown in Fig. 9, where the shaded regions represent uncertainties in the calculation, primarily arising from the estimated frequency noise spectrum of the LO.The scatter points represent single clock instability corresponding to the data in Fig. 4(c) in the main text, divided by √ 2 assuming equal contribution from either atomic servo.The measured instabilities are slightly greater than the numerical calculations using LO frequency noises provided by the manufacturer, which likely because our LO is now in a relatively noisy vibration and thermal environment, and we believe its noise spectrum is now degraded.In order to predict the likelihood of a phase slip and the efficacy of the phase estimation algorithms we employ, we first randomly generate a simulated time series of the LO frequency using its power spectral density (PSD) following Ref.[62], where the LO PSD is sampled with a discretized frequency ∆f .The amplitude A ν (f ) is given by where S ν (f ) is the PSD of frequency noise in units of Hz 2 • Hz −1 and η is a phase randomly drawn from a uniform distribution within [0, 2π] at each f .The amplitude is then converted into a time trace ν(t) via fast Fourier transform

Phase slip probability
To calculate the phase slip probability, we first generate a time trace for up to 1000 cycles with a dead time of 3.5 s and a sampling rate of 10 kHz.We then compute the phase evolution within T through the integral of the instantaneous frequency for each cycle.The distributions of phase evolution over 1000 cycles for various T are shown in Fig. 10.Each distribution is then fitted to a Gaussian function to extract the standard deviation σ.
The phase slip probability P can be computed using the Gaussian error function as where Φ t corresponds to π/2 (π) for standard (quadrature) Ramsey protocol.The computed standard deviations and the associated phase slip probabilities for both the SR and QR protocols are shown in Fig. 11.When bounding the phase slip probability below 10 −4 , we find the maximal interrogation times of 26(3) and 52(3) ms for standard Ramsey and quadrature Ramsey, respectively, consistent with our experimental observations.We note that measurements of clock instability using two interleaved servos in a single clock are blind to both servos simultaneously fringe hopping, e.g., the clock laser occasionally takes sudden frequency steps (as opposed to slower drifts of the LO) that would not be detectable in a self-comparison.In our apparatus, such a large frequency step would primarily come from excessive external vibrations which can rail the active vibration isolation stage holding the clock reference cavity, causing the cavity to experience a large acceleration.We monitor several experiment parameters throughout our data taking to rule this out.First, the AVI stage is monitored throughout the experiment to ensure it stays within its dynamic range.Second, since we are preparing atoms onto desig-  The failure probability for our four-ensemble quadrature Ramsey with phase estimation protocol is given by where P QR,A is the quadrature Ramsey phase slip probability during the shorter interrogation period T A , and P PE is the phase estimation failure probability. .The black curve is a Gaussian fit to the distribution, resulting in a standard deviation of σ = 0.33(1) rad, corresponding to a probability of PPE < 10 −10 using Eq.H3 (Φt = π).This, when accounting for the phase slip probability during TA = 50 ms, results in a failure probability of P QR&PE < 10 −4 , consistent with our experiment observation.
To compute P PE , we perform Monte-Carlo simulations by generating a LO frequency time trace over 1000 cycles.In each cycle, we extract the actual accumulated phase shift θ act,A during T A , which yields a guessed phase shift at T B as θ est,B = θ act,A × (T B /T A ), assuming a constant LO detuning.θ est,B is then compared with the actual phase shift θ est,B during T B to determine the phase estimation deviation as where the phase estimation breaks down when |∆θ PE | > π.
Fig. 12 shows a distribution of the phase estimation deviations computed over 1000 cycles under T A = 50 ms and T B = 85) ms, which is experimentally demonstrated in Fig. 4(c).A Gaussian fit yields in a standard deviation of 0.33(1) rad, corresponding to P PE < 10 −10 calculated using Eq.H3 (with Φ t = π).After accounting for the phase slip probability during T A = 50 ms (P QR,A ≃ 0.8 × 10 −4 ), we find a failure probability of P QR&PE < 1 × 10 −4 .We note that due to the requirement of (2k +1)π/2 phase shifts (k = 0, 1, 2 . . .), the next available T B at T A = 50 ms is 255 ms, which has a phase estimation failure probability of > 10% and is thus not experimentally feasible.To connect our multiple-ensemble protocols to the state-of-the-art optical lattice clocks that make use of higher quality LOs [19][20][21][22][23], we perform Monte Carlo simulations to predict the phase slip probabilities for these systems.In our simulations, we assume LOs with a flicker noise floor of < 4 × 10 −17 (coherence time > 600 ms) and a duty-cycle above 50% [3,6].Similar to Appendix H, we calculate the phase slip probability as a function of interrogation time under a fixed dead time of 600 ms and a flicker noise floor of 4 × 10 −17 .Fig. 13 shows the simulated phase slip probabilities for standard Ramsey (SR) and quadrature Ramsey (QR) approaches, respectively, for state-of-the-art LOs.Bounding the phase slip probabilities below 10 −4 , we predict that a factor of 2 enhancement in coherence times with QR protocol is feasible for state-of-the-art LOs with a duty-cycle > 50%.
We then compute the failure probability four-ensemble quadrature Ramsey and phase estimation (QR & PE) protocol with higher quality LOs, as shown in the black scatter points Fig. 14.We assume a flicker noise floor of 4 × 10 −17 , a fixed dead time of 600 ms for interrogation period T B , and we choose T B = 1.7 × T A .By bounding P QR&PE below 10 −4 , we find a maximal interrogation time T B of around 2.45 s with four-ensemble QR & PE protocol.As the four-ensemble QR & PE protocol doesn't improve over the two-ensemble QR protocol (interrogation time of around 2.7 s) for higher-quality LOs and with > 50% duty-cycles, the simulations confirm the intuition and arguments presented in Section V of the  main text that a third pair of ensembles would be needed to keep track of the phase evolution during T B − T A for a higher duty cycle clock and benefit from multi-ensemble phase estimation with > 2 ensembles.
Appendix J: Potential systematics associated with multiple-ensemble protocols In this section, we consider potential sources of additional systematic shifts arising from the quadrature Ramsey spectroscopy and phase estimation approaches using multiple ensembles, and strategies to mitigate such effects.

Magnetic field gradient
In this work we employ a magnetic field gradient to introduce a π/2 differential phase shift between the ensemble pairs, which is on the order of 20 mG/cm depending on the interrogation time T and magnetic moment of clock transition being used.Under a bias magnetic field of ∼ 5.5 G, the magnitude of first-order Zeeman shift is about 3.1 kHz for the |e, 7/2 → |g, 5/2 transition, and is about 1.87 kHz for the |e, 3/2 → |g, 1/2 transition.While the first-order Zeeman shift is typically cancelled out by averaging transitions between opposite spin states, to support sufficient suppression the magnetic field needs to be stable within ±0.1 mG between experimental cycles for a coherence time of 100 ms, estimated using 10% of the Fourier-limited linewidth.
The addition of a magnetic field gradient also contributes to the second-order Zeeman shift as well as po-tential differential tensor Stark shifts due to spatially varying B-field vectors across the lattice axis.As pointed out in Appendix B, we found a fractional frequency gradient of −1.2(2) × 10 −18 /cm at our operational depth of 15 E rec .In future work, we plan to mitigate this effect by employing a more uniform linear gradient through the addition of Golay coils [61], requiring an upgrade to our experimental apparatus that is currently underway.
For optical lattice clocks with higher quality LOs that enable coherence times approaching second-scale, a smaller magnetic field gradient and a magnetically less sensitivity clock transition can be used to generate the π/2 differential phase shift for quadrature Ramsey spectroscopy.For example, with the 5/2 − 3/2 transition, which has a magnetic field sensitivity of −22.4 Hz/G, only a magnetic field gradient of 3 mG/cm is needed at T = 3 s interrogation time for a pair of ensembles separated by 1 mm [16].In addition, spatially-resolved readout can be carried out with high-resolution imaging to better characterize the gradient, as was demonstrated in Refs.[16,59,63], which can then be used to account for and subtract off the inhomogeneous shifts arising from the gradient.

Systematic shifts from off-resonant clock pulses
The additional clock π/2-pulses employed in the fourensemble phase estimation protocol will lead to additional light shifts on the off-resonant ensembles.Given the measured sensitivity of -13(2) Hz (W cm −2 ) −1 from prior work [64], we estimate a probe light shift of 6×10 −16 caused by the on-resonant clock pulse with a Rabi frequency of 2π × 140 Hz.This sets an upper bound limit below 1 × 10 −15 for the shift arising from the offresonant clock pulses, after accounting for the detuning and Clebsh-Gordan coefficient.While this is a sizable systematic shift for the state-of-the-art clocks, given the large Rabi frequency employed in this work, we note this effect can be reduced by either lengthening the clock pulses, or through the adoption of more sophisticated composite pulse sequences such as hyper-Ramsey spectroscopy [65].

Servo errors from phase decoder
The phase decoder we employ here for quadrature Ramsey spectroscopy assumes a differential phase shift of ∆φ = π/2.Inaccuracy in ∆φ could introduce a frequency offset in the phase decoder (Eq. 4 in the main text).In the limit of δφ ≪ ∆φ, the phase shift error δφ will be converted into a decoder error of δφ/2, and thus result in a frequency offset of δφ/2πT .We note that the phase gradient across the lattice is π/2/cm, which means that the phase difference across an ensemble (1/e 2 radius of 500 µm) is on the order of 100 mrad (between its top and bottom).However, as we per-form a spatial average over each ensemble, to first order the phase gradient across each individual ensemble will contribute to decoherence within each ensemble rather than inaccuracy in the differential phase shift.The dominant source of error will instead be inaccuracy in the differential phase shift between ensemble pairs rather than the phase gradient across each individual ensemble.
In our case, the differential phase error primarily comes from drifts in the magnetic field gradient and any uncertainty in separation between center-of-mass of the ensembles.We note that the magnetic field gradient was characterized with synchronous differential comparison in prior works [27,41], demonstrating long-term differential instability and inaccuracy at 10 −19 level.The ensemble separations are well controlled with the data acquisition card with a timing precision of 100 ns (10 MHz clock reference), corresponding to a separation error less than < 0.1 µm (compared to cm-scale separation).
Assuming a phase error of 10 mrad, this would lead to an error on the order of 1 × 10 −17 at an interrogation time of T = 100 ms, which could be further suppressed well below the 10 −18 level using longer interrogation times with a better LO.In addition, averaging interleaved measurements on opposite hyperfine transitions, as is typically employed to cancel the first order Zeeman shift, will to lowest order also cancel this error.Finally, measurements of δφ can be performed with higher precision if necessary by taking advantage of the ability to probe well beyond the LO coherence time using synchronous differential comparisons [16,27,41].

Cross talk and line pulling
We observe cross talk between the two hyperfine clock transitions used for the four-ensemble phase estimator at the level of roughly 3% (Fig. 3(c)), which could potentially contribute to inaccuracy.We note that line pulling from off-resonant excitation of unintended Zeeman transitions is not unique to our protocol, and the same considerations apply to existing optical lattice clocks [3,4,6,11,66], whether or not other hyperfine clock transitions are employed.We also note that the 3% cross talk in our current demonstration is a technical issue, primarily limited by clock pulse infidelity, rather than a fundamental limit.The maximum off-resonant excitation is given by the ratio of Clebsch-Gordon coefficients squared and the off-resonant Rabi frequency as [4] where γ CG corresponds to the ratio of Clebsch-Gordon coefficient squared, Ω is the Rabi frequency, and ∆ν split is the splitting between nearby transitions.There are multiple ways to mitigate this effect: a) reducing the Rabi frequency of the clock pulse; b) employing a larger bias magnetic field; c) improving the initial state preparation to reduce imperfections in spin-polarization and thus cross talk.With these improvements, we expect it will be feasible to achieve cross talk below 0.1%, and to characterize and control the corresponding systematic uncertainty arising from line-pulling at below the 10 −18 level.

Residual tensor Stark shifts due to multiple mF states
For standard optical lattice clocks operating with transitions with a single m F state, the lattice laser frequency is often detuned such that the scalar component and tensor components null out the lattice Stark shift [67][68][69][70].However, for our four-ensemble phase estimation protocol, which makes use of two hyperfine clock transitions with multiple m F states, such as the |e, 7/2 → |g, 5/2 and |e, 3/2 → |g, 1/2 transitions as used in this work, the lattice stark shift cannot be nulled out for both transitions at a single lattice frequency.Here we estimate the magnitude of this residual shift using the tensor coefficient of -0.0058 (23) mHz/E rec from prior works [71,72], which yields a residual tensor Stark shift of 2.6(1.0)× 10 −17 between the two clock transitions at 15 E rec , where the uncertainty is limited by knowledge of the tensor coefficient.However, we would like to point out that the accuracy of the clock need not be limited by the systematic uncertainty of the transition used to track the phase for shorter interrogation times.The only requirement is that the tensor Stark shift of the transition (|e, 7/2 → |g, 5/2 in our case) that has not been nulled should never be allowed to drift off by enough to result in a phase shift on the order of ±π, and the lattice frequency can be set to null the tensor Stark shift for the clock transition used for the longer interrogation time (|e, 3/2 → |g, 1/2 in our case) and to actually lock the LO.

Running wave contamination due to lattice beam intensity mismatch
To prepare multiple atom ensembles, we employ two double-passing acousto-optical modulators (AOMs) to detune the frequency of the retro-reflected lattice laser beam [27].This results in approximately 50% in an intensity mismatch between the incoming and retro-reflected lattice laser beams due to the diffraction efficiency (90%) of the acousto-optic modulator (AOM), optical path loss and imperfect retro-reflection.The mismatch could lead to a running wave that introduces additional light shifts, which would require precise modeling to constrain at the level of 10 −18 [73].We note that, however, this systematic effect is not new or unique to our system.For instance, McGrew et al. [3] constrains the systematic uncertainty associated with a 15% intensity mismatch to below 1 × 10 −19 .And we expect that this mismatch could be further mitigated with the use of two phasereferenced, counter-propagating lattice laser beams [74] or a bow-tie cavity [14], by doing a handoff of the ensembles from a loading and transportation lattice to a balanced spectroscopy lattice, or by loading multiple ensembles by moving the location of the MOT field zero rather than moving the lattice.

FIG. 1 .
FIG. 1. Quadrature Ramsey spectroscopy with a multi-ensemble clock.(a) Two ensembles of strontium atoms separated by 1 cm are prepared using a movable optical lattice.A magnetic field gradient ∆B(z) is applied to introduce a differential phase shift (∆φ) of π/2 between the two ensembles for a given Ramsey interrogation time T .(b) Theory curves showing the expected Ramsey fringes for the two atom ensembles with a differential phase of π/2 [36].The 4 regions (I, II, III, and IV) correspond to quadratures for decoding the clock laser detuning.(c) Parametric plot of the measured excitation fractions from the two ensembles (P1, P2) for a π/2 phase difference at a short Ramsey free evolution time (T = 7.5 ms) for the full randomly sampled range of detunings corresponding to ∆fi ∈ [−1/2T, 1/2T ], with the light blue curve a fit to the expected ellipse.Each labeled quadrature corresponds to the matching region in Fig. 1(b) (see Appendix C).

FIG. 2 .
FIG. 2. Increasing the interrogation time and reducing clock instability with quadrature Ramsey spectroscopy.(a) Schematic of clock operation using Ramsey spectroscopy with two atom ensembles.The differential phase, ∆φ, is tuned to either 0 or π/2 to compare the clock performance between standard Ramsey spectroscopy and quadrature Ramsey spectroscopy.The clock laser detuning (∆f ) is decoded using the excitation fractions from both clocks, and is then fed back to the LO.(b) Instability as a function of Ramsey free evolution time (T ) for both standard Ramsey (SR) with ∆φ = 0 (red) and quadrature Ramsey (QR) with ∆φ = π/2 (blue).Each point corresponds to the extracted Allan deviation of the self-comparison at 1 s, obtained through linear fit as shown in panel (c).Error bars represent 1σ standard deviation.The shaded region represents a conservative estimate of the region where phase slips do not occur for SR with experiments taken over an hour for at least 2000 runs.Dashed line is a fit to the data based on Dick noise limited instability after accounting for finite pulse length, assuming the LO noise is dominated by flicker frequency noise (see Appendix F).Inset: A representative plot of the measured frequency difference between the two interleaved servos with an interrogation time of T = 40 ms for both SR (red) and QR (blue).A phase slip occurs for SR, leading to a frequency offset of 1/T (dashed line).(c) Lowest measured Allan deviations for clock self-comparisons for SR at T = 20 ms (red) and QR at T = 40 ms (blue).The Allan deviation for SR at T = 40 ms (black), where a phase slip occurred as shown in the inset of (b), is shown for reference.Error bars indicate 1σ standard deviation.The line corresponds to a fit assuming only white frequency noise (1/ √ τ scaling).Due to the finite attack time of the servos we only fit to data after 40 s.Each measurement averages for approximately 1 hour.

FIG. 3 .
FIG.3.Independent control over multiple atomic ensembles with only global addressing.(a) Left: The lattice is loaded with 4 evenly-spaced atomic ensembles with 0.33 cm nearest-neighbor separation, and a linear magnetic field.Right: Relevant energy levels and the pulse sequence used to prepare the ensembles in different hyperfine states.All four ensembles are initially prepared in the stretched hyperfine ground state 1 S0, mF = 9/2 , or |g, 9/2 .A first π-pulse (i) transfers the population into 3 P0, mF = 7/2 (or |e, 7/2 ), followed by a π/2-pulse (ii) on the |e, 7/2 → |g, 5/2 transition that prepares a 50:50 superposition.A second π-pulse (iii) prepares a superposition of two nuclear spins (|e, 3/2 and |e, 7/2 ).During the free evolution time Tp, an applied field gradient (about 45 G/cm) introduces a phase shift of 0, π, 2π, and 3π for each atom ensemble, respectively.Final π-pulse (iv) and π/2-pulse (v) prepare the ensembles into opposite states, with ensembles (1, 3) on |g, 5/2 and ensembles (2, 4) on |e, 3/2 .(b) Global Rabi spectroscopy of the four atom ensembles after state preparation at a bias magnetic field of about 5.5 G. (c) Independent coherent Rabi oscillations for ensembles (1, 3) (orange, red) on the |g, 5/2 → |e, 7/2 transition and for ensembles (2, 4) (blue, light blue) on the |e, 3/2 → |g, 1/2 transition, respectively.Each plot represents the average of 4 independent scans, with each data point indicating the mean.Error bars correspond to 1σ standard deviation.The solid curves correspond to fits to the expected sinusoidal Rabi oscillations, and are used to bound the cross-talk between the two hyperfine clock transitions.

2 and 4 ;
FIG.4.Phase estimation using four atom ensembles with independent control through frequency multiplexing.(a) Basic principle and pulse sequence for clock operation using phase estimation with four ensembles[34][35][36].We interrogate on the |e, 7/2 → |g, 5/2 transition for clock pair(1,3) for a free evolution time of TA, during which the magnetic field gradient is chosen such that the phase difference between ensembles 1 and 3 is π/2.Asynchronously, we interrogate on the |e, 3/2 → |g, 1/2 transition for TB = 1.7 × TA, which is chosen such that the phase difference between ensembles 2 and 4 is also π/2 due to the difference in magnetic sensitivities for the two transitions.An estimator extracts the LO detuning using information from all 4 clocks, which is then used to feed back on the LO.(b) Comparison of inferred phase shifts using various decoders for a long interrogation time.The data is taken at (TA, TB) = (50, 85) ms in a single four-ensemble experiment for 1000 measurement cycles.The inferred phase shifts are post-processed and extracted using various decoders.Top: inferred with standard Ramsey decoder using only ensemble 2; Middle: inferred with quadrature Ramsey decoder using ensembles 2 and 4; Bottom: inferred with quadrature Ramsey with phase estimation decoder using all four ensembles.The dashed lines indicate the [−π/2, π/2] and [−π, π] regions.The numbers above the histograms represent the summed occurrences within each region.The y-axes are shown in log scale for clarity.(c) Achievable interrogation time T and clock instability by probing using two ensembles with ∆φ = 0 (standard Ramsey, SR), two ensembles with ∆φ = π/2 (quadrature Ramsey, QR), and four ensembles with ∆φ = π/2 (quadrature Ramsey with phase estimation as shown in panel (a), QR & PE).Each point corresponds to the extracted Allan deviation of the self-comparison at 1 s, obtained through linear fit.Error bars represent 1σ standard deviation.In the case of QR & PE, the interrogation time T refers to TB.The shaded regions represent conservative bounds on where phase slips do not occur for each case with experiments taken over an hour for at least 2000 runs.The dashed line is a fit to the data based on Dick noise limited instability after accounting for finite pulse length, assuming the LO noise is dominated by flicker frequency noise (see Appendix F).

FIG. 5 .
FIG. 5. Magnetic field gradient (∆B) adjustment by tuning the offset voltage of the Bz bias coils.The B field gradient is measured on an ensemble pair separated by 1 cm via spectroscopy on the 1 S0 → 3 P1 transition.Each data point has a 1σ standard deviation of roughly 0.7 mG/cm.The solid line is a linear fit to the data.The fit residuals are shown in the lower panel.

FIG. 7 .
FIG.7.Point-by-point comparison between known randomly applied phase shifts φ 2nd ∈ [−π, π] (to the phase of the final global π/2-pulse in synchronous Ramsey spectroscopy) and the phase shifts extracted from the decoder using quadrature Ramsey spectroscopy (Eq.4).The red scatter points correspond to the experimental data shown in Fig.1(c) for a short Ramsey free evolution time (T = 7.5 ms).The blue scatter points correspond to Monte-Carlo simulations using the QPN and LO noise spectrum as input parameters and are shown for comparison.

FIG. 8 .
FIG.8.Experimental sequence for the self-clock comparisons.Two independent atomic servos (I) and (II) are interleaved and compared to determine the clock comparison instability.The typical experimental cycle is shown in the bottom (not to scale), which takes about 1.5 s, and thus leads to roughly 3 s dead times in each servo.An additional 250 ms is spent on nuclear spin state initialization when probing four atom ensembles with phase estimation protocol (not shown in the figure).

FIG. 9 .
FIG. 9. Numerical calculations of the expected clock instability using the manufacturer's measurements of the LO instability as inputs into Eq.G5.The shaded region represents upper and lower bound limits from the uncertainty in the LO noise estimates.The experimental data are taken from Fig. 4(c) in the main text, divided by √ 2 to infer the single clock instability.The discrepancy between experimental data and numerical calculations is likely because the noise spectrum of our LO is degraded in a relatively noisy vibration and thermal environment.Error bars correspond to 1σ standard deviation.

FIG. 10 .
FIG. 10.Distribution of LO phase evolution under various interrogation time T .The distribution is computed with a time trace of LO frequency at each T over the course of 1000 measurement cycles.The time trace is generated with a sampling rate of 10 kHz using our experimental parameters, including the LO noise PSD and a dead time of 3.5 s.

3 .
Failure probability of the quadrature Ramsey with phase estimation protocol

FIG. 12 .
FIG.12.Monte-Carlo simulation of phase estimation failure probability.Plot shows the simulated distribution of phase estimation deviations computed over 1000 cycles under TA = 50 ms and TB = 85 ms (the same interrogation times as experimentally demonstrated in Fig.4(c)).The black curve is a Gaussian fit to the distribution, resulting in a standard deviation of σ = 0.33(1) rad, corresponding to a probability of PPE < 10 −10 using Eq.H3 (Φt = π).This, when accounting for the phase slip probability during TA = 50 ms, results in a failure probability of P QR&PE < 10 −4 , consistent with our experiment observation.

TABLE I .
Relevant transitions used in state initialization of four atom ensembles.
Monte Carlo simulations of phase slip probabilities for state-of-the-art LOs.Computed phase slip probabilities PSR (PQR ) under various interrogation times for standard Ramsey (SR, Φt = π/2) and quadrature Ramsey (QR, Φt = π), respectively.A fixed dead time of 600 ms and the flicker noise floor of state-of-the-art LOs (4 × 10 −17 ) are used as input parameters.The shaded area bounds the phase slip probability below 10 −4 .A factor of 2 enhancement in coherence times is feasible with the QR protocol.
FIG. 14. Simulated failure probability for state-of-the-art LOs.We assume a flicker noise floor of 4 × 10 −17 , a fixed dead time of 600 ms for interrogation period TB, and we choose TB = 1.7 × TA.The black scatter points correspond to the failure probability P QR&PE .By bounding P QR&PE below 10 −4 , we find a maximal interrogation time TB of around 2.45 s with four-ensemble QR & PE protocol.