The quantum Allan variance

The instability of an atomic clock is characterized by the Allan variance, a measure widely used to describe the noise of frequency standards. We provide an explicit method to find the ultimate bound on the Allan variance of an atomic clock in the most general scenario where N atoms are prepared in an arbitrarily entangled state and arbitrary measurement and feedback are allowed, including those exploiting coherences between succeeding interrogation steps. While the method is rigorous and general, it becomes numerically challenging for large N and long averaging times.


Introduction
Recent years have seen spectacular improvements in the performance of optical atomic clocks [1,2]. This progress has been enabled by experimental techniques for trapping and manipulating ions [3] and neutral atoms [4] previously developed for quantum information processing and simulation but now fruitfully applied to metrology. In particular, clocks based on a single ion have demonstrated an instability of 7 × 10 −18 within 46 hours of averaging [5], while optical lattice clocks using thousands of trapped neutral atoms can now reach instabilities of 1.6 × 10 −18 in less than 7 hours [6][7][8].
In nearly all ion-based optical clocks, and a number of the latest optical lattice clocks [8,9], the measurement of the accumulated atomic phase is limited by the quantum projection noise of the uncorrelated atoms in the clock. However, unlike in optical interferometry where a phase is itself the measurand, the purpose of an atomic clock is to continuously correct the fluctuating frequency of a classical local oscillator (LO). Quantum projection noise therefore does not, on its own, determine achievable clock stability, which depends on a more complex interplay between phase noise, decoherence, and experimentally adjustable parameters such as the duration of the interrogation pulses applied to the atoms. Thus, while the use of quantum entanglement to enhance clock stability has been studied theoretically for a quarter-century [10][11][12][13][14][15][16][17] and demonstrated in proof-of-concept experiments [18][19][20][21][22], it is not yet clear whether it offers significant benefits in the presence of realistic LO noise and without the artificial constraints on interrogation time imposed in demonstration experiments.
By comparing the performance of an existing clock or known measurement protocol to a bound valid for all possible clocks and protocols, including those using quantum entanglement or adaptive measurement strategies, while properly accounting for the effects of realistic noise, one can determine whether the development and implementation of more sophisticated clock protocols is a worthwhile investment of time and resources. We have already seen a fruitful instance of this decision-making strategy in gravitationalwave detection experiments using quantum-enhanced optical interferometry where it was found that the basic quantum-enhanced interferometer, fed with coherent light and a squeezed vacuum state at the two input ports, saturates the fundamental bounds for the precision of lossy interferometry [23][24][25][26] in the parameter regime relevant for gravitational-wave detection [27]. This finding allows experimentalists to focus on improving technical parameters of their setup, such as increasing laser power and reducing loss, rather than wasting resources to prepare more sophisticated non-classical quantum states of light with no practical benefit. This paper aims to provide analogous decision-making tools for atomic clock operation. We present a general formula for computing, from the LO noise correlations and decoherence properties of a particular system, a general lower bound for the Allan variance, the experimentally accessible measure of clock instability [28]. General scheme for stabilization of a local oscillator (LO) to an atomic reference frequency ω 0 . For the ith interrogation, the atoms are prepared in an initial state ρ 0 , interact with the LO whose (stabilized) frequency is ω (i) (t) and evolve for a time T into ρ (i) . A measurement Π x i is applied which may depend on the outcomes of previous measurements , based on all available data x i , serves to steer the LO frequency ω (i+1) (t) during the next interrogation nearer to ω 0 . This cycle repeats k times within an averaging interval τ . The gray sine link at the bottom indicates that the atomic states used in different interrogation cycles may be entangled, so that the analysis applies to scenarios in which quantum coherences are preserved from one interrogation cycle to the next.

Formulation of the problem
Let ω(t) be the time-dependent frequency of the LO (laser or microwave generator) which is to be stabilized to the reference atomic frequency ω 0 , as in Fig. 1. Each step of the protocol lasts for a time T , during which atoms prepared in a state ρ 0 interact with radiation from the LO. At the ith step, the result x i of a measurement on the final atomic state ρ (i) is combined with all previous measurement results in a record whereω(x i ) denotes the correction and ω LO (t) is the frequency of the free-running, uncorrected LO. The figure of merit quantifying the instability of the clock is the fractional Allan variance (AVAR) [28]: where · represents averaging over frequency fluctuations of the LO as well as the measurement record distribution and ω(t) is the corrected frequency of the clock: The AVAR quantifies the expected discrepancy between two consecutive observations of the clock frequency, each averaged over a time τ . We assume that the averaging time τ is a multiple of the interrogation time τ = kT . Given a spectrum of LO frequency fluctuations, we want to minimize σ 2 (τ ) for a given τ and input state ρ 0 over all measurement and feedback strategies. We call the resulting minimum the quantum Allan variance (QAVAR). This is a horrendous task because the 2τ time period considered when calculating σ 2 (τ ) includes a potentially large number 2k − 1 of coupled measurement and feedback steps to optimize simultaneously. The problem is therefore much more demanding than minimizing estimation variance in the single-shot protocols considered in optical interferometry, and is conceptually much closer to waveform estimation [29].
This problem has been attacked before [14][15][16][30][31][32][33], but to the best of our knowledge none of these attacks have yielded a rigorous and general bound. In [14,33] the AVAR figure of merit was replaced with an instantaneous variance which diverges for several noise processes commonly seen in experiments. Other works such as [31] assumed to operate in the regime where correlations between ω(t) in different interrogation steps can be neglected. Some analyses consider only specific protocols [15,16], or present a numerical approach that may help in improving practical interrogation schemes [30,32] but provides no closed formula for the fundamental bound on the Allan variance.
We consider an atomic system of N two-level atoms whose evolution, during the i-th step reads: where Λ T is a general quantum map [34] describing atomic decoherence processes. Using shorthand notation U[ρ] = UρU † , the unitary operator describes the interrogation dynamics generated by H = N n=0 nP n , where P n is a projection on a subspace containing atomic states with n excitations. Note that we go to a frame rotating at the LO frequency. For clarity we will hereafter replace [ω (i) (t) − ω 0 ] in the above unitary with ω (i) (t), measuring frequencies as detunings from atomic resonance, since a constant frequency shift does not affect the AVAR. We choose a flat window function so the LO detuning from atomic resonance is sensed evenly over the total interrogation time T , as in ideal Ramsey interrogation with two π/2 pulses of negligible duration. A more general window function that weights the LO detuning unevenly, though possible, generically degrades the stability of the clock through the Dick effect [35]. We choose to ignore the Dick effect as the main goal of the present paper is to provide fundamental stability bounds set by the intrinsic LO noise and atomic decoherence, rather than to model in detail the effect of current technical limitations. The ideal Dick-effect-free performance we analyze can be approached by using interleaved interrogation techniques to suppress dead time, as recently demonstrated for microwave atomic clocks [36] and proposed for optical lattice clocks [8].
After the i-th interrogation a measurement is performed, described by a POVM This takes into account the possibility of performing generalized measurements, resulting from non-trivial interaction with ancillary systems, which are not necessarily described by projection operators. The measurement outcomes {x i } occur with probabilities p(x i ) = Tr(ρ (i) Π x i ). We allow the choice of measurement operators to depend on all previous measurement results x j (j < i), so that the analysis encompasses adaptive measurement protocols. The joint probability distribution of the first i measurement results reads: As the feedback correctionω(x j ) depends on earlier measurement results (see Eq. (1)), the density matrices ρ (j) implicitly depend on x j−1 .

Main result
Given a stochastic process ω LO (t) representing frequency fluctuations of the free running LO, an initial atomic state ρ 0 and an interrogation time T , the AVAR σ 2 (τ ) achievable in the most general stabilization scheme is lower bounded by the QAVAR σ 2 Q (τ ): where σ 2 LO (τ ) is the AVAR of the free running LO, is the tensor product of atomic states evolved at each of the interrogation steps according to the uncorrected LO and averaged with respect to the process ω LO . The operator L encapsulates information about possible feedback and measurement strategies and is implicitly defined by in which A derivation of the bound is given in the Appendix A. A more explicit form of the QAVAR reads where |λ r and λ r are eigenvectors and eigenvalues of ρ. The second term in Eqs. (6) and (11) can be understood as a bound on the strongest correlations that can be engineered Figure 2. Fractional Allan variance for atomic clock models whose LO noise (black, dotted) approximates that of Ref. [39]. Bounds for the AVAR are derived for clocks using a single atom (N = 1, left) or two atoms (N = 2, right). In the absence of coherences between subsequent interrogations, the calculated fundamental bound (black, solid) reaches the asymptotic 1/τ scaling regime, where it can be extrapolated to long times (black dashed). In the N = 2 case the modest gains obtainable by allowing entangled states of the two atoms lie within the thick black line. Neglecting noise time-correlations while calculating AVAR leads to visibly looser bounds (gray, dotted) proving that white noise approximations are not justified in this setup. Bounds for the most general scheme, allowing entanglement between the probe states used at different interrogation steps (gray, solid). Simulations of a simple atomic clock with the same LO noise model (open circles) obey the bounds, and seem to approach them asymptotically at long times.
between LO frequency fluctuations and the corrections applied in the two averaging intervals, given the information available from the atoms. The QAVAR is a lower bound on the achievable AVAR, much as the Quantum Fisher Information (QFI) is an upper bound on the classical Fisher information (FI) that imposes a lower bound on estimation variance via the Cramér-Rao inequality [37,38]. Note, however, that our formulation is intrinsically Bayesian due to the averaging with respect to the uncorrected ω LO process, which plays the role of a prior. In the case of the QFI it is known that there always exists a measurement for which the FI saturates the QFI. By inspecting the derivation of the QAVAR one finds that the bound is tight provided that one allows for experimentally implausible collective measurements on probe states at different interrogation steps. Whether such collective measurements are indeed more powerful than adaptive strategies using separate measurements at each step is a question which we leave open. In any case, given a particular measurement and feedback strategy that approaches the QAVAR, one is sure to have reached the ultimate optimum.
The QAVAR allows us to set aside the issue of choosing the optimal measurement and feedback and focus on the choice of optimal probe states and interrogation time. One can carry out brute-force minimization of the QAVAR over input probe states ρ 0 , but a much more effective approach is developed in [33,40]. Starting with a random probe state one obtains the corresponding optimal measurement/feedback strategy described by the operator L. Then, given L, one looks for the state performing optimally under this particular strategy, iterating the procedure until the results converge. For the example presented below we will assume that the initial probe states ρ 0 , while they may involve entanglement among the atoms at a given interrogation step, are uncorrelated from interrogation to interrogation, so that the input is of the form ρ ⊗2k−1 0 . However, we may relax this constraint and assume that the joint state of atomic probes is arbitrary and may even involve entanglement between different interrogation steps. Though most such schemes are hardly realistic, they cover as a special case a recent experimentally feasible proposal where coherence is partially preserved through succeeding interrogation steps [41]. Such an approach, with no restrictions on independence of the atomic probes, provides the most fundamental bounds and sets a benchmark to which realistic protocols may be compared.
These fundamental bounds also apply to protocols where different subensembles of atoms are interrogated for different times [42,43]. Although we have presented our model in terms of a fixed cycle time T , each cycle in the averaging time is included in our model as a separate subsystem and all these subsystems are treated formally as evolving in parallel before being measured in a collective measurement. By allowing entanglement and joint measurements between different subsystems (time-steps), we can model interrogations that last several cycles. For instance, consider the singleatom case and take two successive interrogation steps. If instead of a product state of two subsystems corresponding to two independent interrogation steps we consider an entangled state (|00 + |11 )/ √ 2 and let it evolve over time T , its evolution will be isomorphic to the evolution of the single atom state (|0 + |1 )/ √ 2 evolving over interrogation time 2T . In this way, allowing entanglement between time-steps lets us cover all schemes where effective interrogation times are multiples of some minimum cycle time T .

Example
Assume the LO frequency fluctuations of our free-running LO are modeled by a combination of the Ornstein-Uhlenbeck process with white Gaussian frequency noise so that the autocorrelation function reads: We choose the noise parameters to be α = 2 (rad/s) 2 , β = 0.4 (rad/s) 2 s, γ = 0.5 s −1 , so that the uncorrected AVAR closely resembles that of the laser system used in [39], where ω 0 ≈ 3.25 × 10 15 rad/s, for the important region from 0.1-10 s. For simplicity we neglect any additional atomic decoherence effects, which are negligible compared to LO noise in state-of-the-art atomic clocks. In order to calculate the QAVAR we must first find the operators ρ and ρ ′ defined in Eqs. (7) and (10). Without loss of optimality we may assume that the input state ρ 0 is supported on the fully symmetric subspace of N two-level atoms [14]. We will denote by |n the symmetric state of N atoms with n excitations. Let |n = |n 1 ⊗ . . . ⊗ |n K be a state where n i atoms are excited at the i-th interrogation step. Thanks to the Gaussianity of the process we may use the cumulant expansion [44] to obtain an explicit formula for ρ and ρ ′ in the |n basis (see the Appendix B for details): (14) where ⌈x⌉ is the ceiling function. With explicit formulas for ρ and ρ ′ , the QAVAR can be calculated using (11). Figure 2 depicts the QAVAR as a function of the averaging time for different probe states. For each averaging time τ we look for the optimal interrogation time T = τ /k and choose the sub-multiple k that minimizes the QAVAR. This limits us to relatively short averaging times: for increasing τ the optimal T stabilizes and k increases linearly with τ , making the problem intractable due to the tensor product structure of ρ and ρ ′ . Still, we expect the QAVAR to approach σ 2 Q (τ ) ≈ c/(ω 2 0 τ ) for some constant c when τ is much larger than characteristic noise correlation time scales. Where we observe such behavior, we may extrapolate our results and estimate the best long-term clock stability obtainable with a given atomic reference and LO noise spectrum. From the numerically optimized bounds for the current model we extrapolate rough long-term stability bounds of this form, obtaining c ≈ 1.33 rad 2 s for N = 1 whereas for N = 2 we find 0.78 rad 2 s in the case of product states and 0.73 rad 2 s when allowing optimally entangled states (which were almost identical to the (|00 + |11 )/ √ 2 Bell state). This shows that in this case the use of entangled states provides little benefit.
At larger atom numbers, for which the computational cost of our algorithm currently exceeds our resources, the use of entanglement may produce greater gains. In particular, if no additional atomic decoherence effects are taken into account, other approaches indicate the possibility of Heisenberg scaling of Allan variance with increasing N [16]. Our fundamental bounds, when calculated for the corresponding noise model, must allow such scaling to be consistent with these earlier results. That said, one must keep in mind that from some N onwards even tiny levels of atomic decoherence will play a dominant role, and it follows from general quantum metrological considerations that the scaling must revert to the standard scaling of precision in the asymptotic regime of large N with at most a constant factor improvement from entanglement [25,26].
We have also plotted the QAVAR corresponding to the most general scenario where the probe states used at different interrogation steps are optimally entangled. In this case we were not able to reach the 1/τ regime needed to extrapolate to long-term stability.
We also plot the Allan variance observed in 5000 s of simulated clock operation. The modular simulator we have developed can generate LO noise traces for a range of noise processes including white, flicker, random-walk and the Ornstein-Uhlenbeck process used here, simulate the response of arbitrary references to this LO noise, and derive a frequency trace corresponding to the observed output of the frequency standard when it is stabilized to the reference according to a given feedback policy. For the simulated clock data in Fig. 2, the atomic reference used simple Ramsey interrogation of 1 or 2 atoms in a product state, with a fixed probe time T = 0.5 s. The feedback was derived using the integrator servo algorithm commonly used in optical frequency standards [45]. For N = 2, the simulated clock's stability is little worse than the bound for long averaging times, so we may conclude that a 2-ion clock operating with an LO similar to the one considered here can nearly achieve the best possible long-term stability using conventional Ramsey interrogation, and more exotic measurement strategies could provide only a minor additional improvement.
One further explanation is due in connection with [31] where Allan variance plots, computed for clocks using a hierarchy of atomic subensembles with ever-longer interrogation times, display an initial 1/τ 2 scaling only to revert to 1/τ scaling once atomic decoherence becomes significant. Even in the absence of atomic decoherence the scaling must, for any finite N, eventually revert to 1/τ when there are no longer sufficient atoms to keep adding levels to the hierarchy and thus coherently monitor the LO phase for ever-longer times. In [31] it was assumed that there were sufficient atoms available to track the laser phase up to the atomic lifetime. In our small-N example we encounter the atom-number limit too soon to observe the initial quadratic scaling.

Numerical complexity of the method
The examples provided above were limited to small values of N and relatively short averaging times τ . This was due to a rapid increase in computational cost that prevented us from going beyond the regime of ∼ 1, 2 atoms and averaging times of a few seconds. For the moment, our results are thus directly relevant to atomic clock implementations using few ions and where noise memory effects are weak enough that a few interrogation steps suffice to reach the regime where the Allan variance scales as 1/τ . The computational complexity of obtaining results for larger atomic ensembles and longer noise correlation times can be expressed as a function of two parameters, the number of atoms N and the number of interrogation steps k within the Allan variance averaging time τ .
Although the Hilbert space describing N two-level atoms is 2 N -dimensional in general, for phase or frequency estimation the optimal input probe states are generally supported on the N + 1-dimensional fully symmetric subspace [14,46,47]. Moreover, if only LO noise is considered then the output state will be confined to this space as well, and hence the number of entires in the density matrix of the output state scales as N 2 . The same holds if, in addtition to LO noise, one allows for decoherence processes that act globally on the atoms, such as magnetic field fluctuations that dephase all the atoms collectively. Even for decoherence processes that violate this condition it may sometimes still be possible to describe the output density matrix in a much more efficient way than using the full 2 N -dimensional space. For example, if individual losses are considered, the output state can be written as a direct sum of states supported on fully symmetric subspaces corresponding to the different possible numbers of surviving atoms, so that the space required to describe the output state scales as N 3 . Similarly in the case of individual dephasing processes, permutation symmetry of the output density matrix (not to be confused with permutation symmetry of vectors in the Hilbert space) admits a direct product structure where the number of parameters scales as N 3 [48]. Therefore, the scaling of Hilbert space dimension with the size N of the atomic system, when considering a single interrogation step, is reasonable and without much affort it is possible to analyze problems involving N of the order of 10 2 atoms. A direct analysis of optical lattice clocks, which typically use ∼ 10 3 atoms, could be undertaken using techniques that take advantage of the fact that atoms in such setups are independent or only weakly correlated, such as Gaussian state approximation methods.
The other parameter that dominates the numerical complexity is the averaging time τ . For averaging times much larger than the noise correlation time, we observe as expected that the optimal interrogation time tends to some value that does not depend significantly on the averaging time. This implies that when extending τ the number k of optimally tailored interrogations steps within the averaging period will scale linearly with τ . Since the formalism involves a tensor product of Hilbert spaces corresponding to different interrogation steps the result will be an exponential scaling of the full Hilbert space with k and hence with τ as well. Summing up, the Hilbert space that needs to be considered scales as N 2k , and from this it is clear that the dominant source of numerical inefficiency of the algorithm is the exponential scaling with k.
Physical intuition suggest, however, that when noise correlations die out on certain characteristic time scales, it might be possible to model the problem in a more efficient way by making use of the fact that correlations between subspaces corresponding to different interrogation steps will effectively be of finite range. This immediately suggests the use of well-developed methods from quantum information and quantum manybody physics used to describe systems with finite-range correlations, such as matrix product states or operators [49,50] or quantum Markov chains [51]. Here we outline the steps needed to express our algorithm in matrix product operator formalism. The implementation itself and its application to realistic models will be the subject of a separate publication.
(i) The first step is to choose a bond dimension for matrix product operators, and write appropriate approximations of states ρ, ρ ′ (7,10). One needs to investigate how large the bond dimension must be in order to faithfully represent the state.
When noise correlations die out on a finite time scale, a rough first choice is to take the bond dimension to be of the order of the number of interrogation steps that span the time scale over which noise correlations die out.
(ii) Having written the relevant states in the matrix product operator form one plugs them into equation (9) and finds the operator L, also assumed to be in matrix product operator form. This can be done by performing numerical minimization of a norm of the expression ρ ′ − 1 2 (ρL + Lρ) . One can follow here a standard procedure known from matrix product operator calculations, and perform minimization sequentially over a matrix corresponding to a subsystem representing a given interrogation time and proceed sequentially on different subsystems back and forth until a satisfactory convergence is achieved.
(iii) The obtained operator L is plugged into (6), which is easy as the trace of a product of matrix product operators is directly calculable, to obtain the final bound.
We hope that an efficient implementation of this method will allow the treatment of averaging times long enough to provide rigorous bounds for all typically encountered models in present day atomic clock setups, at least in the regime of small atom number, thanks to the fact that the increase in computational cost due to an increase of the averaging time will now only be linear instead of exponential. It is also important to reiterate that the averaging time needs to be increased only up to the point where one clearly observes the Allan variance taking on the 1/τ scaling characteristic of white noise processes. As observed in the examples of Sec. 4 this occurs for relatively short times in the models we have considered. The prospects for analyzing practical experimental systems are therefore better than one might naïvely expect. In order to model experimentally realistic averaging times of hours or days, it is probably sufficient to compute the QAVAR bound out to the regime of 10 − 100 second averaging times, beyond which the results can be reliably extrapolated. We indeed expect that the matrix product operator approach can enhance the numerical efficiency by the necessary one to two orders of magnitude.

Conclusions
We have provided a rigorous derivation of the fundamental bound on the AVAR in atomic clock operation, taking into account all measurement and feedback strategies allowed by quantum theory. We have provided explicit numerical bounds for atomic clocks using 1 or 2 atoms for a given realistic LO noise model and have shown that, in these scenarios, conventional Ramsey interrogation schemes perform almost as well as these bounds allow in the asymptotic regime of long averaging times. We have also proposed an approach to deriving systematic approximations to the fundamental bound using matrix product operators, which will allow the extension of the algorithm to more atoms and longer averaging times. The computation of such bounds for realistic models will be a valuable guide in choosing which quantum metrological protocols or technical enhancements to pursue in order to obtain real performance enhancements in practical frequency standards. Let us assume that τ = kT , so that there will be K = 2k − 1 measurement/feedback steps within the time period (0, 2τ ) considered when calculating the AVAR. Let Π x i , ω(x i ) describe a general measurement and feedback strategy at the i-th step of the protocol, where x i = (x 1 , . . . , x i ). By definition (2) the AVAR reads: and ⌊x⌋ is the floor function. Λ T describes all the decoherence effects that affect the atoms during the interrogation whereas H is responsible for the unitary sensing part. Their actual form is irrelevant for the proof. We rewrite (A.2) as: where all feedback corrections have been gathered intõ This leads to: Note that in the first term we performed a trivial integral over dx K p(x K ) and what we get is the uncorrected AVAR σ 2 LO . In order to deal with the other two terms we define: so that p(x K ) can be written as: The advantage of this form is that all the dependence on x i is now in Π x i whereas ρ (i) LO represents evolution of the state assuming no corrections to a free-running LO had been applied. As a result we arrive at: Minimization of σ 2 (τ ) in (A.9) over measurement/feedback procedures, which amounts to optimization over L 1 , L 2 operators, has the same formal structure as a standard quantum Bayesian estimation problem with a quadratic cost function [33,37]. Using the same arguments as in the standard Bayesian problem one first proves that standard projective measurements, where Π x i Π x ′ i = δ x i ,x ′ i Π x i , are sufficient to reach the optimal performance and hence: We should now minimize the above formula over L. This is a hard task if we want to keep the constraints imposed on L by the structure seen in (A.11). However, we can obtain a lower bound on σ 2 (τ ) by unconstrained minimization of (A.12) over arbitrary hermitian operators L, just as in the solution to the standard quantum Bayesian problem [33,37]. Physically, this approach is equivalent to relaxing the product structure of measurement operators appearing in (A.11) and allowing for arbitrary collective measurements performed on all subsystems representing different interrogation steps Analogously, using: we can obtain an explicit formula for ρ ′ : (B.6)