Correlation spectroscopy with multi-qubit-enhanced phase estimation

Ramsey interferometry is a widely used tool for precisely measuring transition frequencies between two energy levels of a quantum system, with applications in time-keeping, precision spectroscopy, quantum optics, and quantum information. Often, the coherence time of the quantum system surpasses the one of the oscillator probing the system, thereby limiting the interrogation time and associated spectral resolution. Correlation spectroscopy overcomes this limitation by probing two quantum systems with the same noisy oscillator for a measurement of their transition frequency difference; this technique has enabled very precise comparisons of atomic clocks. Here, we extend correlation spectroscopy to the case of multiple quantum systems undergoing strong correlated dephasing. We model Ramsey correlation spectroscopy with $N$ particles as a multi-parameter phase estimation problem and demonstrate that multiparticle quantum correlations can assist in reducing the measurement uncertainties even in the absence of entanglement. We derive precision limits and optimal sensing techniques for this problem and compare the performance of probe states and measurement with and without entanglement. Using one- and two-dimensional ion Coulomb crystals with up to 91 qubits, we experimentally demonstrate the advantage of measuring multi-particle quantum correlations for reducing phase uncertainties, and apply correlation spectroscopy to measure ion-ion distances, transition frequency shifts, laser-ion detunings, and path-length fluctuations. Our method can be straightforwardly implemented in experimental setups with globally-coherent qubit control and qubit-resolved single-shot read-out and is thus applicable to other physical systems such as neutral atoms in tweezer arrays.


I. INTRODUCTION
The ability to estimate the phase of a wave is key to practical applications such as time keeping with atomic clocks [1], rotation and acceleration sensing [2], gravimetry [3], but also to probing fundamental physics [4] and measuring fundamental constants of nature [5].Sensing techniques such as optical or matter wave interferometry rely on phase comparisons of two light waves or matter waves, respectively.In optical atomic clocks for instance, the phase of an atomic superposition state is compared to the phase of the laser having created the superposition.In most of these applications, a large number of uncorrelated photons or atoms are probed, giving rise to a measurement uncertainty governed by the standard quantum limit according to which the uncertainty decreases inversely with the square root of the number of particles being probed.If, however, quantum correlations exists between the particles, the scenario becomes much more interesting and complex.
In this context, phase estimation based on quantum measurements constitutes a subfield of quantum metrology, which aims at making sensitive measurements of physical quantities by harnessing quantum resources, in particular entanglement [6].To this end, it has been shown that entangled input states can be used to beat the standard quantum limit [7,8] and that entanglement can be a resource for achieving optimal phase sensing over a wider range of phases [9].However, as entangled states easily decohere under environmental noise, the performance gain of entanglement-enhanced metrology protocols can be jeopardized by decoherence processes [10,11]; the achievable precision bounds depend on whether the noise is Markovian or contains temporal or spatial correlations [12,13].
Furthermore, from a practical point of view, entanglement-generating resources are often not readily available in precision experiments.For this reason, it is of interest to consider quantum metrology protocols using quantum correlations other than entanglement that might be easier to implement for carrying out quantumenhanced measurements [14].In this paper, we will focus on correlation spectroscopy [15,16], a phase estimation technique for probing the phase difference of qubits subjected to spatially correlated noise, and extend it to networks of N qubits.In the following, we provide our motivation for studying this measurement scenario.
Coherent probing of ultra-narrow atomic transitions in combination with outstanding characterization of systematic level shifts has led to the development of optical atomic clocks with unprecedented precision [17].To verify a clock's performance, its frequency has to be compared with another clock.The uncertainty with which the frequency difference of the clocks can be determined within a given measurement time is usually not limited by the lifetime of the atomic energy levels but rather by the local oscillator's phase noise that sets an upper bound to the useful probe time.This limitation can be overcome by synchronous probing of the two clocks with the same local oscillator and correlating the measurement outcomes.In the case of ensemble-averaged signals, such as in optical lattice clocks where the excitation probability of a large number of atoms is measured, [18,19], the correlations are purely classical.If, however, the measurements probe the quantum state of individual atoms, the correlations can become non-classical, even in the absence of any entanglement [20].
It is in this context that correlation spectroscopy [15,16] has been developed, a technique for measuring transition frequency differences in the presence of correlated phase noise with probe times that can be significantly longer than the coherence time of each system with respect to the local oscillator [21][22][23][24][25].It is based on a synchronous standard Ramsey-type interrogation of two or more atoms by the same oscillator: a first π/2 pulse rotates the Bloch vector into the equatorial plane where it precesses during the free evolution time with a rate set by the detuning of the oscillator from the atomic transition.The second π/2 pulse in conjunction with a state detection in the energy eigenbasis enables the measurement of a spin projection in the equatorial plane.However, instead of measuring expectation values of individual atoms, a parity measurement is used to correlate the measurement outcomes of pairs of atoms.By this approach, transition frequency differences can be measured by observing parity oscillations as a function of the duration of the free evolution time.While correlation spectroscopy only achieves a maximum parity oscillation contrast of 0.5 and therefore does not achieve the optimum signal-to-noise ratio obtainable by preparing maximally entangled states of the two systems [26][27][28], it is technically much easier to implement.
The detection of phase shifts in the presence of strong correlated phase noise is a common scenario that appears in a wide variety of sensing platforms.Apart from the example of multiple clocks probed by the same oscillator, spatially correlated noise can result from the spatial proximity of the qubits [29][30][31], instabilities of the local oscillator probing them [29,32], or from the coupling of the qubits to a common bosonic mode [33,34].A similar scenario appears also in interferometers and optomechanical sensors where a displacement noise of the mirrors and radiation pressure induce a correlated noise on the different output modes [35][36][37].
In this work, we investigate the parameter phase estimation scenario as sketched in Fig. 1: we consider a set of N qubits all of which were prepared in states with Bloch vectors in the equatorial plane and subjected to correlated phase noise (panel (a)).We want to estimate the angle between the Bloch vectors of a pair of qubits (i, j) by applying a second π/2 pulse, measuring all qubits in (b) Correlations Ci,j between pairs of qubits (i, j) are measured for an estimation of the angle between the respective Bloch vectors.Is it possible to reduce the measurement uncertainty of Ci,j obtained from a finite number of experimental repetitions by taking into account all pair correlations that can be simultaneously measured, or even all N -qubit correlations?
the energy eigenbasis and correlating the measurement outcomes (±1) to obtain the correlation C i,j (panel (b)).
We ask the question whether the measurement uncertainty of one of the correlations, e. g.C 1,2 , obtained from a finite number of experimental repetitions could be reduced by taking into consideration all other pair correlations that were simultaneously recorded instead of analyzing only the measurements of the particular pair, e. g. (1,2).We analyze this simple model and show that this is indeed the case.Moreover, we prove that the uncertainty can be even further reduced by considering the N -qubit correlations between all particles.In this way, we generalize the notion of correlation spectroscopy to a quantum sensor network of N two-level quantum systems [23,25] and demonstrate, in theory and experiment, an improvement compared to the traditional pair-correlations method.We derive precision bounds when all the N 2 pair correlations of the outcomes are used and for the case where all N -qubit correlations are exploited.These methods are implemented in experiments with ion crystals and are used to estimate ion-ion distances and transition frequency shifts.Finally, we propose theoretical schemes for further improvement with entangled measurements and initial states.
The manuscript is structured as follows: in section II we describe the principle of N -qubit correlation spectroscopy and how the analysis of measured correlations can be used for inferring relative phase shifts between the qubits as well as tracking correlated phase shifts on all qubits in the time domain.Section III demonstrates the implementation of the measurement protocol in one-and two-dimensional ion crystals with up to 91 ions.In sec.IV, we discuss lower bounds to the measurement uncertainties when analyzing pair correlations or N -qubit correlations and demonstrate that these bounds are nearly saturated in our experiments.We further dis- 2. (a) Measurement protocol: Ramsey experiments are simultaneously carried out on an ensemble of N qubits subject to correlated dephasing, phase-shifting all qubits by a random phase ϕ, and single-qubit phase shifts φn.The analysis of correlations between measurement outcomes on different qubits taken at the same time (column vector qm) enables the estimation of phase difference between qubits; similarly, the analysis of correlations between measurement outcomes taken at different times on the same qubit (row vector qn) provides information about the temporal evolution of phases (for details see main text).(b) We implemented correlation spectroscopy with ensembles of trapped and laser-cooled ions, such as the two-dimensional 91-ion crystal held in a monolithic ion trap shown in the picture.
cuss general quantum precision limits and show that the input states used in our experiments are near-optimal in terms of the achievable measurement precision.Section V discusses applications of the method in trappedion experiments for the determination of transition frequency differences, ion-ion distances and tracking of local oscillator noise.In sec.VI, we discuss the application of our measurement protocol to other experimental platforms.

II. MODEL: N-QUBIT CORRELATION SPECTROSCOPY
We consider a data set consisting of m = 1, . . ., M realizations of Ramsey experiments with a free evolution time T , each of which is simultaneously carried out on an ensemble of N qubits (see Fig. 2a).Prior to the second π/2 pulse, the state of the N qubits during the m-th realization is with phases where φ i is a qubit-specific phase and ϕ m a random phase that is common to all qubits, i. e. in experiments φ i appears as a spatially varying phase whereas ϕ m encodes temporal changes.To achieve an unambiguous definition of these phases, we define ϕ m to be the phase change between the mth experiment and the first one.The phase ϕ m could result from a stochastic process coupling the qubits to an environment inducing correlated dephasing; alternatively, it could be engineered in the experiment, for example by randomly shifting the phase of the first Ramsey pulse with respect to the second one.The qubitspecific phases φ i = (k 1 − k 2 )r i + ∆ i T arise if the qubits have different detunings ∆ i with respect to the local oscillator frequency or if the qubits are excited from different spatial directions for the first and the second π/2 pulse.Here, k 1 (k 2 ) is the k-vector of the running wave inducing the first (second) π/2 pulse and r i is the qubit position vector.We assign an outcome q im = 1 or −1 to the measurement, depending on whether we observe qubit i in the mth measurement in the state |0 or |1 .The probability of observing the outcome q im is given by p(q im ) = 1 2 (1 + q im sin φ im ).Here, we will study two closely related problems: 1.We want to carry out a multi-parameter estimation of the qubit-dependent phases φ i in experiments where the random phases ϕ m are uniformly distributed over the interval [0, 2π); this situation can arise if, for example, the probe time is much longer than the coherence time of the qubits.We can model this problem by preparing the qubits in a state, which contains no entanglement but quantum correlations in the form of non-zero quantum discord [20,38].Next, the qubits are subjected to the unitary operation U φ = exp( i 2 i φ i σ z i ) followed by a global π/2 pulse around the x-axis, U X = exp(−i π 4 i σ x i ), and finally a projective measurement of all qubits is carried out in the computational basis.Given a set of measurement outcomes stored in the matrix Q = (q im ), the goal is to devise a strategy for estimating all phase differences φ i − φ j with optimal precision.Note that this is a special case of a quantum sensor network [39][40][41][42], where the linear functions we wish to estimate are all the phase differences.2. We are interested in characterizing the stochastic process that gives rise to temporally fluctuating random phases ϕ m .Because of the symmetry of the problem in space and time as showing up in eq. ( 2), a strategy for estimating the single-qubit phases φ i can equally well be applied to an estimation of ϕ m by analyzing the transposed matrix Q t of measurement results.
Let us first understand the fundamental precision limits in estimating the phase differences.Given a pure product state, 2 − N 2 N i=1 (|0 i + e iφi |1 i ), and in the absence of noise, the precision in estimating each phase independently from M measurements is σ φi = 1 √ M .Hence the minimal uncertainty in estimating a phase difference ∆φ = φ 2 − φ 1 is This approach basically amounts to inferring ∆φ from the relative phase shifts of two Ramsey fringes.Since this is the best achievable precision with a product state, we refer to it hereafter as the noiseless precision bound.
Considering N = 2 qubits and correlated dephasing, as in Eq. ( 3), the phase difference ∆φ is estimated using standard correlation spectroscopy.Using error propagation of quantum projection noise, the uncertainty in the estimation of ∆φ from M measurements equals The uncertainty diverges when the phase difference approaches 0 or π, i. e. the points where the parity reaches an extremum.It becomes minimal for ∆φ = π/2, where min the noiseless precision bound.The √ 2 difference in the uncertainty stems from the reduced contrast (< 0.5) in correlation spectroscopy.
For N > 2, one can ask whether the uncertainty of the phase difference estimation can be lowered by employing a more sophisticated analysis.Here, we provide an affirmative answer: we show that the uncertainty can be reduced by estimating the phase differences using all the N 2 pair correlations of measurement outcomes, and that in addition a further reduction is achieved by using all the multi-particle correlations.The intuition behind this improvement is based on the following argument: An estimate of the single-qubit phase differences ∆φ ij = φ i − φ j from the observed correlations makes it possible to estimate the random phases ϕ m of each experimental realization.In the limit of a large number of qubits, the near-perfect estimation of ϕ m enables an "unscrambling" of the Ramsey fringes and in consequence a reconstruction of single-qubit Ramsey fringes with contrast close to 1 instead of 1/2 as for the two-qubit parity fringe.This implies that we should be able to retrieve the noiseless precision bound of σ ∞ ∆φ = 2 M in the limit of N → ∞.Here, and in the remainder of this section, we assume that in the absence of correlated dephasing Ramsey fringes would have the full contrast, i. e. that there is no other source of decoherence.Later, this assumption will be dropped and we will also consider the influence of additional single-qubit dephasing on the measurement uncertainty.In the following, we will discuss different approaches for analyzing the multi-qubit correlations.
Correlation spectroscopy with many qubits: pair correlations.-Ramsey measurements of individual qubits contain no useful information as measuring Ẑi = |0 0| − |1 1| results in Ẑi = Tr(U X U φ ρU φ † U † X Ẑi ) = 0. Yet, information about transition frequency and position differences is obtained from correlation measurements [15], for which the correlated dephasing only reduces the maximum range of correlations by a factor of two.A fit of the correlation matrix C = (C ij ) yields estimates φi of the single-qubit phases φ i up to an irrelevant global offset phase.In experiments where ∆k = k 1 − k 2 = 0, this approach can be used to determine differences in transition frequencies up to a global sign factor.If, on the other hand, T = 0 and |∆k| = 0, information about the spatial arrangement of the qubits is obtained.
Single-qubit phase estimates with N-particle correlations.-For an estimation of the single-qubit phases φ = (φ 1 , φ 2 , . ..), we calculate the likelihood of observing the single-shot measurement outcome q = (q 1 , . . ., q N ), Given a set of measurements Q = (q im ), a maximum likelihood estimation of φ is obtained via evaluation of the log-likelihood function Note that the calculation of the integral in eq. ( 7) can be replaced by an average over N + 1 evenly distributed phases ϕ m = 2πm/N as the highest Fourier component of the integral kernel has a period of 2π/N .N-particle correlations for estimating the collective random phases ϕ m .-Once an estimate φi of the single-qubit phases is available, single-shot Ramsey spectroscopy can be used for estimating the random phase ϕ m of an experimental run from the vector of outcomes q m ≡ (q im ) N i=1 .Towards this end, we calculate the likelihood function and use it for a Bayesian estimate of the random phase φm = arg( 2π 0 This approach allows for tracking the temporal fluctuations of the local oscillator with respect to the qubit transition frequencies.
We note that the Bayesian approach can also be applied to estimating the vector of single-qubit phases φ.Given the estimates φm , the single-qubit phase differences ∆φ ij can be estimated by calculating the likelihood function (1 + q im sin(φ i + φm )), (11) where q i ≡ (q im ) M m=1 , in order to obtain the Bayesian estimate φi = arg( 2π 0 This approach is computationally fast albeit less precise than maximum-likelihood estimation (MLE) of φ.As discussed further below, the resulting uncertainties approach the ones obtained with maximum likelihood estimation only in the limit of large number qubits whereas the performance is unsatisfactory for small numbers of qubits.

III. EXPERIMENTAL IMPLEMENTATION AND MEASUREMENT RESULTS
Measurements on linear and planar 40 Ca + ion crystals are performed in two different experimental setups that will be described in the following.
The centerpiece of the apparatus for trapping planar crystals is a novel microfabricated monolithic linear Paul trap, shown in Fig. 2b, which allows us to create the anisotropic potentials required for trapping 2D ion crystals while simultaneously maintaining sufficient optical access perpendicular to the crystal plane for ion imaging.The trap provides a potential in which the ions are strongly confined in the direction perpendicular to the crystal plane, at an oscillation frequency of 2.196 MHz, and weakly confined along the two other directions, in which the crystal is extended, at oscillation frequencies of about 679.8 kHz and 343.0 kHz.Further details on this new ion trap apparatus can be found in Ref. 43.Ions are loaded into the trap via laser ablation and are Dopplercooled on the S 1/2 ↔ P 1/2 dipole transition.For encoding a qubit in an ion we use the two 4S 1/2 , m = ±1/2 Zeeman ground states, coherently coupled by a magnetic radiofrequency field oscillating at approximately 11.4 MHz.We distinguish the two qubit states by shelving the population of one of them in the long-lived 3D 5/2 Zeeman states, followed by fluorescence detection: The qubits are measured with high fidelity by exciting the ions on the S 1/2 ↔ P 1/2 transition and imaging the ion fluorescence onto an electron-multiplying CCD camera.For the shelving operation, we employ π-pulses induced by a frequency-stable 729 nm laser, coming from a direction perpendicular to the crystal plane.
In contrast to the apparatus for manipulating 2D crystals, long strings of 40 Ca + ions are trapped in a macro-scopic linear Paul trap providing a very anisotropic trapping potential with radial oscillation frequencies of about 2.5-3 MHz and an axial oscillation frequency of about 120 kHz.After Doppler cooling, the radial modes of the ion string are cooled close to the ground state by sideband cooling and the axial modes sub-Doppler cooled by polarization-gradient cooling [44].The qubit is encoded in one of the two 4S 1/2 Zeeman ground states and one of the metastable 3D 5/2 Zeeman states.The ionqubit can be coherently manipulated using 729 nm laser light resonantly exciting the S 1/2 ↔D 5/2 transition.Two laser beams with k-vectors parallel (perpendicular) to the linear ion crystal are available for collectively coupling to the qubits with the same coupling strength.Further details about this experimental setup can be found in Ref. 45.
In a first measurement, we investigate multi-qubit enhanced phase estimation in a 91-ion planar crystal; the results are shown in Fig. 3.We probe the ground-state qubits with a Ramsey probe time of 10 ms; here, magnetic field inhomogeneities gave rise to qubit-dependent phases φ i and correlated dephasing was the result of temporal fluctuations of the magnetic field's magnitude.Fig. 3(a) shows the measured pair correlations C ij used for a first estimate φi of the single-qubit phases shown in panel (b).Panel (c) displays the outcomes of an individual Ramsey experiment plotted against φi together with a single-shot Ramsey fringe obtained from an estimate of the collective random phase ϕ m .In (d), the matrix elements C ij are plotted versus the improved estimate φi − φj obtained by maximum-likelihood estimation based on eq. ( 8), for which we maximized the likelihood by a gradient-based optimization algorithm [46].The plot shows that the contrast of the resulting fringe is close to the maximum possible value.Similarly, averaging over experiments carried out at similar values of ϕ m results in single-qubit Ramsey fringes with contrast close to 1 (panel (e)).By subdividing the data sets into subsets of 200 measurements each, it is possible to measure the uncertainty of the phase difference estimates φi − φj .Pink data points in Fig. 3(f) show the uncertainty based on estimating the phase difference from the pair correlation between two qubits, which becomes minimal for a phase difference of π/2.The measured uncertainties are in agreement with the bound provided by quantum projection noise in the presence of correlated dephasing.The lowest uncertainty is obtained by maximum likelihood estimation using N -qubit correlations (blue data points, blue line: average over all points).The dashed line is the lower bound in the limit of M → ∞ and N → ∞.All data presented in Fig. 3 and subsequent figures is available online [47].
The same data set can also be used for investigating the measurement uncertainties as a function of the number of qubits as shown in Fig. 4. Towards this end, we split the data into subsets, each containing a fixed number of qubits with single-qubit phases that are approximately evenly distributed over the interval [0, 2π).The mea- surements of each of these sets is further split into subsets containing M = 200 experimental realizations from which we reconstruct the single-qubit phases for an estimate of the measurement uncertainties.Dark blue data points represent reconstructions based on MLE (eqs.( 7) and ( 8)), light blue points are the results of the noncompetitive Bayesian approach (eqs.( 11) and ( 12)).As shown in section IV, the uncertainties of the MLE estimates can be fitted by eq. ( 17) with a fringe contrast of C 0 = 0.995, which could result from state-assignment errors and slow drifts of trap parameters over the duration of the measurement.The inset shows that the phase uncertainty decreases inversely proportional to the number of samples in a given set and is thus still projection-noise limited at M = 10 4 samples.Note that we used an unbiased estimator for the determination of the uncertainties displayed in Fig. 3 and Fig. 4 assuming normally distributed measurement results [48].

IV. BOUNDS TO THE ACHIEVABLE PHASE ESTIMATION UNCERTAINTY
In this section, we compare the experimentally measured uncertainties to the theoretically achievable minimum uncertainties for an unbiased estimator and M experimental samples.The noiseless precision bound 2 M of eq. ( 4) cannot be experimentally achieved as it assumes noiseless dynamics, i.e. that the single-qubit Ramsey fringes (cf.Fig. 3e) can be measured with unity contrast.However, this assumption is unrealistic in noisy experiments affected by strong correlated dephasing and small levels of single-qubit dephasing.
The impact of these two noise sources on the precision is different.Uncorrelated dephasing reduces the fringe contrast to C 0 < 1, which unavoidably degrades the precision.In contrast, the effect of strong correlated dephasing can be overcome with a suitable data analysis for a large number of ions.This can be understood as follows: given a large number of ions, the random phase in each shot, ϕ m , can be estimated with an error that goes to zero as N → ∞.Another way to understand this is that for correlated dephasing there are decoherence-free subspaces (unlike the case of uncorrelated dephasing).The density matrix has elements inside and outside the decoherence-free subspaces, and as N → ∞ the contribution of the elements outside the protected subspaces goes to zero.Here, we will take these factors into consideration.We derive heuristic precision bounds that depend on both the number of qubits N and a finite Ramsey contrast C 0 (in the absence of correlated dephasing).We will start with a simple analytical model for an estimator based on N -particle correlations and compare the predictions to numerical simulations based on the classical Fisher information.

A. Measurement uncertainties with N-particle correlations
Here we derive a simple analytic model to approximate the precision bounds given correlated and uncorrelated dephasing.
Before we proceed to the model let us first understand the effect of uncorrelated dephasing and generalize the noiseless precision bound to a finite contrast limit that takes into account this dephasing and serves as a similar benchmark.Given an uncorrelated dephasing, the product state mentioned in section II becomes a mixed state: Measuring the i-th qubit in the Ŷi basis, the probability of 1 is: p = 1 2 (1 + C 0 sin (φ i )) .Then the uncertainty in the determination of p with M repetitions is given by projection noise, where is an effective number of measurements.Assuming φ i is drawn from a uniform distribution, the uncertainty becomes with The uncertainty in estimating ∆φ = φ i − φ j is σ 2 φi + σ 2 φj , and thus equal to Since this is the minimal obtainable uncertainty given a contrast of C 0 and assuming no correlated dephasing we refer to it as the finite contrast precision bound.
We derive now a simple model for precision bounds given also correlated dephasing.The idea is to first find the uncertainty in estimating the common random phase ϕ m and then insert this uncertainty as an uncorrelated dephasing of each qubit.
We first need to find the uncertainty in estimating ϕ m (given that all the other phases are known).Note that this is exactly the same calculation as performed above for φ i , just taking M = N , therefore: . This limits the precision with which the random phases ϕ m can be estimated.We can thus take the distribution of the random phase to be Gaussian with this variance.By averaging the single qubit probability over the Gaussian distribution we observe that the contrast of the unscrambled singlequbit Ramsey fringe (Fig. 3e) gets reduced to ), (16) if C 0 was the Ramsey contrast in the absence of correlated dephasing.We can now apply the same reasoning again to estimate the uncertainty with which the shift of unscrambled Ramsey fringes can be determined in order to estimate the uncertainty of the phase difference φ i −φ j which becomes For the case of large qubit number and high contrast C 0 , this expression can be approximated by eq. ( 15) if the replacement

B. Fisher information based bounds
We use a Fisher information (FI) analysis to calculate the achievable minimum uncertainties.According to the Cramer-Rao bound, the Fisher information matrix sets a bound on the achievable uncertainty with any unbiased estimator, where COV is the covariance matrix of the parameters φ = (φ i ) and I −1 is the inverse of the FI matrix.In case I is singular, i.e. information can be obtained only about a subspace of the parameters, I −1 is the Moore-Penrose pseudoinverse, defined only on this subspace.This implies that the variance of the phase difference φ i − φ j is given by where v ij is a column vector with components (v ij ) n = δ in − δ jn and I −1 the inverse of the relevant FI matrix.The Fisher information matrix I = (I ij ) can be calculated by the following formula: where φ = (φ i ) is the vector of parameters and p k the probability distribution of the observations.As a simple example, observe that for a single-parameter Bernoulli distribution, p (φ), the FI about φ is I = p(1−p) and hence . Given M identical independent Bernoulli trials, the FI about φ is multiplied by a factor of M and This expression coincides with the uncertainty of eqs.( 4)- (5).Furthermore, note that the FI is a generalization of N eff defined in section IV A.

Fisher information bound for the pair correlations
In pair correlation analysis, we estimate the ion phases (φ i ) using the pair correlations of the measurement outcomes, i.e. the correlation matrix C i,j defined in eq. ( 6) and presented in Fig. 3(a).More precisely we take the averages 1 M M m=1 q i,m q j,m i =j and estimate the phases according to it.According to the central limit theorem, the averages converge to a Gaussian random variable N (µ i,j ) i =j , M −1 Σ , where µ i,j = q i q j = 1 2 C 2 0 cos (2 (φ i − φ j )) and Σ is the covariance matrix Σ (i,j),(k,m) = q i q j q k q m − q i q j q k q m .An explicit calculation of the covariance matrix elements is presented in Appendix A.
Since the distribution is normal the FI matrix about (φ i ) is given by [49]: is the information gained due to the change in the mean values, i.e. the signal, and Σ is the covariance matrix of the different correlations representing the noise.Applying Eq. ( 21) for a single pair correlation (i, j) we retrieve the uncertainty in Eq. ( 5): the only linear combination of φ i , φ j that has a non-vanishing FI is φ i − φ j , for which the FI is The minimal uncertainty per measurement is 2, and a divergence occurs for φ i − φ j = nπ (n ∈ Z) due to the vanishing derivative and non-vanishing noise.
Since information about φ i − φ j is encoded not only in the (i,j) correlations but in other pairs as well, using all pairs improves the uncertainty, and removes the divergence around nπ.We use eq.( 21) to perform an exact numerical calculation of the FI.The behavior of the FI is presented in Figs. 10, 11 in Appendix A. It can be seen from the figures that as N → ∞ the FI with pair correlations does not saturate the noiseless precision bound.The reason for this is the information encoded only in higher moments.Using an analytical approximation we show in Appendix A that the variance for large N converges to , while the finite contrast precision bound to the variance is . As C 0 gets smaller the variance with pair correlations gets closer to the this bound since the information from higher moments becomes smaller.

Fisher information bound with N-particle correlations
When using the full counting statistics, the probability distribution entering into the calculation of the Fisher information matrix is given by eq. ( 7) with the replacement q i → C 0 q i in order to account for a Ramsey contrast C 0 < 1.In Appendix A, we numerically calculate the Fisher information matrix for finding the lower limit to the achievable uncertainty as a function of qubit number N and contrast C 0 (see Fig. 11).When N becomes large, an exact evaluation of the Fisher information matrix by eq. ( 20) becomes impractical as a summation over 2 N terms would have to be carried out.For N > 24, we sampled bit strings from the underlying probability distribution for a Monte-Carlo calculation of the empirical Fisher information matrix.The uncertainty achievable in experiments with a finite number of repetitions are numerically investigated in Appendix B.

Improving precision limits using entanglement
In our experiments, the qubits are initialized to a product state and measured in a local X basis.Hence, no entanglement occurs in these experiments, and an analysis based on classical Fisher information suffices.This raises the question of whether non-classical protocols that involve entangled states or different measurement bases can yield an advantage.It turns out that this is indeed the case: more general quantum protocols can obtain the noiseless precision bound of 2 M with an initial product state for every N and with an entangled initial state we can further reduce the uncertainty to We prove in Appendix C that this uncertainty is optimal.
To obtain these results we use the quantum Fisher information (QFI), which is the FI optimized over all possible measurement strategies [50,51].After averaging the quantum state over the random phase (Eq.( 3)), we show that the noiseless precision bound can be achieved for every N with a suitable measurement strategy (see appendix C 1).To gain intuition, let us examine the case of N = 2: when measuring in the local X basis, eq. ( 5) predicts σ ∆φ ≥ 2 √ M .However if we first measure the total number of excitations, i.e.Z 1 + Z 2 , and then measure in the local X basis, the noiseless precision bound of 2 M is achieved.Optimizing over both initial states and measurement strategies, we prove in appendix C 2 that the ultimate precision limit is N −1 N 2 M .Several initialization strategies saturate this bound, in particular any initial pure state with Z j = 0, Z j Z k = − 1 N −1 for all j = k achieves it.The reason for this improvement is the minimal j =k Z j Z k , which guarantees minimal uncertainty .
The reason this limit grows with N is frustration: one cannot make all pairs of spins anti-parallel.While the number of pairs is N 2 , the minimal and thus the optimal Z j Z k is − 1 N −1 .It can be immediately observed that the symmetric Dicke state with N/2 excitations satisfies these conditions and thus is optimal.Another optimal strategy is to employ a probabilistic initialization to products of antiparallel Bell states, i.e. in each experiment different pairs are being entangled to form an anti-parallel Bell state.With these two initialization strategies, the optimal sensitivity can be achieved with local measurements in the X or Y basis.This bound is plotted as a red curve in Fig. 4 along with other theoretical limits.A detailed derivation of the bound and the required initial states and measurements is presented in Appendix C 2.
These theoretical quantum limits imply that some improvement can indeed be obtained using entangled states or non-local measurements, however this improvement becomes negligible in the limit of large number of ions.This potential improvement and a comparison between the different precision limits is presented in Fig. 4.

V. APPLICATIONS IN TRAPPED-ION EXPERIMENTS
In the following, different applications of correlation spectroscopy in trapped-ion experiments will be presented.

A. Measurement of ion positions
Very anisotropic potentials are required for confining many ions in the form of a linear string.As a consequence of the weak axial confinement, these strings have lengths that are no longer small as compared to the distance between the ions and the nearest trap electrode.Therefore, the trapping potential can no longer be modeled as being purely harmonic and anharmonicities, which might affect the ion string's normal modes of motion, have to be considered.
We reconstruct the trapping potential in the axial direction by Ramsey experiments probing an optical qubit on the S 1/2 ↔D 5/2 transition, in which the first (second) π/2 pulse is realized by a laser beam impinging on the ions from the axial (perpendicular) direction.This setting results in qubit-specific phases φ i = kx i where k is the wave number and x i denotes the coordinate of the i th ion along the direction of the ion string.To suppress energy-dependent phase contributions, we use short π/2 pulses without any free-evolution time in between.Following the previously outlined procedure, we first reconstruct the qubit phase φ i and the measurement contrast by fitting the correlation matrix.Next we use these phases for reconstructing the time-dependent random phases ϕ m .Using N-qubit correlations, we finally use the Bayesian approach of eq. ( 12) for an improved phase estimate of φ i , shown as open symbols in Figure 5 (a) for a string of 62 ions.
In a second step, we extract the trapping potential from the measured correlations.We approximate the potential by Taylor-expanding it up to fourth order, , where ω 0 is the oscillation frequency of a single ion and l 3 (l 4 ) account for the cubic (quartic) anharmonicity of the potential.By calculating the ion positions in this potential, we fit the measured φ i and find ω 0 = (2π) 109.728(3)kHz, l 3 = 2.1(7) mm, l 4 = 0.8(6) mm, where the error bars are obtained from nonlinear regression assuming quantum projection noise as the only source of errors.We compare the measured phases φ i to the ones obtained from fitting the potential (φ f it i ) by calculating the residual position errors, δx i = (φ − φ f it i )/k.Fig. 5 (b) shows that these residuals have a standard deviation of 6.0 nm, barely above the theoretically expected error σ ∆φ = 5.1 nm.Moreover, the absence of spatial correlations in the residuals demonstrates that Taylor-expanding the potential up to the fourth order is an adequate approximation to the exact potential.
To further test the method, we carried out the reconstruction of the potential for a fixed set of trap parameters but different number of ions (10 ≤ N ≤ 62) and obtained consistent results.Fig. 5 (c) shows the inferred oscillation frequency ω 0 (light blue points) and the lowest collective mode frequency ω c (dark blue points).For an independent cross-check, the latter was also measured via sideband spectroscopy on the S 1/2 to D 5/2 transition (red squares).We observe that the correlation measurement systematically underestimates the mode frequency by about 220 Hz (Fig. 5 (d)).This discrepancy could be explained by the perpendicular laser beam being misaligned by about 1 mrad.Apart from this systematic error, the match between the two methods is quite good for N > 10 ions: the inset shows the difference of the predicted mode frequencies, which have a standard deviation of only 14 Hz if the N=10 data point is excluded on the basis of the rather uneven distribution of the phases φ i over the interval from 0 to 2π.Systematic effects in the measured frequencies by imperfect laser beam misalignment could be further reduced by replacing the perpendicular beam by another axial beam that is counterpropagating to the axial beam in place, because small alignment errors of the beams with the direction of the ion string would affect the measurement outcomes only in second order.

B. Measurement of transition frequency differences
Correlation spectroscopy with long probe times provides a tool for precisely measuring spatial transition frequency variations, which are relevant for frequency standards and quantum simulation experiments.For 40 Ca + ions, the dominant frequency shifts are Zeeman and electric quadrupole shifts.We measure the spatial dependence of these shifts by probing the stretched S 1/2 , m = ±1/2 ↔D 5/2 , m = ±5/2 transitions with a Ramsey time of τ = 40 ms duration with a 51-ion string.In contrast to the experiments of subsection V A, both Ramsey pulses are realized by the same laser beam.Writing the spatially resolved shifts ∆ ± i as and ∆ B i = (∆ + i − ∆ − i )/2 enables a separation of electric quadrupole and magnetic field shifts.
Figure 6(a,b) displays the measured quadrupole shift together with a calculated shift obtained from a measurement of the ion positions and the known quadrupole moment θ(3d, 5/2) of the D 5/2 level [26].The systematic variation of the residuals on the scale of 0.5 Hz could be explained by a 1.5σ error in the determination of θ(3d, 5/2) or by a misalignment of the perpendicular laser beam by 3 mrad.Figure 6(c,d) display the level shifts by the inhomogeneous magnetic field produced by the permanent magnets defining the quantization axis.We fit the Zeeman shifts with a third-order polynomial of the ion positions in order to extract the residuals.The latter have a standard deviation of 109 mHz, approaching the minimal uncertainty of 103 mHz predicted by the noiseless limit.
Figure 7(a) shows the measured single-ion phases for a two-dimensional 91-ion crystal in the presence of a spatially-varying magnetic field.We probe the groundstate transition, S 1/2 , m = −1/2 ↔ S 1/2 , m = +1/2, by a Ramsey experiment of τ = 5 ms duration.We fit a linear function to the measured phases and show the contour lines of constant phases from the fit in Fig. 7(a).We extract a magnetic-field gradient of 0.85(1) G/m from the linear fit with an angle of 38.6(4) degrees with respect to the horizontal direction.The maximal measured transition frequency difference between the ions is 218.2(8)Hz.The spatial distribution of the residuals from the linear fit, shown in Fig. 7(b), reveal that the magnetic field contains higher-order terms in addition to the linear gradient.We further fit a quadratic function to the residuals, and show the remaining residuals in Fig. 7(c,d) along the two orthogonal directions.The majority of the spatial structure in the magnetic-field can be explained with linear and quadratic terms, as the remaining residuals show almost no systematic structure.Experimentally, it is straightforward to cancel linear variations of the magnetic field across the ion crystal with permanent magnets or coils placed outside the vacuum system.

C. Single-shot Ramsey interferometry
The data taken for probing the spatial dependence of phase shifts can also be analyzed in the time-domain: we probe temporal fluctuations of the local oscillator's phase at the locations of the ions by single-shot Ramsey interferometry.Fig. 8 shows examples of such temporal phase changes that are caused by magnetic field fluc- tuations, laser frequency noise and optical path length fluctuations, respectively.Panel (a) shows a magneticfield change of about 3 µG at the location of the ions induced by the arrival of an elevator at the lab floor.The magnetic field was sensed by a 49-ion string probed by a 40 ms Ramsey experiment on the Zeeman ground state qubit transition.For the data shown in (b), the S 1/2 , m = 1/2 ↔D 5/2 , m = 3/2 transition was probed for τ = 20 ms.Here, 371 data point were acquired, each containing 50 experiments that were recorded at a repetition rate of 25 Hz.Laser phase noise gave rise to phase fluctuations for which an autocorrelation was calculated.The spectral density of the autocorrelation function reveals distinct components at low frequencies contributing to the laser noise.The dominant component at ∼ 8 Hz introduces a frequency excursion on the order of 1 Hz. Figure 9 shows differential path length fluctuations in the time domain, measured with short Ramsey experiments using two different laser beam paths for the two Ramsey pulses.For durations below 2 s, the data shows phase fluctuations (ϕ(t + τ ) − ϕ(t)) 2  t between experiments separated by a time τ that increase in proportion to τ as shown in the inset.The phase diffusion is predominantly caused by path length fluctuations in the two optical fibers delivering the light to the ion trap.

VI. DISCUSSION AND OUTLOOK
We have investigated many-qubit correlation spectroscopy for probing qubits subjected to spatially correlated noise.The technique enables phase comparisons between any pair of qubits, provided that the qubit states are about evenly distributed over the equatorial plane of the Bloch sphere.The latter condition does not impose a strong restriction, as in most experimental setups it should be possible to deliberately imprint spatial phase gradients on the qubit array to satisfy this requirement.In the limit of large qubit number and perfectly correlated noise, the quantum correlations induced by the noise enable a nearly complete restoration of the Ramsey contrast, which in the case of two-qubit correlation spectroscopy is upper-bounded to 50%.The increased contrast gives rise to a fourfold reduction in measurement time needed to achieve the targeted measurement uncertainty.
Many-qubit correlation spectroscopy is easy to implement as it requires only standard Ramsey spectroscopy enhanced by single-qubit read-out.The technique is therefore not limited to trapped-ion experiments but could be used in any multi-qubit physical system with high-fidelity single-shot read-out of individual qubits.In particular, it might be applicable to atomic clock experiments in tweezer arrays.Recently, experiments applying Ramsey correlation spectroscopy to subensembles of atoms held in optical lattices or tweezer arrays have demonstrated very small frequency gradients and impressive optical atomic coherence times reaching tens of seconds [25,52,53].In one-or two-dimensional tweezer arrays, which feature single-atom detection of tens to hundreds of atoms [25,31,54], our method is directly applicable and could assist in reducing the measurement time required for characterizing spatially varying transition frequency shifts across the atomic array.With the further development of atomic clocks networks connected by phase-stable photonic links [55,56], multi-qubit correlation spectroscopy could be applied for mutual frequency comparisons of the clocks, too.Another application of the technique might be found in quantum information processing experiments where spatially correlated noise can degrade the device performance.For example, in the atomic tweezer experiments reported in Ref.31, an auxialliary atomic species was employed for sensing and in-sequence correction of correlated phase noise.Here, an application of our protocol to the sensing species might increase the maximum noise level for which the correction can still be applied.
In the context of trapped-ion experiments, many-qubit correlation spectroscopy proves to be a valuable tool for characterizing various aspects of the experimental setup with high precision.Our experiments demonstrate that experimentally observed uncertainties come close to the theoretically predicted ones.The resulting reduction in measurement time for achieving a desired uncertainty could be of interest for tracking the frequency of an unstable laser and providing feedback from individual measurements for improving its stability.Another application of multi-qubit correlation spectroscopy, which we did not explore in this paper, is to use it for thermometry and detection of structural phase transitions in ion crystals [57,58].In this context, insufficiently cooled (lowfrequency) motional modes could give rise to a reduction of fringe contrast that could be detected in correlation spectroscopy experiments.
While in our experiments both the read-out and the initial states are non-entangled, we have theoretically shown the possibility of improving the precision using an entangled initial state or non-local measurements.While being negligible for large N , this improvement can be considerable assuming a relatively small N .A possible experimental realization requires initialization to a symmetric Dicke state.In trapped-ion experiments, these Dicke states could be engineered [59] by preparing the ions' center-of-mass mode in a Fock state with N/2 quanta, followed by a rapid adiabatic passage on its red-sideband transition [60], which converts motional quanta into collective electronic excitations [61].Finding simple protocols for generating optimal initial states and implementing these protocols in a sensing experiment is an interesting challenge for future work.
the European Union's Horizon 2020 research and innovation programme under grant agreement No 817482.Furthermore, we acknowledge support by the Austrian Science Fund through the SFB BeyondC (F7110) and funding by the Institut für Quanteninformation GmbH.This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 801110 and the Austrian Federal Ministry of Education, Science and Research (BMBWF).TG acknowledges funding provided by the Institute for Quantum Information and Matter and the Quantum Science and Technology Scholarship of the Israel Council for Higher Education.Pair correlations-We analyze the Fisher information obtained using only pair correlations.As mentioned in the main text the relevant random variables are the pair correlations which, according to the central limit theorem, converge to a Gaussian distribution, where µ i,j is the average of q i q j and Σ is the covariance matrix of the {q i q j } i,j>i .
Hence the problem boils down to calculating the FI matrix for this Gaussian distribution.The FI matrix about − → φ given this Gaussian distribution is presented in eq. ( 21) in the main text.We write it here as where: Hence we need to calculate {µ i,j } i,j>i and Σ in order to get the FI matrix.Let us first assume only correlated dephasing (no uncorrelated dephasing).For the mean values, we have Let us now calculate Σ.The diagonal terms of Σ read: Regarding the non-diagonal terms, let us begin with nonoverlapping pairs (i, j) , (k, n): Hence: For overlapping pairs, such as (i, j) , (i, n) we have: The derivatives matrix, D, is: Inserting equations (A3),(A4),(A5), (A6) into (A2), we can perform exact numerical calculations of the FI.Given an uncorrelated dephasing in addition to the correlated dephasing, the probabilities p ± of observing outcomes q i = ±1 are modified to It can be immediately observed that µ i,j = 1 2 C 2 0 cos (φ i − φ j ) .The covariance matrix is modified as follows: Let us analyze the uncertainty using pair correlations as N → ∞.The behavior in the limit of large N is presented in Fig. 10.It can be observed that for C 0 = 1 this uncertainty does not converge to the noiseless precision bound of 2 M , but approximately to 3 M .
The limit of 3 M can be derived analytically, assuming that the phases are distributed evenly in kπ (k ∈ Z) .To obtain this result, we use the following approximation of the variance: where u 1,2 is the vector that corresponds to φ 2 − φ 1 , i.e. (−1, 1, 0, ..., 0) .This approximation is obtained by using a Cauchy-Schwarz inequality twice: , and it can be understood as a single parameter estimation bound where the derivative of the mean is ): note that since the denominator goes as N 2 , we can omit in the calculation of this term any contributions that are smaller than N 2 .Since Du 1,2 is a real vector, this term is given by the sum (summation convention is used) (Du 1,2 ) (i,j) Σ (i,j)(k,n) (Du 1,2 ) (k,n) .We can neglect the N terms of identical pairs.From overlapping pairs (i, j) , (i, n) the contribution is: The contribution from the non-overlapping pairs is: 4 .Therefore altogether: which matches the numerical results.
For general contrast C 0 these expressions are modified to: Hence we obtain that for a general contrast var (φ 1 − φ 2 ) from pair correlations converges to ≈ .This implies that as C 0 becomes smaller the uncertainty using pair correlations converges to the finite contrast bound of the variance .The reason for this convergence is that the Fisher information obtained from higher moments goes with higher powers of C 0 , in general the Fisher information obtained from the 2k−th moments goes as C 2k 0 and thus the contribution from the higher moment gets smaller for smaller C 0 .In fact the FI with pair correlations coincides with the finite contrast bound up to a second order in C 2 0 : . This raises a natural question: is the variance with all m ≤ k particle correlations equal to N-qubit correlations-The information about the phase differences when taking all correlations into account is the information contained in the full distribution averaged over the random phase: (1 + q i sin(φ i + ϕ)).

(A9)
Hence the precision bound is given by the FI matrix about φ with this distribution.Since the FI matrix involves summation over all 2 N possible q vectors, I i,j = q (∂φ i p(q))(∂ φ j p(q)) p(q) , an exact calculation becomes intractable for large N. Hence to make an efficient calculation of the FI we use the fact that I i,j = (∂φ i p)(∂ φ j p)  simple least-squares estimation, i.e. minimizing V † V where , and maximumlikelihood estimation.Note that since the relevant distribution, eq.(A1), is Gaussian, the maximum-likelihood estimation becomes a weighted least-squares estimation [62]: The difference between the two estimation methods is thus rooted in the weights given by the inverse of the covariance matrix, Σ −1 .The maximum-likelihood is in general asymptotically efficient, i.e. saturates the Fisher information, whereas the simple least square is more straightforward as it does not require evaluation of the covariance matrix.
A comparison of both approaches is presented in Fig. 12.It can be observed that for a large number of samples (here M = 10 4 ) the maximum likelihood indeed saturates the FI, while the simple least-squares method does not saturate it.Interestingly, for a smaller number of samples (here M = 200) maximum likelihood does not saturate the FI and a simple least-squares approximation outperforms it.In fact, for some phases simple least squares even outperform the FI (due to its bias for small number of samples).
Estimating the phases from N-qubit correlations-A maximum-likelihood estimation of the phases by analysis of N-qubit correlations via eq.( 8) satisfies the FI-based bound in the limit of infinite sample size.We carried out numerical simulations of the phase estimation process to investigate the influence of a finite number of samples on the phase uncertainties.The simulations showed the uncertainty increasing over the FI-bound with decreasing sample size M ; however, the effect was not very pronounced: in simulations with 20 and 100 qubits, we observed an increase by about 10% for M = 50 and by about 1% for M = 500.12. Phase estimation errors using simple least squares and maximum likelihood estimation for 200 and 10 4 samples and 20 qubits (data from numerical simulation).Green lines and blue lines correspond to the distribution of the estimation errors with simple least squares and maximum likelihood respectively.The green circles and blue diamonds correspond to the average estimation error.The red solid line correspond to the Cramer-Rao limit.(a) For 200 samples the estimation errors are above the limit and simple least squares performs better than maximum likelihood.This behavior is due to the small number of samples.(b) For 10 4 samples the behavior matches the expectations: maximum likelihood coincides with the Fisher information and outperforms simple least squares.

Appendix C: Improving the precision with entangled states and non-local measurements
Let us inquire about the optimality of our scheme by asking the following question: what is the optimal precision optimized over all possible initial states and measurement strategies?The figure of merit is Var (φ j − φ k ) .We first show that the noiseless precision bound can be obtained by modifying the measurement basis to a non-local one.In a second step, we consider also entangled input states and find an optimal initial state for this sensing task.Let us introduce the following notation for this part: The state | z is the product state |z 1 |z 2 ...|z N where |z i is an eigenstate of Z i with eigenvalue z i ∈ {±1} .A different notation to the same state would be | q where q i = 1−zi 2 , thus z i = 1 → q i = 0, z i = −1 → q i = 1.

Improving precision with non-local measurements
Let us first write the quantum state after the time evolution.The initial state is a pure product state: |+ N , with |+ = 1 √ 2 (|0 + |1 ) (an eigenstate of X with eigenvalue +1).After a free evolution the state evolves into the state of Eq. ( 1) with phases φ im = φ i + ϕ m .ϕ m is a random phase that is distributed uniformly in [0, 2π] .This random phase induces a correlated dephasing, i.e. the final state, after averaging out ϕ m , becomes a mixture of Dicke states: where ρ j is a Dicke state with j excitation: Given ρ f φ , the fundamental precision limit is set by the quantum Fisher information matrix (QFIM) about the parameters φ 1 , φ 2 , ..., φ N ; this is the Fisher information matrix optimized over all possible measurement strategies.Hence the covariance matrix of the estimators, Σ, satisfies: where I is the QFIM and thus for any j, k: where u (j,k) is the parameter vector that corresponds to φ j − φ k .For a general mixed state ρ, given its spectral decomposition ρ = k p k |k k|, the QFIM is given by [50]: where {p k } k are the eigenvalues of ρ and the matrix elements,(•) k,l = k| • |l , are with respect to the eigenbasis of ρ.
For pure states this expression is reduced to: Let us calculate the QFIM of our ρ f , which we denote as I.It can be observed that in this special case I is a weighted sum of the QFIM of each |ψ j [63]: where I (j) is the QFIM of |ψ j .For every |ψ j we have ∂ φi |ψ j = −i1 2 (Z i + I) |ψ j , inserting this into equation (C1) we get that the QFIM of each |ψ j is: Now |ψ j is a symmetric superposition of all states with j excitations, from symmetry we get: and for k = l: Hence all the non-diagonal terms of I (j) are: .
The diagonal terms of I (j) read: Inserting these terms into equation (C2), we get that I reads: It is now simple to see that for any k = m the vector that corresponds to φ k − φ m is an eigenvector of I with an eigenvalue of 1.The variance per measurement is thus and this is exactly the noiseless precision bound.Since the strong commutativity condition is satisfied (all Hamiltonian terms commute with each other), we know that there exists a basis that saturates this QFI [64].This implies that there exists a measurement strategy such that the noiseless precision bound is obtained.
As a simple example we examine the case of two qubits.The density matrix for two qubits reads hence the bound is saturated.A general optimal measurement strategy would be: 1. first measure i Z i (this measurement collapses the density matrix into one of the Dicke states, |ψ j ). 2. Measure the Dicke state, |ψ j , in its optimal measurement basis.
The optimal measurement basis of |ψ j can be written implicitly as proven in [64]: projecting into a (Gram-Schmidt) orthogonalization of {|ψ j , |∂ φ k ψ j } φ k would be optimal.For example, for N = 3, given The construction for |ψ 2 is equivalent.

Optimal initial states
We show that the average variance of phase difference | z , and any other state that is an equal superposition of states with i z i = 0 saturate this optimal precision.Another strategy is a probabilistic initialization from an ensemble of products of anti-parallel Bell pairs.It can be shown that all optimal strategies involve eigenstates of i Z i with eigenvalue 0. Since these states are robust against correlated dephasing, this optimal precision is achieved irrespective of whether there is correlated dephasing or not.In the following derivation we use techniques similar to those used in reference [40].
Using the QFIM we prove that the optimal achievable variance, 1 Therefore we seek to lower bound j<k 2 (1− Zj Z k ) .The minimal possible var (φ j − φ k ) is therefore obtained when Z j Z k = −1, however there is no state that satisfies Z j Z k = −1 for all j, k.To lower bound This basically proves that 1 ( N 2 ) j<k Var (φ j − φ k ) ≥ 2 N −1 N .To show that this lower bound is saturable we need to find an initial state |ψ , for which all these inequalities are saturated, namely: where I is the QFIM given this |ψ .We observed that a necessary condition is Z i = 0 and identical Z j Z k = − 1 N −1 .Let us show that this is also a sufficient condition: given that this condition is satisfied the QFIM is It can be now observed that any u (j,k) is an eigenvector of this matrix with eigenvalue N N −1 and thus for any j, k : u † (j,k) I −1 u (j,k) = 2 N −1 N .Hence any initial pure state that satisfies the conditions: saturates this QFIM.We can immediately observe that the symmetric Dicke state: | z satisfies these conditions and thus saturates this bound.Other strategies exist, such as preparing a classical ensemble of products of anti-parallel Bell states, and they will be discussed later.For now Let us focus on the symmetric Dicke state.
(other than the symmetric Dicke state) is then basically solving a system of linear equations.First observe that the conditions imply ( i Z i ) 2 = 0, hence any optimal pure state is an eigenstate of i Z i with eigenvalue 0. We can therefore write the states as p z z i = 0 and l p l = 1, with a constraint of 0 ≤ p l ≤ 1 for all l.The number of equations, N 2 +N +1, is smaller than the number of variables, N N/2 , which implies that there exist solutions other than the symmetric one.Finding solutions that are simple to prepare, in terms of entanglement or circuit complexity, is an interesting problem and we leave it as an open question.

FIG. 1 .
FIG. 1. Measurement scenario.(a) In a network of quantum sensors comprised of qubits Qi, the qubits are prepared in Bloch states lying in the equatorial plane and subjected to correlated dephasing that randomly rotates all Bloch vectors by the same angle, as indicated by solid and dashed arrows.(b)Correlations Ci,j between pairs of qubits (i, j) are measured for an estimation of the angle between the respective Bloch vectors.Is it possible to reduce the measurement uncertainty of Ci,j obtained from a finite number of experimental repetitions by taking into account all pair correlations that can be simultaneously measured, or even all N -qubit correlations?

FIG. 3 .
FIG. 3. Many-qubit correlation spectroscopy of a 91-ion planar crystal based on M = 26852 experimental repetitions.(a) Measured correlation matrix with correlations |Cij| ≤ 1/2 limited by correlated dephasing.(b) Single-qubit phases φi estimated by fitting the correlation matrix.(c) Measurement outcomes qim of an individual experiment (black dots) used for a Bayesian estimate of the common random phase φm (blue curve: fitted Ramsey fringe).(d) Correlation matrix elements Cij plotted as a function of the phase difference φi − φj obtained by analyzing N-qubit correlations.(e) Single-qubit Ramsey fringes with nearly full contrast obtained from binning into sets of similar common random phase ϕm.The red and the blue curve are just two out of 91 measured fringes.(f) Measurement uncertainties inferred from subdividing the data into 134 data sets with 200repetitions each.The pink dots do not cover the entire range of 0 to π as we omit those qubit pairs (i,j) for which we measure |Cij| > 0.5 for one or several subsets.Uncertainties obtained from individual elements Cij (pink dots) and the analysis of N-qubit correlations (blue dots, solid light blue line: average over the data points).The dashed red line indicates the noiseless precision bound achievable in the limit of N → ∞, the solid black line the two-qubit limit σ N =2 ∆φ .

FIG. 4 .
FIG. 4. Phase uncertainties vs number of qubits for M = 200 samples.Dark blue dots represent uncertainties estimated from experimental data by MLE, light blue dots the uncertainties of the Bayesian estimation.The prediction of eq.(17) is shown as the solid blue curve for a contrast C0 = 1 and as a dashed blue curve for C0 = 0.995.The latter is obtained by fitting the experimental data.The black curve represents the noiseless precision bound of

FIG. 6 .
FIG. 6. Transition frequency shift measurement obtained by probing the quadrupole transitions between stretched states with 1500 experimental repetitions each.(a) Quadrupole shift of the D 5/2 , m = ±5/2 states (red circles: measured shift, black line: predicted shift).The frequency shift was measured with respect to the first ion; in the figure, a constant offset was added so that the averaged shift equaled the calculated average quadrupole shift.(b) Measurement residuals.(c) Differential Zeeman shift of the S 1/2 , m = 1/2 ↔D 5/2 , m = 5/2 transition frequency.An offset was added so that the average shift became equal to zero.The black line is a fit to the data by a third-order polynomial.(d) Residuals.The gray rectangle indicates the measurement uncertainty (1σ) predicted for quantum projection noise.

FIG. 7 .
FIG. 7. Transition frequency shift measurement in a twodimensional crystal obtained by probing the ground-state transition with 19736 experimental repetitions.(a) Singleion phases together with contour lines of constant phases obtained from a linear fit.(b) Residuals from the linear fit.(c,d) Residuals from the quadratic fit along (c) x and (d) z axis.The gray rectangle indicates the 1σ measurement uncertainty predicted for quantum projection noise.

FIG. 8 .
FIG. 8. Single-shot Ramsey interferometry tracking temporal phase fluctuations.(a) Phase fluctuations induced by a time-varying magnetic field probed with a 49-ion string.(b,c) Tracking laser frequency variations with a 51-ion string: Autocorrelation function A(T ) of laser phase fluctuations (left plot) together with its spectral density |F[A(T )]| 2 (right plot).

FIG. 9 .
FIG. 9. Measurement of temporal relative optical path length fluctuations in two beam path delivering the laser pulses to a string of 40 ions.
Appendix A: Bounds to the achievable phase estimation uncertainty

C 2 0 4 l?
We leave it as an open question (and conjecture) as we do not have a proof to this.

p 2 =
∂ φi ln (p) ∂ φj ln (p) .This allows us to make a Monte-Carlo calculation of the FI matrix by sampling ∂ φi ln (p) ∂ φj ln (p) .Simulation results are shown in Fig. 11 for the case of evenly distributed single-qubit phases, φ j = 2πj/N .

FIG. 10 .
FIG. 10.Sensitivity (per measurement) with pair correlations compared to precision bounds (numerical analysis).(a) Uncertainty in estimating ∆φ using all pair correlations as a function of ∆φ for different number of ions.The dashed red line is the noiseless precision bound, pair correlations of N = 20, 50, 150 correspond to light blue (top), red (middle), dark blue (bottom) points respectively.Inset: standard deviation averaged over all phases as a function of N , the top (bottom) lines correspond to pair correlations (noiseless precision bound).(b) Uncertainty as a function of the contrast C0.The red dashed line corresponds to the finite contrast precision bound.Light blue (top) and blue (bottom) points correspond to pair correlations for N = 30, 120 respectively.

5 FIG. 11 .
FIG. 11.Uncertainty with N-qubit correlations, pair correlations and precision limits (all values are per measurement).(a) precision bounds as a function of N : grey (upper) curve and blue dots correspond to pair correlations and N -qubit correlations respectively.The light blue curve correspond to the analytical approximation and the red dashed line is the fundamental noiseless limit of √ 2. (b) Precision bounds as a function of C0: Pink dots and grey (upper) curve correspond to all N correlations and pair correlations respectively for N = 30.Same for blue dots and grey lower curve for N = 100.The red dashed line correspond to the finite contrast precision bound.
FIG.12.Phase estimation errors using simple least squares and maximum likelihood estimation for 200 and 10 4 samples and 20 qubits (data from numerical simulation).Green lines and blue lines correspond to the distribution of the estimation errors with simple least squares and maximum likelihood respectively.The green circles and blue diamonds correspond to the average estimation error.The red solid line correspond to the Cramer-Rao limit.(a) For 200 samples the estimation errors are above the limit and simple least squares performs better than maximum likelihood.This behavior is due to the small number of samples.(b) For 10 4 samples the behavior matches the expectations: maximum likelihood coincides with the Fisher information and outperforms simple least squares.

i zi=0 √p z z j z k = − 1 N − 1
p z | z .The conditions then become a system of N 2 + N + 1 linear equa-tions for the distribution {p z } : ∀j, k, i i zi=0 , i zi=0 |01 + e i(φ1−φ2) |10 .Measuring the local X basis, we project onto the states | x = |x 1 |x 2 ...|x N where |x i is an eigenstate of X i with eigenvalues x i ∈ {±1} and obtain the probabilities: clearly this is exactly the variance with a single pair correlation and thus does not saturate the QFI.It can be observed that optimizing over all local measurement bases is equivalent to optimizing over φ 1 , φ 2 and thus no local measurement saturates the QFI.There exists however a non-local measurement strategy that saturates the QFI: consider first measuring Z 1 + Z 2 and then measuring the local X basis.With probability 1/2 we get |00 , |11 in the first measurement and thus no information and with probability 1/2 we collapse into |ψ 1 which yields a Fisher information of 1. Therefore the total Fisher information is: initialization strategies, all involve entanglement, that saturate this bound.In particular the symmetric Dicke state: 1 ( N 2 ) j>k Var (φ j − φ k ) is lower-bounded by 2 N −1 Nand we find several i zi=0