# Correlation spectroscopy with multi-qubit-enhanced phase estimation

H. Hainzer  Institut für Quantenoptik und Quanteninformation, Österreichische
Akademie der Wissenschaften, Technikerstraße 21a, 6020 Innsbruck, Austria
Institut für Experimentalphysik, Universität Innsbruck, Technikerstraße 25,
6020 Innsbruck, Austria  D. Kiesenhofer  Institut für Quantenoptik und
Quanteninformation, Österreichische Akademie der Wissenschaften,
Technikerstraße 21a, 6020 Innsbruck, Austria  Institut für Experimentalphysik,
Universität Innsbruck, Technikerstraße 25, 6020 Innsbruck, Austria  T.
Ollikainen  Institut für Quantenoptik und Quanteninformation, Österreichische
Akademie der Wissenschaften, Technikerstraße 21a, 6020 Innsbruck, Austria  M.
Bock  Institut für Quantenoptik und Quanteninformation, Österreichische
Akademie der Wissenschaften, Technikerstraße 21a, 6020 Innsbruck, Austria  F.
Kranzl  Institut für Quantenoptik und Quanteninformation, Österreichische
Akademie der Wissenschaften, Technikerstraße 21a, 6020 Innsbruck, Austria
Institut für Experimentalphysik, Universität Innsbruck, Technikerstraße 25,
6020 Innsbruck, Austria  M. K. Joshi  Institut für Quantenoptik und
Quanteninformation, Österreichische Akademie der Wissenschaften,
Technikerstraße 21a, 6020 Innsbruck, Austria  G. Yoeli  Racah Institute of
Physics, The Hebrew University of Jerusalem, Givat Ram, Jerusalem, 91904,
Israel  R. Blatt  Institut für Quantenoptik und Quanteninformation,
Österreichische Akademie der Wissenschaften, Technikerstraße 21a, 6020
Innsbruck, Austria  Institut für Experimentalphysik, Universität Innsbruck,
Technikerstraße 25, 6020 Innsbruck, Austria  T. Gefen<EMAIL_ADDRESS>Institute for Quantum Information and Matter, Caltech, Pasadena, CA, USA  C.
F. Roos<EMAIL_ADDRESS>Institut für Quantenoptik und
Quanteninformation, Österreichische Akademie der Wissenschaften,
Technikerstraße 21a, 6020 Innsbruck, Austria  Institut für Experimentalphysik,
Universität Innsbruck, Technikerstraße 25, 6020 Innsbruck, Austria

###### Abstract

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 quantum-enhanced 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
$\pi/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 $\pi/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.

Figure 1: Measurement scenario. (a) In a network of quantum sensors comprised
of qubits $Q_{i}$, 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 $C_{i,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 $C_{i,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 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 $\pi/2$ pulse,
measuring all qubits in the energy eigenbasis and correlating the measurement
outcomes ($\pm 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\choose 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 discuss 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 trapped-ion 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.

Figure 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 $\varphi$, and single-qubit phase
shifts $\phi_{n}$. The analysis of correlations between measurement outcomes
on different qubits taken at the same time (column vector $\mathbf{q}_{m}$)
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 $\mathbf{q}_{n}$) 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.

##  II Model: N-qubit correlation spectroscopy

We consider a data set consisting of $m=1,\ldots,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
$\pi/2$ pulse, the state of the $N$ qubits during the m-th realization is

$2^{-\frac{N}{2}}\underset{i=1}{\overset{N}{\Pi}}(|0\rangle+e^{i\phi_{im}}|1\rangle),$
(1)

with phases

$\phi_{im}=\phi_{i}+\varphi_{m}$ (2)

where $\phi_{i}$ is a qubit-specific phase and $\varphi_{m}$ a random phase
that is common to all qubits, i. e. in experiments $\phi_{i}$ appears as a
spatially varying phase whereas $\varphi_{m}$ encodes temporal changes. To
achieve an unambiguous definition of these phases, we define $\varphi_{m}$ to
be the phase change between the $m$th experiment and the first one. The phase
$\varphi_{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 qubit-specific
phases $\phi_{i}=(\mathbf{k}_{1}-\mathbf{k}_{2})\mathbf{r}_{i}+\Delta_{i}T$
arise if the qubits have different detunings $\Delta_{i}$ with respect to the
local oscillator frequency or if the qubits are excited from different spatial
directions for the first and the second $\pi/2$ pulse. Here,
$\mathbf{k}_{1}(\mathbf{k}_{2})$ is the k-vector of the running wave inducing
the first (second) $\pi/2$ pulse and $\mathbf{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 $m$th measurement in the state
$|0\rangle$ or $|1\rangle$. The probability of observing the outcome $q_{im}$
is given by $p(q_{im})=\frac{1}{2}(1+q_{im}\sin\phi_{im})$.

Here, we will study two closely related problems:

1\. We want to carry out a multi-parameter estimation of the qubit-dependent
phases $\phi_{i}$ in experiments where the random phases $\varphi_{m}$ are
uniformly distributed over the interval $[0,2\pi)$; 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

$\displaystyle\begin{split}&\rho=\frac{1}{2\pi}\int_{0}^{2\pi}d\varphi|\Psi(\varphi)\rangle\langle\Psi(\varphi)|\;\;\;\mbox{with}\\\
&|\Psi(\varphi)\rangle=2^{-\frac{N}{2}}\prod_{i=1}^{N}(|0\rangle_{i}+e^{i\varphi}|1\rangle_{i}),\end{split}$
(3)

a state, which contains no entanglement but quantum correlations in the form
of non-zero quantum discord [38, 20]. Next, the qubits are subjected to the
unitary operation
$U_{\bm{\phi}}=\exp(\frac{i}{2}\sum_{i}\phi_{i}\sigma_{i}^{z})$ followed by a
global $\pi/2$ pulse around the x-axis,
$U_{X}=\exp(-i\frac{\pi}{4}\sum_{i}\sigma_{i}^{x})$, 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 $\phi_{i}-\phi_{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 $\varphi_{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 $\phi_{i}$ can equally well be applied to
an estimation of $\varphi_{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^{-\frac{N}{2}}\prod_{i=1}^{N}(|0\rangle_{i}+e^{i\phi_{i}}|1\rangle_{i}),$
and in the absence of noise, the precision in estimating each phase
independently from $M$ measurements is $\sigma_{\phi_{i}}=\frac{1}{\sqrt{M}}.$
Hence the minimal uncertainty in estimating a phase difference
$\Delta\phi=\phi_{2}-\phi_{1}$ is

$\sigma_{\Delta\phi}=\sqrt{\sigma_{\phi_{1}}^{2}+\sigma_{\phi_{2}}^{2}}=\sqrt{\frac{2}{M}}.$
(4)

This approach basically amounts to inferring $\Delta\phi$ 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 $\Delta\phi$ is estimated using standard correlation spectroscopy.
Using error propagation of quantum projection noise, the uncertainty in the
estimation of $\Delta\phi$ from $M$ measurements equals

$\sigma_{\Delta\phi}^{N=2}=\frac{\sqrt{4-\cos^{2}\left(\phi_{i}-\phi_{j}\right)}}{\sqrt{M}|\sin\left(\phi_{i}-\phi_{j}\right)|}.$
(5)

The uncertainty diverges when the phase difference approaches 0 or $\pi$, i.
e. the points where the parity reaches an extremum. It becomes minimal for
$\Delta\phi=\pi/2$, where $\min\sigma_{\Delta\phi}^{N=2}=\frac{2}{\sqrt{M}}$,
which is larger by a factor of $\sqrt{2}$ than the noiseless precision bound.
The $\sqrt{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 $\binom{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 $\Delta\phi_{ij}=\phi_{i}-\phi_{j}$ from the
observed correlations makes it possible to estimate the random phases
$\varphi_{m}$ of each experimental realization. In the limit of a large number
of qubits, the near-perfect estimation of $\varphi_{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
$\sigma_{\Delta\phi}^{\infty}=\sqrt{\frac{2}{M}}$ in the limit of
$N\rightarrow\infty$.

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
$\hat{Z}_{i}=|0\rangle\langle 0|-|1\rangle\langle 1|$ results in
$\langle\hat{Z}_{i}\rangle=\mbox{Tr}(U_{X}U_{\bm{\phi}}\rho{U_{\bm{\phi}}}^{\dagger}U_{X}^{\dagger}\hat{Z}_{i})=0$.
Yet, information about transition frequency and position differences is
obtained from correlation measurements [15],

$C_{ij}\equiv\langle\hat{Z}_{i}\hat{Z}_{j}\rangle=\frac{1}{2}\cos(\phi_{i}-\phi_{j}),$
(6)

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 $\hat{\phi}_{i}$ of the single-qubit phases $\phi_{i}$ up to
an irrelevant global offset phase. In experiments where $\mathbf{\Delta
k}=\mathbf{k}_{1}-\mathbf{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 $|\mathbf{\Delta k}|\neq 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 $\bm{\phi}=(\phi_{1},\phi_{2},\ldots)$,
we calculate the likelihood of observing the single-shot measurement outcome
$\mathbf{q}=(q_{1},\ldots,q_{N})$,

$P(\mathbf{q}|\bm{\phi})=\frac{2^{-N}}{2\pi}\int_{0}^{2\pi}d\varphi\prod_{i=1}^{N}(1+q_{i}\sin(\phi_{i}+\varphi)).$
(7)

Given a set of measurements $Q=(q_{im})$, a maximum likelihood estimation of
$\bm{\phi}$ is obtained via evaluation of the log-likelihood function

${\cal L}(Q|\bm{\phi})=\log\prod_{m=1}^{M}P(\mathbf{q}_{m}|\bm{\phi}).$ (8)

Note that the calculation of the integral in eq. (7) can be replaced by an
average over $N+1$ evenly distributed phases $\varphi_{m}=2\pi m/N$ as the
highest Fourier component of the integral kernel has a period of $2\pi/N$.

_N-particle correlations for estimating the collective random phases
$\varphi_{m}$._ — Once an estimate $\hat{\phi}_{i}$ of the single-qubit phases
is available, single-shot Ramsey spectroscopy can be used for estimating the
random phase $\varphi_{m}$ of an experimental run from the vector of outcomes
$\mathbf{q}_{m}\equiv(q_{im})_{i=1}^{N}$. Towards this end, we calculate the
likelihood function

$P(\mathbf{q}_{m}|\\{\hat{\phi}_{i}\\},\varphi)=2^{-N}\prod_{i=1}^{N}(1+q_{im}\sin(\hat{\phi}_{i}+\varphi))$
(9)

and use it for a Bayesian estimate of the random phase

$\hat{\varphi}_{m}=\arg(\int_{0}^{2\pi}d\varphi
e^{i\varphi}P(\mathbf{q}_{m}|\\{\hat{\phi}_{i}\\},\varphi)).$ (10)

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 $\bm{\phi}$. Given the estimates
$\hat{\varphi}_{m}$, the single-qubit phase differences $\Delta\phi_{ij}$ can
be estimated by calculating the likelihood function

$P(\mathbf{q}_{i}|\phi_{i},\\{\hat{\varphi}_{m}\\})=2^{-N}\prod_{m=1}^{M}(1+q_{im}\sin(\phi_{i}+\hat{\varphi}_{m})),$
(11)

where $\mathbf{q}_{i}\equiv(q_{im})_{m=1}^{M}$, in order to obtain the
Bayesian estimate

$\hat{\phi}_{i}=\arg(\int_{0}^{2\pi}d\phi
e^{i\phi}P(\mathbf{q}_{i}|\\{\hat{\varphi}_{m}\\},\phi)).$ (12)

This approach is computationally fast albeit less precise than maximum-
likelihood estimation (MLE) of $\bm{\phi}$. 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 40Ca+ ion crystals are performed in two
different experimental setups that will be described in the following.

Figure 3: Many-qubit correlation spectroscopy of a 91-ion planar crystal based
on $M=26852$ experimental repetitions. (a) Measured correlation matrix with
correlations $|C_{ij}|\leq 1/2$ limited by correlated dephasing. (b) Single-
qubit phases $\hat{\phi}_{i}$ estimated by fitting the correlation matrix. (c)
Measurement outcomes $q_{im}$ of an individual experiment (black dots) used
for a Bayesian estimate of the common random phase $\hat{\varphi}_{m}$ (blue
curve: fitted Ramsey fringe). (d) Correlation matrix elements $C_{ij}$ plotted
as a function of the phase difference $\phi_{i}-\phi_{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
$\varphi_{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 200 repetitions each. The pink dots do not cover the entire
range of 0 to $\pi$ as we omit those qubit pairs (i,j) for which we measure
$|C_{ij}|>0.5$ for one or several subsets. Uncertainties obtained from
individual elements $C_{ij}$ (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\rightarrow\infty$, the solid black line the two-qubit limit
$\sigma_{\Delta\phi}^{N=2}$.

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 Doppler-cooled on the
$\mathrm{S}_{1/2}\leftrightarrow\mathrm{P}_{1/2}$ dipole transition. For
encoding a qubit in an ion we use the two 4S${}_{1/2},m=\pm 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 3D5/2 Zeeman states, followed by
fluorescence detection: The qubits are measured with high fidelity by exciting
the ions on the $\mathrm{S}_{1/2}\leftrightarrow\mathrm{P}_{1/2}$ transition
and imaging the ion fluorescence onto an electron-multiplying CCD camera. For
the shelving operation, we employ $\pi$-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
40Ca+ ions are trapped in a macroscopic 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
4S1/2 Zeeman ground states and one of the metastable 3D5/2 Zeeman states. The
ion-qubit can be coherently manipulated using 729 nm laser light resonantly
exciting the S${}_{1/2}\leftrightarrow$D5/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 $\phi_{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 $\hat{\phi}_{i}$ of the single-qubit phases shown in panel (b).
Panel (c) displays the outcomes of an individual Ramsey experiment plotted
against $\hat{\phi}_{i}$ together with a single-shot Ramsey fringe obtained
from an estimate of the collective random phase $\varphi_{m}$. In (d), the
matrix elements $C_{ij}$ are plotted versus the improved estimate
$\hat{\phi}_{i}-\hat{\phi}_{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 $\varphi_{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 $\hat{\phi}_{i}-\hat{\phi}_{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 $\pi/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\rightarrow\infty$ and $N\rightarrow\infty$. All data presented in
Fig. 3 and subsequent figures is available online [47].

Figure 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 $C_{0}=1$ and as a
dashed blue curve for $C_{0}=0.995$. The latter is obtained by fitting the
experimental data. The black curve represents the noiseless precision bound of
$\sqrt{\frac{2}{M}}.$ The red curve represents the optimum uncertainty that is
obtainable with entangled input states (see Appendix C). The reduction in
measurement uncertainty provided by preparation of entangled input states
rapidly shrinks with increasing $N$. The inset displays the measured
uncertainty for $N=91$ on the number samples, with the black curve
representing the noiseless precision bound.

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\pi)$. The measurements 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 non-competitive 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 $\sqrt{\frac{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, $\varphi_{m}$, can be estimated with an
error that goes to zero as $N\rightarrow\infty.$ 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\rightarrow\infty$
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.

###  IV.1 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:
$\frac{1}{2^{N}}\underset{i}{\prod}\left(\mathrm{I}+C_{0}\hat{R}_{\phi_{i}}\right)$,
with
$\hat{R}_{\phi_{i}}=\cos\left(\phi_{i}\right)\hat{X}_{i}+\sin\left(\phi_{i}\right)\hat{Y}_{i}.$
Measuring the $i$-th qubit in the $\hat{Y}_{i}$ basis, the probability of $1$
is: $p=\frac{1}{2}\left(1+C_{0}\sin\left(\phi_{i}\right)\right).$ Then the
uncertainty in the determination of $p$ with $M$ repetitions is given by
projection noise, $\sigma_{p}=\sqrt{\frac{p(1-p)}{M}}$. Because of
$\sigma_{p}=\left|\frac{dp}{d\phi_{i}}\right|\sigma_{\phi_{i}}$ and
$\frac{dp}{d\phi_{i}}=\frac{C_{0}}{2}\cos\left(\phi_{i}\right)$, we have
$\sigma_{\phi_{i}}=\frac{1}{\sqrt{N_{\mathrm{eff}}(C_{0},\phi_{i})}}$ where

$N_{\mathrm{eff}}(C_{0},\phi_{i})=MC_{0}^{2}\frac{(1-\sin^{2}\phi_{i})}{1-C_{0}^{2}\sin^{2}\phi_{i}}$
(13)

is an effective number of measurements. Assuming $\phi_{i}$ is drawn from a
uniform distribution, the uncertainty becomes
$\sigma_{\phi_{i}}=\frac{1}{\sqrt{N_{\mathrm{eff}}(C_{0})}}$ with

$\displaystyle N_{\mathrm{eff}}(C_{0})$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}d\phi_{i}\,N_{\mathrm{eff}}(C_{0},\phi_{i})$
$\displaystyle=M(1-\sqrt{1-C_{0}^{2}}).$ (14)

The uncertainty in estimating $\Delta\phi=\phi_{i}-\phi_{j}$ is
$\sqrt{\sigma_{\phi_{i}}^{2}+\sigma_{\phi_{j}}^{2}},$ and thus equal to

$\sigma_{\Delta\phi}=\frac{\sqrt{2}}{\sqrt{M}\sqrt{1-\sqrt{1-C_{0}^{2}}}}.$
(15)

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 $\varphi_{m}$ and then insert this uncertainty as an uncorrelated
dephasing of each qubit.

We first need to find the uncertainty in estimating $\varphi_{m}$ (given that
all the other phases are known). Note that this is exactly the same
calculation as performed above for $\phi_{i}$, just taking $M=N$, therefore:
$\sigma_{\varphi_{m}}=\frac{1}{\sqrt{N}\sqrt{1-\sqrt{1-C_{0}^{2}}}}$. This
limits the precision with which the random phases $\varphi_{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 ${\cal N}\left(\varphi_{m},\sigma^{2}\right)$,

$\langle\cos\left(\phi_{i}+\varphi\right)\rangle_{{\cal
N}\left(\varphi_{m},\sigma^{2}\right)}=\cos\left(\phi_{i}+\varphi_{m}\right)\exp\left(-\frac{\sigma^{2}}{2}\right),$

we observe that the contrast of the unscrambled single-qubit Ramsey fringe
(Fig. 3e) gets reduced to

$C_{unscr}=C_{0}\exp(-\frac{1}{2N(1-\sqrt{1-C_{0}^{2}})}),$ (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 $\phi_{i}-\phi_{j}$ which becomes

$\sigma_{\Delta\phi}(N,C_{0})=\frac{\sqrt{2}}{\sqrt{M(1-\sqrt{1-C_{unscr}^{2}})}}.$
(17)

For the case of large qubit number and high contrast $C_{0}$, this expression
can be approximated by eq. (15) if the replacement $C_{0}\rightarrow
C_{0}\exp(-\frac{1}{N})$ is made.

###  IV.2 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,

$\text{COV}\left(\bm{\phi}\right)\geq I^{-1},$ (18)

where COV is the covariance matrix of the parameters $\bm{\phi}=(\phi_{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 $\phi_{i}-\phi_{j}$ is given
by

$\mbox{Var}(\phi_{i}-\phi_{j})\geq\mathbf{v}_{ij}^{t}I^{-1}\mathbf{v}_{ij},$
(19)

where $\mathbf{v}_{ij}$ is a column vector with components
$(\mathbf{v}_{ij})_{n}=\delta_{in}-\delta_{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:

$I_{ij}=\sum_{k}p_{k}^{-1}(\bm{\phi})\frac{\partial
p_{k}(\bm{\phi})}{\partial\phi_{i}}\frac{\partial
p_{k}(\bm{\phi})}{\partial\phi_{j}},$ (20)

where $\bm{\phi}=(\phi_{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\left(\phi\right)$, the
FI about $\phi$ is
$I=\frac{\left(\partial_{\phi}p\right)^{2}}{p\left(1-p\right)}$ and hence
$\sigma_{\phi}=\frac{\sqrt{p\left(1-p\right)}}{\partial_{\phi}p}.$ Given $M$
identical independent Bernoulli trials, the FI about $\phi$ is multiplied by a
factor of $M$ and thus
$\sigma_{\phi}=\frac{\sqrt{p\left(1-p\right)}}{\sqrt{M}\partial_{\phi}p}.$
This expression coincides with the uncertainty of eqs. (4)-(5). Furthermore,
note that the FI is a generalization of $N_{\mathrm{eff}}$ defined in section
IV.1.

####  IV.2.1 Fisher information bound for the pair correlations

In pair correlation analysis, we estimate the ion phases
$\left(\phi_{i}\right)$ 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
$\left\\{\frac{1}{M}\underset{m=1}{\overset{M}{\sum}}q_{i,m}q_{j,m}\right\\}_{i\neq
j}$ and estimate the phases according to it. According to the central limit
theorem, the averages converge to a Gaussian random variable ${\cal
N}\left(\left(\mathbf{\mu}_{i,j}\right)_{i\neq j},M^{-1}\Sigma\right),$ where
$\mathbf{\mu}_{i,j}=\langle
q_{i}q_{j}\rangle=\frac{1}{2}C_{0}^{2}\cos\left(2\left(\phi_{i}-\phi_{j}\right)\right)$
and $\Sigma$ is the covariance matrix
$\Sigma_{\left(i,j\right),\left(k,m\right)}=\langle
q_{i}q_{j}q_{k}q_{m}\rangle-\langle q_{i}q_{j}\rangle\langle
q_{k}q_{m}\rangle.$ An explicit calculation of the covariance matrix elements
is presented in Appendix A.

Since the distribution is normal the FI matrix about $\left(\phi_{i}\right)$
is given by [49]:

$I=\left(\partial_{\phi}\bm{\mu}\right)^{\dagger}\Sigma^{-1}\left(\partial_{\phi}\bm{\mu}\right).$
(21)

$\left(\partial_{\phi}\bm{\mu}\right)_{i,\left(k,m\right)}=\partial_{\phi_{i}}\mu_{k,m}$
is the information gained due to the change in the mean values, i.e. the
signal, and $\Sigma$ is the covariance matrix of the different correlations
representing the noise.

Applying Eq. (21) for a single pair correlation $\left(i,j\right)$ we retrieve
the uncertainty in Eq. (5): the only linear combination of $\phi_{i},\phi_{j}$
that has a non-vanishing FI is $\phi_{i}-\phi_{j},$ for which the FI is
$\frac{\sin^{2}\left(\phi_{i}-\phi_{j}\right)}{4-\cos^{2}\left(\phi_{i}-\phi_{j}\right)},$
i.e.
$\sigma_{\Delta\phi}=\frac{\sqrt{4-\cos^{2}\left(\phi_{i}-\phi_{j}\right)}}{|\sin\left(\phi_{i}-\phi_{j}\right)|}.$
The minimal uncertainty per measurement is $2$, and a divergence occurs for
$\phi_{i}-\phi_{j}=n\pi\;\left(n\in\mathbb{Z}\right)$ due to the vanishing
derivative and non-vanishing noise.

Since information about $\phi_{i}-\phi_{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\pi.$ 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\rightarrow\infty$ 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 $\frac{4-C_{0}^{2}}{C_{0}^{2}},$
while the finite contrast precision bound to the variance is
$\frac{2}{1-\sqrt{1-C_{0}^{2}}}.$ As $C_{0}$ gets smaller the variance with
pair correlations gets closer to the this bound since the information from
higher moments becomes smaller.

####  IV.2.2 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}\rightarrow 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.

####  IV.2.3 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 $\sqrt{\frac{2}{M}}$ with an initial product state for
every $N$ and with an entangled initial state we can further reduce the
uncertainty to $\sqrt{\frac{N-1}{N}}\sqrt{\frac{2}{M}}.$ 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
$\sigma_{\Delta\phi}\geq\frac{2}{\sqrt{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 $\sqrt{\frac{2}{M}}$ is achieved.

Optimizing over both initial states and measurement strategies, we prove in
appendix C.2 that the ultimate precision limit is
$\sqrt{\frac{N-1}{N}}\sqrt{\frac{2}{M}}.$ Several initialization strategies
saturate this bound, in particular any initial pure state with $\langle
Z_{j}\rangle=0,\langle Z_{j}Z_{k}\rangle=-\frac{1}{N-1}$ for all $j\neq k$
achieves it. The reason for this improvement is the minimal $\underset{j\neq
k}{\sum}\langle Z_{j}Z_{k}\rangle,$ 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\choose 2}$, the minimal
$\underset{j<k}{\sum}\langle Z_{j}Z_{k}\rangle$ is $-N/2$ and thus the optimal
$\langle Z_{j}Z_{k}\rangle$ is $-\frac{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 anti-
parallel 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.

###  V.1 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}\leftrightarrow$D5/2
transition, in which the first (second) $\pi/2$ pulse is realized by a laser
beam impinging on the ions from the axial (perpendicular) direction. This
setting results in qubit-specific phases $\phi_{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 $\pi/2$ pulses without any free-evolution time in between.
Following the previously outlined procedure, we first reconstruct the qubit
phase $\phi_{i}$ and the measurement contrast by fitting the correlation
matrix. Next we use these phases for reconstructing the time-dependent random
phases $\varphi_{m}$. Using N-qubit correlations, we finally use the Bayesian
approach of eq. (12) for an improved phase estimate of $\phi_{i}$, shown as
open symbols in Figure 5 (a) for a string of 62 ions.

Figure 5: (a) Reconstruction of the axial trapping potential. (a) Measured
phases $\phi_{i}$ (open symbols) and nearest-neighbour distances (full
symbols) obtained by fitting a model potential. (b) Residuals. Shaded area:
theoretical minimum measurement uncertainty. (c) Center-of-mass mode frequency
measured by correlation spectroscopy, $\omega_{c}/(2\pi)$ (dark blue circles)
and by sideband spectroscopy, $\omega_{sb}/(2\pi)$, (squares) vs number of
ions, together with the fitted value of $\omega_{0}/(2\pi)$ (light blue
circles). (d) Measured $(\omega_{sb}-\omega_{c})/(2\pi)$.

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, $V(z)=\frac{1}{2}m\omega_{0}^{2}z^{2}(1+z/l_{3}+(z/l_{4})^{2})$, where
$\omega_{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
$\phi_{i}$ and find $\omega_{0}=(2\pi)~{}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 $\phi_{i}$ to the ones obtained from fitting the potential
($\phi_{i}^{fit}$) by calculating the residual position errors, $\delta
x_{i}=(\phi-\phi_{i}^{fit})/k$. Fig. 5 (b) shows that these residuals have a
standard deviation of 6.0 nm, barely above the theoretically expected error
$\sigma_{\Delta\phi}=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\leq N\leq
62$) and obtained consistent results. Fig. 5 (c) shows the inferred
oscillation frequency $\omega_{0}$ (light blue points) and the lowest
collective mode frequency $\omega_{c}$ (dark blue points). For an independent
cross-check, the latter was also measured via sideband spectroscopy on the
S1/2 to D5/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 $\phi_{i}$ over the interval from 0
to $2\pi$. 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 counter-propagating 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.

###  V.2 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 40Ca+ 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=\pm 1/2\leftrightarrow$D${}_{5/2},m=\pm 5/2$ transitions with a
Ramsey time of $\tau=40$ ms duration with a 51-ion string. In contrast to the
experiments of subsection V.1, both Ramsey pulses are realized by the same
laser beam. Writing the spatially resolved shifts $\Delta_{i}^{\pm}$ as
$\Delta_{i}^{Q}=(\Delta_{i}^{+}+\Delta_{i}^{-})/2$ and
$\Delta_{i}^{B}=(\Delta_{i}^{+}-\Delta_{i}^{-})/2$ enables a separation of
electric quadrupole and magnetic field shifts.

Figure 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=\pm 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\leftrightarrow$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$\sigma$) predicted for
quantum projection noise.

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 $\theta(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$\sigma$ error in the determination of $\theta(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: Transition frequency shift measurement in a two-dimensional crystal
obtained by probing the ground-state transition with 19736 experimental
repetitions. (a) Single-ion 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$\sigma$ measurement uncertainty predicted for
quantum projection noise.

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
ground-state transition,
$\mathrm{S}_{1/2},m=-1/2\leftrightarrow\mathrm{S}_{1/2},m=+1/2$, by a Ramsey
experiment of $\tau=5~{}\mathrm{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)~{}\mathrm{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)~{}\mathrm{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.

###  V.3 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 fluctuations, laser frequency noise and optical path
length fluctuations, respectively. Panel (a) shows a magnetic-field change of
about $3\,\mu$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\leftrightarrow$D${}_{5/2},m=3/2$ transition was probed for
$\tau=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 $\sim 8~{}\mathrm{Hz}$ introduces a frequency excursion
on the order of $1~{}\mathrm{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
$\langle(\varphi(t+\tau)-\varphi(t))^{2}\rangle_{t}$ between experiments
separated by a time $\tau$ that increase in proportion to $\tau$ 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.

Figure 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 $|\mathcal{F}[A(T)]|^{2}$
(right plot). Figure 9: Measurement of temporal relative optical path length
fluctuations in two beam path delivering the laser pulses to a string of 40
ions.

##  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 [54, 25, 31], 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 (low-frequency) 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.

###### Acknowledgements.

We acknowledge useful discussions with Alex Retzker. The project leading to
this application has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement No 741541), and from 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łodowska‐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.

##  Appendix A Bounds to the achievable phase estimation uncertainty

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

$\left\\{\frac{1}{M}\underset{m=1}{\overset{M}{\sum}}q_{i,m}q_{j,m}\right\\}_{i,j>i}\rightarrow{\cal
N}\left(\left\\{\mu_{i,j}\right\\}_{i,j>i},\Sigma\right),$ (22)

which, according to the central limit theorem, converge to a Gaussian
distribution, where $\mu_{i,j}$ is the average of $q_{i}q_{j}$ and $\Sigma$ is
the covariance matrix of the $\left\\{q_{i}q_{j}\right\\}_{i,j>i}$.

Hence the problem boils down to calculating the FI matrix for this Gaussian
distribution. The FI matrix about $\overrightarrow{\phi}$ given this Gaussian
distribution is presented in eq. (21) in the main text. We write it here as

$I=D^{\dagger}\Sigma^{-1}D,$ (23)

where:

$D=\left(\begin{array}[]{cccc}\partial_{\phi_{1}}\mu_{1,2}&\partial_{\phi_{2}}\mu_{1,2}&\cdots&\partial_{\phi_{N}}\mu_{1,2}\\\
\vdots&\vdots&\vdots&\vdots\\\
\partial_{\phi_{1}}\mu_{N-1,N}&\partial_{\phi_{2}}\mu_{N-1,N}&\cdots&\partial_{\phi_{N}}\mu_{N-1,N}\end{array}\right).$

Hence we need to calculate $\left\\{\mu_{i,j}\right\\}_{i,j>i}$ and $\Sigma$
in order to get the FI matrix. Let us first assume only correlated dephasing
(no uncorrelated dephasing). For the mean values, we have

$\displaystyle\begin{split}&\mu_{i,j}=2\left[\frac{1}{2\pi}\underset{0}{\overset{2\pi}{\int}}\cos^{2}\left(\frac{1}{2}\left(\phi_{i}+\varphi_{m}\right)\right)\cos^{2}\left(\frac{1}{2}\left(\phi_{j}+\varphi_{m}\right)\right)+\right.\\\
&\left.\sin^{2}\left(\frac{1}{2}\left(\phi_{i}+\varphi_{m}\right)\right)\sin^{2}\left(\frac{1}{2}\left(\phi_{j}+\varphi_{m}\right)\right)\;d\varphi_{m}\vphantom{\frac{1}{2\pi}\underset{0}{\overset{2\pi}{\int}}\cos^{2}\left(\phi_{i}+\phi_{r}\right)\cos^{2}\left(\phi_{j}+\phi_{r}\right)}\right]-1\\\
&=\frac{1}{2}\cos\left(\phi_{i}-\phi_{j}\right).\end{split}$ (24)

Let us now calculate $\Sigma$. The diagonal terms of $\Sigma$ read:

$\displaystyle\begin{split}&\Sigma_{\left(i,j\right),\left(i,j\right)}=\langle
q_{i}^{2}q_{j}^{2}\rangle-\langle
q_{i}q_{j}\rangle^{2}=1-\frac{1}{4}\cos^{2}\left(\phi_{i}-\phi_{j}\right)\\\
&=\frac{7}{8}-\frac{\cos\left(2\left(\phi_{i}-\phi_{j}\right)\right)}{8}\end{split}$
(25)

Regarding the non-diagonal terms, let us begin with non-overlapping pairs
$\left(i,j\right),\left(k,n\right)$:

$\displaystyle\begin{split}&\langle
q_{i}q_{j}q_{k}q_{n}\rangle=\frac{1}{8}\cos(\text{$\phi_{i}$}+\text{$\phi_{j}$}-\text{$\phi_{k}$}-\text{$\phi_{n}$})\\\
&+\frac{1}{4}\cos(\text{$\phi_{i}$}-\text{$\phi_{j}$})\cos(\text{$\phi_{k}$}-\phi_{n}),\end{split}$
$\langle q_{i}q_{j}\rangle\langle
q_{k}q_{n}\rangle=\frac{1}{4}\cos\left(\phi_{i}-\phi_{j}\right)\cos\left(\phi_{k}-\phi_{n}\right)$

Hence:

$\Sigma_{\left(i,j\right),\left(k,n\right)}=\frac{1}{8}\cos\left(\text{{\hbox{\phi_{i}}}}+\text{{\hbox{\phi_{j}}}}-\text{{\hbox{\phi_{k}}}}-\text{{\hbox{\phi_{n}}}}\right)$
(26)

For overlapping pairs, such as $\left(i,j\right),\left(i,n\right)$ we have:

$\displaystyle\begin{split}&\Sigma_{\left(i,j\right),\left(i,n\right)}=\langle
q_{i}^{2}q_{j}q_{n}\rangle-\langle q_{i}q_{j}\rangle\langle
q_{i}q_{n}\rangle=\\\ &\langle q_{j}q_{n}\rangle-\langle
q_{j}q_{i}\rangle\langle q_{i}q_{n}\rangle=\\\
&\frac{3}{8}\cos\left(\phi_{j}-\phi_{n}\right)-\frac{1}{8}\cos\left(\phi_{j}+\phi_{n}-2\phi_{i}\right)\end{split}$
(27)

The derivatives matrix, $D$, is:

$D_{\left(i,j\right),m}=\begin{cases}-\frac{1}{2}\sin\left(\phi_{i}-\phi_{j}\right)&m=i\\\
\frac{1}{2}\sin\left(\phi_{i}-\phi_{j}\right)&m=j\\\ 0&m\neq i,j\end{cases}$
(28)

Inserting equations (24),(25),(26), (27) into (23), we can perform exact
numerical calculations of the FI.

Given an uncorrelated dephasing in addition to the correlated dephasing, the
probabilities $p_{\pm}$ of observing outcomes $q_{i}=\pm 1$ are modified to
$p_{\pm}=\frac{1}{2}(1\pm C_{0}\sin\left(\phi_{i}+\varphi_{m}\right)),$ i.e. a
finite contrast of $0\leq C_{0}\leq 1$. It can be immediately observed that
$\mu_{i,j}=\frac{1}{2}C_{0}^{2}\cos\left(\phi_{i}-\phi_{j}\right).$ The
covariance matrix is modified as follows:

$\displaystyle\begin{split}&\Sigma_{\left(i,j\right),\left(i,j\right)}=\left(1-\frac{C_{0}^{4}}{8}\right)-C_{0}^{4}\frac{\cos\left(2\left(\phi_{i}-\phi_{j}\right)\right)}{8},\\\
&\Sigma_{\left(i,j\right),\left(k,n\right)}=\frac{1}{8}C_{0}^{4}\cos(\text{$\phi_{i}$}+\text{$\phi_{j}$}-\text{$\phi_{k}$}-\text{$\phi_{n}$}),\\\
&\Sigma_{\left(i,j\right),\left(i,n\right)}=\left(\frac{1}{2}C_{0}^{2}-\frac{1}{8}C_{0}^{4}\right)\cos\left(\phi_{j}-\phi_{n}\right)-\\\
&\frac{1}{8}C_{0}^{4}\cos\left(\phi_{j}+\phi_{n}-2\phi_{i}\right).\end{split}$

Let us analyze the uncertainty using pair correlations as
$N\rightarrow\infty.$ 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 $\sqrt{\frac{2}{M}},$ but
approximately to $\sqrt{\frac{3}{M}}.$

The limit of $\sqrt{\frac{3}{M}}$ can be derived analytically, assuming that
the phases are distributed evenly in $k\pi\;\left(k\in\mathbb{Z}\right).$ To
obtain this result, we use the following approximation of the variance:

$\displaystyle\begin{split}&\text{var}\left(\phi_{1}-\phi_{2}\right)=u_{1,2}^{\dagger}\left(D^{\dagger}\Sigma^{-1}D\right)^{-1}u_{1,2}\\\
&\approx\frac{||u_{1,2}||^{4}}{||Du_{1,2}||^{4}}\left(Du_{1,2}\right)^{\dagger}\Sigma\left(Du_{1,2}\right),\end{split}$

where $u_{1,2}$ is the vector that corresponds to $\phi_{2}-\phi_{1},$ i.e.
$\left(-1,1,0,...,0\right).$ This approximation is obtained by using a Cauchy-
Schwarz inequality twice:
$v^{\dagger}M^{-1}v\geq\frac{||v||^{4}}{v^{\dagger}Mv},$ and it can be
understood as a single parameter estimation bound where the derivative of the
mean is $\frac{1}{||u_{1,2}||^{2}}Du_{1,2}$ and the variance is
$\frac{1}{||Du_{1,2}||^{4}}\left(Du_{1,2}\right)^{\dagger}\Sigma\left(Du_{1,2}\right).$
We now calculate this approximation to show that in the limit of large $N$ it
converges to $3.$

Clearly $||u_{12}||^{4}=4$, and
$||Du_{1,2}||^{4}\approx\left(2\underset{k=1}{\overset{N}{\sum}}\frac{1}{4}\sin\left(2k\phi\right)^{2}\right)^{2}\approx\frac{N^{2}}{16}.$
We now need to calculate
$\left(Du_{1,2}\right)^{\dagger}\Sigma\left(Du_{1,2}\right)$: 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)
$\left(Du_{1,2}\right)_{\left(i,j\right)}\Sigma_{\left(i,j\right)\left(k,n\right)}\left(Du_{1,2}\right)_{\left(k,n\right)}.$
We can neglect the $N$ terms of identical pairs. From overlapping pairs
$\left(i,j\right),\left(i,n\right)$ the contribution is:

$\displaystyle\begin{split}&2\langle
q_{i}^{2}q_{j}q_{n}\rangle\left(Du_{1,2}\right)_{\left(i,j\right)}\left(Du_{1,2}\right)_{\left(i,n\right)}\approx\\\
&\frac{2}{8}\cos\left(\phi_{j}-\phi_{n}\right)\sin\left(\phi_{1}-\phi_{j}\right)\sin\left(\phi_{1}-\phi_{n}\right)+\\\
&\frac{2}{8}\cos\left(\phi_{j}-\phi_{n}\right)\sin\left(\phi_{2}-\phi_{j}\right)\sin\left(\phi_{2}-\phi_{n}\right)\approx\\\
&N^{2}/16.\end{split}$

The contribution from the non-overlapping pairs is:

$\displaystyle\begin{split}&2\langle
q_{i}q_{j}q_{k}q_{n}\rangle\left(Du_{1,2}\right)_{\left(i,j\right)}\left(Du_{1,2}\right)_{\left(k,n\right)}\approx\\\
&-\frac{2}{8}\cos(\text{$\phi_{1}$}+\text{$\phi_{j}$}-\text{$\phi_{2}$}-\text{$\phi_{n}$})\frac{1}{2}\sin\left(\phi_{1}-\phi_{j}\right)\frac{1}{2}\sin\left(\phi_{2}-\phi_{n}\right)\\\
&\approx-N^{2}/48\end{split}$

Hence
$\left(Du_{1,2}\right)^{\dagger}\Sigma\left(Du_{1,2}\right)\approx\frac{N^{2}}{16}\left(1-\frac{1}{4}\right)=\frac{N^{2}}{16}\frac{3}{4}.$
Therefore altogether:

$\frac{||u_{1,2}||^{4}}{||Du_{1,2}||^{4}}\left(Du_{1,2}\right){}^{\dagger}\Sigma\left(Du_{1,2}\right)\approx
3,$ (29)

which matches the numerical results.

For general contrast $C_{0}$ these expressions are modified to:

$\displaystyle\begin{split}&||Du_{1,2}||^{4}\approx\frac{N^{2}}{16}C_{0}^{8},\\\
&2\langle
q_{i}^{2}q_{j}q_{n}\rangle\left(Du_{1,2}\right)_{\left(i,j\right)}\left(Du_{1,2}\right)_{\left(i,n\right)}\approx
C_{0}^{8}N^{2}/16\\\ &2\langle
q_{i}q_{j}q_{k}q_{n}\rangle\left(Du_{1,2}\right)_{\left(i,j\right)}\left(Du_{1,2}\right)_{\left(k,n\right)}\approx-
C_{0}^{6}N^{2}/48.\end{split}$

Hence we obtain that for a general contrast
$\text{var}\left(\phi_{1}-\phi_{2}\right)$ from pair correlations converges to
$\approx\frac{4-C_{0}^{2}}{C_{0}^{2}}.$ This implies that as $C_{0}$ becomes
smaller the uncertainty using pair correlations converges to the finite
contrast bound of the variance $\frac{2}{1-\sqrt{1-C_{0}^{2}}}$. 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_{0}^{2k}$ 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_{0}^{2}$:
$\frac{2}{1-\sqrt{1-C_{0}^{2}}}=\frac{4-C_{0}^{2}+\mathcal{O}\left(C_{0}^{4}\right)}{C_{0}^{2}}.$
This raises a natural question: is the variance with all $m\leq k$ particle
correlations equal to
$\frac{2}{C_{0}^{2}}+\frac{2}{C_{0}^{2}}\left(1-\underset{l=1}{\overset{k}{\sum}}\frac{2}{l}\binom{2l-2}{l-1}\left(\frac{C_{0}^{2}}{4}\right)^{l}\right)$?
We leave it as an open question (and conjecture) as we do not have a proof to
this.

Figure 10: Sensitivity (per measurement) with pair correlations compared to
precision bounds (numerical analysis). (a) Uncertainty in estimating
$\Delta\phi$ using all pair correlations as a function of $\Delta\phi$ 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 $C_{0}.$ 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.

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:

$P(\mathbf{q})=\frac{2^{-N}}{2\pi}\int_{0}^{2\pi}d\varphi\prod_{i=1}^{N}(1+q_{i}\sin(\phi_{i}+\varphi)).$
(30)

Hence the precision bound is given by the FI matrix about $\mathbf{\phi}$ with
this distribution. Since the FI matrix involves summation over all $2^{N}$
possible $\mathbf{q}$ vectors,
$I_{i,j}=\underset{\mathbf{q}}{\sum}\frac{\left(\partial_{\phi_{i}}p\left(\mathbf{q}\right)\right)\left(\partial_{\phi_{j}}p\left(\mathbf{q}\right)\right)}{p\left(\mathbf{q}\right)}$,
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}=\langle\frac{\left(\partial_{\phi_{i}}p\right)\left(\partial_{\phi_{j}}p\right)}{p^{2}}\rangle=\langle\partial_{\phi_{i}}\ln\left(p\right)\partial_{\phi_{j}}\ln\left(p\right)\rangle.$
This allows us to make a Monte-Carlo calculation of the FI matrix by sampling
$\partial_{\phi_{i}}\ln\left(p\right)\partial_{\phi_{j}}\ln\left(p\right).$
Simulation results are shown in Fig. 11 for the case of evenly distributed
single-qubit phases, $\phi_{j}=2\pi j/N$.

Figure 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 $\sqrt{2}.$ (b) Precision bounds as a function
of $C_{0}$: 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.

##  Appendix B Numerical simulations

Estimating the phases with pair correlations: Maximum likelihood and least-
squares estimation— We calculated precision bounds of the phases given the
pair correlations; in this part we discuss estimation methods using pair
correlations and the saturability of these precision limits. We compare
between two estimation methods: simple least-squares estimation, i.e.
minimizing $\mathbf{V^{\dagger}V}$ where
$\mathbf{V}=\left\\{\frac{1}{M}\underset{m=1}{\overset{M}{\sum}}q_{i,m}q_{j,m}-\mu_{i,j}\right\\}_{i,j},$
and maximum-likelihood estimation. Note that since the relevant distribution,
eq. (22), is Gaussian, the maximum-likelihood estimation becomes a weighted
least-squares estimation [62]:

$\underset{\mathbf{\phi}}{\text{max}}\mathcal{L}\left(\mathbf{V}|\mathbf{\phi}\right)=\underset{\mathbf{\phi}}{\text{min}}\mathbf{V}^{\dagger}\Sigma^{-1}\mathbf{V}.$

The difference between the two estimation methods is thus rooted in the
weights given by the inverse of the covariance matrix, $\Sigma^{-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).

Figure 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.

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$.

##  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
$\frac{1}{\binom{N}{2}}\underset{j>k}{\sum}\text{Var}\left(\phi_{j}-\phi_{k}\right).$
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
$|\vec{z}\rangle$ is the product state
$|z_{1}\rangle|z_{2}\rangle...|z_{N}\rangle$ where $|z_{i}\rangle$ is an
eigenstate of $Z_{i}$ with eigenvalue $z_{i}\in\left\\{\pm 1\right\\}.$ A
different notation to the same state would be $|\vec{q}\rangle$ where
$q_{i}=\frac{1-z_{i}}{2},$ thus $z_{i}=1\rightarrow q_{i}=0,$
$z_{i}=-1\rightarrow q_{i}=1$.

###  C.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: $|+\rangle^{N}$, with
$|+\rangle=\frac{1}{\sqrt{2}}\left(|0\rangle+|1\rangle\right)$ (an eigenstate
of $X$ with eigenvalue $+1$). After a free evolution the state evolves into
the state of Eq. (1) with phases $\phi_{im}=\phi_{i}+\varphi_{m}.$
$\varphi_{m}$ is a random phase that is distributed uniformly in
$\left[0,2\pi\right].$ This random phase induces a correlated dephasing, i.e.
the final state, after averaging out $\varphi_{m}$, becomes a mixture of Dicke
states:

$\rho_{f}=\underset{j=0}{\overset{N}{\oplus}}\rho_{j},$

where $\rho_{j}$ is a Dicke state with $j$ excitation:

$\rho_{j}=\frac{\binom{N}{j}}{2^{N}}|\psi_{j}\rangle\langle\psi_{j}|,\;|\psi_{j}\rangle=\frac{1}{\sqrt{\binom{N}{j}}}\underset{\vec{q},\,\underset{k}{\sum}q_{k}=j}{\sum}e^{i\vec{\phi}\cdot\vec{q}}|\vec{q}\rangle.$

Given $\rho_{f}\left(\vec{\phi}\right),$ the fundamental precision limit is
set by the quantum Fisher information matrix (QFIM) about the parameters
$\phi_{1},\phi_{2},...,\phi_{N}$; this is the Fisher information matrix
optimized over all possible measurement strategies. Hence the covariance
matrix of the estimators, $\Sigma$, satisfies:

$\Sigma\geq I^{-1},$

where $I$ is the QFIM and thus for any $j,k$:

$\text{var}\left(\phi_{j}-\phi_{k}\right)\geq
u_{(j,k)}^{\dagger}I^{-1}u_{\left(j,k\right)},$

where $u_{\left(j,k\right)}$ is the parameter vector that corresponds to
$\phi_{j}-\phi_{k}$.

For a general mixed state $\rho,$ given its spectral decomposition
$\rho=\underset{k}{\sum}p_{k}|k\rangle\langle k|,$ the QFIM is given by [50]:

$I_{i,j}=2\underset{k,l}{\sum}\frac{\left(\frac{\partial\rho}{\partial\phi_{i}}\right)_{k,l}\left(\frac{\partial\rho}{\partial\phi_{j}}\right)_{l,k}}{\left(p_{l}+p_{k}\right)},$

where $\left\\{p_{k}\right\\}_{k}$ are the eigenvalues of $\rho$ and the
matrix elements,$\left(\bullet\right)_{k,l}=\langle k|\bullet|l\rangle$, are
with respect to the eigenbasis of $\rho.$

For pure states this expression is reduced to:

$I_{i,j}=4\left(\langle\partial_{\phi_{i}}\psi|\partial_{\phi_{j}}\psi\rangle-\langle\partial_{\phi_{i}}\psi|\psi\rangle\langle\psi|\partial_{\phi_{j}}\psi\rangle\right).$
(31)

Let us calculate the QFIM of our $\rho_{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
$|\psi_{j}\rangle$ [63]:

$I=\underset{j=0}{\overset{N}{\sum}}\frac{\binom{N}{j}}{2^{N}}I^{\left(j\right)},$
(32)

where $I^{\left(j\right)}$ is the QFIM of $|\psi_{j}\rangle.$ For every
$|\psi_{j}\rangle$ we have
$\partial_{\phi_{i}}|\psi_{j}\rangle=-i\frac{1}{2}\left(Z_{i}+\mathds{I}\right)|\psi_{j}\rangle,$
inserting this into equation (31) we get that the QFIM of each
$|\psi_{j}\rangle$ is:

$I_{k,l}^{\left(j\right)}=\left(\langle\psi_{j}|Z_{k}Z_{l}|\psi_{j}\rangle-\langle\psi_{j}|Z_{k}|\psi_{j}\rangle\langle\psi_{j}|Z_{l}|\psi_{j}\rangle\right).$

Now $|\psi_{j}\rangle$ is a symmetric superposition of all states with $j$
excitations, from symmetry we get:

$\langle Z_{k}\rangle=\frac{1}{N}\left(N-j-j\right)=\frac{N-2j}{N}$

and for $k\neq l$:

$\langle Z_{k}Z_{l}\rangle=\frac{\left(N-2j\right)^{2}-N}{N\left(N-1\right)}$

Hence all the non-diagonal terms of $I^{\left(j\right)}$ are:

$\displaystyle I_{k,l}^{\left(j\right)}$
$\displaystyle=\frac{\left(N-2j\right)^{2}-N}{N\left(N-1\right)}-\frac{\left(N-2j\right)^{2}}{N^{2}}$
$\displaystyle=\frac{4j\left(j-N\right)}{N^{2}\left(N-1\right)}.$

The diagonal terms of $I^{\left(j\right)}$ read:

$I_{k,k}^{\left(j\right)}=1-\left(\frac{N-2j}{N}\right)^{2}=\frac{4j\left(N-j\right)}{N^{2}}$

Inserting these terms into equation (32), we get that $I$ reads:

$I=\begin{cases}\underset{j=0}{\overset{N}{\sum}}\frac{\binom{N}{j}}{2^{N}}\frac{4j\left(j-N\right)}{N^{2}\left(N-1\right)}=-\frac{1}{N}&k\neq
l\\\
\underset{j=0}{\overset{N}{\sum}}\frac{\binom{N}{j}}{2^{N}}\frac{4j\left(N-j\right)}{N^{2}}=\frac{N-1}{N}.&k=l\end{cases}$

It is now simple to see that for any $k\neq m$ the vector that corresponds to
$\phi_{k}-\phi_{m}$ is an eigenvector of $I$ with an eigenvalue of $1.$ The
variance per measurement is thus

$\text{var}\left(\phi_{k}-\phi_{m}\right)=2$

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

$\rho_{f}=\frac{1}{4}\left(|11\rangle\langle 11|+|00\rangle\langle
00|\right)+\frac{1}{2}|\psi_{1}\rangle\langle\psi_{1}|,$

with
$|\psi_{1}\rangle=\frac{1}{\sqrt{2}}\left(|01\rangle+e^{i\left(\phi_{1}-\phi_{2}\right)}|10\rangle\right)$.
Measuring the local $X$ basis, we project onto the states
$|\vec{x}\rangle=|x_{1}\rangle|x_{2}\rangle...|x_{N}\rangle$ where
$|x_{i}\rangle$ is an eigenstate of $X_{i}$ with eigenvalues
$x_{i}\in\left\\{\pm 1\right\\}$ and obtain the probabilities:

$\displaystyle\begin{split}&\text{even
}\vec{x}:\frac{1}{8}+\frac{1}{4}\cos\left(\frac{\phi_{1}-\phi_{2}}{2}\right)^{2},\\\
&\text{odd
}\vec{x}:\frac{1}{8}+\frac{1}{4}\sin\left(\frac{\phi_{1}-\phi_{2}}{2}\right)^{2},\end{split}$

where $\vec{x}$ odd (even) stands for $\\#\left(x_{i}=-1\right)$ odd (even).
This leads to:

$\text{var}\left(\phi_{2}-\phi_{1}\right)=\frac{4-\cos\left(\phi_{1}-\phi_{2}\right)^{2}}{\sin\left(\phi_{1}-\phi_{2}\right)^{2}}\geq
4,$

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 $\phi_{1},\phi_{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\rangle,|11\rangle$ in the first measurement and thus no information
and with probability $1/2$ we collapse into $|\psi_{1}\rangle$ which yields a
Fisher information of $1$. Therefore the total Fisher information is:

$\frac{1}{2}\cdot 0+\frac{1}{2}\cdot
1=\frac{1}{2}\rightarrow\text{var}\left(\phi_{2}-\phi_{1}\right)=2,$

hence the bound is saturated.

A general optimal measurement strategy would be:
1\. first measure $\underset{i}{\sum}Z_{i}$ (this measurement collapses the
density matrix into one of the Dicke states, $|\psi_{j}\rangle$).
2\. Measure the Dicke state, $|\psi_{j}\rangle$, in its optimal measurement
basis.

The optimal measurement basis of $|\psi_{j}\rangle$ can be written implicitly
as proven in [64]: projecting into a (Gram-Schmidt) orthogonalization of
$\left\\{|\psi_{j}\rangle,|\partial_{\phi_{k}}\psi_{j}\rangle\right\\}_{\phi_{k}}$
would be optimal. For example, for $N=3$, given
$|\psi_{1}\rangle=\frac{1}{\sqrt{3}}\left(|011\rangle+|101\rangle+|110\rangle\right)$
(we can assume for convenience $\phi_{1}=\phi_{2}=\phi_{3}=0$, this can be
achieved adaptively by local operations), an optimal measurement basis would
be:

$\begin{split}&\frac{1}{\sqrt{3}}\left(|011\rangle+|101\rangle+|110\rangle\right),\frac{1}{\sqrt{6}}(-2|011\rangle+|101\rangle+|110\rangle),\\\
&\frac{1}{\sqrt{2}}\left(|101\rangle-|110\rangle\right).\end{split}$

The construction for $|\psi_{2}\rangle$ is equivalent.

###  C.2 Optimal initial states

We show that the average variance of phase difference
$\frac{1}{\binom{N}{2}}\underset{j>k}{\sum}\text{Var}\left(\phi_{j}-\phi_{k}\right)$
is lower-bounded by $2\frac{N-1}{N}$ and we find several initialization
strategies, all involve entanglement, that saturate this bound. In particular
the symmetric Dicke state:

$|\psi\rangle=\frac{1}{\sqrt{\binom{N}{N/2}}}\underset{\vec{z}\;\text{with
}\underset{i}{\sum}z_{i}=0}{\sum}|\vec{z}\rangle,$

and any other state that is an equal superposition of states with
$\underset{i}{\sum}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 $\underset{i}{\sum}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].

We denote the final and initial state as $|\psi_{f}\rangle,|\psi\rangle$
respectively, where $|\psi_{f}\rangle=U|\psi\rangle,$ with
$U=\exp\left(-i\frac{1}{2}\vec{\phi}\cdot\vec{Z}\right).$ Since
$|\partial_{\phi_{j}}\psi_{f}\rangle=-\frac{i}{2}Z_{j}U|\psi\rangle$, the QFIM
of $|\psi_{f}\rangle$ reads:

$I_{i,j}=\left(\langle\psi|Z_{i}Z_{j}|\psi\rangle-\langle\psi|Z_{i}|\psi\rangle\langle\psi|Z_{j}|\psi\rangle\right).$

Using the QFIM we prove that the optimal achievable variance,
$\frac{1}{\binom{N}{2}}\underset{j>k}{\sum}\text{Var}\left(\phi_{j}-\phi_{k}\right)$
, is $2\frac{N-1}{N}$.
Proof: By definition of QFIM

$\text{var}\left(\phi_{j}-\phi_{k}\right)\geq
u_{(j,k)}^{\dagger}I^{-1}u_{\left(j,k\right)}.$

The Cauchy-Schwarz inequality
$\left(u^{\dagger}I^{-1}u\right)\left(u^{\dagger}Iu\right)\geq|u^{\dagger}\sqrt{I^{-1}}\sqrt{I}u|^{2}$
implies:

$\left(u_{\left(j,k\right)}^{\dagger}I^{-1}u_{j,k}\right)\geq\frac{||u_{\left(j,k\right)}||^{4}}{\left(u_{\left(j,k\right)}^{\dagger}Iu_{\left(j,k\right)}\right)}$

Hence:

$\displaystyle\begin{split}&\text{var}\left(\phi_{j}-\phi_{k}\right)\geq\frac{||u_{\left(j,k\right)}||^{4}}{\left(u_{\left(j,k\right)}^{\dagger}Iu_{\left(j,k\right)}\right)}=\frac{4}{I_{jj}+I_{kk}-2I_{jk}}\\\
&\geq\frac{2}{\left(1-\langle Z_{j}Z_{k}\rangle\right)}.\end{split}$

Therefore we seek to lower bound $\underset{j<k}{\sum}\frac{2}{\left(1-\langle
Z_{j}Z_{k}\rangle\right)}$. The minimal possible
$\text{var}\left(\phi_{j}-\phi_{k}\right)$ is therefore obtained when $\langle
Z_{j}Z_{k}\rangle=-1$, however there is no state that satisfies $\langle
Z_{j}Z_{k}\rangle=-1$ for all $j,k$. To lower bound
$\underset{j<k}{\sum}\frac{2}{\left(1-\langle Z_{j}Z_{k}\rangle\right)}$ we
use the Cauchy-Schwarz inequality:

$\displaystyle\begin{split}&2\left(\underset{j<k}{\sum}\frac{1}{\left(1-\langle
Z_{j}Z_{k}\rangle\right)}\right)\left(\underset{j<k}{\sum}\left(1-\langle
Z_{j}Z_{k}\rangle\right)\right)\\\ &\geq
2\left(\underset{j<k}{\sum}1\right)^{2}=2\underset{j<k}{\sum}{\binom{N}{2}}\end{split}$

The first inequality is due to
$\left(\underset{i}{\sum}\frac{1}{x_{i}}\right)\left(\underset{i}{\sum}x_{i}\right)\geq\left(\underset{i}{\sum}\frac{1}{\sqrt{x_{i}}}\cdot\sqrt{x_{i}}\right)^{2}=\left(\underset{i}{\sum}1\right)^{2}$
, which is just the Cauchy-Schwarz inequality.
Hence:

$\underset{j<k}{\sum}\text{Var}\left(\phi_{j}-\phi_{k}\right)\geq
2\underset{j<k}{\sum}\frac{\binom{N}{2}}{\underset{j<k}{\sum}\left(1-\langle
Z_{j}Z_{k}\rangle\right)}.$

Note that
$\langle\underset{i}{\sum}Z_{i}\rangle^{2}=2\underset{j<k}{\sum}\langle
Z_{j}Z_{k}\rangle+N$, therefore $2\underset{j<k}{\sum}\langle
Z_{j}Z_{k}\rangle\geq-N$. Hence:

$\displaystyle\begin{split}&\underset{j<k}{\sum}\text{Var}\left(\phi_{j}-\phi_{k}\right)\geq\underset{j<k}{\sum}\frac{4\binom{N}{2}}{2\binom{N}{2}+N}\\\
&=2\underset{j<k}{\sum}\frac{N-1}{N}.\end{split}$

This basically proves that
$\frac{1}{\binom{N}{2}}\underset{j<k}{\sum}\text{Var}\left(\phi_{j}-\phi_{k}\right)\geq
2\frac{N-1}{N}$. To show that this lower bound is saturable we need to find an
initial state $|\psi\rangle$, for which all these inequalities are saturated,
namely:

$\underset{j<k}{\sum}u_{(j,k)}^{\dagger}I^{-1}u_{\left(j,k\right)}=2\underset{j<k}{\sum}\frac{N-1}{N},$

where $I$ is the QFIM given this $|\psi\rangle$. We observed that a necessary
condition is $\langle Z_{i}\rangle=0$ and identical $\langle
Z_{j}Z_{k}\rangle=-\frac{1}{N-1}$. Let us show that this is also a sufficient
condition: given that this condition is satisfied the QFIM is

$I_{i,j}=\begin{cases}1&i=j\\\ -\frac{1}{N-1}&i\neq j\end{cases}.$

It can be now observed that any $u_{\left(j,k\right)}$ is an eigenvector of
this matrix with eigenvalue $\frac{N}{N-1}$ and thus for any $j,k$ :
$u_{(j,k)}^{\dagger}I^{-1}u_{\left(j,k\right)}=2\frac{N-1}{N}.$

Hence any initial pure state that satisfies the conditions:

$\displaystyle\forall i\;\langle Z_{i}\rangle=0\;\text{and}\;\forall
j,k\;\langle Z_{j}Z_{k}\rangle=-\frac{1}{N-1},$ (33)

saturates this QFIM. We can immediately observe that the symmetric Dicke
state:

$|\psi\rangle=\frac{1}{\sqrt{\binom{N}{N/2}}}\underset{\vec{z}\;\text{with
}\underset{i}{\sum}z_{i}=0}{\sum}|\vec{z}\rangle$

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.
To show that indeed $\text{Var}\left(\phi_{j}-\phi_{k}\right)=2\frac{N-1}{N}$
can be achieved with $|\psi\rangle$, we need to find a read-out strategy that
achieves this bound, i.e. a measurement with a classical FI matrix that equals
the QFIM. We show that local measurements in $X$ saturate this optimal
variance.

To show this let us first write $|\psi_{f}\rangle,$ the final probe state
given the initial symmetric Dicke state:

$\displaystyle|\psi_{f}\rangle=\frac{1}{\sqrt{\binom{N}{N/2}}}\underset{\underset{i}{\sum}z_{i}=0,z_{1}=1}{\sum}\left(e^{-i\frac{1}{2}\vec{\phi}\cdot\vec{z}}|\vec{z}\rangle+e^{i\frac{1}{2}\vec{\phi}\cdot\vec{z}}|-\vec{z}\rangle\right)=$
(34)
$\displaystyle\frac{\sqrt{2}}{\sqrt{\binom{N}{N/2}}}\underset{\underset{i}{\sum}z_{i}=0,z_{1}=1}{\sum}\cos\left(\frac{1}{2}\vec{\phi}\cdot\vec{z}\right)|+_{\vec{z}}\rangle-i\sin\left(\frac{1}{2}\vec{\phi}\cdot\vec{z}\right)|-_{\vec{z}}\rangle,$
(35)

where
$|\pm_{\vec{z}}\rangle=\frac{1}{\sqrt{2}}\left(|\vec{z}\rangle\pm|-\vec{z}\rangle\right).$

Let us now use theorem $2$ of reference [64]: Given a pure probe state
$|\Psi\left(\vec{\phi}\right)\rangle,$ then a projective measurement that
consists of rank $1$ projectors $\left\\{\Pi_{k}\right\\}_{k}$ saturates the
QFIM if and only if for every $k$, $j$:

$\text{Im}\left(\langle\Psi|\Pi_{k}|\partial_{\phi_{j}}\Psi_{\perp}\rangle\right)=0,$
(36)

where
$|\partial_{\phi_{j}}\Psi_{\perp}\rangle:=|\partial_{\phi_{j}}\Psi\rangle-|\Psi\rangle\langle\Psi|\partial_{\phi_{j}}\Psi\rangle,$
i.e. it is the projection of $|\partial_{\phi_{j}}\Psi\rangle$ onto the
orthogonal subspace of $|\Psi\rangle.$ The full proof of this theorem is
presented in ref. [64]. Let us briefly explain the intuition behind this
theorem: given $|\Psi\rangle$ the probability of detecting the $k$-th result
is $p_{k}=\langle\Psi|\Pi_{k}|\Psi\rangle.$ The derivative of this probability
with respect to $\phi_{j}$ is
$\partial_{\phi_{j}}p_{k}=2\text{Re}\langle\Psi|\Pi_{k}|\partial_{\phi_{j}}\Psi\rangle.$
The parallel part of $|\partial_{\phi_{j}}\Psi\rangle,$ i.e.
$|\partial_{\phi_{j}}\Psi_{\parallel}\rangle:=\langle\Psi|\partial_{\phi_{j}}\Psi\rangle|\Psi\rangle,$
does not contribute to the derivative because
$\text{Re}\langle\Psi|\Pi_{k}|\partial_{\phi_{j}}\Psi_{\parallel}\rangle=0.$
The derivative can therefore be written as
$\partial_{\phi_{j}}p_{k}=2\text{Re}\langle\Psi|\Pi_{k}|\partial_{\phi_{j}}\Psi_{\perp}\rangle.$
Hence if
$\text{Im}\langle\Psi|\Pi_{k}|\partial_{\phi_{j}}\Psi_{\perp}\rangle\neq 0$
then some of the information about $\phi_{j}$ is being lost when measuring in
this basis, i.e. a change $\phi_{j}$ is being translated to a change in the
phase and not the probability. If
$\text{Im}\langle\Psi|\Pi_{k}|\partial_{\phi_{j}}\Psi_{\perp}\rangle=0$ for
every $j,k$ then no information about $\vec{\phi}$ is lost and thus the QFIM
is being saturated. We remark that this intuitive argument is correct only for
pure states.

Let us apply this theorem to our case: we need to show that the condition of
eq. (36) is satisfied for our $|\psi_{f}\rangle$ and local $X$ measurements.
The rank $1$ projectors in our case are thus
$\left\\{\Pi_{\vec{x}}=|\vec{x}\rangle\langle\vec{x}|\right\\}_{\vec{x}}$.
Hence we need to show that for every $\vec{x}$:
$\text{Im}\left(\langle\psi_{f}|\Pi_{\vec{x}}|\partial_{\phi_{j}}\psi_{f\perp}\rangle\right)=0.$

Let us use the following identity:

$\langle\vec{x}|\frac{1}{\sqrt{2}}\left(|\vec{z}\rangle+|-\vec{z}\rangle\right)=\frac{1}{\sqrt{2^{N+1}}}\left[\left(-1\right)^{\vec{q}_{\vec{x}}\cdot\vec{q}_{\vec{z}}}+\left(-1\right)^{\vec{q}_{\vec{x}}\cdot\vec{q}_{-\vec{z}}}\right],$

with $\left(\vec{q}_{\vec{z}}\right)_{i}=\frac{1}{2}\left(1-z_{i}\right)$ and
$\vec{q}_{\vec{x}}$ and analogously
$\left(\vec{q}_{\vec{x}}\right)_{i}=\frac{1}{2}\left(1-x_{i}\right)$. Note
that
$\vec{q}_{\vec{x}}\cdot\vec{q}_{\vec{z}}+\vec{q}_{\vec{x}}\cdot\vec{q}_{-\vec{z}}=\\#\left(x_{i}=-1\right)$,
hence if $\\#\left(x_{i}=-1\right)$ is even then
$\left(-1\right)^{\vec{q}_{\vec{x}}\cdot\vec{q}_{\vec{z}}}=\left(-1\right)^{\vec{q}_{\vec{x}}\cdot\vec{q}_{-\vec{z}}}$,
and if it is odd
$\left(-1\right)^{\vec{q}_{\vec{x}}\cdot\vec{q}_{\vec{z}}}=-\left(-1\right)^{\vec{q}_{\vec{x}}\cdot\vec{q}_{-\vec{z}}}$.
Therefore:

$\langle\vec{x}|\frac{1}{\sqrt{2}}\left(|\vec{z}\rangle+|-\vec{z}\rangle\right)=\begin{cases}0&\vec{x}\;\text{odd}\\\
\pm\frac{1}{\sqrt{2^{N}}}&\vec{x}\;\text{even},\end{cases}$ (37)

and:

$\langle\vec{x}|\frac{1}{\sqrt{2}}\left(|\vec{z}\rangle-|-\vec{z}\rangle\right)=\begin{cases}\pm\frac{1}{\sqrt{2^{N}}}&\vec{x}\;\text{odd}\\\
0&\vec{x}\;\text{even},\end{cases}$ (38)

Inserting equations (37)-(38) into equation (35) we can observe that:

$\displaystyle\langle\vec{x}|\psi_{f}\rangle=\frac{1}{\sqrt{2^{N-1}{\binom{N}{N/2}}}}\cdot$
$\displaystyle\begin{cases}\left(-i\right)\underset{\underset{i}{\sum}z_{i}=0,z_{1}=1}{\sum}\left(-1\right)^{\underset{i}{\sum}\left(\vec{q}_{\vec{x}}\right)_{i}\left(\vec{q}_{\vec{z}}\right)_{i}}\sin\left(\frac{1}{2}\vec{\phi}\cdot\vec{z}\right)&\vec{x}\;\text{odd}\\\
\underset{\underset{i}{\sum}z_{i}=0,z_{1}=1}{\sum}\left(-1\right)^{\underset{i}{\sum}\left(\vec{q}_{\vec{x}}\right)_{i}\left(\vec{q}_{\vec{z}}\right)_{i}}\cos\left(\frac{1}{2}\vec{\phi}\cdot\vec{z}\right)&\vec{x}\;\text{even}.\end{cases}$

Therefore for any value of $\vec{\phi},$ $\langle\vec{x}|\psi_{f}\rangle$ is
either real (for even $\vec{x}$) or imaginary (for odd $\vec{x}$). Similarly
it is simple to observe that
$\langle\vec{x}|\partial_{\phi_{j}}\psi_{f\perp}\rangle$ is real (imaginary)
for even $\vec{x}$ (odd $\vec{x}$). To sum up:

$\displaystyle\begin{cases}\vec{x}\;\text{odd}:\\\
\langle\vec{x}|\psi_{f}\rangle,\langle\vec{x}|\partial_{\phi_{j}}\psi_{f\perp}\rangle\;\text{imaginary}\Rightarrow\langle\psi_{f}|\Pi_{\vec{x}}|\partial_{\phi_{j}}\psi_{f\perp}\rangle\;\text{real}\\\
\vec{x}\;\text{even}:\\\
\langle\vec{x}|\psi_{f}\rangle,\langle\vec{x}|\partial_{\phi_{j}}\psi_{f\perp}\rangle\;\text{real}\Rightarrow\langle\psi_{f}|\Pi_{\vec{x}}|\partial_{\phi_{j}}\psi_{f\perp}\rangle\;\text{real}\;\;.\end{cases}$

Hence the condition in eq. (36) is satisfied, and thus the local $X$ basis
indeed saturates the QFIM.

We remark that another strategy to saturate the QFIM is to choose the initial
state of each experiment from a classical ensemble of products of anti-
parallel Bell states. An anti-parallel Bell state is defined as
$|i,j\rangle=|0\rangle_{i}|1\rangle_{j}+|0\rangle_{j}|1\rangle_{i}.$ A Bell-
product state is then a product of N/2 such anti-parallel Bell pairs, we
denote any such state as
$\underset{k=1}{\overset{N/2}{\prod}}|i_{k},j_{k}\rangle.$ It can be observed
that the QFIM of any $\underset{k=1}{\overset{N/2}{\prod}}|i_{k},j_{k}\rangle$
is $\mathds{1}-\underset{k}{\prod}X_{i_{k},j_{k}}$ , where
$X_{i_{k},j_{k}}=|i_{k}\rangle\langle j_{k}|+h.c.$. The total number of these
states is $\left(N-1\right)\left(N-3\right)...1$, the total number of states
with a specific pair $\left(i_{k},j_{k}\right)$ is $\left(N-3\right)...1$.
Hence by sampling from a uniform distribution of these Bell-product states the
QFIM becomes the optimal one: $I_{i,j}=\begin{cases}-\frac{1}{N-1}&i\neq j\\\
1&i=j\end{cases}.$ This QFIM is saturated with local $X$ (or $Y$) measurements
since these measurements saturate the QFIM of each Bell-product state
individually.

Finally, we remark about the set of optimal initial pure states. These are the
states that satisfy the conditions in equation 33. The problem of finding
optimal states (other than the symmetric Dicke state) is then basically
solving a system of linear equations. First observe that the conditions imply
$\langle(\underset{i}{\sum}Z_{i})^{2}\rangle=0,$ hence any optimal pure state
is an eigenstate of $\underset{i}{\sum}Z_{i}$ with eigenvalue $0$. We can
therefore write the states as
$\underset{\underset{i}{\sum}z_{i}=0}{\sum}\sqrt{p}_{\vec{z}}|\vec{z}\rangle.$
The conditions then become a system of ${N\choose 2}+N+1$ linear equations for
the distribution $\left\\{p_{\vec{z}}\right\\}$ : $\forall
j,k,i\underset{\underset{i}{\sum}z_{i}=0}{\sum}p_{\vec{z}}z_{j}z_{k}=-\frac{1}{N-1},\underset{\underset{i}{\sum}z_{i}=0}{\sum}p_{\vec{z}}z_{i}=0$
and $\underset{l}{\sum}p_{l}=1,$ with a constraint of $0\leq p_{l}\leq 1$ for
all $l$. The number of equations, ${N\choose 2}+N+1,$ is smaller than the
number of variables, ${N\choose 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.

## References

  * Ludlow  _et al._ [2015] A. D. Ludlow, M. M. Boyd, J. Ye, E. Peik, and P. O. Schmidt, Optical atomic clocks, Rev. Mod. Phys. 87, 637 (2015).
  * Canuel  _et al._ [2006] B. Canuel, F. Leduc, D. Holleville, A. Gauguet, J. Fils, A. Virdis, A. Clairon, N. Dimarcq, C. J. Bordé, A. Landragin, and P. Bouyer, Six-axis inertial sensor using cold-atom interferometry, Phys. Rev. Lett. 97, 010402 (2006).
  * Peters  _et al._ [1999] A. Peters, K. Y. Chung, and S. Chu, Measurement of gravitational acceleration by dropping atoms, Nature 400, 849 (1999).
  * Hamilton  _et al._ [2015] P. Hamilton, M. Jaffe, P. Haslinger, Q. Simmons, H. Müller, and J. Khoury, Atom-interferometry constraints on dark energy, Science 349, 849 (2015).
  * Bouchendira  _et al._ [2011] R. Bouchendira, P. Cladé, S. Guellati-Khélifa, F. m. c. Nez, and F. m. c. Biraben, New determination of the fine structure constant and test of the quantum electrodynamics, Phys. Rev. Lett. 106, 080801 (2011).
  * Giovannetti  _et al._ [2011] V. Giovannetti, S. Lloyd, and L. Maccone, Advances in quantum metrology, Nat. Photonics 5, 222 (2011).
  * Wineland  _et al._ [1992] D. J. Wineland, J. J. Bollinger, W. M. Itano, F. L. Moore, and D. J. Heinzen, Spin squeezing and reduced quantum noise in spectroscopy, Phys. Rev. A 46, R6797 (1992).
  * Bollinger  _et al._ [1996] J. J. . Bollinger, W. M. Itano, D. J. Wineland, and D. J. Heinzen, Optimal frequency measurements with maximally correlated states, Phys. Rev. A 54, R4649 (1996).
  * Marciniak  _et al._ [2022] C. D. Marciniak, T. Feldker, I. Pogorelov, R. Kaubruegger, D. V. Vasilyev, R. van Bijnen, P. Schindler, P. Zoller, R. Blatt, and T. Monz, Optimal metrology with programmable quantum sensors, Nature 603, 604 (2022).
  * Huelga  _et al._ [1997] S. F. Huelga, C. Macchiavello, T. Pellizzari, A. K. Ekert, M. B. Plenio, and J. I. Cirac, Improvement of frequency standards with quantum entanglement, Phys. Rev. Lett. 79, 3865 (1997).
  * Demkowicz-Dobrzański  _et al._ [2012] R. Demkowicz-Dobrzański, J. Kołodyński, and M. Guţă, The elusive Heisenberg limit in quantum-enhanced metrology, Nat. Commun. 3, 1063 (2012).
  * Chin  _et al._ [2012] A. W. Chin, S. F. Huelga, and M. B. Plenio, Quantum metrology in non-markovian environments, Phys. Rev. Lett. 109, 233601 (2012).
  * Jeske  _et al._ [2014] J. Jeske, J. H. Cole, and S. F. Huelga, Quantum metrology subject to spatially correlated Markovian noise: restoring the Heisenberg limit, New. J. Phys. 16, 073039 (2014).
  * Braun  _et al._ [2018] D. Braun, G. Adesso, F. Benatti, R. Floreanini, U. Marzolino, M. W. Mitchell, and S. Pirandola, Quantum-enhanced measurements without entanglement, Rev. Mod. Phys. 90, 035006 (2018).
  * Chwalla  _et al._ [2007] M. Chwalla, K. Kim, T. Monz, P. Schindler, M. Riebe, C. F. Roos, and R. Blatt, Precision spectroscopy with two correlated atoms, Appl. Phys. B 89, 483 (2007).
  * Olmschenk  _et al._ [2007] S. Olmschenk, K. C. Younge, D. L. Moehring, D. N. Matsukevich, P. Maunz, and C. Monroe, Manipulation and detection of a trapped Yb+ hyperfine qubit, Phys. Rev. A 76, 052314 (2007).
  * Boulder Atomic Clock Optical Network Collaboration  _et al._(2021)Boulder Atomic Clock Optical Network (BACON) Collaboration, Beloy, and Bodine et al. [BACON] Boulder Atomic Clock Optical Network (BACON) Collaboration, K. Beloy, and M. Bodine et al., Frequency ratio measurements at 18-digit accuracy using an optical clock network, Nature 591, 564 (2021).
  * Takamoto  _et al._ [2011] M. Takamoto, T. Takano, and H. Katori, Frequency comparison of optical lattice clocks beyond the Dick limit, Nat. Photonics 5, 288 (2011).
  * Nicholson  _et al._ [2012] T. L. Nicholson, M. J. Martin, J. R. Williams, B. J. Bloom, M. Bishof, M. D. Swallows, S. L. Campbell, and J. Ye, Comparison of two independent sr optical clocks with $1\mathbf{\times{}}{10}^{-17}$ stability at ${10}^{3}\text{ }\text{ }\mathbf{s}$, Phys. Rev. Lett. 109, 230801 (2012).
  * Lanyon  _et al._ [2013] B. P. Lanyon, P. Jurcevic, C. Hempel, M. Gessner, V. Vedral, R. Blatt, and C. F. Roos, Experimental generation of quantum discord via noisy processes, Phys. Rev. Lett. 111, 100504 (2013).
  * Chou  _et al._ [2011] C. W. Chou, D. B. Hume, M. J. Thorpe, D. J. Wineland, and T. Rosenband, Quantum coherence between two atoms beyond Q$={10}^{15}$, Phys. Rev. Lett. 106, 160801 (2011).
  * Marti  _et al._ [2018] G. E. Marti, R. B. Hutson, A. Goban, S. L. Campbell, N. Poli, and J. Ye, Imaging optical frequencies with $100\text{ }\text{ }\mu\mathrm{Hz}$ precision and $1.1\text{ }\text{ }\mu\mathrm{m}$ resolution, Phys. Rev. Lett. 120, 103201 (2018).
  * Shaniv  _et al._ [2019] R. Shaniv, N. Akerman, T. Manovitz, Y. Shapira, and R. Ozeri, Quadrupole shift cancellation using dynamic decoupling, Phys. Rev. Lett. 122, 223204 (2019).
  * Clements  _et al._ [2020] E. R. Clements, M. E. Kim, K. Cui, A. M. Hankin, S. M. Brewer, J. Valencia, J.-S. Chen, C.-W. Chou, D. R. Leibrandt, and D. B. Hume, Lifetime-limited interrogation of two independent 27Al+ clocks using correlation spectroscopy, Phys. Rev. Lett. 125, 243602 (2020).
  * Young  _et al._ [2020] A. W. Young, W. J. Eckner, W. R. Milner, D. Kedar, M. A. Norcia, E. Oelker, N. Schine, J. Ye, and A. M. Kaufman, Half-minute-scale atomic coherence and high relative stability in a tweezer clock, Nature 588, 408 (2020).
  * Roos  _et al._ [2006] C. F. Roos, M. Chwalla, K. Kim, M. Riebe, and R. Blatt, ’Designer atoms’ for quantum metrology, Nature 443, 316 (2006).
  * Megidish  _et al._ [2019] E. Megidish, J. Broz, N. Greene, and H. Häffner, Improved test of local Lorentz invariance from a deterministic preparation of entangled states, Phys. Rev. Lett. 122, 123605 (2019).
  * Manovitz  _et al._ [2019] T. Manovitz, R. Shaniv, Y. Shapira, R. Ozeri, and N. Akerman, Precision measurement of atomic isotope shifts using a two-isotope entangled state, Phys. Rev. Lett. 123, 203001 (2019).
  * Monz  _et al._ [2011] T. Monz, P. Schindler, J. T. Barreiro, M. Chwalla, D. Nigg, W. A. Coish, M. Harlander, W. Hänsel, M. Hennrich, and R. Blatt, 14-qubit entanglement: creation and coherence, Phys. Rev. Lett. 106, 130506 (2011).
  * Bradley  _et al._ [2019] C. E. Bradley, J. Randall, M. H. Abobeih, R. C. Berrevoets, M. J. Degen, M. A. Bakker, M. Markham, D. J. Twitchen, and T. H. Taminiau, A ten-qubit solid-state spin register with quantum memory up to one minute, Phys. Rev. X 9, 031045 (2019).
  * Singh  _et al._ [2022] K. Singh, C. E. Bradley, S. Anand, V. Ramesh, R. White, and H. Bernien, Mid-circuit correction of correlated phase errors using an array of spectator qubits, arXiv:2208.11716 10.48550/ARXIV.2208.11716 (2022).
  * Braverman  _et al._ [2018] B. Braverman, A. Kawasaki, and V. Vuletić, Impact of non-unitary spin squeezing on atomic clock performance, New J. Phys. 20, 103019 (2018).
  * Sawyer  _et al._ [2012] B. C. Sawyer, J. W. Britton, A. C. Keith, C.-C. J. Wang, J. K. Freericks, H. Uys, M. J. Biercuk, and J. J. Bollinger, Spectroscopy and thermometry of drumhead modes in a mesoscopic trapped-ion crystal using entanglement, Phys. Rev. Lett. 108, 213003 (2012).
  * von Lüpke  _et al._ [2020] U. von Lüpke, F. Beaudoin, L. M. Norris, Y. Sung, R. Winik, J. Y. Qiu, M. Kjaergaard, D. Kim, J. Yoder, S. Gustavsson, L. Viola, and W. D. Oliver, Two-qubit spectroscopy of spatiotemporally correlated quantum noise in superconducting qubits, PRX Quantum 1, 010305 (2020).
  * Kawamura and Chen [2004] S. Kawamura and Y. Chen, Displacement-noise-free gravitational-wave detection, Phys. Rev. Lett. 93, 211103 (2004).
  * Chen  _et al._ [2006] Y. Chen, A. Pai, K. Somiya, S. Kawamura, S. Sato, K. Kokeyama, R. L. Ward, K. Goda, and E. E. Mikhailov, Interferometers for displacement-noise-free gravitational-wave detection, Phys. Rev. Lett. 97, 151103 (2006).
  * Gefen  _et al._ [2022] T. Gefen, R. Tarafder, R. X. Adhikari, and Y. Chen, Quantum precision limits of displacement noise free interferometers, arXiv:2209.02998 10.48550/arXiv.2209.02998 (2022).
  * Modi  _et al._ [2012] K. Modi, A. Brodutch, H. Cable, T. Paterek, and V. Vedral, The classical-quantum boundary for correlations: Discord and related measures, Rev. Mod. Phys. 84, 1655 (2012).
  * Proctor  _et al._ [2018] T. J. Proctor, P. A. Knott, and J. A. Dunningham, Multiparameter estimation in networked quantum sensors, Phys. Rev. Lett. 120, 080501 (2018).
  * Ge  _et al._ [2018] W. Ge, K. Jacobs, Z. Eldredge, A. V. Gorshkov, and M. Foss-Feig, Distributed quantum metrology with linear networks and separable inputs, Phys. Rev Lett. 121, 043604 (2018).
  * Eldredge  _et al._ [2018] Z. Eldredge, M. Foss-Feig, J. A. Gross, S. L. Rolston, and A. V. Gorshkov, Optimal and secure measurement protocols for quantum sensor networks, Phys. Rev. A 97, 042337 (2018).
  * Rubio  _et al._ [2020] J. Rubio, P. A. Knott, T. J. Proctor, and J. A. Dunningham, Quantum sensing networks for the estimation of linear functions, J. Phys. A Math. Theor. 53, 344001 (2020).
  * Kiesenhofer  _et al._ [2023] D. Kiesenhofer, H. Hainzer, A. Zhdanov, P. C. Holz, M. Bock, T. Ollikainen, and C. F. Roos, Controlling two-dimensional Coulomb crystals of more than 100 ions in a monolithic radio-frequency trap, arXiv:2302.00565 10.48550/arXiv.2302.00565 (2023).
  * Joshi  _et al._ [2020] M. K. Joshi, A. Fabre, C. Maier, T. Brydges, D. Kiesenhofer, H. Hainzer, R. Blatt, and C. F. Roos, Polarization-gradient cooling of 1D and 2D ion Coulomb crystals, New. J. Phys. 22, 103013 (2020).
  * Kranzl  _et al._ [2022] F. Kranzl, M. K. Joshi, C. Maier, T. Brydges, J. Franke, R. Blatt, and C. F. Roos, Controlling long ion strings for quantum simulation and precision measurements, Phys. Rev. A 105, 052426 (2022).
  * Schmidt  _et al._ [2009] M. Schmidt, E. van den Berg, M. Friedlander, and K. Murphy, Optimizing costly functions with simple constraints: A limited-memory projected quasi-newton algorithm, in  _Proceedings of the 12 th International Conference on Articial Intelligence and Statistics (AISTATS)_ (2009).
  * Hainzer  _et al._ [2022] H. Hainzer, D. Kiesenhofer, T. Ollikainen, M. Bock, F. Kranzl, M. K. Joshi, G. Yoeli, R. Blatt, T. Gefen, and C. F. Roos, Correlation spectroscopy with multi-qubit-enhanced phase estimation: data. Zenodo (2022).
  * Holtzman [1950] W. H. Holtzman, The unbiased estimate of the population variance and standard deviation, The American Journal of Psychology 63, 615 (1950).
  * Cover [1999] T. M. Cover,  _Elements of information theory_ (John Wiley & Sons, 1999).
  * Braunstein and Caves [1994] S. L. Braunstein and C. M. Caves, Statistical distance and the geometry of quantum states, Phys. Rev. Lett. 72, 3439 (1994).
  * Liu  _et al._ [2019] J. Liu, H. Yuan, X.-M. Lu, and X. Wang, Quantum fisher information matrix and multiparameter estimation, J. Phys. A Math. Theor. 53, 023001 (2019).
  * Bothwell  _et al._ [2022] T. Bothwell, C. J. Kennedy, A. Aeppli, D. Kedar, J. M. Robinson, E. Oelker, A. Staron, and J. Ye, Resolving the gravitational redshift across a millimetre-scale atomic sample, Nature 602, 420 (2022).
  * Zheng  _et al._ [2022] X. Zheng, J. Dolde, V. Lochab, B. N. Merriman, H. Li, and S. Kolkowitz, Differential clock comparisons with a multiplexed optical lattice clock, Nature 602, 425 (2022).
  * Madjarov  _et al._ [2019] I. S. Madjarov, A. Cooper, A. L. Shaw, J. P. Covey, V. Schkolnik, T. H. Yoon, J. R. Williams, and M. Endres, An atomic-array optical clock with single-atom readout, Phys. Rev. X 9, 041052 (2019).
  * Kómár  _et al._ [2014] P. Kómár, E. M. Kessler, M. Bishof, L. Jiang, A. S. Sørensen, J. Ye, and M. D. Lukin, A quantum network of clocks, Nat. Phys. 10, 582 (2014).
  * Nichol  _et al._ [2022] B. C. Nichol, R. Srinivas, D. P. Nadlinger, P. Drmota, D. Main, G. Araneda, C. J. Ballance, and D. M. Lucas, An elementary quantum network of entangled optical atomic clocks, Nature 609, 689 (2022).
  * Dubin [1993] D. H. E. Dubin, Theory of structural phase transitions in a trapped Coulomb crystal, Phys. Rev. Lett. 71, 2753 (1993).
  * Kiethe  _et al._ [2021] J. Kiethe, L. Timm, H. Landa, D. Kalincev, G. Morigi, and T. E. Mehlstäubler, Finite-temperature spectrum at the symmetry-breaking linear to zigzag transition, Phys. Rev. B 103, 104106 (2021).
  * Linington and Vitanov [2008] I. E. Linington and N. V. Vitanov, Robust creation of arbitrary-sized Dicke states of trapped ions by global addressing, Phys. Rev. A 77, 010302(R) (2008).
  * Um  _et al._ [2016] M. Um, J. Zhang, D. Lv, Y. Lu, S. An, J.-N. Zhang, H. Nha, M. S. Kim, and K. Kim, Phonon arithmetic in a trapped ion system, Nat. Commun. 7, 11410 (2016).
  * Lechner  _et al._ [2016] R. Lechner, C. Maier, C. Hempel, P. Jurcevic, B. P. Lanyon, T. Monz, M. Brownnutt, R. Blatt, and C. F. Roos, Electromagnetically-induced-transparency ground-state cooling of long ion strings, Phys. Rev. A 93, 053401 (2016).
  * Kay [1993] S. M. Kay,  _Fundamentals of statistical signal processing: estimation theory_ (Prentice-Hall, Inc., 1993).
  * Demkowicz-Dobrzanski  _et al._ [2009] R. Demkowicz-Dobrzanski, U. Dorner, B. J. Smith, J. S. Lundeen, W. Wasilewski, K. Banaszek, and I. A. Walmsley, Quantum phase estimation with lossy interferometers, Phys. Rev. A 80, 013825 (2009).
  * Pezze  _et al._ [2017] L. Pezze, M. A. Ciampini, N. Spagnolo, P. C. Humphreys, A. Datta, I. A. Walmsley, M. Barbieri, F. Sciarrino, and A. Smerzi, Optimal measurements for simultaneous quantum estimation of multiple phases, Phys. Rev. Lett. 119, 130504 (2017).