# Nearly Heisenberg-limited noise-unbiased frequency estimation by tailored
sensor design

Francisco Riberi  Department of Physics and Astronomy, Dartmouth College, 6127
Wilder Laboratory, Hanover, New Hampshire 03755, USA Gerardo A. Paz-Silva
Centre for Quantum Computation and Communication Technology (Australian
Research Council),
Centre for Quantum Dynamics, Griffith University, Brisbane, Queensland 4111,
Australia  Lorenza Viola  Department of Physics and Astronomy, Dartmouth
College, 6127 Wilder Laboratory, Hanover, New Hampshire 03755, USA

###### Abstract

We consider entanglement-assisted frequency estimation by Ramsey
interferometry, in the presence of dephasing noise from general
spatiotemporally correlated environments. By working in the widely employed
local estimation regime, we show that even for infinite measurement
statistics, noise renders standard estimators biased or ill-defined. We
introduce ratio estimators which, at the cost of doubling the required
resources, are insensitive to noise and retain the asymptotic precision
scaling of standard ones. While ratio estimators are applicable also in the
limit of Markovian noise, we focus on non-Markovian dephasing from a bosonic
bath and show how knowledge about the noise spectrum may be used to maximize
metrological advantage, by tailoring the sensor’s geometry. Notably,
Heisenberg scaling is attained up to a logarithmic prefactor by maximally
entangled states, while optimal Zeno scaling is afforded by one-axis twisted
spin-squeezed states.

##  I Introduction

High-precision estimation of transition frequencies (or energy splittings) is
a fundamental task in quantum metrology, with implications ranging from atomic
spectroscopy [1, 2, 3] to time-keeping with atomic clocks [4, 5, 6]. In the
context of Ramsey interferometry [7, 8] with a quantum sensor comprising $N$
probes, the use of initial entangled states can yield, in principle,
asymptotic precision bounds which surpass the optimal $N^{-1/2}$ standard
quantum limit (SQL) achievable classically. Under ideal conditions, and
assuming that no interaction between the probes is permitted after
preparation, the ultimate $N^{-1}$ precision bound is set by the Heisenberg
limit (HL), and is saturated by maximally entangled Greenberger-Horne-
Zeilinger (GHZ) states [9, 10, 11].

In practice, noise inevitably degrades the attainable precision, to an extent
that depends on the model specifics [12]. In many quantum sensors, dephasing
noise which couples through the same operator as the signal provides the
dominant noise mechanism. While no gain over the SQL can be achieved for
independent Markovian noise (that is, noise with no spatial and temporal
correlations) [9, 13], temporal correlations alone can be exploited to achieve
a superclassical Zeno limit $\propto N^{-3/4}$ [14, 15, 16]. Strong temporal
correlations have been detected in a variety of systems via quantum noise
spectroscopy experiments, see e.g., [17, 18, 19, 20, 21, 22, 23], including
for non-classical noise environments [24, 25, 26]. Spatial noise correlations
also tend to emerge due to probe proximity [27, 28, 29, 30], making noise
substantially more harmful than uncorrelated one. For perfectly correlated
(collective) Markovian noise, GHZ states are the most fragile, resulting in an
$N$-independent precision scaling [28, 31], and sub-SQL scaling is also
precluded in the non-Markovian regime [32, 33]. For noise with partial spatial
correlations, superclassical precision scaling can be restored by tailored
error-correcting codes in the Markovian case [34], or by means of a randomized
protocol for general temporally correlated scenarios [33].

While, as the above shows, the impact noise on the scaling of precision has
been extensively studied, far less attention has been devoted thus far to the
fact that noise may also introduce unwanted _bias_ , compromising accuracy if
unaccounted for. Most research on bias in quantum metrology has focused on
estimation in the regime of limited data, using Bayesian approaches [35, 36,
37]. More recently, in the context of Markovian noise, bias due to finite-
frequency error-corrected sensing was addressed in [38] through post-
processing, while a purification-based protocol was proposed in [39] to
mitigate bias due to imperfect knowledge of the noise model.

Here, we focus on general spatiotemporally correlated dephasing noise,
potentially arising from nonclassical sources, and address its impact on both
the mean and the variance of the possible measurement outcomes. We first show
how, even in the ideal limit of infinite measurement statistics and perfect
knowledge of the noise spectral properties, standard estimators for both GHZ
and one-axis twisted (OAT) states become systematically biased and possibly
ill-defined. By suitably modifying the estimation protocol, we introduce ratio
estimators, which are insensitive to dephasing by construction and match the
asymptotic $N$-scaling of standard estimators. Importantly, no detailed noise
knowledge is needed for constructing ratio estimators, which remain applicable
in the limiting case of noise with no temporal or spatial correlations.

As our second main result, we further show that, if spatiotemporal
correlations are present, access to the noise spectral features allows for the
achievable precision scaling to be optimized, for a given initial state. We
make our analysis concrete by examining a setting where $N$ qubits are placed
on a one-dimensional (1D) regular lattice with tunable separation and couple
to a bosonic bath. We show that engineering _negative_ noise correlations
yields a far greater scaling advantage than achievable via randomized
protocols that spatially decorrelate the probes on average [33]. Remarkably,
OAT states saturate the $N^{-3/4}$ Zeno limit which is optimal for independent
non-Markovian dephasing [16], whereas a novel nearly-Heisenberg
$N^{-1}\sqrt{\log(N)}$ scaling emerges for GHZ states.

The content is organized as follows. Section II summarizes the main tools
needed to understand our work, including the spatiotemporally correlated
dephasing setting we consider and some basic estimation-theoretic notions used
in standard metrology protocols. In Sec. III, we demonstrate the emergence of
noise-induced bias, and show how the latter can be countered by introducing
ratio estimators, at the cost of doubling the total time resources needed to
carry out the experiment. The performance of these new estimators is
illustrated in the simplest instance of perfectly correlated (collective) non-
Markovian dephasing. In Sec. IV, we specialize our analysis to a 1D lattice
geometry, and show explicitly how to optimize the sensor’s performance by
tailoring the spatial separation of the qubit probes, while preparing them in
a GHZ and OAT state, respectively. Some considerations about the robustness of
our predicted superclassical precision scaling against variations of the
underlying noise model are also included. We briefly summarize our conclusions
in Sec. V, and provide complete detail about our scaling derivations for both
GHZ and OAT states in a separate appendix.

##  II Background

###  II.1 Open quantum system setting

We consider $N$ qubit probes, each coupled to the target frequency $b$, and a
“parallel” dephasing noise environment (or bath). In the interaction picture
with respect to the free bath Hamiltonian, $H_{\rm B}$, the joint evolution is
generated by

$\displaystyle H_{\rm SB}(t)=bJ_{z}\otimes
I_{B}+\frac{1}{2}\sum_{n=1}^{N}\,\sigma^{z}_{n}\otimes B_{n}(t)\,,$ (1)

where $\sigma^{u}_{n}$, $J_{u}\equiv\tfrac{1}{2}\sum_{n=1}^{N}\sigma_{n}^{u}$,
$u\in\\{x,y,z\\}$, are Pauli matrices and collective spin operators,
respectively, and the bath operators $\\{B_{n}(t)\\}$ describe a noise process
which we take to be zero-mean, stationary, and Gaussian. Under the assumption
that the initial joint state is factorized, $\rho_{\rm
SB}(0)\equiv\rho_{0}\otimes\rho_{\rm B}$, with $[\rho_{\rm B},H_{\rm B}]=0$,
the statistical properties of the noise are fully captured by the two-point
correlation functions,

$C_{nm}(t)\equiv\langle B_{n}(t)B_{m}(0)\rangle=\text{Tr}_{\rm
B}\\{B_{n}(t)B_{m}(0)\rho_{\rm B}\\}.$

For a classical noise environment, $B_{n}(t)$ and $\langle\bullet\rangle$
denote a stochastic commuting process and an ensemble average, respectively.
In general, the noise is non-trivially correlated both in space, across
different qubits, and in time (non-Markovian); “$\delta$-correlated”
(Markovian) noise is included as a special case, $C_{nm}(t)=c_{nm}\delta(t)$,
for some Hermitian matrix $c_{nm}$ that encodes the spatial noise
correlations.

In the frequency domain, the presence of temporal correlations translates into
“colored” classical (symmetrized, $+$) and quantum (anti-symmetrized, $-$)
noise spectra, given by

$S^{\pm}_{nm}(\omega)\equiv\int_{-\infty}^{\infty}\\!dt\,e^{-i\omega
t}\langle[B_{n}(t),B_{m}(0)]_{\pm}\rangle,$

with $[B_{n}(t),B_{m}(0)]_{-}$ and $[B_{n}(t),B_{m}(0)]_{+}$ denoting the
commutator and anti-commutator, respectively, of the relevant noise operator
[40]. Let $|\vec{\alpha}\rangle\equiv\bigotimes_{n=1}^{N}|\alpha_{n}\rangle$,
with $\alpha_{n}=\pm 1$ corresponding to the eigenstates
$\\{|\\!\\!\\!\uparrow\rangle,|\\!\\!\\!\downarrow\rangle\\}$ of
$\sigma_{n}^{z}$, denote the $z$ basis. The time-evolved state of the system
may then be generally represented as

$\langle\vec{\alpha}|\rho(t)|\vec{\beta}\rangle=e^{ibt\sum_{n=1}^{N}(\beta_{n}-\alpha_{n})}e^{-\gamma(t)+i\varphi_{0}(t)+i\varphi_{1}(t)}\langle\vec{\alpha}|\rho_{0}|\vec{\beta}\rangle,$
(2)

in terms of real functions $\gamma(t)$, $\varphi_{0}(t)$, $\varphi_{1}(t)$
which involve products of a state-dependent component and a corresponding
time-dependent dynamic coefficient [33]. Irrespective of the classical or
quantum nature of the noise, $\gamma(t)$ governs the decay of off-diagonal
coherence elements,

$\gamma(t)\equiv\sum_{n,m=1}^{N}(\alpha_{n}-\beta_{n})(\alpha_{m}-\beta_{m})\kappa_{nm}(t)$
(3)

where the decay dynamic coefficient $\kappa_{nm}(t)$ is expressible in terms
of a frequency overlap integral:

$\displaystyle\kappa_{nm}(t)=\frac{1}{32\pi}\;\\!\int_{-\infty}^{\infty}\\!\\!d\omega\,\frac{\sin^{2}(\omega
t/2)}{\omega^{2}}\,S_{nm}^{+}(\omega).$ (4)

On the other hand, phase evolution is distinctive of a quantum, non-commuting
environment. One may show that $\varphi_{0}(t)$ arises from a unitary “Lamb-
shift” contribution due to bath-mediated entanglement between the qubits,
whereas $\varphi_{1}(t)$ is linked to whether the dephasing can be described
as “random unitary,” hence as classical in nature [33, 41]. In particular, we
can write

$\varphi_{0}(t)\equiv\sum_{n,m=1}^{N}(\beta_{n}\beta_{m}-\alpha_{n}\alpha_{m})\xi_{nm}(t),$
(5)

with a corresponding phase dynamic coefficient $\xi_{nm}(t)$ which can be
expressed as

$\displaystyle\xi_{nm}(t)=\frac{1}{32\pi}\;\\!\int_{-\infty}^{\infty}\\!\\!d\omega\,\frac{\omega
t-\sin(\omega t)}{\omega^{2}}\,S_{nm}^{-}(\omega),$ (6)

and depends explicitly on the quantum spectra $S_{nm}^{-}(\omega)$.
Structurally similar relationships hold for $\varphi_{1}(t)$ although, for
symmetry reasons, this phase will not play a role in the model systems we
shall focus on [33]. Under the assumption that the spectra have vanishing
support above a high-frequency cutoff, say, $S^{\pm}_{nm}(\omega)\approx 0$
for $|\omega|\gtrsim\omega_{c}$, Eq. (4) implies a quadratic dependence upon
time of $\kappa_{nm}(t)$, hence of $\gamma(t)$, in the short-time limit
$\omega_{c}t\ll 1$. This contrasts with the linear behavior
($\dot{\gamma}(t)=\,$const) that arises in the formal limit
$\omega_{c}\rightarrow\infty$ of Markovian noise, described by semigroup
dynamics [12].

###  II.2 Noiseless frequency estimation

####  II.2.1 Basic estimation notions

A typical estimation protocol consists of tree distinct steps. On a first
stage the $N$ qubit probes are initialized in a (possibly entangled) input
state $\rho_{0}$. In the absence of noise, the system is then left to evolve
unitarily under the Hamiltonian in Eq. (1) with $B_{n}(t)=0$, for an
“encoding” period of duration $\tau$, after which a suitable observable, say,
${\cal O}$, is measured. This process is repeated $\nu\equiv T/\tau\gg 1$
times over a total duration $T$, resulting in a vector of independently
distributed measurement outcomes
$\vec{\mu}\equiv\\{\mu_{1},\ldots,\mu_{\nu}\\}$, from which information about
$b$ is extracted. An estimator $\hat{b}(\vec{\mu})$ is a function of random
variable that associates each set of outcomes with an estimate of $b$. Let the
mean and variance of $\hat{b}$ be denoted, respectively, by
$\langle\hat{b}(\vec{\mu})\rangle$ and
$\Delta\hat{b}^{2}(\vec{\mu})=\langle(\hat{b}(\vec{\mu})-\langle\hat{b}(\vec{\mu})\rangle)^{2}\rangle$,
where expectations $\langle\cdot\rangle$ are taken over all possible
measurement outcomes. The estimator is unbiased if
$\langle\hat{b}(\vec{\mu})\rangle=b$ and is asymptotically unbiased if the
bias vanishes as $\nu\rightarrow\infty$. More precisely, we say that
$\hat{b}(\vec{\mu})$ is consistent if it converges in probability to the true
value for an infinitely large sample, that is, $\forall\varepsilon>0$,
$\lim_{\nu\rightarrow\infty}{\mathbb{P}}(|\hat{b}(\vec{\mu})-b|>\varepsilon)=0$,
implying $\lim_{\nu\rightarrow\infty}\hat{b}(\vec{\mu})=b$. A consistent
estimator with vanishing variance as $\nu\rightarrow\infty$ is asymptotically
unbiased. The estimators we are interested in in this work fall under this
category.

In a setting where the relationship between the expectation value of
$\mathcal{O}$ and $b$ is well-known, $\langle\mathcal{O}(\tau)\rangle\equiv
f(b)$, an estimator $\hat{b}(\vec{\mu})$ is customarily constructed as
follows: First, the frequency can be written as a function of
$\langle\mathcal{O}(\tau)\rangle$ by inverting the above relationship,
$b\equiv f^{-1}(\langle\mathcal{O}(\tau)\rangle)$. Then, after performing
$\nu\gg 1$ measurements, we estimate $\langle\mathcal{O}(\tau)\rangle$ as the
sample mean of the measurements outcomes,

$\langle\mathcal{O}(\tau)\rangle\approx\langle\hat{\mathcal{O}}(\tau)\rangle=\frac{1}{\nu}\,\sum_{i=1}^{\nu}\mu_{i}\equiv\hat{f}(b).$

From this, the estimator readily follows:
$\hat{b}(\vec{\mu})=f^{-1}(\hat{f}(b))$. An analytic expression for
$\Delta\hat{b}^{2}$ can be obtained if we further assume that the target
frequency is confined within a small neighborhood of a known value $b_{0}$,
with a corresponding mean value $\langle\mathcal{O}(\tau)\rangle_{0}$, such
that $f(b)$ is injective. To do so, we first linearize around $b_{0}$,

$\displaystyle\hat{b}\approx f^{-1}(f(b_{0}))+\frac{\partial}{\partial
f}\,f^{-1}(f(b))\big{|}_{f(b_{0})}(\hat{f}(b)-f(b_{0})),$

whereby it follows, using the expression for the derivative of an inverse
function, that

$\displaystyle\hat{b}-b_{0}\approx(\partial_{b}f(\hat{b})\big{|}_{b_{0}})^{-1}(\langle\hat{\mathcal{O}}(\tau)\rangle-\langle\mathcal{O}(\tau)\rangle_{0}).$

Because the variance is independent of additive constants ($b_{0}$ and
$\langle\mathcal{O}(\tau)\rangle_{0}$ in our case), it is then clear that the
variance of $\hat{b}$ is proportional to
$\Delta\hat{\mathcal{O}}^{2}(\tau)=\langle\hat{\mathcal{O}}^{2}(\tau)\rangle-\langle\hat{\mathcal{O}}(\tau)\rangle^{2}$
in this limit, leading to the well-known error propagation formula

$\displaystyle\Delta\hat{b}^{2}(\tau)\approx(\partial_{b}\,\langle\mathcal{O}(\tau)\rangle\big{|}_{b_{0}})^{-2}\,\nu^{-1}\Delta\hat{\mathcal{O}}(\tau)^{2}.$
(7)

This procedure, known as the method of moments [3], has been widely used for
quantifying estimation precision in metrology settings. The approximate
equality in Eq. (7) allows us to obtain an analytic expression for the
uncertainty in the asymptotic limit of $\nu\gg 1$ and high prior knowledge,
and is commonly treated as an exact measure of the frequency uncertainty. We
henceforth follow this standard practice, and write the relationship in Eq.
(7) as an equality.

For the particular case of an initial $N$-qubit GHZ state, $|{\rm
GHZ}\rangle\equiv\tfrac{1}{\sqrt{2}}(|\\!\\!\uparrow\rangle^{\otimes
N}\\!+|\\!\\!\downarrow\rangle^{\otimes N})$, we will explicitly verify that
the analytic expression in Eq. (7) accurately reproduces the behavior of the
average over all measurement outcomes of an appropriate survival probability
(see Sec. II.2.2). In addition, we are also interested on frequency estimation
through Ramsey interferometry for spin-squeezed initial states, which exhibit
a reduced variance of the collective spin $J_{u}$ along a particular
direction, and are easier to generate and measure experimentally than GHZ
states are [42]. It can be shown that their uncertainty $\Delta\hat{b}(\tau)$
can scale superclassically when measuring an angular momentum component,
$\mathcal{O}=J_{\hat{n}}$. Here, we consider OAT states obtained from product
coherent spin states (CSS) as follows [43, 44, 45],

$|\text{OAT}\rangle=e^{-i\beta J_{x}}e^{-i\theta
J_{z}^{2}}|\text{CSS}\rangle_{x},\qquad|\text{CSS}\rangle_{x}\equiv|+\rangle^{\otimes
N},$ (8)

with $\beta,\theta\in[0,2\pi]$ being rotation and squeezing angles set to
minimize the initial variance along $y$, and $\sigma_{x}|+\rangle=|+\rangle$.
For such an OAT state, the precision $\Delta\hat{b}$ resulting from measuring
$J_{y}$ scales like $N^{-5/6}$ in the noiseless scenario [33].

Lastly, we remark that if $\hat{b}$ is an unbiased estimator, the variance
associated with measuring any ${\cal O}$ is lower-bounded by the quantum
Cramér-Rao bound (QCR) [46, 47],

$\displaystyle\Delta\hat{b}^{2}\geq\Delta\hat{b}^{2}_{\mathrm{QCR}}=(\nu
F_{\mathrm{Q}}[\rho_{b}(\tau)])^{-1},$ (9)

where $F_{\mathrm{Q}}[\rho_{b}(\tau)]$ is the quantum Fisher information (QFI)
of the output state $\rho_{b}(\tau)$. By construction, the QFI is derived from
the maximization of its classical, operator-dependent counterpart, the _Fisher
Information_ over the set of all positive operator valued measurements
(POVMs), which includes all possible observables $\mathcal{O}$ [48].

####  II.2.2 Case study: GHZ state

In the absence of noise, an input GHZ state saturates the HL [3], as can be
shown by performing a measurement of the survival probability $p_{b,0}(\tau)$.
Since the dynamics is unitary in this case, at a given interrogation time
$\tau$ the latter is given by $p_{b,0}(\tau)=(1+\cos(Nb\tau))/2$. Let us apply
the method of moments by defining the relevant observable

$\mathcal{O}\equiv\Pi_{|\text{GHZ}\rangle}-(I-\Pi_{|\text{GHZ}\rangle}),\quad\Pi_{|\text{GHZ}\rangle}=|\text{GHZ}\rangle\langle\text{GHZ}|,$
(10)

in terms of the projector $\Pi_{|\text{GHZ}\rangle}$ onto the initial state.
In this way, $\langle\mathcal{O}(\tau)\rangle=\cos(Nb\tau)\equiv f(b)$. The
frequency can then be written as a function of
$\langle\mathcal{O}(\tau)\rangle$ by inverting this relationship,

$b=\arccos(\langle\mathcal{O}(\tau)\rangle)/(N\tau)\equiv
f^{-1}(\langle\mathcal{O}(\tau)\rangle).$

After performing $\nu\gg 1$ measurements, we can estimate
$\langle\mathcal{O}(\tau)\rangle$ in terms of the ratio of the number of
detections of the initial state, $\nu_{+}$, over the total number of trials:
$\langle\hat{\mathcal{O}}(\tau)\rangle=\hat{f}(b)=2(\nu_{+}/\nu)-1$. This
leads to to the estimator

$\hat{b}(\tau)=f^{-1}(\hat{f}(b))=\frac{1}{N\tau}\arccos(2(\nu_{+}/\nu)-1).$

Thanks to the fact that the relevant POVM has only two outcomes, the resulting
measurement probability distribution is binomial, and the estimator’s variance
can thus be computed numerically by averaging over all measurement outcomes,

$\displaystyle\Delta\hat{b}^{2}(\tau)=\langle(\hat{b}(\vec{\mu})-\langle\hat{b}(\vec{\mu})\rangle)^{2}\rangle,$
(11)
$\displaystyle\langle\bullet\rangle=\sum_{\nu_{+}=0}^{\nu}{\nu\choose\nu_{+}}\;p_{b,0}(\tau)^{\nu_{+}}\,(1-p_{b,0}(\tau))^{\nu-\nu_{+}}\;(\bullet)\;.$

On the other hand, $\Delta\hat{b}^{2}(\tau)$ may be calculated by means of Eq.
(7). Using the second moment of the binomial distribution,
$\Delta\hat{\mathcal{O}}(\tau)^{2}=\nu^{-1}4p_{b,0}(\tau)(1-p_{b,0}(\tau))$,
Eq. (7) reduces to $\Delta\hat{b}^{2}(\tau)\approx(N^{2}\,T\,\tau)^{-1}$, that
is, the HL. It follows that this particular choice of measurement operator
$\mathcal{O}$ is optimal, as it saturates the QCR in Eq. (9). We can verify
that such an approximate analytic expression accurately reproduces the
behavior of the average over all measurement outcomes: as one may see from
Fig. 1 (top), the agreement is excellent already for $\nu=400$ measurements.

Figure 1: Standard estimator uncertainty in the noiseless setting and noise-
induced bias. Top: Standard estimator uncertainty as a function of time in the
ideal limit of no dephasing. Sample variance (blue, dashed), and analytic
expression $\Delta\hat{b}(t)\approx(N\sqrt{Tt})^{-1}$ derived from the method
of moments (grey, solid). Bottom: Performance of the standard estimator in the
presence of collective spin-boson dephasing from a 1D zero-temperature
environment (see Sec. IV.1). The spectral density [Eq. (22)] is supra-Ohmic
with an exponential cutoff, with $s=3$, $\omega_{c}=1$, and different noise
strengths $\alpha$, as shown in the legend. The limiting analytic expression
in the absence of noise (grey, solid) is also included. The measurement time
$\tau=0.2067$ is optimal in the lattice setting. We consider an initial GHZ
state with $N=100$ qubits, and $\nu=400$ measurement shots.

##  III Noisy frequency estimation

###  III.1 Noise-induced bias

Crucially, the method of moments we described above relies on the assumption
that the function $f$ to be one-to-one, at least in a neighborhood of $b_{0}$.
Despite being ubiquitous in the literature, this inversion procedure on which
the construction of $\hat{b}(\vec{\mu})$ hinges upon generally becomes
problematic in the presence of noise. To illustrate the issue, consider as
before an initial GHZ state, but this time evolving under the dephasing
Hamiltonian of Eq. (1) which, to simplify our argument, we take to be
collective (permutation-invariant), that is, $B_{n}(t)\equiv B(t),\forall n$.

In the presence of dephasing, the survival probability becomes
$p_{b}(\tau)=(1+\cos(Nb\tau)e^{-\gamma_{\text{GHZ}}(t)})/2$ and,
correspondingly, the mean value is multiplied by a decay factor,
$\langle\mathcal{O}(\tau)\rangle=e^{-\gamma_{\text{GHZ}}(t)}\cos(Nb\tau)\equiv
g(b)$. Thus, using the same inversion formula to estimate $b$ as in the
noiseless scenario,
$\hat{b}(\tau)=\arccos(2(\nu_{+}/\nu)-1)/(N\tau)=f^{-1}(g(b))$, leads to a
systematic _noise-induced bias_ as $\gamma_{\text{GHZ}}(t)$ is not accounted
for. This is illustrated in Fig. 1(bottom), which shows how the performance of
the standard estimator increasingly degrades as the dephasing rate increases.
Alternatively, computation of $\hat{b}$ using the correct inversion formula
would require perfect knowledge of the decay parameter
$\gamma_{\text{GHZ}}(\tau)$,

$\hat{b}(\tau)=\frac{1}{N\tau}\arccos\big{[}e^{\gamma_{\text{GHZ}}(\tau)}(2(\nu_{+}/\nu)-1)\big{]},$

which would imply previous noise characterization to the same level of
precision desired for the estimation of $b$. Not only would this add to the
resource count, but there is no known noise characterization protocol capable
of achieving sub-SQL scaling. Even more importantly from a conceptual
standpoint, the above recipe returns imaginary values for $\hat{b}$ whenever
$e^{\gamma_{\text{GHZ}}(\tau)}(2(\nu_{+}/\nu)-1)\geq 1$. Note that this
happens for a non-vanishing set of outcomes whenever
$\gamma_{\text{GHZ}}(\tau)>0$, rendering quantities such as
$\langle\hat{b}(\vec{\mu})\rangle$ ill-defined.

Despite the pitfalls involved in this procedure, the analytical expression for
the uncertainty which follows from the method of moments has found widespread
application in the literature. For a fixed total time $T$, we have [12]

$\displaystyle\Delta\hat{b}(\tau)^{2}=\frac{1}{T\tau
N^{2}}\bigg{[}\frac{1-e^{-2\gamma_{\text{GHZ}}(\tau)}\cos^{2}(Nb\tau)}{e^{-2\gamma_{\text{GHZ}}(\tau)}\sin^{2}(Nb\tau)}\bigg{]}.$
(12)

The performance of Eq. (12) can be quantitatively assessed once a specific
noise model is chosen. For collective spin-boson dephasing [32], which we will
also recover as a limiting case in Sec. IV.1, we have
$\gamma_{\text{GHZ}}(\tau)=N^{2}\kappa(\tau)$. Optimizing with respect to the
phase argument $\varphi\equiv N\tau b$ by demanding that
$\varphi=(2k+1)\pi/2$, $k\in\mathbb{N}$, and minimizing with respect to $\tau$
in the short-time regime $\omega_{c}\tau\ll 1$ where
$\gamma_{\text{GHZ}}(\tau)\approx N^{2}\kappa_{0}^{2}(\omega_{c}\tau)^{2}$ for
some dimensionless constant $\kappa_{0}^{2}$, the optimal performance can be
shown to be SQL-limited. Specifically,
$\Delta\hat{b}^{\text{GHZ}}_{\text{opt,coll}}=e^{1/4}(\kappa_{0}\omega_{c}/T)^{1/2}N^{-1/2}$,
for an optimal measurement time $\tau^{\text{GHZ}}_{\text{opt,
coll}}=1/2(\kappa_{0}\omega_{c})^{-1}N^{-1}$.

###  III.2 Ratio estimator for GHZ states

To circumvent the inversion issues plaguing the estimator linked to Eq. (12),
we propose to modify the estimation protocol in such a way that two suitable
observables are independently measured, and their outcomes combined to obtain
a noise-robust estimator. Concretely, for the GHZ state we are focusing on, we
propose to not only measure the observable $\mathcal{O}$ given in Eq. (10),
but also
$\mathcal{O}^{\prime}\equiv\Pi_{|{\rm{GHZ}^{\prime}}\rangle}-(I-\Pi_{|{\rm{GHZ}^{\prime}}\rangle})$,
where we now project into a related state within the GHZ class:

$|{\rm{GHZ}^{\prime}}\rangle\equiv\frac{1}{\sqrt{2}}\big{(}|\\!\\!\uparrow\rangle^{\otimes
N}+i|\\!\\!\downarrow\rangle^{\otimes N}\big{)}.$

Accordingly, the probability of measuring $\Pi_{|{\rm{GHZ}^{\prime}}\rangle}$
is given by
$p_{b}^{\prime}(\tau)=(1+\sin(Nb\tau)e^{-\gamma_{\text{GHZ}}(t)})/2$. With
each observable being measured $\nu=T/\tau$ times, our protocol implies a
total of $2\nu$ measurements: that is, we double the time resources needed to
carry out the experiment from $T$ to $2T$. Since
$\langle\mathcal{O}^{\prime}(\tau)\rangle=e^{-\gamma_{\text{GHZ}}(t)}\sin(Nb\tau)\equiv
h(b)$, we can solve for the frequency by taking the ratio between the two mean
values,
$b=\arctan[{\langle\mathcal{O}^{\prime}(\tau)\rangle}/\langle\mathcal{O}(\tau)\rangle]/(N\tau)$.
It follows that a noise-robust _ratio_ estimator for $b$ can be constructed as

$\hat{b}_{\text{R}}^{\text{GHZ}}(\tau)\equiv\frac{1}{N\tau}\arctan\big{[}\langle\hat{\mathcal{O}}^{\prime}(\tau)\rangle/\langle\hat{O}(\tau)\rangle\big{]},$
(13)

where
$\langle\hat{\mathcal{O}}^{\prime}(\tau)\rangle=\hat{h}(b)=2(\nu^{\prime}_{+}/\nu)-1$
and $\langle\hat{\mathcal{O}}(\tau)\rangle=\hat{g}(b)=2(\nu_{+}/\nu)-1$ are
the finite-sample estimators to the respective mean values, obtained from
$\nu^{\prime}_{+}$ (respectively, $\nu_{+}$) detections in $\nu$ trials. In
Fig. 2 (top), the limiting estimator function for $\hat{b}_{\text{R}}(\tau)$,
$\arctan[\tan(Nb\tau)]/(N\tau)$, is compared to an average over all
measurement outcomes which now involves two independent probability
distributions corresponding to $\mathcal{O}$ and $\mathcal{O}^{\prime}$, for
two different numbers of repetitions $\nu$ and $\nu^{\prime}$, showing
excellent asymptotic convergence.

In a similar way, for the variance of the ratio estimator in Eq. (13), we have

$\displaystyle\Delta\hat{b}_{\text{R}}^{\text{GHZ}}(\tau)^{2}=\langle(\hat{b}_{\text{R}}(\vec{\mu},\vec{\mu^{\prime}})-\langle\hat{b}(\vec{\mu},\vec{\mu}^{\prime})\rangle)^{2}\rangle,$
(14) $\displaystyle\langle\bullet\rangle$
$\displaystyle=\bigg{\\{}\sum_{\nu_{+},\nu^{\prime}_{+}=0}^{\nu}{\nu\choose\nu_{+}}\;p_{b}(\tau)^{\nu_{+}}\,(1-p_{b}(\tau))^{\nu-\nu_{+}}$
$\displaystyle\hskip
56.9055pt{\nu\choose\nu^{\prime}_{+}}\;p^{\prime}_{b}(\tau)^{\nu^{\prime}_{+}}\,(1-p^{\prime}_{b}(\tau))^{\nu-\nu^{\prime}_{+}}\;(\bullet)\bigg{\\}}.$

Assuming that our prior knowledge confines the frequency within a neighborhood
of $b_{0}$ in which both $g(b)$ and $h(b)$ are injective, we can linearize to
obtain an analytic expression for the uncertainty in terms of the dispersions
$\Delta\hat{O}^{2}(\tau)$ and
$\Delta\hat{O}^{\prime\,2}(\tau)=\nu^{-1}4\,p_{b}^{\prime}(\tau)(1-p_{b}^{\prime}(\tau))$.
This leads to

$\displaystyle\Delta\hat{b}_{\text{R}}^{\text{GHZ}}(\tau)^{2}$
$\displaystyle\\!\\!=\\!$
$\displaystyle\\!\bigg{(}\\!\frac{\partial\hat{b}_{\text{R}}}{\partial\langle\mathcal{O}^{\prime}(\tau)\rangle}\\!\bigg{)}^{2}\\!\Delta\hat{\mathcal{O}}^{\prime
2}(\tau)+\\!\bigg{(}\\!\frac{\partial\hat{b}_{\text{R}}}{\partial\langle\mathcal{O}(\tau)\rangle}\\!\bigg{)}^{2}\\!\Delta\hat{\mathcal{O}}^{2}(\tau)$
(15) $\displaystyle\\!\\!=\\!$ $\displaystyle\frac{1}{T\tau
N^{2}}\bigg{[}e^{2\gamma_{\text{GHZ}}(\tau)}-\frac{1}{2}\sin^{2}(2Nb\tau)\bigg{]},$

to be contrasted with Eq. (12). Interestingly, in this case the optimal
performance is achieved when the phase $\varphi^{\prime}\equiv 2Nb\tau$ is
$\varphi=(2k+1)\pi/4$. In Fig. 2 (bottom), the analytic expression
$\Delta\hat{b}_{\text{R}}^{2}(\tau)$ in Eq. (15) is compared against the
sample variance for finite $\nu$, showing remarkable agreement in the region
where the estimator is linear. Note that the sample variance diverges around
the points where the estimator is discontinuous, forcing our prior knowledge
of the frequency to be confined to an interval of length $\pi/(N\tau)$.

Figure 2: (Color online) Performance of ratio estimator for collective
dephasing. Top: Ratio estimator’s sample mean for $\nu=30$ (red, solid),
$\nu=400$ (blue, dashed), and analytic $\nu\rightarrow\infty$ mean value
(grey, solid). Bottom: Ratio estimator’s sample variance for $\nu=400$ and
limiting analytic expression (grey, solid). In both cases, a GHZ state of
$N=100$ qubits is considered, subject to collective spin-boson dephasing from
a 1D zero-temperature noise environment; a supra-Ohmic spectral density with
$s=3,\omega_{c}=1$ and noise strength $\alpha=1$ is used, and the measurement
time $\tau\approx 0.0028$ is chosen to minimize uncertainty.

In order to compare the performance of $\Delta\hat{b}_{\text{R}}(\tau)^{2}$
against the standard figure of merit of Eq. (12), for consistency it is
necessary to let $T\mapsto 2T$ in the latter, so that the same fixed resources
are ensured. It can be checked, by plotting both uncertainties as a function
of time, that $\Delta\hat{b}_{\text{R}}(\tau)$ is slightly bigger than its
standard counterpart at all times (see inset of Fig. 3(left)). The optimal
interrogation time for the ratio estimator is harder to compute analytically.
However, evaluating the ratio estimator uncertainty at the time
$\tau^{\text{GHZ}}_{\text{opt, coll}}$ that minimizes the standard one, we can
verify that it reaches the same SQL scaling as its standard counterpart,
$\Delta\hat{b}_{\text{R}}(\tau^{\text{GHZ}}_{\text{opt,
coll}})=(2\sqrt{e}-1)^{1/2}(\kappa_{0}\omega_{c}/T)^{1/2}N^{-1/2}$, albeit
with a slightly larger constant prefactor ($(2\sqrt{e}-1)^{1/2}\approx 1.52$
as opposed to $e^{1/4}\approx 1.28$).

A similar analysis for the more general case of non-collective spin-boson
dephasing will be carried out in the context of the noise-tailored lattice
sensor of Sec. IV.

###  III.3 Ratio estimator for general Ramsey interferometry

While the emergence of noise-induced bias was illustrated for GHZ states in
the above discussion, standard estimators used in Ramsey interferometry with
different initial states are plagued by similar issues in the presence of
dephasing noise. Here, we focus on dephasing noise from a bosonic bath and a
class of OAT states, as defined in Eq. (8), which includes CSS product states
as a limiting case ($\theta=\beta=0$) as well as OATs with squeezing and
rotation angles set to minimize the initial variance along $y$. In the $N\gg
1$ limit, we have [43, 33]:

$\displaystyle\theta_{\text{opt}}(N)$ $\displaystyle\approx
12^{1/6}\,2^{2/3}\,N^{-2/3}$ (16) $\displaystyle\beta_{\text{opt}}(N)$
$\displaystyle\approx\pi/2-3^{-1/6}N^{-1/3}-3^{1/6}2^{-1}\;N^{-2/3}.$

After an encoding time $\tau$, the mean value $\langle J_{y}(\tau)\rangle$ can
be shown to be given by an expression of the form (see Appendix A.2.1 for more
detail)

$\displaystyle\langle
J_{y}(\tau)\rangle=e^{-\gamma_{\text{OAT}}(\tau)}\mathcal{Q}(\tau)\sin(b\tau),$
(17)

where $\gamma_{\text{OAT}}(\tau)$ depends on decay coefficients
$\kappa_{nm}(\tau)$, $\mathcal{Q}(\tau)$ is a function of both the angles
$\theta,\beta$ and, for a nonclassical environment, the coefficients
$\xi_{nm}(\tau)$ determine the phase evolution $\varphi_{0}(\tau)$.

For a product state, it is possible to obtain the exact expression
$\mathcal{Q}(\tau)=\sum_{n}\prod_{m\neq n}\cos(\,\xi_{nm}(\tau))$, whereas for
arbitrary squeezing and twisting angles the reduced dynamics become highly
entangled and an exact analytic approach is unfeasible. However, in the short-
time regime $\omega_{c}\tau\ll 1$ where the relevant evolution takes place,
the above mean values can be computed to great accuracy by performing a
second-order cumulant expansion over the qubit operators [32, 33]. In this
same limit, any quantum noise contribution due to $\varphi_{0}(t)$ enters as a
correction to the decoherence dynamics determined by the decay terms
$\kappa_{nm}(t)$. Thus, it may be disregarded to first approximation, making
$\mathcal{Q}(\tau)$ approximately time-independent, say,
$\mathcal{Q}_{0}\equiv\mathcal{Q}_{0}(\theta,\beta)$. Inverting this
relationship yields

$\hat{b}(\tau)=\frac{1}{\tau}\arcsin[e^{\gamma_{\text{OAT}}}\langle\hat{J}_{y}(\tau)\rangle\mathcal{Q}_{0}^{-1}],$

with $\langle\hat{J}_{y}(\tau)\rangle=\sum_{m=-N/2}^{N/2}(\nu_{m}/\nu)$ being
the sample average of $\langle J_{y}(\tau)\rangle$ after $\nu$ detections.
Similar to the GHZ case, evaluating this estimator would require precise
knowledge of the noise parameters $\gamma_{\text{OAT}}(\tau)$,
$\mathcal{Q}_{0}$, and may yield an imaginary number for the subset of
outcomes where
$e^{\gamma_{\text{OAT}}(\tau)}\langle\hat{J}_{y}(\tau)\rangle\mathcal{Q}_{0}>1$
in the argument of the arcsin. We now show how this issue can again be
circumvented by constructing a ratio estimator after measuring two distinct
observables – in this case, the two orthogonal angular momentum components
$J_{x}$ and $J_{y}$.

Since, from Eq. (17), the target frequency $b$ enters the reduced dynamics
unitarily, say, $U(\varphi)=\exp(-i\varphi J_{z})$, with $\varphi=bt$, we can
derive the equations of motion with respect to $\varphi$ for $\langle
J^{\ell}_{v}(\varphi)\rangle$, $v\in\\{x,y\\}$, for the two relevant moments
$\ell\in\\{1,2\\}$. For $\ell=1$, taking the derivatives in $\langle
J_{v}(\varphi)\rangle\equiv\text{Tr}\\{J_{v}U(\varphi)\rho_{0}(t)U(\varphi)^{\dagger}\\}$,
with $\rho_{0}(t)$ depending parametrically on the noisy dynamical
coefficients, we obtain:

$\displaystyle\partial_{\varphi}\langle J_{x}(\varphi)\rangle$
$\displaystyle=-i\langle[J_{x},J_{z}](\varphi)\rangle=\langle
J_{y}(\varphi)\rangle,$ $\displaystyle\partial_{\varphi}\langle
J_{y}(\varphi)\rangle$
$\displaystyle=-i\langle[J_{y},J_{z}](\varphi)\rangle=-\langle
J_{x}(\varphi)\rangle,$

which immediately leads to $\partial^{2}_{\varphi}\langle
J_{v}(\varphi)\rangle=-\langle J_{v}(\varphi)\rangle$. Thus,

$\displaystyle\langle J_{x}(t)\rangle$
$\displaystyle=\cos(\varphi)\,\mathcal{F}(t)+\sin(\varphi)\,\mathcal{G}(t),$
$\displaystyle\langle J_{y}(t)\rangle$
$\displaystyle=-\sin(\varphi)\,\mathcal{F}(t)+\cos(\varphi)\,\mathcal{G}(t),$

for some state-dependent functions $\mathcal{F}(t),\mathcal{G}(t)$ which fully
capture the noise effects. If we impose the initial conditions $\langle
J_{x}(0)\rangle\equiv\mathcal{J}_{x}$, $\langle J_{y}(0)\rangle=0$, which hold
for the CSS and OAT states we are interested in, we are led to the following
ratio estimator:

$\hat{b}_{\text{R}}^{\text{OAT}}(\tau)=\frac{1}{\tau}\arctan[\langle\hat{J}_{y}(\tau)\rangle/\langle\hat{J}_{x}(\tau)\rangle].$
(18)

Again, implementing the same number of measurements $\nu=T/\tau$ for each
observable doubles the time resources from $T$ to $2T$. Linearizing the
expressions in the asymptotic limit $\nu\gg 1$ and high prior knowledge yields
the desired expression for the squared uncertainty:

$\displaystyle\Delta\hat{b}_{\text{R}}^{\text{OAT}}(\tau)^{2}$
$\displaystyle\\!=\\!\bigg{(}\\!\frac{\partial\hat{b}_{\text{R}}}{\partial\langle
J_{x}(\tau)\rangle}\\!\bigg{)}^{2}\\!\Delta\hat{J}_{x}^{2}(\tau)+\\!\bigg{(}\\!\frac{\partial\hat{b}_{\text{R}}}{\partial\langle
J_{y}(\tau)\rangle}\\!\bigg{)}^{2}\\!\Delta\hat{J}_{y}^{2}(\tau)$
$\displaystyle=\frac{\langle J_{x}(\tau)\rangle^{2}\Delta
J_{y}^{2}(\tau)+\langle J_{y}(\tau)\rangle^{2}\Delta
J_{x}^{2}(\tau)}{T\tau[\langle J_{x}(\tau)\rangle^{2}+\langle
J_{y}(\tau)\rangle^{2}]^{2}}.$ (19)

This provides a ratio estimator useful to carry out asymptotically unbiased
quantum sensing in the presence of dephasing, starting from a class of
squeezed states with $\langle J_{x}(0)\rangle\neq 0$ and $\langle
J_{y}(0)\rangle=0$.

Figure 3: Ratio estimator comparative performance for collective spin-boson
dephasing.  Left: GHZ optimal uncertainty vs. qubit number for the standard
(grey circles) and ratio estimator (blue triangles). Both saturate the SQL,
with the standard uncertainty outperforming
$\Delta\hat{b}^{\text{GHZ}}_{\text{R opt}}$ by a constant $\text{C}\approx
0.85$ factor. Inset: GHZ $N=100$ qubit uncertainty as a function of time for
standard (grey, solid) and ratio estimator (blue, solid). Right: CSS optimal
uncertainty vs. qubit number for the standard (grey, circles) and ratio
estimator (blue, triangles). Both estimators show $N^{-1/4}$ scaling with the
same constant factor when evaluated at their respective optimal phases. Inset:
CSS $N=100$ qubit uncertainty as a function of time for standard (grey, solid)
and ratio estimator (blue, dashed). The limit of no quantum noise,
$\xi(\tau)=0$, for $\Delta\hat{b}(\tau)$ in Eq. (21) (orange, dotted) is
included for comparison as a lower bound. Collective spin-boson dephasing
noise from a 1D zero-temperature environment, and a supra-Ohmic spectral
density as in Fig. 2 is assumed.

The expression in Eq. (19) can be simplified by obtaining the $\varphi$
dependence of $\langle J_{x}^{2}(\varphi)\rangle$ and $\langle
J_{y}^{2}(\varphi)\rangle$. Taking first- and second-order derivatives with
respect to $\varphi$ and solving the resulting partial differential equation
system leads to the following mean values:

$\displaystyle\langle J_{x}^{2}(\varphi)\rangle$
$\displaystyle=\mathcal{K}(t)+\mathcal{A}(t)\cos(2\varphi)+\mathcal{B}(t)\sin(2\varphi),$
$\displaystyle\langle J_{y}^{2}(\varphi)\rangle$
$\displaystyle=\mathcal{K}(t)-\mathcal{A}(t)\cos(2\varphi)-\mathcal{B}(t)\sin(2\varphi).$

Here, $\mathcal{K}(t)$, $\mathcal{A}(t)$ and $\mathcal{B}(t)$ are state-
dependent quantities that capture the noise effects. Enforcing the initial
conditions $\langle J_{x}^{2}(0)\rangle=\mathcal{J}_{x}^{2}$ , $\langle
J_{y}^{2}(0)\rangle=\mathcal{J}_{y}^{2}$ and
$\langle[J_{x},J_{y}]_{+}(0)\rangle=0$, again valid for the CSS and OAT states
considered here, we get $\mathcal{B}(t)=0$. Substituting in Eq. (19) yields

$\displaystyle\Delta\hat{b}_{\text{R}}^{\text{OAT}}(\tau)^{2}\\!=\\!\frac{2\mathcal{K}(\tau)-\mathcal{A}[\tau)(1+\cos(4\varphi)]-\mathcal{F}(\tau)^{2}\sin(2\varphi)^{2}}{2T\tau\mathcal{F}(\tau)^{2}}.$

On the other hand, the precision attained by using the standard method of
moments for a set of measurements taking the same fixed total time $2T$ is
given by:

$\displaystyle\Delta\hat{b}(\tau)^{2}$ $\displaystyle=\frac{\Delta
J_{y}^{2}(\tau)}{2T\tau(\partial_{\varphi}\langle J_{y}(\tau)\rangle)^{2}}$
$\displaystyle=\frac{\mathcal{K}(\tau)+\mathcal{A}(\tau)\cos(2\varphi)-(\mathcal{F}(\tau)\sin(\varphi))^{2}}{2T\tau\mathcal{F}(\tau)^{2}\cos(\varphi)^{2}}.$
(20)

Let us now compare the performance of both uncertainties for a product CSS
state along $x$ undergoing collective spin-boson dephasing, as considered for
the GHZ state. In this scenario the relevant mean values can be computed
exactly [32, 33], from which the quantities
$\mathcal{F}(t),\mathcal{K}(t),\mathcal{A}(t)$ are easily deduced.
Substituting into Eqs. (19)-(20) at their respective optimal phases
$\varphi=\pi(2n+1/4)$ and $\varphi=n\pi,n\in\mathbb{N}$, for fixed total time
$2T$, we find for the CSS squared uncertainties arising from the ratio
estimator and the method of moments:

$\displaystyle\Delta\hat{b}_{\text{R}}^{\text{CSS}}(\tau)^{2}$
$\displaystyle=\frac{1}{2T\tau}\bigg{[}\frac{e^{\kappa(\tau)}(N+1)}{2N\,\cos(\xi(\tau))^{2N-2)}}-1\bigg{]},$
(21) $\displaystyle\Delta\hat{b}^{\text{CSS}}(\tau)^{2}$
$\displaystyle=\frac{(N+1)e^{\kappa(t)}-(N-1)e^{-\kappa(t)}\cos(2\xi(\tau))^{N-2}}{2N\,(2T)\,\tau\,\cos(\xi(\tau))^{2N-2}}.$

For the standard estimation setting, the optimal measurement time and
precision for $\Delta\hat{b}^{\text{CSS}}(\tau)$ have been computed to be [32]
$\tau^{\text{CSS}}_{\text{opt, coll}}=(\kappa_{0}\,\omega_{c})^{-1}N^{-1/2}$
and $\Delta\hat{b}^{\text{CSS}}_{\text{opt,
coll}}=(\kappa_{0}\,\omega_{c}/T)^{1/2}N^{-1/4}$, respectively. For the ratio
estimator, we can carry out an expansion in the short-time limit to derive the
optimal measurement time, $\tau^{\text{CSS}}_{\text{R opt,
coll}}=(\kappa_{0}\,\omega_{c})^{-1}N^{-1/2}=\tau^{\text{CSS}}_{\text{opt,
coll}}$, and performance $\Delta\hat{b}^{\text{CSS}}_{\text{ R opt,
coll}}=(\kappa_{0}\,\omega_{c}/T)^{1/2}N^{-1/4}$. In line with our findings
for the GHZ state, the ratio estimator optimal precision follows the same
$N$-scaling as its standard counterpart; furthermore, for a CSS it is
multiplied by the same constant factor and hence attains identical precision:
$\Delta\hat{b}^{\text{CSS}}_{\text{ R opt
coll}}=\Delta\hat{b}^{\text{CSS}}_{\text{opt coll}}$. This is shown in Fig. 3,
where we plot the optimal uncertainty obtained by numerical optimization as a
function of qubit number for both estimators and CSS as well as GHZ initial
states evaluated at their respective ideal phases. The results confirm that
the analytic predictions for both scaling and constant factors are valid to
great accuracy.

##  IV Noise-tailored lattice sensor

###  IV.1 Spin-boson model

As a concrete illustrative setting, we consider a 1D dephasing spin-boson
model, whereby the qubits interact with a collection of oscillators vibrating
at frequencies $\Omega_{k}$, in thermal equilibrium at inverse temperature
$\beta$. Then,

$B_{n}(t)=\sum_{k}(g_{k}\,e^{ikr_{n}}e^{i\Omega_{k}t}\,b_{k}^{\dagger}+\text{H.c.}),$

in terms of bosonic operators $b_{k}$, $b^{\dagger}_{k}$, with
$g_{k}\in{\mathbb{R}}$ being a coupling strength for momentum mode $k>0$. We
assume a linear dispersion $\Omega_{k}\equiv vk$, with $v>0$ a speed
parameter, and envision the qubits to be arranged in a regular lattice with
tunable spacing $x_{0}$, that is, the position of qubit $n$ obeys
$r_{n}=nx_{0}$, $1\leq n\leq N$. By contrast, the collective limit analyzed in
the previous section is obtained when all qubits share the same position,
$r_{n}=x_{0},\,\forall n$, leading to _perfect spatial correlations_ which
cannot be controlled. For a qubit pair $n,m$, a _transit time_ proportional to
their distance may then be defined by $t_{nm}\equiv|n-m|x_{0}/v$.

Figure 4: Noise-optimized superclassical precision scaling. Left: GHZ optimal
uncertainty vs. qubit number. Grey circles: Exact numerical optimization.
Blue, solid line: Analytic expression for
$\Delta\hat{b}^{\text{GHZ}}_{\text{opt}}$. Orange, dashed: Asymptotic scaling
limit, Eq. (25). Inset: GHZ uncertainty vs. time for $N=100$ qubits and
optimal lattice spacing $x_{0}=0.4296$. Grey, solid: Exact expression. Blue,
dashed: Short-time approximate uncertainty. Right: Optimal uncertainty vs.
qubit number for an OATS ideally squeezed along $y$ before evolution. Blue
triangles: Numerical optimization for purely classical noise. Gray circles:
Numerical optimization including quantum noise. Orange, dashed: Asymptotic
scaling limit, Eq. (26). Inset: OATS uncertainty vs. time for $N=30$ qubits
and optimal spacing $x_{0}=0.46$. Blue, dashed: Purely classical dephasing.
Grey, solid: Classical and quantum dephasing. Orange, dotted: Short-time
approximate uncertainty. All other parameters are as in Fig. 2.

Notably, for the above setting we have $\varphi_{1}(t)\equiv 0$ under standard
physical assumptions on the coupling [49, 33]. The dynamic coefficients
determining $\gamma(t)$ and $\varphi_{0}(t)$ are given in Eq. (3) and Eq. (5),
respectively. To compute $\kappa_{nm}(t)$ and $\xi_{nm}(t)$, we further assume
a continuum of bosonic modes, characterized by a spectral density of the form

$J(\omega)\equiv\alpha\omega_{c}(\omega/\omega_{c})^{s}e^{-\omega/\omega_{c}},\quad\alpha>0,\>s>0,$
(22)

with $s$ being the so-called Ohmicity parameter. In a low-temperature regime
where thermal effects may be taken to be negligible
($\coth(\beta|\omega|/2\approx 1$), the noise classical and quantum spectra
are then calculated as

$S^{+}_{nm}(\omega)=4\pi J(\omega)\cos(\omega
t_{nm}),\;\;S^{-}_{nm}(\omega)=S^{+}_{nm}(\omega)\text{sgn}(\omega),$

where $\text{sign}(\omega)$ is the sign function.

The cutoff frequency $\omega_{c}$ defines a metrologically relevant short-time
limit, given by $\omega_{c}t\ll 1$, where the quantum advantage in estimation
precision may be maximized [15, 32, 33]. In this regime, the dynamic
coefficients entering $\gamma(t)$ and $\varphi_{0}(t)$ can be approximated as

$\displaystyle\kappa_{nm}(t)$ $\displaystyle\approx$
$\displaystyle\kappa_{0}^{2}(\omega_{c}t)^{2}\delta_{1}(|n-m|x_{0}),\quad\kappa_{0}^{2}\equiv\alpha\,\Gamma({s}+1),$
$\displaystyle\xi_{nm}(t)$ $\displaystyle\approx$
$\displaystyle\xi_{0}^{3}(\omega_{c}t)^{3}\delta_{2}(|n-m|x_{0}),\quad\xi_{0}^{3}\equiv({\alpha}/{6})\,\Gamma(s+2),$

with $\Gamma$ being the Euler Gamma function. Their temporal dependence is
thus factored from the relevant spatial correlations, which are encoded in the
functions

$\displaystyle\delta_{1}(|n-m|x_{0})$ $\displaystyle\equiv$
$\displaystyle\;u^{{s}+1}\,T_{{s}+1}(u),$ (23)
$\displaystyle\delta_{2}(|n-m|x_{0})$ $\displaystyle\equiv$
$\displaystyle\;u^{s+2}\,T_{s+2}(u),$ (24)

where $u\equiv[1+(|n-m|x_{0})^{2}]^{-1/2}$ and one can show that $T_{n}(x)$
denotes the $n$-th order Chebyshev polynomial of the first kind. Provided that
knowledge of the spectral density $J(\omega)$, hence of the spatial
correlation functional forms $\delta_{\ell}(x)$ is available, and that
$J(\omega)$ exhibits supra-Ohmic low-frequency behavior ($s>1$), spatial
correlations with $n\neq m$ can be made _negative_ by tuning the distance
$x_{0}$. We now show how this may engender a drastic improvement in the
scaling of estimation precision with respect to both collective noise [32] and
our previous randomized Ramsey protocol [33].

###  IV.2 Initial GHZ state

Figure 5: Ratio estimator uncertainty optimal performance for different
Ohmicity parameters.  Minimal uncertainty as function of $N$ of a GHZ (left)
and OAT state (right) for different spectral density with a high-frequency
exponential cutoff and different Ohmicity parameters $s=0$, (solid lilac
stars) $s=1$ (solid green diamonds), $s=2$ (hollow orange squares), $s=3$
(hollow grey circles), $s=4$ (hollow blue triangles), and $s=5$ (hollow brown
rectangles). The solid black line describes the asymptotic lower bound of
$3^{1/2}\Gamma(3)^{1/4}(\omega_{c}/T)^{1/2}N^{1/2}$ and
$3/4N^{-1}\sqrt{\log(N)}$ for OATS and GHZ, respectively. Here,
$\alpha=1=\omega_{c}$, as in previous figures.

Since, despite being non-collective, the noise remains pure dephasing, the
same analysis carried out in Sec. III.2 holds here, thus we may estimate the
target frequency without incurring in noise-induced bias through Eq. (13). The
resulting precision, however, is different, with the dephasing rate
$\gamma_{\text{GHZ}}(t)$ no longer being that of the collective setting. To
determine the GHZ optimal precision, we not only minimize the uncertainty with
respect to time, but also the lattice parameter $x_{0}$, for optimal phase
$Nbt=(2k+1)\,\pi/4$, $k\in\mathbb{N}$. Since $\varphi_{0}(t)\equiv 0$ for
$|\text{GHZ}\rangle$ by symmetry, quantum correlations do not play any role in
this case ($\delta_{2}\equiv 0$ in Eq. (24)). Taking advantage of the fact
that, in the short-time limit,

$\gamma_{\text{GHZ}}(t)\\!\approx\\!(\omega_{c}t)^{2}F_{N}(x_{0}),\;\>F_{N}(x_{0})\equiv\kappa_{0}^{2}\sum_{n,m}^{N}\delta_{1}(x_{0}|n-m|),$

the minimization may be carried out with respect to each variable separately.
Substituting in Eq. (15) and optimizing with respect to $\tau$ leads to the
optimal measurement time, $\tau_{\text{R
opt}}^{\text{GHZ}}=\frac{1}{2}\omega_{c}^{-1}F_{N}(x_{0})^{-1/2}$. It follows
that the best sensing performance is achieved by minimizing the spatial
function $F_{N}(x_{0})$ for the time-optimized uncertainty, leading to

$\Delta\hat{b}_{\text{R}}(\tau^{\text{GHZ}}_{\text{R opt}})\approx
2.96\;{\omega_{c}/T}^{1/2}\,F_{N}(x_{0})^{1/4}\,N^{-1}.$

While the details are provided in Appendix A.1 (see in particular Eq. (28) for
the origin of the above numerical prefactor), the resulting approximate
(analytically) minimized spatial function, $F_{N}(x_{0\text{
opt}}^{\text{GHZ}})$, can be shown to scale logarithmically in the large-$N$
limit: $F_{N}(x_{0\text{ opt}}^{\text{GHZ}})\propto\mathcal{O}(\log(N)^{2})$.
Accordingly, the optimal asymptotic sensing performance is

$\Delta\hat{b}^{\text{GHZ}}_{\text{R opt }}\approx
2.96\,(\omega_{c}/T)^{1/2}\sqrt{\log(N)}\,N^{-1},$ (25)

which is _closer to Heisenberg scaling than any power law_.

Figure 4(left) demonstrates that the agreement between the analytic expression
in Eq. (25) and the exact numerical optimization over $x_{0},\tau$ is
excellent even at finite $N\gtrsim 20$, despite the fact that discrepancies
between the lattice parameter that numerically minimizes $\delta_{1}$ and its
analytic approximation $x_{0\text{ an}}^{\text{GHZ}}$ vanish only
asymptotically in $N$. This indicates a high degree of robustness against
deviations from the exact optimal lattice parameter and, in turn, against
uncertainty in the characterization of the underlying noise spectral density.

###  IV.3 Initial OAT state

Unlike GHZ states, OATS are not immune to the genuinely quantum contribution
of the noise. For clarity, however, let us first consider the classical decay
contribution, described as before by Eqs. (3)-(4), and assess the impact of
the phase $\varphi_{0}(t)$ at a later stage. As detailed in Appendix A.2.2, by
performing a cumulant expansion over appropriate qubit operators [32], we may
derive approximate expressions for the moments $\langle J_{v}(t)\rangle$,
$\langle J_{v}^{2}(t)\rangle$, $v=x,y$, which are remarkably accurate in the
short-time limit $\omega_{c}t\ll 1$. Substituting in Eq. (19), it is then
possible to expand the uncertainty in powers of $t$,

$\Delta\hat{b}_{\text{R}}(t)\approx\big{(}\sqrt{Tt}\,h_{0}(N)\big{)}^{-1}\\!\bigg{(}\sum_{k=0,1,2}\\!a_{2k}(N,x_{0})\,(\omega_{c}t)^{2k}\bigg{)}^{\\!1/2},$

where the coefficients $h_{0}(N)$ and $a_{2k}(N,x_{0})$ can be evaluated
analytically by the same approach used to compute the spatial function
$F_{N}(x_{0})$ for the GHZ state. This expansion is accurate in the region
where the minimum occurs (see inset of Fig. 4(right)), and once again it
allows one to separate the spatial and temporal dependence of the uncertainty
into two distinct contributions. The best performance can then be extracted by
minimizing the expansion with respect to time and, subsequently, with respect
to the lattice parameter. For $s>1$, as detailed in Appendix A.2.2, we find
that the noise-tailored sensor reaches the Zeno limit,

$\Delta\hat{b}_{\text{R
opt}}^{\text{OATS}}\approx{{\omega_{c}}/{(3T)}}^{1/2}\,\Gamma(s+1)^{1/4}N^{-3/4}.$
(26)

This is a mere $N^{-1/12}$ scaling loss with respect to the noiseless OAT
scenario, and is found to be in full agreement with numerical minimization,
see Fig. 4(right).

Accounting for the effect of quantum noise due to $\varphi_{0}(t)$ makes the
expression for the mean values significantly more involved. This prevents us
from deriving a tractable short-time expansion for the uncertainty and
lengthens the computational time required to evaluate
$\Delta\hat{b}_{\text{R}}^{\text{OAT}}(\tau)$ numerically, limiting in turn
the accessible values of $N$. Still, as one may see from the inset in Fig.
4(right), the uncertainty as a function of time has a nearly identical
behavior as in the presence of classical noise alone, with the contribution of
$\varphi_{0}(t)$ entering as a small correction. In line with the above,
numerical minimization in the range of $N$ we were able to test leads to
similar values for the optimal uncertainty as when quantum noise was
disregarded, see Fig. 4(right). Such a behavior is plausible in light of the
fact that, as in the GHZ case, $\xi_{nm}(t)\propto(\omega_{c}t)^{3}$, as
compared to $\kappa_{nm}(t)\propto(\omega_{c}t)^{2}$, which causes the decay
contributions to dominate in the short-time limit. Altogether, this strongly
suggests that a Zeno limit as in Eq. (26) is realized asymptotically, also
when quantum noise is fully included.

###  IV.4 Parameter robustness

While the calculations reported thus far for the spin-boson setting have been
carried out by assuming a spectral density with Ohmicity parameter $s=3$ and
an exponential high-frequency decay, that is,
$J(\omega)=\alpha\,\omega_{c}^{-2}\,\,\omega^{3}\,e^{-\omega/\omega_{c}}$, it
is important to assess how sensitively the resulting precision scaling depend
upon the details of the noise spectra. Numerical optimization indicates that
the same nearly-Heisenberg uncertainty scaling and Zeno limit derived for GHZ
and OAT states, respectively, hold for other supra-ohmic densities, with the
$N$-independent prefactor increasing with $s$, see Fig. 5. Furthermore, a
linear fit shows that for OAT states, the prefactor is $\Gamma(s+1)/\sqrt{3}$,
which was only rigorously derived for $s=3$.

Figure 5 also shows that for an Ohmic ($s=1$) spectral density, however, the
optimal ratio uncertainty experiences a small scaling loss for both GHZ (left)
and OAT (right) states, which manifests as a slight slope decrease with
respect to the supra-Ohmic cases in the log-log plot. Finally, for a sub-Ohmic
spectral density ($s=0$), such a change in slope and scaling loss are more
pronounced, especially for OAT states. This is in line with the expectation
that negative spatial correlations play a diminished role in the presence of
non-vanishing zero-frequency noise. Understanding in more detail the precision
achievable in the Ohmic and sub-Ohmic scenarios, as well as addressing
spectral densities with non-exponential (e.g., Gaussian) high-frequency
cutoffs or non-trivial spectral features (e.g., high-frequency peaks), remains
a task worth of additional investigation.

##  V Conclusion

We have shown how to construct ratio estimators for entanglement-assisted
frequency estimation by Ramsey interferometry which, at the cost of doubling
the required measurement resources, overcome noise-induce bias and retain the
same precision scaling of standard local estimators in the asymptotic
large-$N$ regime. While our construction is applicable to both Markovian and
non-Markovian dephasing settings, and requires no prior detailed noise
knowledge, we have additionally shown that asymptotically unbiased frequency
estimation with noise-optimized superclassical precision is possible by
suitable sensor’s engineering and tuning of the noise spatial correlations. In
particular, by focusing on spin-boson dephasing in a one-dimensional lattice
geometry with tunable qubit-probe separation, we established a novel scaling
which is a $\log(N)^{1/2}$ factor away from the Heisenberg limit for a GHZ
state, and found that Zeno scaling is reachable with a properly squeezed OAT
state.

An outstanding question for future work is to explore potential experimental
realizations of our tunable lattice-sensor model, for instance in trapped ion
settings [50, 51] or NV centers in diamond [52]. From a theoretical
standpoint, our results hint at the possibility that the Heisenberg limit may
be reachable by further optimizing the initial entangled state, in addition to
the sensor’s geometry. Likewise, it would be interesting to extend the
applicability of our noise-unbiased estimators to generalized Ramsey echo-
based protocols, for which the symmetry condition $\langle J_{y}(0)\rangle=0$
of the input states is relaxed, and pre-measurement anti-squeezing operations
are allowed [53, 54]. Finally, assessing the usefulness of our ideas in a non-
asymptotic metrology setting with a limited amount of data [37, 55], where a
Bayesian estimation paradigm is appropriate, and noise-induced bias competes
with finite-sample effects, may be especially important to expand the
practical reach of quantum metrology applications.

## Acknowledgements

It is a pleasure to thank Augusto Smerzi for valuable input and
clarifications, and Félix Beaudoin for a critical reading of the manuscript.
L.V. also acknowledges early discussions with Félix Beaudoin on the issue of
noise-induced bias. Work at Dartmouth was supported by the US NSF through
Grant No. PHY-2013974, and the Constance and Walter Burke Special Projects
Fund in QIS. Work at Griffith was supported (partially) by the Australian
Government through the Australian Research Council’s Discovery Projects
funding scheme (project No. DP210102291).

*

##  Appendix A Additional detail on technical derivations

###  A.1 Optimal noise-tailored precision scaling for initial GHZ states

Here we derive the GHZ optimal uncertainty achievable by our ratio estimator
in the presence of spin-boson noise with tunable spatial correlations. Our
starting point is Eq. (15), evaluated for optimal phase $\varphi=\pi/4$ in the
short-time regime:

$\Delta\hat{b}_{\text{R}}^{\text{GHZ}}(\tau)^{2}=\frac{1}{N^{2}T\tau}\bigg{(}e^{2\gamma_{\text{GHZ}}(\tau)}-\frac{1}{2}\bigg{)}\approx\frac{1}{N^{2}T\tau}\bigg{(}e^{2\kappa_{0}^{2}(\omega_{c}\tau)^{2}\sum_{n,m=1}^{N}\delta_{1}(x_{0}|n-m|)}-\frac{1}{2}\bigg{)}.$
(27)

One may first minimize with respect to time, leading to the optimal
measurement time and the corresponding time-optimized uncertainty. The result
may be cast in the following form:

$\tau_{\text{R
opt}}^{\text{GHZ}}=\frac{1}{2}\omega_{c}^{-1}F_{N}(x_{0})^{-1/2}\sqrt{1+\frac{1}{2}W(e^{-1/2}/4)},\quad\>\Delta\hat{b}_{\text{R}}(\tau_{\text{R
opt}}^{\text{GHZ}})=\sqrt{2}\,\frac{(1+2\,W(e^{-1/2}/4))}{\sqrt{-W(e^{-1/2}/4)}}\;\sqrt{\frac{\omega_{c}}{T}}\;F_{N}(x_{0})^{1/4}\;N^{-1},$
(28)

where $F_{N}(x_{0})=\kappa_{0}^{2}\,\sum_{n,m=1}^{N}\delta_{1}(x_{0}|n-m|)$ is
the function of the lattice spacing $x_{0}$ also defined in the main text, and
$W(x)$ is the principal branch of the Lambert’s $W$ function [56]. Recall that
the Lambert’s $W$ function $W_{k}(z)$ is defined, given complex numbers $w$
and $z$, by the equation $we^{w}=z$, which holds if and only if $w=W_{k}(z)$
for some integer $k$. For simplicity, we have dropped the subscript $k$ for
the principal branch $k=0$. The numerical prefactor in Eq. (28) can be
computed to be $\sqrt{2}\;(1+2\,W(e^{-1/2}/4))(-W(e^{-1/2}/4))^{-1/2}\approx
2.96$, as quoted in the main text above Eq. (25).

In order to obtain the best precision, we must now minimize $F_{N}(x_{0})$
with respect to $x_{0}$. Using the symmetry of a 1D regular lattice, we may
write the spatial function $F_{N}(x_{0})$ in terms of a single index sum, and
a term including the spatial “self-correlations” $\delta_{1}(0)$ corresponding
to the $n=m$ case, that is,

$F_{N}(x_{0})=\kappa_{0}^{2}\bigg{[}N\delta_{1}(0)+2\sum_{j=1}^{N-1}(N-j)\delta_{1}(j\,x_{0})\bigg{]}.$
(29)

Replacing the explicit functional form of the spatial correlations
$\delta_{1}(|n-m|x_{0})$ for $s=3$ into $F_{N}(x_{0})$ yields

$F_{N}(x_{0})=\kappa_{0}^{2}\bigg{[}N+2\sum_{j=1}^{N-1}(N-j)\,\frac{1-6\,(j\,x_{0})^{2}+(j\,x_{0})^{4}}{(1+(j\,x_{0})^{2})^{4}}\bigg{]}=\kappa_{0}^{2}\left[N+S_{N}(x_{0})\right],$
(30)

with $\kappa_{0}^{2}=\alpha\,\Gamma(4)$ a dimensionless constant and
$S_{N}(x_{0})$ consisting of two finite sums:

$S_{N}(x_{0})=2N\sum_{j=1}^{N-1}\,\frac{1-6\,(j\,x_{0})^{2}+(j\,x_{0})^{4}}{(1+(j\,x_{0})^{2})^{4}}-2\sum_{j=1}^{N-1}j\,\frac{1-6\,(j\,x_{0})^{2}+(j\,x_{0})^{4}}{(1+(j\,x_{0})^{2})^{4}}.$

To evaluate these series, we now show that they can be written as linear
combinations of Polygamma functions, defined by
$\psi^{(m)}(z)\equiv(-1)^{m+1}\int_{0}^{\infty}\frac{t^{m}e^{-zt}}{1-e^{-t}}\,dt$.
The key step is to expand $(1-e^{-t})^{-1}$ in the integrand that defines
$\psi^{(m)}(t)$ by using the geometric series, and summing each contribution
term-by-term. For example, for

$\text{Re}\,\psi^{(3)}\left(\frac{i+x_{0}}{x_{0}}\right)=\frac{1}{2}\left[\psi^{(3)}\\!\left(\frac{i+x_{0}}{x_{0}}\right)+\psi^{(3)}\\!\left(\frac{-i+x_{0}}{x_{0}}\right)\right],$

we have

$\displaystyle\text{Re}\,\psi^{(3)}\bigg{(}\frac{i+x_{0}}{x_{0}}\bigg{)}\\!=\\!\int_{0}^{\infty}\\!\frac{t^{3}e^{-t}}{1-e^{-t}}\cos\Big{(}\frac{t}{x_{0}}\Big{)}dt=\sum_{j=1}^{\infty}\int_{0}^{\infty}\\!\\!t^{3}\;e^{-jt}\cos\Big{(}\frac{t}{x_{0}}\Big{)}dt=6x_{0}^{4}\;\sum_{j=1}^{\infty}\frac{1-6\,(j\,x_{0})^{2}+(j\,x_{0})^{4}}{(1+(j\,x_{0})^{2})^{4}},\qquad$
(31)

where integrating each series coefficient independently is allowed by Fubini’s
theorem. Similarly,

$\displaystyle\text{Re}\;\psi^{(3)}\bigg{(}\frac{i}{x_{0}}+N\bigg{)}=6x_{0}^{4}\sum_{j=N+1}^{\infty}\frac{1-6\,(j\,x_{0})^{2}+(j\,x_{0})^{4}}{(1+(j\,x_{0})^{2})^{4}}.$
(32)

We can then subtract Eq. (32) from Eq. (31) and, multiplying by
$N\,\kappa_{0}^{2}\,x_{0}^{-4}$, we get:

$\displaystyle
N\,\kappa_{0}^{2}\,x_{0}^{-4}\left[\text{Re}\;\psi^{(3)}\left(\frac{i+x_{0}}{x_{0}}\right)-\text{Re}\;\psi^{(3)}\left(\frac{i}{x_{0}}+N\right)\right]=2N\sum_{j=1}^{N-1}\,\frac{1-6\,(j\,x_{0})^{2}+(j\,x_{0})^{4}}{(1+(j\,x_{0})^{2})^{4}},$
(33)

which is the first of the desired sums in $S_{N}(x_{0})$. The second sum can
be computed in a similar fashion by combining the appropriate Polygamma
functions. We obtain, carrying out the algebra:

$\displaystyle S_{N}(x_{0})$ $\displaystyle=$
$\displaystyle\frac{N}{6x_{0}^{5}}\bigg{\\{}3x_{0}\,\text{Re}\;\psi^{(2)}\left(\frac{i+x_{0}}{x_{0}}\right)-3x_{0}\,\text{Re}\;\psi^{(2)}\left(\frac{i}{x_{0}}+N\right)-2\,\text{Im}\;\psi^{(3)}\left(\frac{i+x_{0}}{x_{0}}\right)$
(34) $\displaystyle+$ $\displaystyle
x_{0}\,N\,\text{Re}\;\psi^{(2)}\left(\frac{i+x_{0}}{x_{0}}\right)+\,\text{Im}\;\psi^{(3)}\left(\frac{i}{x_{0}}+N\right)+\,x_{0}N\,\text{Re}\;\psi^{(3)}\left(\frac{i}{x_{0}}+N\right)\bigg{\\}}.$

Eq. (34) can be further simplified by invoking two key properties of the
Polygamma functions:

$\displaystyle\psi^{(m)}(z+1)=\psi^{(m)}(z)-(-1)^{m}m!\;z^{-m-1},\qquad\psi^{(m)}(z)-(-1)^{m}\psi^{(m)}(1-z)=-\pi\frac{d^{m}}{dz^{m}}\cot(\pi
z).$

This allows us to write some of the terms in Eq. (34) in terms of elementary
functions, by using

$\displaystyle\psi^{(m)}\\!\left(\frac{i+x_{0}}{x_{0}}\right)-(-1)^{m}\psi^{(m)}\\!\left(\frac{-i+x_{0}}{x_{0}}\right)=(-1)^{m}m!\left(\frac{i}{x_{0}}\right)^{-m-1}-\pi\frac{d^{m}}{du^{m}}\cot(\pi
u),$

with $u=x_{0}^{-1}$. For example, recalling the above expression for
$\text{Re}\;\psi^{(3)}\left(\frac{i+x_{0}}{x_{0}}\right)$, we have

$\displaystyle\text{Re}\;\psi^{(3)}\\!\left(\frac{i+x_{0}}{x_{0}}\right)=-3+\left(\frac{\pi}{x_{0}}\right)^{4}\left[2+\cosh\left(2\frac{\pi}{x_{0}}\right)\right]\sinh\left(\frac{\pi}{x_{0}}\right)^{-4}\approx-3+8\left(\frac{\pi}{x_{0}}\right)^{4}e^{-2{\pi}/{x_{0}}},$

where the last approximate equality comes from assuming a small-lattice
parameter regime, $x_{0}\lesssim 0.5$, which is where the optimal qubit
separation lies, as shown by numerical optimization. Furthermore, all terms
carrying Polygamma functions with $({i}/{x_{0}}+N)$ in their argument give a
negligible contribution to the series $S_{N}(x_{0})$ in the region where the
minimum of $F_{N}(x_{0})$ occurs. We can then safely drop them in Eq. (34), if
our task is to find the optimal lattice separation. Lastly, we expand the
remaining terms in Eq. (34) in powers of $x_{0}$, for $x_{0}\lesssim 0.5$.
Substituting in Eq. (29) after performing these steps finally leads to the
following analytic approximation to the spatial function:

$\displaystyle
F_{N}(x_{0})\approx\frac{1}{6}\bigg{[}8N\left(\frac{\pi}{x_{0}}\right)^{4}\exp\left(-2\frac{\pi}{x_{0}}\right)+2\,\frac{1}{\pi^{2}}\,\left(\frac{\pi}{x_{0}}\right)^{2}\bigg{]}\equiv
F^{\text{an}}_{N}(x_{0}).$ (35)

Figure 6: Spatial function and optimal lattice parameter. Left: Spatial
function $F_{N}(x_{0})$ given by Eq. (29) (grey, solid) and analytic
approximation $F^{\text{an}}_{N}(x_{0})$ given by Eq. (35) (blue, dashed) vs.
lattice parameter $x_{0}$ for $N=100$ qubits. Right: Optimal lattice parameter
as a function of qubit number. Grey dots: Exact numerical optimization,
yielding the optimal lattice parameter $x^{\text{GHZ}}_{0\text{ opt}}$. Blue,
dashed: Approximate analytic solution $x^{\text{GHZ}}_{0\text{an}}$, given by
Eq. (38). Noise parameters as in Fig. 2, except that now dephasing is clearly
not collective.

In Fig. 6(left) we have plotted $F_{N}(x_{0})$ and $F^{\text{an}}_{N}(x_{0})$
as a function of the lattice parameter $x_{0}$ for $N=100$ qubits. One can see
that the analytical expression $F^{\text{an}}_{N}(x_{0})$ reproduces the
behavior of the spatial function to great accuracy, which only increases for
large $N$. Bearing that in mind, in what follows we drop the superscript “an”
and denote it by the same label as its exact counterpart, $F_{N}(x_{0})$. The
optimal lattice parameter minimizing $F_{N}(x_{0})$ lies around $x_{0\text{
opt}}^{\text{ GHZ}}\approx 0.44$; however, there is a neighborhood around this
point in which the spatial function stays approximately constant. Importantly,
since small departures from the optimal $x_{0}$ value do not drastically alter
$F_{N}(x_{0})$, no prohibitively precise noise characterization is needed.
This grants our noise-tailored sensor a degree of parameter robustness with
respect to deviations from the ideal lattice separation.

To obtain an approximate expression for $x_{0\text{ opt}}^{\text{GHZ}}$, it is
convenient to let $v\equiv{\pi}/{x_{0}}$ and write
$F_{N}(v)=\frac{1}{6}\left(8Nv^{4}e^{-2v}+2\frac{1}{\pi^{2}}\,v^{2}\right)$.
Taking the derivative with respect to $v$ and equating to zero then yields the
following condition:

$\displaystyle
8Ne^{-2v}v^{2}(2-v)+\frac{1}{\pi^{2}}=0\quad\Rightarrow\quad-8\pi^{2}Ne^{-2v}v^{3}\approx
1.$ (36)

The second approximate equality is fulfilled in the $N\gg 1$ limit, where the
optimal lattice parameter obeys $x_{0\text{ opt}}^{\text{ GHZ}}\ll 1$, and
hence $v^{\text{GHZ}}_{\text{ an}}\gg 2$. This can again be solved with the
help of the Lambert $W$ function. The relevant solution is given by

$\displaystyle v^{\text{GHZ}}_{\text{ an}}=\frac{\pi}{x_{0\text{
opt}}^{\text{GHZ}}}\approx-\frac{3}{2}\;W_{-1}\\!\left(-(3\;\pi^{2/3}\;N^{1/3})^{-1}\right).$
(37)

Since the argument in the Lambert function vanishes for $N\gg 1$, we can use
an expansion for $W_{-1}(x)$, valid when $x<0$ and $|x|\ll 1$ [56], to cast
Eq. (37) into a simpler form. That is, $W_{-1}(x)\approx L1-L2+\frac{L1}{L2}$,
with $L1\equiv\log(|x|)$, and $L2\equiv\log(-\log(|x|))$. In practice, keeping
the first two terms of the expansion, $W_{-1}(x)\approx L1-L2$, greatly
simplifies the ensuing algebra and still provides an accurate enough
characterization for our purposes. Substituting in Eq. (37) and solving for
$x_{0\text{ an}}^{\text{GHZ}}$, we finally obtain

$\displaystyle x_{0\text{ opt}}^{\text{GHZ}}\approx x_{0\text{
an}}^{\text{GHZ}}=\frac{2\pi}{\log(3\pi^{2/3}N^{1/3})+\log(\log(3\pi^{2/3}N^{1/3}))}.$
(38)

We have plotted the analytic approximation $x_{0\text{ an}}^{\text{GHZ}}$ and
exact optimal lattice parameter $x_{0\text{ opt}}^{\text{GHZ}}$, obtained by
numerical optimization, as function of qubit number $N$ in Fig. 6(right). We
can see that Eq. (38) has an $N$-dependent offset. This is the result of
having discarded the constant factor, $(v-2)\approx v$ in Eq. (36) when
minimizing the spatial function, which only vanishes in the asymptotic limit.
Despite this, the analytic value lies close enough to the true value of the
optimal lattice parameter to yield an excellent approximation to the optimal
uncertainty, see Fig. 4 in main text. Upon substituting Eq. (38) into Eq.
(35), we arrive at an expression for the minimum of the spatial function,
namely,

$\displaystyle F_{N}(x_{0\text{ an}}^{\text{
GHZ}})=\frac{1}{{6\pi^{2}}}\Big{[}\log\left(\pi^{2}N\right)+\log\left(\log^{3}\left(27\pi^{2}N\right)\right)\Big{]}^{2}\bigg{[}\frac{\left(\log\left(\log^{3}\left(27\pi^{2}N\right)\right)+\log\left(\pi^{2}N\right)\right)^{2}}{\log\left(27\pi^{2}N\right)}+\frac{1}{2}\bigg{]},$
(39)

which is of $\mathcal{O}(\log(N)^{2})$ in the $N\gg 1$ limit. Finally,
replacing $F_{N}(x^{\text{GHZ}}_{0\text{ an}})$ in the time-optimized
uncertainty, Eq. (28), we obtain an analytic formula for the overall minimun
uncertainty,

$\displaystyle\Delta\hat{b}^{\text{GHZ}}_{\text{R
opt}}(N)\approx\sqrt{2}\;\frac{(1+2\,W(e^{-1/2}/4))^{1/4}}{\sqrt{-W(e^{-1/2}/4)}}\;\sqrt{\frac{\omega_{c}}{T}}\;\left[F_{N}(x^{\text{GHZ}}_{0\text{
an}})\right]^{1/4}N^{-1}\propto\sqrt{\log(N)}\,N^{-1},$ (40)

with $\sqrt{2}\;(1+2\,W(e^{-1/2}/4))^{1/4}(\sqrt{-W(e^{-1/2}/4)})^{-1}\approx
2.96$ as noted above, and with the last proportionality relation holding for
$N\gg 1$. Remarkably, the agreement between the above analytic expression and
the exact numerical result is excellent even at finite number of qubits,
$N\gtrsim 20$, despite the discrepancies between $x^{\text{GHZ}}_{0\text{
opt}}$ and $x_{0\text{ an}}^{\text{GHZ}}$; see again Fig. 4 in the main text.

###  A.2 Initial OAT states

####  A.2.1 Expectation values of collective angular momentum operators

Generally, the dynamics of a state $|\text{OATS}\rangle=e^{-i\beta
J_{x}}e^{-i\theta J_{z}^{2}}|\text{CSS}\rangle_{x}$ subject to Gaussian
dephasing noise are highly entangled, and exact expressions for the relevant
mean values are not available. To evaluate the mean values, a cumulant
expansion over the qubit operators, which matches numerical calculations to
great accuracy, was devised in Ref. [32]. In Ref. [33] we have further
provided the relevant expressions $\langle J_{y}(\tau)\rangle,\langle
J_{y}^{2}(\tau)\rangle$ for non-collective spin-boson dephasing, where qubit
$n$ is assumed to be at position $\vec{r}_{n}$; analogous formulas for
$\langle J_{x}(\tau)\rangle,\langle J_{x}^{2}(\tau)\rangle$ may be readily
deduced. Here, we work in 1D, with $r_{n}\equiv nx_{0}$ in the limit of short
interrogation time $\omega_{c}\tau\ll 1$. In this regime, the quantum noise
(phase) contributions to the dynamics may be disregarded to first
approximation, as we argued in the main text. The relevant cumulant expansion
for the required expectation values can then be greatly simplified:

$\displaystyle\langle J_{x}(\tau)\rangle$ $\displaystyle\approx
e^{-\kappa_{0}^{2}(\omega_{c}\tau)^{2}}\mathcal{Q}_{0}(\theta,\beta)\cos(\varphi),$
(41) $\displaystyle\langle J_{y}(\tau)\rangle$ $\displaystyle\approx
e^{-\kappa_{0}^{2}(\omega_{c}\tau)^{2}}\mathcal{Q}_{0}(\theta,\beta)\sin(\varphi),$
(42) $\displaystyle\langle J_{v}^{2}(\tau)\rangle$
$\displaystyle\approx\frac{N}{4}+e^{-2\kappa_{0}^{2}(\omega_{c}\tau)^{2}}\left[C_{1}(\theta,\beta)\,G_{+}(x_{0})\pm\cos(2\varphi)\,C_{2}(\theta,\beta)\,G_{-}(x_{0})\right],\quad
v\in\\{x,y\\},$ (43)

with
$G_{\pm}(x_{0})\equiv\sum_{j=1}^{N-1}\,(N-j)\,e^{\pm\,2\kappa_{0}^{2}(\omega_{c}\tau)^{2}\tilde{g}_{1}(jx_{0})}$
capturing the lattice spatial correlations, and
$\mathcal{Q}_{0}(\theta,\beta),C_{1}(\theta,\beta)$, $C_{2}(\theta,\beta)$
being functions of the angles $\theta,\beta$. In the $N\gg 1$ limit we are
interested in, the latter take the form:

$\displaystyle Q_{0}(\theta,\beta)$
$\displaystyle=\frac{N}{2}e^{N\theta^{2}/8},$ (44) $\displaystyle
C_{1}(\theta,\beta)$
$\displaystyle=\frac{1}{8}\,\left[2\,\cos(\beta)^{2}-(1+e^{-N\theta^{2}/2})\sin(\beta)^{2}-2e^{-N\theta^{2}/2}\sin(2\beta)\sin\left(\theta/2\right)\right],$
(45) $\displaystyle C_{2}(\theta,\beta)$
$\displaystyle=\frac{1}{16}\left[2\sin(\beta)^{2}+4\sin(2\beta)\sin\left(\theta/2\right)e^{-N\theta^{2}/8}+(\cos(2\beta)+3)\,e^{-N\theta^{2}/2}\right].$
(46)

We emphasize that numerical calculations show that, in the time regime of
interest, for as low as $N=30$ qubits, the ratio uncertainty which follows
from these approximate expressions already closely reproduces its counterpart
including quantum noise contributions, whose expressions are much more
involved; see inset in Fig. 4(right) of the main text.

####  A.2.2 Optimal noise-tailored precision scaling

Let us now derive the optimal sensor performance for an OATS under classical
dephasing (that is, ignoring contributions from the quantum spectra, hence
$\varphi_{0}(t)$). As remarked in the main text, the squeezing and twisting
angles were chosen to minimize the initial variance along the $y$ axis,
according to Eq. (16). The numerator and denominator in
$\Delta\hat{b}^{\text{OAT}}_{R}(\tau)$, obtained by substituting the OATS mean
values from Eqs. (41)-(43) into Eq. (19) must first be minimized with respect
to their phase argument $\varphi$. Since an optimally squeezed initial
variance along $y$ has a maximally anti-squeezed component along $x$, the
ideal phase is such that this latter term vanishes:
$\varphi=k\pi,k\in\mathbb{N}$. In this case, Eq. (19) becomes proportional to
the precision derived from method of moments when measuring the observable
$\mathcal{O}=J_{y}$:

$\Delta\hat{b}_{\text{R}}^{\text{OAT}}(\tau)^{2}=\frac{\Delta
J_{y}^{2}(\tau)}{T\tau\,\langle J_{x}(\tau)\rangle^{2}},$ (47)

where we used $\partial_{\varphi}\langle J_{y}(\tau)\rangle=\langle
J_{x}(\tau)\rangle$. For the same fixed total time $2T$, however, when we are
only measuring $J_{y}$ and computing the precision through the method of
moments, we have $\Delta\hat{b}^{\text{OAT}}(\tau)^{2}=\Delta
J_{y}^{2}(\tau)/(2T\tau\,[\partial_{\varphi}\langle
J_{y}(\tau)\rangle]^{2}\,)$: that is, the variance of the ratio estimator is
larger by a factor 2.

The variance $\Delta J_{y}^{2}(\tau)$ in the numerator of Eq. (47) can now be
expanded in powers of $(\omega_{c}\tau)$ up to the fourth order, whereas we
only need to keep the denominators’ leading-order contribution; replacing in
$\Delta\hat{b}_{R}(\tau)$ yields

$\displaystyle\Delta\hat{b}_{\text{R}}(\tau)\approx\frac{\left[\,a_{0}(N,x_{0})+a_{2}(N,x_{0})(\omega_{c}\tau)^{2}+a_{4}(N,x_{0})\,(\omega_{c}\tau)^{4}\,\right]^{1/2}}{\sqrt{T\tau}\,h_{0}(N)},$
(48)

where the coefficients have a well-defined $N$-dependence through the ideal
squeezing and twisting angles of Eq. (16), as well as a spatial dependence
coming from the $G_{\pm}(x_{0})$ in the short-time expansion. Analytic
formulas for the coefficients can be derived by means of the Polygamma
functions approach worked out for the GHZ state in Sec. A.1. This yields
$h_{0}(N)\approx N/2$ and

$\displaystyle a_{0}(N,x_{0})$ $\displaystyle\approx$
$\displaystyle\frac{3^{2/3}}{8}N^{1/3},$ $\displaystyle a_{2}(N,x_{0})$
$\displaystyle\approx$
$\displaystyle\frac{1}{2}\Gamma(s+1)\bigg{\\{}\frac{1}{6x_{0}^{2}}+4\bigg{(}\frac{\pi}{x_{0}}\bigg{)}^{4}e^{-\frac{2\pi}{x_{0}}}\bigg{[}\bigg{(}\frac{1}{3}\bigg{)}^{1/3}N^{1/3}-\bigg{(}\frac{1}{3}\bigg{)}^{2/3}N^{2/3}+\frac{1}{3}N\bigg{]}+\frac{1}{2}\big{(}3^{1/3}N^{2/3}-3^{2/3}N^{1/3}\big{)}\bigg{\\}},$
$\displaystyle a_{4}(N,x_{0})$ $\displaystyle\approx$
$\displaystyle\bigg{\\{}\left[\Gamma(s+3)\left(\frac{1}{16\,3^{1/3}}N^{1/3}-\frac{1}{16\,3^{2/3}}N^{2/3}\right)+\frac{1}{2}\,N\,\Gamma(s+1)^{2}\right]-\frac{\Gamma(s+3)}{480}\frac{1}{x_{0}^{2}}$
$\displaystyle+\,\,\Gamma(s+3)e^{-\frac{2\pi}{x_{0}}}\left(\frac{\pi}{x_{0}}\right)^{6}\left(-N^{1/3}\frac{1}{30\,3^{1/3}}+\frac{1}{30\,3^{2/3}}N^{2/3}+N\frac{1}{90}\right)\bigg{\\}}.$

Similar to the GHZ state, Eq. (48) allows us to minimize with respect to
$\tau$ and $x_{0}$ independently. We can obtain the optimal interrogation time
by setting the first derivative of Eq. (48) with respect to $\tau$ equal to
zero:

$\displaystyle\tau^{\text{OATS}}_{\text{R
opt}}=\frac{1}{\sqrt{6}a_{4}(N,x_{0})}\left[\left(12a_{0}(N,x_{0})a_{4}(N,x_{0})+a_{2}(N,x_{0})^{2}\right)^{1/2}-a_{2}(N,x_{0})\right]^{1/2}.$
(49)

Evaluating Eq. (48) at $\tau^{\text{OATS}}_{\text{R opt}}$ leads, in turn, to
the following approximate time-optimized frequency uncertainty:

$\displaystyle\Delta\hat{b}_{\text{R}}(\tau^{\text{OATS}}_{\text{R opt}})$
$\displaystyle=$
$\displaystyle\frac{2^{5/4}}{3^{3/4}}\,\frac{\,[a_{2}(x_{0},N)\Delta+12a_{0}(N,x_{0})a_{4}(N,x_{0})]^{1/2}}{\,h_{0}(N)N\left[a_{4}(N,x_{0})\Delta\right]^{1/4}},$
(50) $\displaystyle\Delta$ $\displaystyle\equiv$
$\displaystyle\left(12a_{0}(N,x_{0})a_{4}(N,x_{0})+a_{2}(N,x_{0})^{2}\right)^{1/2}\\!\\!-a_{2}(N,x_{0}).$

Analytic minimization with respect to $x_{0}$, however, is involved. Numerical
evidence shows that the value of the ideal lattice parameter is remarkably
close to the value $x_{0\text{ an}}^{\text{OATS}}$ that minimizes
$a_{2}(x_{0},N)$. Following a similar procedure as in deriving Eq. (37), we
can then provide an asymptotic solution to the equation
$\partial_{x_{0}}a_{2}(x_{0},N)=0$:

$\displaystyle x_{0\text{
an}}^{\text{OATS}}=-\frac{2\pi}{3}\log\left[\frac{\log(3\pi^{2/3}N^{1/3})}{3\pi^{2/3}N^{1/3}}\right].$

Figure 7: Optimal lattice separation for OAT states under classical noise.
Grey dots: Exact numerical optimization, yielding $x_{\text{
opt}}^{\text{OATS}}$. Blue dashed line: Approximate analytic solution
$x_{\text{opt an}}^{\text{OATS}}$, Eq. (38). A supra-Ohmic spectral density
with an exponential cutoff, and $\alpha=1,s=3,\omega_{c}=1$ is assumed.

Much like for the GHZ state, the values of $x_{0\text{ an}}^{\text{OATS}}$
differ from the numerical optimization by a small offset, see Fig. 7;
nonetheless, parameter robustness leads to the analytic value obtained by
computing $\Delta\hat{b}_{\text{R}}(\tau^{\text{OATS}}_{\text{R
opt}},x_{0\text{ an}}^{\text{OATS}})$ to be extremely close to the exact
minimum. Evaluating the expansion coefficients at $x_{0\text{
an}}^{\text{OATS}}$, and disregarding all non-leading orders in $N$, we get:

$\displaystyle a_{0}(N,x_{\text{
an}}^{\text{OATS}})\approx\frac{3^{2/3}}{2^{7/3}}N^{1/3},\;\;\;\;a_{2}(N,x_{\text{
an}}^{\text{OATS}})\approx\frac{3^{1/3}}{2^{5/3}}\,\Gamma(s+1)\,N^{2/3},\;\;\;\;\;a_{4}(N,x_{\text{an}}^{\text{OATS}})\approx\frac{1}{4}\;\Gamma(s+1)^{2}\,N.$

For supra-Ohmic spectral densities with $s=3$, in particular, Zeno-like
scaling is then found upon substituting the above approximate expressions in
Eq. (50),

$\displaystyle\Delta\hat{b}_{\text{R
opt}}^{\text{OATS}}(N)\approx\frac{1}{\sqrt{3}}\,\Gamma(s+1)^{1/4}\sqrt{\frac{\omega_{c}}{T}}N^{-3/4},$
(51)

attained at an optimal time $\tau_{\text{R opt}}^{\text{OATS}}\approx
3^{-1/3}\,\Gamma(s+1)^{-1/2}N^{-1/6}$, in excellent agreement with numerical
optimization. We remark that, based on numerical fittings, Eq. (51) has been
found to hold for all $s>1$ spectral densities we investigated, although it
has only been rigorously derived for $s=3$.

## References

  * 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).
  * Wineland  _et al._ [1994] D. J. Wineland, J. J. Bollinger, W. M. Itano, and D. J. Heinzen, Squeezed atomic states and projection noise in spectroscopy, Phys. Rev. A 50, 67 (1994).
  * Pezzè  _et al._ [2018] L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, Quantum metrology with nonclassical states of atomic ensembles, Rev. Mod. Phys. 90, 035005 (2018).
  * Diddams  _et al._ [2004] S. A. Diddams, J. C. Bergquist, S. R. Jefferts, and C. W. Oates, Standards of time and frequency at the outset of the 21st century, Science 306, 1318 (2004).
  * 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).
  * Colombo  _et al._ [2022] S. Colombo, E. Pedrozo-Penafiel, and V. Vuletić, Entanglement-enhanced optical atomic clocks, Appl. Phys. Lett. 121, 210502 (2022).
  * Lee  _et al._ [2002] H. Lee, P. Kok, and J. P. Dowling, A quantum Rosetta stone for interferometry, J. Mod. Opt. 49, 2325 (2002).
  * Degen  _et al._ [2017] C. L. Degen, F. Reinhard, and P. Cappellaro, Quantum sensing, Rev. Mod. Phys. 89, 035002 (2017).
  * 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).
  * Giovannetti  _et al._ [2004] V. Giovannetti, S. Lloyd, and L. Maccone, Quantum-enhanced measurements: Beating the standard quantum limit, Science 306, 1330 (2004).
  * [11] Super-Heisenberg scaling is achievable in non-linear, interaction-based quantum metrology using unentangled states, see S. Boixo, A. Datta, M. J. Davis, S. T. Flammia, A. Shaji, and C. M. Caves, “Quantum Metrology: Dynamics versus Entanglement,” Phys. Rev. Lett. 101, 040403 (2008).
  * Haase  _et al._ [2018] J. F. Haase, A. Smirne, J. Kolodynski, R. Demkowicz-Dobrzanski, and S. F. Huelga, Fundamental limits to frequency estimation: a comprehensive microscopic perspective, New J. Phys. 20, 053009 (2018).
  * Escher  _et al._ [2011] B. M. Escher, R. L. de Matos Filho, and L. Davidovich, General framework for estimating the ultimate precision limit in noisy quantum-enhanced metrology, Nat. Phys. 7, 406 (2011).
  * Matsuzaki  _et al._ [2011] Y. Matsuzaki, S. C. Benjamin, and J. Fitzsimons, Magnetic field sensing beyond the standard quantum limit under the effect of decoherence, Phys. Rev. A 84, 012103 (2011).
  * 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).
  * Macieszczak [2015] K. Macieszczak, Zeno limit in frequency estimation with non-Markovian environments, Phys. Rev. A 92, 010102 (2015).
  * Bylander  _et al._ [2011] J. Bylander, S. Gustavsson, F. Yan, F. Yoshihara, K. Harrabi, G. Fitch, D. G. Cory, Y. Nakamura, J.-S. Tsai, and W. D. Oliver, Noise spectroscopy through dynamical decoupling with a superconducting flux qubit, Nat. Phys. 7, 565 (2011).
  * Álvarez and Suter [2011] G. A. Álvarez and D. Suter, Measuring the spectrum of colored noise by dynamical decoupling, Phys. Rev. Lett. 107, 230501 (2011).
  * Romach  _et al._ [2015] Y. Romach, C. Müller, T. Unden, L. J. Rogers, T. Isoda, K. M. Itoh, M. Markham, A. Stacey, J. Meijer, S. Pezzagna, B. Naydenov, L. P. McGuinness, N. Bar-Gill, and F. Jelezko, Spectroscopy of surface-induced noise using shallow spins in diamond, Phys. Rev. Lett. 114, 017601 (2015).
  * Malinowski  _et al._ [2017] F. K. Malinowski, F. Martins, L. Cywiński, M. S. Rudner, P. D. Nissen, S. Fallahi, G. C. Gardner, M. J. Manfra, C. M. Marcus, and F. Kuemmeth, Spectrum of the nuclear environment for gaas spin qubits, Phys. Rev. Lett. 118, 177702 (2017).
  * Wang  _et al._ [2017] Y. Wang, M. Um, J. Zhang, S. An, M. Lyu, J.-N. Zhang, L.-M. Duan, D. Yum, and K. Kim, Single-qubit quantum memory exceeding ten-minute coherence time, Nat. Photon. 11, 646 (2017).
  * Chan  _et al._ [2018] K. W. Chan, W. Huang, C. H. Yang, J. C. C. Hwang, B. Hensen, T. Tanttu, F. E. Hudson, K. M. Itoh, A. Laucht, A. Morello, and A. S. Dzurak, Assessment of a Silicon quantum dot spin qubit environment via noise spectroscopy, Phys. Rev. Appl. 10, 044017 (2018).
  * [23] G. Wang, Y. Zhu, B. Li, C. Li, L. Viola, A. Cooper, and P. Cappellaro, Digital noise spectroscopy with a quantum sensor, arXiv:2212.09216 https://doi.org/10.48550/arXiv.2212.09216.
  * Quintana  _et al._ [2017] C. M. Quintana, Y. Chen, D. Sank, A. G. Petukhov, T. C. White, D. Kafri, B. Chiaro, A. Megrant, R. Barends, B. Campbell, Z. Chen, A. Dunsworth, A. G. Fowler, R. Graff, E. Jeffrey, J. Kelly, E. Lucero, J. Y. Mutus, M. Neeley, C. Neill, P. J. J. O’Malley, P. Roushan, A. Shabani, V. N. Smelyanskiy, A. Vainsencher, J. Wenner, H. Neven, and J. M. Martinis, Observation of classical-quantum crossover of $1/f$ flux noise and its paramagnetic temperature dependence, Phys. Rev. Lett. 118, 057702 (2017).
  * Yan  _et al._ [2018] F. Yan, D. Campbell, P. Krantz, M. Kjaergaard, D. Kim, J. L. Yoder, D. Hover, A. Sears, A. J. Kerman, T. P. Orlando, S. Gustavsson, and W. D. Oliver, Distinguishing coherent and thermal photon noise in a circuit quantum electrodynamical system, Phys. Rev. Lett. 120, 260504 (2018).
  * 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).
  * 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).
  * Dorner [2012] U. Dorner, Quantum frequency estimation with trapped ions and atoms, New J. Phys. 14, 043011 (2012).
  * Boter  _et al._ [2020] J. M. Boter, X. Xue, T. Krähenmann, T. F. Watson, V. N. Premakumar, D. R. Ward, D. E. Savage, M. G. Lagally, M. Friesen, S. N. Coppersmith, M. A. Eriksson, R. Joynt, and L. M. K. Vandersypen, Spatial noise correlations in a Si/Sige two-qubit device from Bell state coherences, Phys. Rev. B 101, 235133 (2020).
  * [30] J. Yoneda, J. S. Rojas-Arias, P. Stano, K. Takeda, A. Noiri, T. Nakajima, D. Loss, and S. Tarucha, Noise-correlation spectrum for a pair of spin qubits in silicon, arXiv:2208.14150 https://doi.org/10.48550/arXiv.2208.14150.
  * [31] Generalized Ramsey protocols which leverage the existence of a decoherence free subspace and attain HL are possible for estimation of frequency differences or sums [28].
  * Beaudoin  _et al._ [2018] F. Beaudoin, L. M. Norris, and L. Viola, Ramsey interferometry in correlated quantum noise environments, Phys. Rev. A 98, (R)020102 (2018).
  * Riberi  _et al._ [2022] F. Riberi, L. M. Norris, F. Beaudoin, and L. Viola, Frequency estimation under non-Markovian spatially correlated quantum noise, New J. Phys. 24, 103011 (2022).
  * Layden and Cappellaro [2011] D. Layden and P. Cappellaro, Spatial noise filtering through error correction for quantum sensing, npj Quantum Inf. 4, 30 (2011).
  * Liu and Yuan [2016] J. Liu and H. Yuan, Valid lower bound for all estimators in quantum parameter estimation, New J. Phys. 18, 093009 (2016).
  * Chabuda  _et al._ [2016] K. Chabuda, I. D. Leroux, and R. Demkowicz-Dobrzański, The quantum Allan variance, New J. Phys. 18, 083035 (2016).
  * Rubio and Dunningham [2019] J. Rubio and J. Dunningham, Quantum metrology in the presence of limited data, New J. Phys. 21, 043347 (2019).
  * Rojkov  _et al._ [2022] I. Rojkov, D. Layden, P. Cappellaro, J. Home, and F. Reiter, Bias in error-corrected quantum sensing, Phys. Rev. Lett. 128, 140503 (2022).
  * Yamamoto  _et al._ [2022] K. Yamamoto, S. Endo, H. Hakoshima, Y. Matsuzaki, and Y. Tokunaga, Error-mitigated quantum metrology via virtual purification, Phys. Rev. Lett. 129, 250503 (2022).
  * Paz-Silva  _et al._ [2017] G. A. Paz-Silva, L. M. Norris, and L. Viola, Multiqubit spectroscopy of Gaussian quantum noise, Phys. Rev. A 95, 022121 (2017).
  * Fagnola  _et al._ [2019] F. Fagnola, J. E. Gough, H. I. Nurdin, and L. Viola, Mathematical models of Markovian dephasing, J. Phys. A 52, 385301 (2019).
  * Ma  _et al._ [2011] J. Ma, X. Wang, C. P. Sun, and F. Nori, Quantum spin squeezing, Phys. Rep. 509, 89 (2011).
  * Kitagawa and Ueda [1993] M. Kitagawa and M. Ueda, Squeezed spin states, Phys. Rev. A 47, 5138 (1993).
  * Leroux  _et al._ [2010] I. D. Leroux, M. H. Schleier-Smith, and V. Vuletić, Implementation of cavity squeezing of a collective atomic spin, Phys. Rev. Lett. 104, 073602 (2010).
  * Brask  _et al._ [2015] J. B. Brask, R. Chaves, and J. Kołodyński, Improved quantum magnetometry beyond the standard quantum limit, Phys. Rev. X 5, 031010 (2015).
  * Helstrom [1976] C. W. Helstrom,  _Quantum Detection and Estimation Theory_ (Academic Press, New York, NY, 1976).
  * Kay [1993] S. M. Kay,  _Fundamentals of Statistical Signal Processing: Estimation Theory_ (Prentice Hall, Englewood Cliffs, NJ, 1993).
  * Pezzè and Smerzi [2014] L. Pezzè and A. Smerzi, Quantum theory of phase estimation, in  _Atom Interferometry, Proc. Int. School Phys. “Enrico Fermi”, Course 188, Varenna_ , edited by G. M. Tino and M. A. Kasevich (IOP Press, Amsterdam, 2014) p. 691.
  * Ischi  _et al._ [2005] B. Ischi, M. Hilke, and M. Dubé, Decoherence in a ${N}$-qubit solid-state quantum register, Phys. Rev. B 71, 195325 (2005).
  * Brown  _et al._ [2016] K. R. Brown, J. Kim, and C. Monroe, Co-designing a scalable quantum computer with trapped atomic ions, npj Quantum Inf. 2, 16034 (2016).
  * Postler  _et al._ [2018] L. Postler, Á. Rivas, P. Schindler, A. Erhard, R. Stricker, D. Nigg, T. Monz, R. Blatt, and M. Müller, Experimental quantification of spatial correlations in quantum dynamics, Quantum 2, 90 (2018).
  * Zhou  _et al._ [2020] H. Zhou, J. Choi, S. Choi, R. Landig, A. M. Douglas, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, P. Cappellaro, _et al._ , Quantum metrology with strongly interacting spin systems, Phys. Rev. X 10, 031003 (2020).
  * Schulte  _et al._ [2020] M. Schulte, V. J. Martinez-Lahuerta, M. S. Scharnagl, and K. Hammerer, Ramsey interferometry with generalized one-axis twisting echoes, Quantum 4, 268 (2020).
  * [54] M. S. Scharnagl, T. Kielinski, and K. Hammerer, Optimal Ramsey interferometry with echo protocols based on one-axis twisting, arXiv:2307.08510 https://doi.org/10.48550/arXiv.2307.08510.
  * [55] J. J. Meyer, S. Khatri, D. S. Franca, J. Eisert, and P. Faist, Quantum metrology in the finite-sample regime, arXiv:2307.06370 https://doi.org/10.48550/arXiv.2307.06370.
  * Corless  _et al._ [1996] R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. J. Jeffrey, and D. E. Knuth, On the Lambert W function, Adv. Comput. Math. 5, 329 (1996).