#  $B\to K^{*}M_{X}$ vs $B\to KM_{X}$ as a probe of a scalar-mediator dark
matter scenario

Alexander Berezhnoya and Dmitri Melikhova,b,c aD. V. Skobeltsyn Institute of
Nuclear Physics, M. V. Lomonosov Moscow State University, 119991, Moscow,
Russia
bJoint Institute for Nuclear Research, 141980 Dubna, Russia
cFaculty of Physics, University of Vienna, Boltzmanngasse 5, A-1090 Vienna,
Austria

###### Abstract

Recently, Belle II reported the observation of the decay $B\to KM_{X}$,
$M_{X}$ the missing mass, with the branching ratio much exceeding ${\cal
B}(B\to K\nu\bar{\nu})$ which is the only Standard Model (SM) process
contributing to this reaction. If confirmed, this might be an indication of
new nonSM particles produced in this decay. One of the possible explanations
of the observed effect could be light dark-matter (DM) particles produced via
a scalar mediator field. We give simple arguments, that a combined analysis of
the $B\to KM_{X}$ and $B\to K^{*}M_{X}$ reactions would be a clean probe of
the scalar mediator scenario: (i) making use of an observed value ${\cal
B}(B\to KM_{X})\simeq 5.4\,{\cal B}(B\to K\nu\bar{\nu})_{\rm SM}$ and (ii)
assuming that the effect is due to the light dark matter coupling to the top
quark via a scalar mediator field, one finds an upper limit ${\cal B}(B\to
K^{*}M_{X})<2.8\,{\cal B}(B\to K^{*}\nu\bar{\nu})_{\rm SM}$. The measurement
of the ${\cal B}(B\to K^{*}M_{X})$ will allow one to determine the scalar
mediator mass within this scenario.

1\. Introduction. A recent Belle II observation belle  of $B\to KM_{X}$ (for
which also the notation $B\to K\not{E}$, $\not{E}$ the missing energy, is
used) at the level much exceeding the Standard model (SM) prediction for $B\to
K\nu\bar{\nu}$,

$\displaystyle{\cal B}(B^{+}\to K^{+}M_{X})=(2.4\pm 0.7)\times
10^{-5}\simeq(5.4\pm 1.5)\,{\cal B}(B^{+}\to K^{+}\nu\bar{\nu})_{\rm SM},$ (1)

opened a window for immediate discussions of possible new physics effects
capable to explain this result (see e.g. recent publications 2023hiller ;
2023damir ; 2023a ; 2023dm ; 2023b ; 2023c ). One of the discussed options is
the decay into Dark Matter (DM) particles 2023dm  with multiple scenarios for
the content of these particles and the possible mediators. The goal of this
short communication is to remark that combining the present Belle II result
for the $B\to KM_{X}$ with the hypothesis of a DM origin of the enhancement of
$\Gamma(B\to KM_{X}$ coupled to the SM particles via a scalar mediator field
leads to interesting and rather robust constraints independent of further
details of the DM model.

2\. The $B\to(K,K^{*})\phi\to(K,K^{*})\bar{\chi}\chi$ decays. As an example we
consider a simple but rather general Lagrangian describing the iteraction of
dark-matter particles $\chi$ with the top quark via a scalar mediator field
$\phi$ Leffa ; Leffb :

$\displaystyle\mathcal{L_{\rm
int}}=-\frac{ym_{f}}{v}\phi\bar{f}f-\frac{1}{2}\kappa\phi\bar{\chi}\chi\ ,$
(2)

The corresponding effective Lagrangian describing the flavour-changing neutral
current (FCNC) $b\to s\phi$ vertex has the form Leffa ; Leffb

$\displaystyle\mathcal{L}_{b\rightarrow
s\phi}=\frac{ym_{b}}{v}\frac{3\sqrt{2}G_{F}m_{q}^{2}V^{*}_{qs}V_{qb}}{16\pi^{2}}\left(\phi\bar{s}_{L}b_{R}+{\rm
h.c.}\right).$ (3)

Essential for our discussion is merely the operator structure of this vertex,
$\phi\,\bar{s}_{L}b_{R}$. Using the QCD equations of motion, it is
straigtforward to calculate the amplitudes

$\displaystyle\langle K|\bar{s}_{L}b_{R}|B\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2}\langle
K|\bar{s}(1-\gamma_{5})b|B\rangle=\frac{1}{2}\langle
K|\bar{s}b|B\rangle=\frac{1}{2}\frac{M_{B}^{2}-M_{K}^{2}}{m_{b}-m_{s}}f_{0}^{B\to
K}(q^{2}),$ $\displaystyle\langle K^{*}|\bar{s}_{L}b_{R}|B\rangle$
$\displaystyle=$ $\displaystyle\frac{1}{2}\langle
K^{*}|\bar{s}(1-\gamma_{5})b|B\rangle=-\frac{1}{2}\langle
K^{*}|\bar{s}\gamma_{5}b|B\rangle=-i(\epsilon
q)\frac{M_{K^{*}}}{m_{b}+m_{s}}A_{0}^{B\to K^{*}}(q^{2}),$ (4)

where the dimensionless form factors $f_{0}$ and $A_{0}$ are well-known
quantities parametrizing the $\langle K|\bar{s}\gamma_{\mu}b|B\rangle$ and
$\langle K^{*}|\bar{s}\gamma_{\mu}\gamma_{5}b|B\rangle$ amplitudes wsb  and
the corresponding decay rates have the form

$\displaystyle\Gamma_{B\rightarrow K\phi}$ $\displaystyle=$
$\displaystyle\left(\frac{ym_{b}}{v}\frac{3\sqrt{2}G_{F}m_{t}^{2}V^{*}_{ts}V_{tb}}{16\pi^{2}}\right)^{2}\frac{\lambda^{1/2}(M_{B}^{2},M_{K}^{2},M_{\phi}^{2})}{16\pi
M_{B}^{3}}|\langle K|\bar{s}_{L}b_{R}|B\rangle|^{2},$ (5)
$\displaystyle\Gamma_{B\rightarrow K^{*}\phi}$ $\displaystyle=$
$\displaystyle\left(\frac{ym_{b}}{v}\frac{3\sqrt{2}G_{F}m_{t}^{2}V^{*}_{ts}V_{tb}}{16\pi^{2}}\right)^{2}\frac{\lambda^{1/2}(M_{B}^{2},M_{K}^{2},M_{\phi}^{2})}{16\pi
M_{B}^{3}}\sum_{K^{*}\mbox{polar}}|\langle
K^{*}|\bar{s}_{L}b_{R}|B\rangle|^{2}.$ (6)

Performing summation over the $K^{*}$ polarizations, we find for the ratio of
the differential distributions $\Gamma_{B\rightarrow K^{*}M_{X}}$ and
$\Gamma_{B\rightarrow K^{*}M_{X}}$ proceeding via the virtual scalar $\phi$
the following expression:

$\displaystyle\frac{\Gamma_{B\rightarrow K^{*}M_{X}}}{\Gamma_{B\rightarrow
KM_{X}}}=\frac{\lambda^{3/2}(M_{B}^{2},M_{K^{*}}^{2},M_{X}^{2})}{\lambda^{1/2}(M_{B}^{2},M_{K}^{2},M_{X}^{2})(M_{B}^{2}-M_{K}^{2})^{2}}\left[\frac{m_{b}-m_{s}}{m_{b}+m_{s}}\right]^{2}\left|\frac{A_{0}(M_{X}^{2})}{f_{0}(M_{X}^{2})}\right|^{2}$
(7)

We point out that the $\phi$-propagator (for which the finite-width effects
should be taken into account), and all factors related to the transition of
the virtual $\phi$ into the $\chi\chi$-pair cancel in the ratio.

3\. Numerical estimates. For the form factors, we make use of the results from
ballP  and ballV  in the form of convenient parametrizations of ms2000 :

$\displaystyle f_{0}(q^{2})$ $\displaystyle=$
$\displaystyle\frac{0.33}{1-0.7r_{V}+0.27r_{V}^{2}},\quad
r_{V}=q^{2}/m_{B_{s}^{*}}^{,}$ (8) $\displaystyle A_{0}(q^{2})$
$\displaystyle=$ $\displaystyle\frac{0.37}{(1-0.46r_{P})(1-r_{P})},\quad
r_{P}=q^{2}/m_{B_{s}}^{2}.$ (9)

We then come to the following numerical prediction for the ratio of the
differential distribitions in the reactions $B\to K\phi\to K\chi\chi$ and
$B\to K^{*}\phi\to K^{*}\chi\chi$. Notice that using more recent
parametrizations for the form factors, e.g., those from damir2023epjc ;
hpqcd2023 , does not lead to any visible changes of the results in Fig. 1.

Fig. 1: The ratio $d\Gamma(B\to K^{*}M_{X})/d\Gamma(B\to KM_{X})$ vs $M_{X}$,
the missing mass of the outgoing $\chi\bar{\chi}$ pair of DM particles
interacting with the top quark via scalar mediator field. Phase-space
restriction related to the finite masses of the DM $\chi$ particles should be
taken into account (i.e. not all values of $M_{X}$ shown in the plot are
accessible in the decay). We point out that all factors related to the
$\phi$-propagator (including its possible finite-width effects) and to the
subsequent decay $\phi\to\chi\chi$ precisely cancel in the ratio independently
of the masses of $\chi$ and the details of the $\phi\chi\chi$ vertex.

4\. Discussion. As clear from the plot, one can obtain an upper limit on the
ratio of the partial widths due to the contribution of the diagrams with a
scalar mediator $\phi$:

$\displaystyle\delta\Gamma(B\to K^{*}M_{X})_{\phi}<1.05\,\delta\Gamma(B\to
KM_{X})_{\phi}.$ (10)

This upper limit is obtained from Fig. 1 at $M_{X}\to 0$. Noteworthy, this
upper limit is independent of the coupling of the mediator to the DM
particles. The upper limit of Eq. (10) is a rather conservative estimate and
in practice one may observe a sizeable suppression of the $B\to K^{*}M_{X}$
compared to $B\to KM_{X}$: The $M_{X}^{2}$-region providing the dominant
contribution to the $B\to(K,K^{*})\bar{\chi}\chi$ amplitudes is determined by
the resonant structure of the $\phi$-propagator and corresponds to nonzero
$M_{X}^{2}$ leading to a suppression $\delta\Gamma(B\to
K^{*}M_{X})_{\phi}<\delta\Gamma(B\to KM_{X})_{\phi}$. The accurate
experimental determination of the ratio $\delta\Gamma(B\to
K^{*}M_{X})_{\phi}/\delta\Gamma(B\to KM_{X})_{\phi}$ opens the way to
measuring the scalar mediator mass.

Now, assuming that the observed enhancement of the SM $B\to K\nu\bar{\nu}$
cross section is due to the DM particles in the final state interacting via a
scalar mediator we can write

$\displaystyle\Gamma(B\to KM_{X})$ $\displaystyle=$ $\displaystyle\Gamma(B\to
K\nu\bar{\nu})_{\rm SM}+\delta\Gamma(B\to KM_{X})_{\phi},$ (11)
$\displaystyle\Gamma(B\to K^{*}M_{X})$ $\displaystyle=$
$\displaystyle\Gamma(B\to K^{*}\nu\bar{\nu})_{\rm SM}+\delta\Gamma(B\to
K^{*}M_{X})_{\phi}.$ (12)

For ${\cal B}(B\to K^{*}\nu\bar{\nu})/{\cal B}(B\to K\nu\bar{\nu})$ one
obtains $2.5\pm 0.7$ making use of old estimates mns1997a ; mns1997b ; mns1998
or $2.23\pm 0.40$ using the recent analysis of damir2023epjc . Within the
quoted uncertainties these two values are in excellent agreement with each
other.

Collecting together the central values of the discussed constraints:

$\displaystyle\Gamma(B\to K^{*}\nu\bar{\nu})_{\rm SM}$ $\displaystyle\simeq$
$\displaystyle 2.5\,\Gamma(B\to K\nu\bar{\nu})_{\rm SM}\hskip
19.91684pt\mbox{[Theoretical estimates]}$ (13) $\displaystyle\delta\Gamma(B\to
KM_{X})_{\phi}$ $\displaystyle\simeq$ $\displaystyle 4.4\,\Gamma(B\to
K\nu\bar{\nu})_{\rm SM}\hskip 19.91684pt\mbox{[Belle measurement]}$ (14)
$\displaystyle\delta\Gamma(B\to K^{*}M_{X})_{\phi}$ $\displaystyle<$
$\displaystyle 1.05\,\delta\Gamma(B\to KM_{X})_{\phi}\hskip
11.38092pt\mbox{[Our result here]},$ (15)

we then obtain the following upper bound for the $\Gamma(B\to K^{*}M_{X})$
within the scalar mediator scenario:

$\displaystyle\Gamma(B\to K^{*}M_{X})<2.8\,\Gamma(B\to K^{*}\nu\bar{\nu})_{\rm
SM}.$ (16)

Therefore, for the case of scalar mediator, an expected enhancement of the
branching ratio $\Gamma(B\to K^{*}M_{X})$ compared to the SM value is
approximately two times smaller than the enhancement for $\Gamma(B\to
KM_{X})$, Eq. (1). This constraint is well compatible with the present
experimental limit belle2017

$\displaystyle{\cal B}(B\to K^{*}M_{X})\leq 2.7\times 10^{-5}\leq 2.7\,{\cal
B}(B\to K^{*}\nu\bar{\nu})_{\rm SM},$ (17)

where at the second step we have used a recently updated estimate ${\cal
B}(B\to K^{*}\nu\bar{\nu})_{\rm SM}=(1.03\pm 0.15)\times 10^{-5}$
damir2023epjc .

On the other hand, if the “mediator” is a spin-1 field, it is natural to
expect a similar ehnancement for both $K$ and $K^{*}$ mesons in the final
state, ${\cal B}(B\to K^{*}M_{X})/{\cal B}(B\to K^{*}\nu\bar{\nu})_{\rm
SM}\simeq{\cal B}(B\to KM_{X})/{\cal B}(B\to K\nu\bar{\nu})_{\rm SM}$.

We therefore conclude that the combined analysis of ${\cal B}(B\to
K^{*}M_{X})$ and ${\cal B}(B\to K\,M_{X})$ may provide a clear probe of the
scalar mediator scenario and if this scenario is confirmed eventually allow
one to extract the mass of the scalar mediator making use of the results
presented in Fig. 1.

Acknowledgments. The research was carried out within the framework of the
program “Particle Physics and Cosmology” of the National Center for Physics
and Mathematics.

## References

  * (1) A. Glazov (Belle II), News from Belle II, BELLE2-TALK-CONF-2023-123.
  * (2) R. Bause, H. Gisbert, and G. Hiller, Implications of an enhanced $B\to K\nu\bar{\nu}$ branching ratio, arXiv:2309.00075.
  * (3) L. Allwicher, D. Becirevic, G. Piazza, S. Rosauro-Alcaraz, O. Sumensari, Understanding the first measurement of $Br(B\to K\nu\bar{\nu})$, arXiv:2309.02246.
  * (4) T. Felkl, A. Giri, R. Mohanta, M. A. Schmidt, When Energy Goes Missing: New Physics in $b\to s\nu\bar{\nu}$ with Sterile Neutrinos, arXiv:2309.02940.
  * (5) M. Abdughani and Y. Reyimuaji, Constraining light dark matter and mediator with $B^{+}\to K^{+}\nu\bar{\nu}$ data, arXiv:2309.03706.
  * (6) H. K. Dreiner, J. Y. Guenther, and Z. .S. Wang, The Decay $B\to K\nu\bar{\nu}$ at Belle II and a Massless Bino in R-parity-violating Supersymmetry, arXiv:2309.03727.
  * (7) X.-G. He, X.-D. Ma, and G. Valencia, Revisiting models that enhance $B\to K\nu\bar{\nu}$ in light of the new Belle II measurement, arXiv:2309.12741.
  * (8) B. Batell, M. Pospelov, and A. Ritz, Multi-lepton Signatures of a Hidden Sector in Rare B Decays, Phys. Rev. D83, 054005 (2011).
  * (9) K. Schmidt-Hoberg, F. Staub, and M. W. Winkler, Constraints on light mediators: confronting dark matter searches with B physics, Phys. Lett. B727, 506 (2013).
  * (10) M. Wirbel, B. Stech, and M. Bauer, Exclusive Semileptonic Decays of Heavy Mesons, Z. Phys. C29, 637 (1985).
  * (11) P. Ball and R. Zwicky, New results on $B\to\pi,K,\eta$ decay formfactors from light-cone sum rules, Phys. Rev. D71, 014015 (2005).
  * (12) P. Ball and R. Zwicky, $B_{d,s}\to\rho,\omega,K^{*},\phi$ decay form-factors from light-cone sum rules revisited, Phys. Rev. D71, 014029 (2005).
  * (13) D. Melikhov and B.Stech, Weak form-factors for heavy meson decays: an update, Phys. Rev. D62, 014006 (2000).
  * (14) D. Becirevic, G. Piazza, and O. Sumensari, Revisiting $B\to K^{(*)}\nu\bar{\nu}$ decays in the Standard Model and beyond, Eur. Phys. J. C83, 252 (2023).
  * (15) W. G. Parrott, C. Bouchard, and C. T. H. Davies (HPQCD collaboration), Standard Model predictions for $B\to Kl^{+}l^{-}$, $B\to Kl_{1}^{-}l_{2}^{+}$, and $B\to K\nu\bar{\nu})$ using form factors from $N_{f}=2+1+1$ lattice QCD, Phys. Rev. D107, 014511 (2023); Phys. Rev. D107, 119903(E) (2023).
  * (16) J. Grygier et al. [Belle], Search for $\bm{B\to h\nu\bar{\nu}}$ decays with semileptonic tagging at Belle, Phys. Rev. D96, 091101 (2017).
  * (17) P. Colangelo, F. De Fazio, P. Santorelli, and E. Scrimieri, Rare $B\to(K.K^{*})\bar{\nu}\bar{\nu}$ decays at $B$-factories, Phys. Lett. B395, 339 (1997).
  * (18) D. Melikhov, N. Nikitin, and S. Simula, Rare decays $B\to(K,K^{*})(l^{+}l^{-},\nu\bar{\nu})$ in the quark model, Phys. Lett. B410, 290 (1997).
  * (19) D. Melikhov, N. Nikitin, and S. Simula, Right-handed currents in rare exclusive $B\to(K,K^{*})\nu\bar{\nu}$ decays, Phys. Lett. B428, 171 (1997).
  * (20) D. Melikhov, N. Nikitin, and S. Simula, Rare exclusive semileptonic $b\to s$ transitions in the standard model, Phys. Rev. D57, 6814 (1998).


# Supplementary Online Material: Quantum-enhanced sensing of displacements and
electric fields with large trapped-ion crystals

K. A. Gilmore  National Institute of Standards and Technology, Boulder,
Colorado 80305, USA  Department of Physics, University of Colorado, Boulder,
Colorado, 80309, USA  M. Affolter  National Institute of Standards and
Technology, Boulder, Colorado 80305, USA  R. J. Lewis-Swan  Homer L. Dodge
Department of Physics and Astronomy, The University of Oklahoma, Norman,
Oklahoma 73019, USA  Center for Quantum Research and Technology, The
University of Oklahoma, Norman, Oklahoma 73019, USA  D. Barberena  JILA, NIST
and Department of Physics, University of Colorado, 440 UCB, Boulder, Colorado
80309, USA  Center for Theory of Quantum Matter, University of Colorado,
Boulder, CO 80309, USA  E. Jordan  National Institute of Standards and
Technology, Boulder, Colorado 80305, USA  A. M. Rey  JILA, NIST and Department
of Physics, University of Colorado, 440 UCB, Boulder, Colorado 80309, USA
Center for Theory of Quantum Matter, University of Colorado, Boulder, CO
80309, USA  J. J. Bollinger  National Institute of Standards and Technology,
Boulder, Colorado 80305, USA

##  S1 Ideal model of displacement sensor

The ideal displacement sensing protocol is realized by the implementation of
the ODF Hamiltonian,

$\hat{H}_{\text{ODF}}=\frac{g}{\sqrt{N}}(\hat{a}+\hat{a}^{\dagger})\hat{J}_{z}.$
(S1)

Here, $\hat{a}$ ($\hat{a}^{\dagger}$) is an annihilation (creation) operator
for the center of mass (COM) bosonic mode that couples uniformly to all spins
with strength $g$ and
$\hat{J}_{x,y,z}=\sum_{i=1}^{N}\hat{\sigma}_{x,y,z}^{j}/2$ are collective spin
operators, with $\hat{\sigma}_{x,y,z}^{j}$ Pauli operators acting on ion $j$.
As described in the main text, the protocol involves three steps.

  1. 1.

The application of $\hat{H}_{\text{ODF}}$ for a time $\tau$ on an initially
unentangled product state $\ket{\psi(0)}$ of the spins and phonons, which
creates an entangled state with quantum-enhanced sensitivity:

$\ket{\psi_{\text{in}}}=e^{-i\tau\hat{H}_{\text{ODF}}}\ket{\psi(0)}.$ (S2)

This entangled state can be understood by considering the ideal initial
condition of the main text,
$\ket{\psi(0)}=\ket{\left(N/2\right)_{x}}\otimes\ket{0}_{b}$ where
$\ket{0}_{b}$ is the bosonic vacuum and coherent spin state, defined by
$\hat{J}_{\alpha}\ket{m_{\alpha}}=m_{\alpha}\ket{m_{\alpha}}$. Then, from Eq.
(S2) it is straightforward to obtain

$\ket{\psi_{\text{in}}}=\sum_{m_{z}=-N/2}^{N/2}c_{m_{z}}\ket{\alpha_{m_{z}}}_{b}\otimes\ket{m_{z}},$
(S3)

Here, $\ket{\alpha_{m_{z}}}_{b}$ is a coherent state with amplitude
$\alpha_{m_{z}}=-im_{z}g\tau/\sqrt{N}$, while $c_{m_{z}}$ are expansion
coefficients defined by
$\ket{\left(N/2\right)_{x}}=\sum_{m_{z}}c_{m_{z}}\ket{m_{z}}$. This is a type
of generalized spin-phonon cat state, known to be useful for metrology [1, 2,
3].

  2. 2.

The instantaneous application of a displacement $\beta\in\mathbb{R}$. This is
equivalent to the application of the unitary transformation
$\hat{U}_{\beta}=e^{\beta(\hat{a}-\hat{a}^{\dagger})}$ on
$\ket{\psi_{\text{in}}}$, to give the intermediate state

$\ket{\psi^{\beta}_{\text{in}}}=e^{\beta(\hat{a}^{\dagger}-\hat{a})}\ket{\psi_{\text{in}}}.$
(S4)

  3. 3.

Application of $-\hat{H}_{\text{ODF}}$ for a time $\tau$, a process that is
equivalent to time-reversal of the initial entangling operation,

$\displaystyle\begin{split}\ket{\psi_{f}}&=e^{i\tau\hat{H}_{\text{ODF}}}\ket{\psi_{\beta}}=\sum_{m_{z}=-N/2}^{N/2}c_{m_{z}}e^{2i\beta\alpha_{m_{z}}}\ket{\beta}_{b}\otimes\ket{m_{z}}.\end{split}$
(S5)

The final state (S5) can be written in an equivalent form that better
elucidates how the time-reversal maps the displacement to an overall rotation
of the spin degree of freedom. This result is actually general, independent of
the initial spin-boson state. To show this, consider the final state at the
end of the protocol

$\ket{\psi_{f}}=e^{i\tau\hat{H}_{\text{ODF}}}e^{\beta(\hat{a}-\hat{a}^{\dagger})}e^{-i\tau\hat{H}_{\text{ODF}}}\ket{\psi(0)}.$
(S6)

Using that the ODF interaction is equivalent to an $J_{z}$-dependent
displacement of the COM mode,
$\exp(-i\tau\hat{H}_{\text{ODF}})=\mathcal{D}(-ig\tau\hat{J}_{z})$ where
$\mathcal{D}(\alpha)=e^{\alpha\hat{a}^{\dagger}-\alpha^{*}\hat{a}}$ is the
usual bosonic displacement operator, combined with the identity
$\mathcal{D}(\beta)^{\dagger}\hat{a}\mathcal{D}(\beta)=\hat{a}+\beta$ we can
rewrite the last equation as

$\ket{\psi_{f}}=e^{\frac{2ig\tau\beta}{\sqrt{N}}\hat{J}_{z}}e^{\beta(\hat{a}^{\dagger}-\hat{a})}\ket{\psi(0)}.$
(S7)

In this form it is clear that the spins and bosons will be completely
disentangled at the conclusion of the protocol so long as the initial state
$\ket{\psi(0)}$ is also unentangled. Moreover, the initial state of the bosons
is irrelevant and the displacement manifests independently as: i) a
displacement of the boson spin state and ii) a rotation of the initial spin
state about $\hat{z}$ by an angle $\varphi=2g\tau\beta/\sqrt{N}$.

By initializing the spins along the $\hat{x}$-axis, e.g., in the state
$\ket{\left(N/2\right)_{x}}$, we are thus able to observe this rotation
through a measurement of simple spin observables. For example, a small
displacement $\beta$ will lead to a small rotation that generates a nonzero
value of
$\braket{\hat{J}_{y}}_{f}\equiv\bra{\psi_{f}}\hat{J}_{y}\ket{\psi_{f}}$.

The usefulness of this observable is determined by computing the expected
signal $\braket{\hat{J}_{y}}_{f}$ at the end of the protocol and comparing it
to its associated quantum variance $\text{var}(\hat{J}_{y})_{f}$ computed also
at the end of the protocol (note that both quantities are functions of
$\beta$). This is quantified by the signal to noise ratio

$\text{SNR}=\frac{\braket{\hat{J}_{y}}_{f}}{\sqrt{\text{var}(\hat{J}_{y})_{f}}}\approx\frac{\frac{\partial\braket{\hat{J}_{y}}_{f}}{\partial\beta}}{\sqrt{\text{var}(\hat{J}_{y})_{f}}}\Bigg{|}_{\beta=0}\beta,$
(S8)

where $\beta=0$ is the working point of our protocol, having the property that
$\braket{\hat{J}_{y}}_{f}|_{\beta=0}=0$. The inverse of the slope (with
respect to $\beta$) of the SNR is the sensitivity $\Delta\beta$ with which
$\beta$ can be estimated using $\hat{J}_{y}$ as a measurement observable,
given by

$(\Delta\beta)^{2}\equiv\frac{\text{Var}(\hat{J}_{y})_{f}}{\Big{(}\frac{\partial\braket{\hat{J}_{y}}_{f}}{\partial\beta}\Big{)}^{2}}\bigg{|}_{\beta=0}=\frac{1}{4g^{2}\tau^{2}}.$
(S9)

This expression quantifies the ratio of: i) how fast the signal varies as a
function of $\beta$ ($\partial_{\beta}\langle\hat{J}_{y}\rangle_{f}$) and ii)
the fundamental noise present in any measurement of $\hat{J}_{y}$ due to
quantum fluctuations [$\text{Var}(\hat{J}_{y})_{f}$]. For interaction times
$g\tau\gg 1$ the change in the signal,
$\partial_{\beta}\langle\hat{J}_{y}\rangle_{f}=-g\tau\sqrt{N}$, is much larger
than the intrinsic quantum noise
$\sqrt{\text{Var}(\hat{J}_{y})_{f}}\sim\sqrt{N}$ and therefore sensing
displacements with $\beta<1/4$ (SQL) is possible.

##  S2 Role of entanglement

In this section we analyze the precise role that entanglement plays in the
time reversal protocol. To keep the analysis simple, we work in the limit of
large $N$, which we can do because of the many body nature of the system.
Then, since the initial spin state is polarized along $+x$, we can approximate
the spin degree of freedom by a canonically conjugate pair
$(\hat{x}_{s},\hat{p}_{s})$ according to

$\hat{J}_{x}\approx N/2\hskip
56.9055pt\hat{J}_{z}\approx\sqrt{\frac{N}{2}}\,\hat{x}_{s}\hskip
56.9055pt\hat{J}_{y}\approx-\sqrt{\frac{N}{2}}\,\hat{p}_{s},$ (S10)

with the initial state $\ket{\left(N/2\right)_{x}}$ corresponging to the
bosonic vacuum $\ket{0}_{s}$. It is also convenient to analyze the COM using
quadratures instead of creation and annihilation operators:

$\hat{x}_{b}=\frac{\hat{a}+\hat{a}^{\dagger}}{\sqrt{2}}\hskip
56.9055pt\hat{p}_{b}=\frac{\hat{a}-\hat{a}^{\dagger}}{i\sqrt{2}}.$ (S11)

###  S2.1 Properties of the state $\ket{\psi_{\text{in}}}$

Let’s analyze the state after the first step [Eq. (S2)], which is
approximately given by

$\ket{\psi_{\text{in}}}=e^{-ig\tau\hat{x}_{s}\hat{x}_{b}}\ket{0}_{b}\otimes\ket{0}_{s}.$
(S12)

While this evolution certainly mixes spins and bosons, it has a very simple
description in terms of the hybrid quadratures

$\hat{x}_{+}=\frac{\hat{x}_{s}+\hat{x}_{b}}{\sqrt{2}}\hskip
56.9055pt\hat{x}_{-}=\frac{\hat{x}_{s}-\hat{x}_{b}}{\sqrt{2}},$ (S13)

whereby

$\ket{\psi_{\text{in}}}=e^{-i\frac{g\tau}{2}(\hat{x}_{+}^{2}-\hat{x}_{-}^{2})}\ket{0}_{+}\otimes\ket{0}_{-},$
(S14)

and $\ket{0}_{+}\otimes\ket{0}_{-}=\ket{0}_{s}\otimes\ket{0}_{b}$ is the
vacuum in either pair of variables, with wavefunction $\psi(0)\propto
e^{-(x_{+}^{2}+x_{-}^{2})/2}=e^{-(x_{s}^{2}+x_{b}^{2})/2}$. Consider the
$x_{+}$ quadrature. Since everything is gaussian, we can compute the Wigner
function [4] at time $\tau$:

$W(x_{+},p_{+})=\frac{1}{\pi}\exp\big{[}-x_{+}^{2}-(p_{+}-g\tau
x_{+})^{2}\big{]}=\frac{1}{\pi}\exp\bigg{[}-\frac{\bar{x}^{2}}{\xi^{2}}-\xi^{2}\bar{p}^{2}\bigg{]},$
(S15)

where

$\displaystyle\begin{split}\bar{x}=x_{+}\cos\phi+p_{+}\sin\phi\hskip
56.9055pt\bar{p}=p_{+}\cos\phi-x_{+}\sin\phi,\end{split}$ (S16)

are rotated quadratures about an angle $\phi>0$:

$\tan\phi=\frac{\sqrt{4+(g\tau)^{2}}-g\tau}{2},$ (S17)

and $\xi^{2}$ is the squeezing parameter, given by

$\xi^{2}=\frac{1}{1+g^{2}\tau^{2}+g\tau\sqrt{2+g^{2}\tau^{2}}}\approx\frac{1}{2g^{2}\tau^{2}}\hskip
28.45274ptg\tau\gg 1.$ (S18)

This state is represented in Fig. S1 (left) and shows squeezing along the
$\bar{x}$ direction. We also show the Wigner function of the $(x_{-},p_{-})$
degree of freedom (right) which is just a mirror reflection of the previous
one.

Figure S1: We show the Wigner function in the hybrid phase spaces
$(\hat{x}_{+},\hat{p}_{+})$ (left) and $(\hat{x}_{-},\hat{p}_{-})$ (right) for
the state $\ket{\psi_{\text{in}}}$, after the application of the first ODF
pulse. In these variables the state is clearly squeezed, but the relation to
the original spin and boson variables is not straightforward.

A displacement of $x_{b}$, which is what we are trying to measure, manifests
as mirror shifts of the Wigner functions in the horizontal direction. Taking
advantage of this squeezing would require us to measure the $\bar{x}$
quadrature, which in terms of the original variables is very complicated
indeed:

$\hat{\bar{x}}\approx\frac{\hat{x}_{b}\cos\phi+\hat{p}_{b}\sin\phi}{\sqrt{2}}+\frac{\hat{J}_{z}\cos\phi-\hat{J}_{y}\sin\phi}{\sqrt{N}}.$
(S19)

Not only is this a time dependent observable ($\phi$ is a function of $\tau$),
but also we would need to measure both bosons and spins, which is not possible
with the current technical constraints of this experiment, which can only
access the spins. Measuring only the spins immediately after the displacement
of the bosons (implemented by $e^{\beta(\hat{a}^{\dagger}-\hat{a})}$) is
useless because
$e^{-\beta(\hat{a}^{\dagger}-\hat{a})}\hat{O}e^{\beta(\hat{a}^{\dagger}-\hat{a})}=\hat{O}$
for any spin observable $\hat{O}$. Measuring the bosons only is also useless,
as we can see by computing their reduced density matrix $\hat{\rho}_{B}$:

$\braket{x_{b}}{\hat{\rho}_{B}}{x_{b}^{\prime}}=\rho_{B}(x_{b},x_{b}^{\prime})=\frac{1}{\sqrt{\pi}}e^{-(x_{b}^{2}+x_{b}^{\prime
2})/2}e^{-(x_{b}-x_{b}^{\prime})^{2}g^{2}\tau^{2}/4},$ (S20)

whose Wigner function

$W_{S}(x_{b},p_{b})=\frac{1}{\pi\sqrt{1+(g\tau)^{2}}}e^{-x_{b}^{2}}e^{-p_{b}^{2}/(1+g^{2}\tau^{2})},$
(S21)

shows only antisqueezing in the $p_{b}$ quadrature.

###  S2.2 Entanglement dynamics

Given the crucial role that entanglement plays in the protocol, let us compute
an entanglement measure as a function of the time $t$ in a sequence of
duration $2\tau$. The state is given by

$\ket{\psi_{t}}=\begin{cases}e^{-it\hat{H}_{\text{ODF}}}\ket{0}_{b}\otimes\ket{\left(N/2\right)_{x}}&\text{if
}0<t<\tau\\\\[5.0pt]
e^{i(t-\tau)\hat{H}_{\text{ODF}}}\,e^{\beta(\hat{a}^{\dagger}-\hat{a})}\,e^{-i\tau\hat{H}_{\text{ODF}}}\ket{0}_{b}\otimes\ket{\left(N/2\right)_{x}}&\text{if
}\tau<t<2\tau\end{cases},$ (S22)

and we choose the second Renyi entropy

$S_{2}=-\log[\mathrm{Tr}(\hat{\rho}_{B}^{2})]$ (S23)

as our entanglement proxy. The state during the first half ($t<\tau$) is of
the form of $\ket{\psi_{\text{in}}}$ with the replacement $\tau\to t$, so the
Renyi entropy can be directly calculated from the reduced density matrix [Eq.
(S20)] as

$S_{2}(t)=-\log\bigg{[}\int\,dx_{b}\,dx_{b}^{\prime}\,\rho_{B}(x_{b},x_{b}^{\prime})^{2}\bigg{]}=\frac{1}{2}\log(1+g^{2}t^{2}).$
(S24)

For the second half ($\tau<t<2\tau$)

$\ket{\psi_{t}}=e^{i(t-\tau)\hat{H}_{\text{ODF}}}\,e^{\beta(\hat{a}^{\dagger}-\hat{a})}\,e^{-i\tau\hat{H}_{\text{ODF}}}\ket{0}_{b}\otimes\ket{\left(N/2\right)_{x}}=e^{\beta(\hat{a}^{\dagger}-\hat{a})}\,e^{\frac{2ig\beta\hat{J}_{x}}{\sqrt{N}}}\,e^{-i(2\tau-t)\hat{H}_{\text{ODF}}}\ket{0}_{b}\otimes\ket{\left(N/2\right)_{x}},$
(S25)

where we have used again the fact $e^{-i\tau\hat{H}_{\text{ODF}}}$ is an
$\hat{J}_{z}$ dependent displacement. Since the first two exponentials on the
left act on independent degrees of freedom, they do not modify entanglement
measures. Thus $S_{2}$ depends only on the rightmost exponential, whose effect
we have already computed. Thus, we have that

$S_{2}(t)=\begin{cases}\frac{1}{2}\log\big{[}1+g^{2}t^{2}\big{]}^{-1/2}&\text{if
}0<t<\tau\\\ \frac{1}{2}\log\big{[}1+g^{2}(t-2\tau)^{2}\big{]}^{-1/2}&\text{if
}\tau<t<2\tau\end{cases}.$ (S26)

At then end ($t=2\tau$) there is no spin-boson entanglement, in agreement with
Eq. (S7). The initial build up of entanglement ($t<\tau$) is crucial to keep
the fluctuations $\text{Var}(\hat{J}_{y})_{f}$ small at the end of the
protocol. We explore the consequences of omitting this first step in the next
subsection.

###  S2.3 Read-out protocol

The protocol we will study in this subsection is characterized by omitting the
first step of the full time reversal protocol described in the main text.
Mathematically

$\ket{\psi_{f}^{\prime}}=\underbrace{e^{i\tau\hat{H}_{\text{ODF}}}}_{\text{ODF
interaction}}\hskip
5.69046pt\underbrace{e^{\beta(\hat{a}^{\dagger}-\hat{a})}}_{\text{Displacement}}\hskip
5.69046pt\underbrace{\ket{0}_{b}\otimes\ket{\left(N/2\right)_{x}}}_{\text{Initial
state}}.$ (S27)

This has been implemented before for off-resonant electric field sensing [5,
6], and the logic justifying this particular sequence of steps is
semiclassical. The displacement operation transforms the bosonic vacuum into a
coherent state: $\ket{0}_{b}\to\ket{\beta}_{b}$. Naively doing the replacement
$\hat{a},\hat{a}^{\dagger}\to\beta$ in $\hat{H}_{\text{ODF}}$ then makes
$\hat{H}_{\text{ODF}}\to 2g\beta\hat{J}_{z}/\sqrt{N}$ which generates a
rotation of the spin degree of freedom. Measuring $\hat{J}_{y}$ should then
give information about $\beta$. While this argument captures correctly the
result of the measurement signal, it omits the fact that the entanglement
created by $\hat{H}_{\text{ODF}}$ also increases the variance of
$\hat{J}_{y}$, which we now calculate. Working again in the large $N$ limit,
the evolution of the operator $\hat{J}_{y}$ in the Heisenberg picture is given
by

$\displaystyle\begin{split}e^{-\beta(\hat{a}^{\dagger}-\hat{a})}e^{-i\tau\hat{H}_{\text{ODF}}}\hat{J}_{y}e^{i\tau\hat{H}_{\text{ODF}}}e^{\beta(\hat{a}^{\dagger}-\hat{a})}&\longrightarrow-\sqrt{\frac{N}{2}}\Big{(}e^{i\beta\sqrt{2}\hat{p}_{b}}\Big{)}\Big{(}e^{-ig\tau\hat{x}_{s}\hat{x}_{b}}\Big{)}\,\hat{p}_{s}\Big{(}e^{ig\tau\hat{x}_{s}\hat{x}_{b}}\Big{)}\Big{(}e^{-i\beta\sqrt{2}\hat{p}_{b}}\Big{)}\\\
&=-\sqrt{\frac{N}{2}}e^{i\beta\sqrt{2}\hat{p}_{b}}(\hat{p}_{s}+g\tau\hat{x}_{b})e^{-i\beta\sqrt{2}\hat{p}_{b}}\\\
&=-\sqrt{\frac{N}{2}}(\hat{p}_{s}+g\tau\hat{x}_{b}+g\tau\beta\sqrt{2}).\end{split}$
(S28)

Since the initial state in this approximation is
$\ket{0}_{b}\otimes\ket{0}_{s}$, it follows that

$\displaystyle\begin{split}\braket{\hat{J}_{y}}_{f}&\approx-\sqrt{N}\beta
g\tau\\\
\text{var}(\hat{J}_{y})_{f}&\approx\frac{N}{4}(1+g^{2}\tau^{2}).\end{split}$
(S29)

We see that while the slope of the signal increases with $g\tau$, the variance
of $\hat{J}_{y}$ also increases with $g\tau$ in a way that completely
nullifies the effect of the increased signal. In fact, the sensitivity is

$(\Delta\beta)^{2}=\frac{1}{4g^{2}\tau^{2}}+\frac{1}{4}>\frac{1}{4},$ (S30)

and it never gets below the SQL. As can be seen from the last line of Eq.
(S28), this increase in variance is a result of the vacuum noise from the
bosons being fed into the spins, which in turn is a consequence of the build
up of entanglement.

##  S3 Generalized treatment of sensing with technical noise

To incorporate the effects of various sources of technical noise and
adequately describe both the displacement and electric field sensing protocols
of the main text, we adopt a model based on the more general time-dependent
Hamiltonian:

$\hat{H}_{\eta}(t)=-\delta\hat{a}^{\dagger}\hat{a}+\frac{g(t)}{\sqrt{N}}(\hat{a}^{\dagger}+\hat{a})\hat{J}_{z}+i\eta(t)(\hat{a}^{\dagger}-\hat{a}).$
(S31)

Here, we have: i) included a time-dependent ODF interaction strength [$g(t)$]
which can include, e.g., the sign change necessary for time-reversal, and ii)
a time-dependent coherent drive of strength $\eta(t)$ to account for the
concurrent application of the displacement and ODF drive in the electric field
sensing protocol. Moreover, we assume that $\eta(t)$ is sufficiently small so
that it can be accounted for within perturbation theory. Within Eq. (S31),
different sensing protocols simply correspond to different time profiles for
$\eta(t)$ and $g(t)$. Thus, we will first work with a generic time dependence
before specializing the obtained results to particular scenarios.

The inclusion of the static detuning $\delta=\mu-\omega_{z}$, where $\mu$ is
the ODF drive frequency and $\omega_{z}$ is the COM mode frequency, accounts
for the observed deviations of this quantity away from the ideal $\delta=0$ in
the experiment. Specifically, we assume that $\delta$ fluctuates randomly
between experimental trials but is fixed over the duration of a single
experimental trial. We model $\delta$ as a Gaussian random variable with mean
zero and variance $\sigma^{2}$.

Experimentally, we also find the ions are not perfectly initialized in their
motional ground-state. Limited by Doppler cooling, the COM mode is initially
well described as a thermal bosonic state with an occupation $\bar{n}$.
Combined with the spin-state after the initial microwave pulse, the initial
state of the spin-phonon system is thus the density matrix
$\hat{\rho}(0)=\rho_{\text{th}}(\bar{n})\otimes|\left(N/2\right)_{x}\rangle\langle\left(N/2\right)_{x}|$.
The final state after application of $\hat{H}_{\eta}$ for a time $t$ is then

$\hat{\rho}_{\eta}(t)=\hat{U}_{\eta}(t)\hat{\rho}(0)U_{\eta}^{\dagger}(t),$
(S32)

where $\hat{U}_{\eta}(t)$ is the time evolution operator associated to
$\hat{H}_{\eta}(t)$. The achievable sensitivity $\Delta\eta$ to the applied
drive is, for the specific observable $\hat{J}_{y}$ and at time $t$,
constructed from two quantities:

  1. 1.

$\text{Var}\big{[}\hat{J}_{y}(t)\big{]}$, the variance of the observable at
the end of the evolution when $\eta(t)=0$ (equivalent to $\beta=0$ in the
previous section).

  2. 2.

$\delta_{\eta}\braket{\hat{J}_{y}(t)}=\braket{\hat{J}_{y}(t)}_{\eta}-\braket{\hat{J}_{y}(t)}_{0}$,
the change in the observable at the end of the protocol with respect to the
unperturbed evolution to first order in $\eta(t)$ (equivalent to
$\beta\partial_{\beta}\braket{\hat{J}_{y}}$ in the previous section).

However, for calculational purposes it is more convenient to work with
$\braket{\hat{J}_{+}(t)}$, $\braket{\hat{J}_{+}^{2}(t)}$ and
$\delta_{\eta}\braket{\hat{J}_{+}(t)}$.

###  S3.1 Derivation

The general results for $\text{Var}(\hat{J}_{y}(t))$ and
$\delta\braket{\hat{J}_{y}(t)}$ are better expressed in terms of the
Heisenberg picture operators

$\hat{J}_{+}(t)\equiv\hat{U}_{0}^{\dagger}(t)\hat{J}_{+}\hat{U}_{0}(t)\hskip
56.9055pt\hat{a}(t)\equiv\hat{U}_{0}^{\dagger}(t)\hat{a}\,\hat{U}_{0}(t).$
(S33)

Note that this Heisenberg picture is defined in the absence of the
perturbation ($\eta=0$). The relevant quantities can then be calculated using
time-dependent perturbation theory with $\eta$ as a perturbation parameter:

$\displaystyle\begin{split}\braket{\hat{J}_{+}^{n}(t)}&=\mathrm{Tr}\Big{\\{}\hat{J}_{+}^{n}(t)\hat{\rho}(0)\Big{\\}}\\\
\delta_{\eta}\braket{\hat{J}_{+}(t)}&=-i\,\mathrm{Tr}\bigg{\\{}\hat{\rho}(0)\int_{0}^{t}i\eta(s)\,ds\,\Big{[}\hat{J}_{+}(t),\hat{a}(s)^{\dagger}-\hat{a}(s)\Big{]}\bigg{\\}},\end{split}$
(S34)

with $n=1,2$. The evolution of the operators $\hat{J}_{+}(t)$ and $\hat{a}(t)$
is given by the following Heisenberg equations of motion:

$\displaystyle i\dot{\hat{a}}(t)$
$\displaystyle=i\delta\hat{a}(t)-\frac{ig(t)}{\sqrt{N}}\hat{J}_{z}$ (S35)
$\displaystyle i\dot{\hat{J}}_{+}(t)$
$\displaystyle=-\frac{g(t)}{\sqrt{N}}\big{[}\hat{a}(t)^{\dagger}+\hat{a}(t)\big{]}\hat{J}_{+}(t).$
(S36)

As $\hat{J}_{z}$ is a conserved quantity the equation for $\hat{a}(t)$ is
closed. After formal solution of $\hat{a}(t)$ we subsequently solve for
$\hat{J}_{+}(t)$. The final result is

$\displaystyle\begin{split}\hat{a}(t)&=\hat{a}e^{i\delta
t}-\frac{ih(t)e^{i\delta t}}{\sqrt{N}}\hat{J}_{z}\\\
\hat{J}_{+}(t)&=\exp\bigg{[}\frac{i}{\sqrt{N}}(h(t)\hat{a}^{\dagger}+\bar{h}(t)\hat{a})\bigg{]}\exp\bigg{[}-\frac{2ip(t)}{\sqrt{N}}\bigg{(}\hat{J}_{z}+\frac{1}{2}\bigg{)}\bigg{]}\hat{J}_{+},\end{split}$
(S37)

where we have introduced the functions

$\displaystyle\begin{split}h(t)&=\int_{0}^{t}g(s)e^{-i\delta s}\,ds,\hskip
56.9055ptp(t)=\frac{i}{2}\int_{0}^{t}(\dot{\bar{h}}(s)h(s)-\bar{h}(s)\dot{h}(s))\,ds.\end{split}$
(S38)

Expectation values are taken with respect to the initial state, yielding,

$\displaystyle\begin{split}\braket{\hat{J}_{+}(t)}&=\frac{N}{2}\exp\bigg{[}-\frac{|h(t)|^{2}}{N}\bigg{(}\bar{n}+\frac{1}{2}\bigg{)}\bigg{]}\cos\bigg{(}\frac{p(t)}{N}\bigg{)}^{N-1},\\\\[10.0pt]
\braket{\hat{J}_{+}^{2}(t)}&=\frac{N(N-1)}{4}\exp\bigg{[}-\frac{4|h(t)|^{2}}{N}\bigg{(}\bar{n}+\frac{1}{2}\bigg{)}\bigg{]}\cos\bigg{(}\frac{2p(t)}{N}\bigg{)}^{N-2}.\end{split}$
(S39)

The last thing we need is $\delta\braket{\hat{J}_{+}(t)}$, which can be
calculated using Eq. (S34), the solutions in Eq. (S37) and systematic
application of commutation relations:

$\delta_{\eta}\braket{\hat{J}_{+}(t)}=i\sqrt{N}q(t)\exp\bigg{[}-\frac{|h(t)|^{2}}{N}\bigg{(}\bar{n}+\frac{1}{2}\bigg{)}\bigg{]}\cos\bigg{(}\frac{p(t)}{N}\bigg{)}^{N-1},$
(S40)

where

$q(t)=\int_{0}^{t}\int_{0}^{s}g(s)\eta(u)\cos[\delta(u-s)]\,du\,ds.$ (S41)

The previous formulas include all the Hamiltonian effects we will take into
account. However, the system is also subject to light scattering whenever the
ODF interaction is on, and is modelled by the following dissipative
Lindbladian:

$\mathcal{L}\hat{\rho}=\frac{\Gamma}{4}\sum_{i=1}^{N}\Big{(}\hat{\sigma}_{z}^{i}\hat{\rho}\hat{\sigma}_{z}^{i}-\hat{\rho}\Big{)}.$
(S42)

Since the Hamiltonian evolution commutes with $\mathcal{L}$ (in the
superoperator sense), the consequences of this dissipative process can be
calculated at the end and independently from the effects of $\hat{H}$ [7].
Thus, we need to replace
$\braket{\hat{J}_{+}^{n}(t)}\to\braket{\hat{J}_{+}^{n}(t)}e^{-n\Gamma
t_{\text{odf}}/2}$, where $t_{\text{odf}}$ is the amount of time the ODF is
on. In addition, the observable
$\braket{\hat{J}_{+}\hat{J}_{-}+\hat{J}_{-}\hat{J}_{+}}/2=N(N+1)/4$, which was
conserved throughout the protocol in the absence of $\Gamma$, becomes

$\frac{1}{2}\braket{\hat{J}_{+}\hat{J}_{-}+\hat{J}_{-}\hat{J}_{+}}\to\frac{N}{2}+\frac{N(N-1)}{4}e^{-\Gamma
t_{\text{odf}}}$ (S43)

###  S3.2 Summary of results

Using the definition of $\hat{J}_{y}$ in terms of $\hat{J}_{+}$ and
$\hat{J}_{-}$ we thus find that

$\displaystyle\begin{split}\braket{\hat{J}_{y}(t)}&=0\\\
\braket{\hat{J}_{y}^{2}(t)}&=\frac{N}{4}+\frac{N(N-1)e^{-\Gamma
t_{\text{odf}}}}{8}\Bigg{\\{}1-\exp\bigg{[}-\frac{4|h(t)|^{2}}{N}\bigg{(}\bar{n}+\frac{1}{2}\bigg{)}\bigg{]}\cos\bigg{(}\frac{2p(t)}{N}\bigg{)}^{N-2}\Bigg{\\}}\\\
\delta_{\eta}\braket{\hat{J}_{y}(t)}&=q(t)\sqrt{N}e^{-\Gamma
t_{\text{odf}}/2}\exp\bigg{[}-\frac{|h(t)|^{2}}{N}\bigg{(}\bar{n}+\frac{1}{2}\bigg{)}\bigg{]}\cos\bigg{(}\frac{p(t)}{N}\bigg{)}^{N-1}\\\
\braket{\hat{J}_{x}(t)}&=\frac{Ne^{-\Gamma
t_{\text{odf}}/2}}{2}\exp\bigg{[}-\frac{|h(t)|^{2}}{N}\bigg{(}\bar{n}+\frac{1}{2}\bigg{)}\bigg{]}\cos\bigg{(}\frac{p(t)}{N}\bigg{)}^{N-1}\end{split}$
(S44)

We see that $h(t)$ is associated to thermal and zero-point bosonic motion,
$p(t)$ is a squeezing contribution that increases the variance and $q(t)$ is
some kind of detuning averaged total displacement. We also show
$\braket{\hat{J}_{x}(t)}$ because it is necessary for the experimental
calibration of $\sigma$, the COM frequency variance.

###  S3.3 Displacement sensing

The scheme is shown in Fig. S2, with total protocol duration $2\tau$. The time
profiles are $g(t)=g\,\text{sign}(t-\tau)$, where the sign function accounts
for the time reversal step in the latter half of the protocol
($\tau<t<2\tau$), and $\eta(t)=\beta\delta(t-\tau)$ indicates that the drive
instantaneously generates a displacement $\beta$ just before the time reversal
step ($t=\tau$). Plugging this into the definitions of $h(t)$, $p(t)$ and
$q(t)$ we obtain

$\displaystyle\begin{split}|h{(2\tau)}|^{2}=\frac{16g^{2}}{\delta^{2}}\sin\bigg{(}\frac{\delta\tau}{2}\bigg{)}^{4}\hskip
42.67912ptp{(2\tau)}=\frac{g^{2}}{\delta^{2}}\big{[}4\sin(\delta\tau)-\sin(2\delta\tau)-2\delta\tau\big{]}\hskip
42.67912ptq{(2\tau)}=-\frac{\beta g}{\delta}\sin(\delta\tau).\end{split}$
(S45)

We then numerically average $\braket{\hat{J}_{y}(2\tau)}$ and
$\delta\braket{\hat{J}_{y}(2\tau)}$ over $\delta$ and calculate the
sensitivity according to

$(\Delta\beta)^{2}=\frac{\braket{\hat{J}_{y}^{2}(2\tau)}^{\text{av}}}{\big{(}\frac{\partial\braket{\hat{J}_{y}(2\tau)}^{\text{av}}}{\partial\beta}\big{)}^{2}},$
(S46)

where $\braket{...}^{\text{av}}$ indicates the average over detunings. For
short protocol times, $2\sigma\tau\ll 1$, we can calculate perturbatively the
effect of COM fluctuations. Keeping terms up to second order in $\delta$ and
lowest order in $1/N$, we obtain after averaging

$\displaystyle\begin{split}\braket{\hat{J}_{y}^{2}(2\tau)}^{\text{av}}&\approx\frac{N}{4}+Ng^{2}\sigma^{2}\tau^{4}e^{-2\Gamma\tau}\bigg{(}\frac{\bar{n}+1/2}{2}+\frac{g^{2}\tau^{2}}{9}\bigg{)}\\\
\frac{\delta_{\eta}\braket{\hat{J}_{y}(2\tau)}^{\text{av}}}{\beta}&\approx-\sqrt{N}g\tau
e^{-\Gamma\tau}\bigg{(}1-\frac{\sigma^{2}\tau^{2}}{6}\bigg{)}.\end{split}$
(S47)

Hence the sensitivity is

$(\Delta\beta)^{2}\approx\underbrace{\frac{e^{2\Gamma\tau}}{4g^{2}\tau^{2}}}_{\text{Ideal
result}}+\underbrace{\frac{\sigma^{2}e^{2\Gamma\tau}}{12g^{2}}}_{\text{Signal
reduction}}+\underbrace{\frac{\sigma^{2}\tau^{2}(\bar{n}+1/2)}{2}}_{\text{Spin-
phonon
correlations}}+\underbrace{\frac{g^{2}\sigma^{2}\tau^{4}}{9}}_{\text{Spin-spin
correlations}}.$ (S48)

Figure S2: (Left) Schematic of the time-reversal protocol for displacement
sensing. The ODF is turned on for a time $\tau$, after which a displacement
$\beta$ is quickly applied. Finally, the second ODF pulse is applied for a
time $\tau$. (Right) Sensitivity $(\delta\beta)^{2}$ for experimentally
realistic parameters $N=150$, $g=2\pi\times 3.91$ kHz, $\bar{n}=5$,
$\Gamma_{\text{el}}=610\text{ s}^{-1}$ and $\sigma=2\pi\times 40$ Hz. We
compare the exact result obtained from numerics (light blue) versus the
perturbative expression Eq. (S48) (gray) and see that the optimal sensitivity
is captured correctly by the perturbative result.

We have indicated the nature of each of the terms involved in the sensitivity.
The first term corresponds (when setting $\Gamma=0$) to the ideal result
obtained in Eq. (S9). The next term arises from a reduction in the signal due
to the presence of $\sigma$. Nevertheless, it usually is the smallest
contribution. More relevant is the third term, which is a manifestation of the
remaining spin-phonon correlations at the end of the protocol as a consequence
of the imperfect echo. It has the typical $\bar{n}+1/2$ factor associated to
the size of the boson wavefunction. Finally, the imperfect echo also results
in spin squeezing, which creates spin-spin correlations and an increase in the
variance of $\hat{J}_{y}$. This is the content of the last term. We compare
Eq. (S48) against exact numerics for the sensitivity in Fig. S2 (right panel).
There is good agreement at short times and, importantly, the optimal
sensitivity is captured correctly. Discrepancies appear only at longer times.

From Eq. (S48) we obtain Eq. (5) in the main text by neglecting the signal
reduction term, which we find does not affect any of the results presented.
The application of the perturbative formulas is justified for Fig. 3(a) in the
main text since the signal-to-noise ratio is computed at the time of optimal
sensitivity.

###  S3.4 Electric Field Sensing

Here, we are interested in sensing the drive strength itself [$\eta(t)$ in Eq.
(S31)] since it is proportional to the applied electric field. We want the
drive to act for as long as possible, so we keep $\eta(t)=\eta$ constant
throughout the full evolution and add ODF evolutions for short spans of time
at the end and/or the beginning of the sensing protocol. We will consider two
schemes here, which we dub as “Classical Protocol” and “Quantum Protocol”. In
the former, the ODF interaction is only used as a readout mechanism at the end
of the sequence, closely following the discussion in Sec. S2.3, thus building
spin-phonon entanglment even under ideal conditions and limiting the
achievable sensitivity. We improve this situation in the “Quantum Protocol” by
adding an extra ODF pulse at the beginning that creates an initially entangled
state, and allows for an enhancement in sensitivity.

####  S3.4.1 Classical Protocol

The scheme is shown in Fig. S3 (left) with total protocol duration $T$ and the
ODF interaction is only applied at the end of the protocol for a time $\tau$.
Again, we calculate $h(T)$, $p(T)$ and $q(T)$:

$|h(T)|^{2}=\frac{4g^{2}}{\delta^{2}}\sin\bigg{(}\frac{\delta\tau}{2}\bigg{)}^{2}\hskip
34.14322ptp(T)=\frac{g^{2}}{\delta^{2}}[\sin(\delta\tau)-\delta\tau]\hskip
34.14322ptq(T)=\frac{\eta g}{\delta^{2}}\big{\\{}\cos(\delta
T)-\cos\big{[}\delta(T-\tau)\big{]}\big{\\}}$ (S49)

As in the previous case, we numerically average the expectation values over
$\delta$ and calculate the sensitivity according to

$(\Delta\eta)^{2}=\frac{\braket{\hat{J}_{y}^{2}(T)}^{\text{av}}}{\big{(}\partial\braket{\hat{J}_{y}(T)}^{\text{av}}/\partial\eta\big{)}^{2}}.$
(S50)

To second order in $\sigma$, we find after averaging that

$\displaystyle\begin{split}\braket{\hat{J}_{y}^{2}(T)}^{\text{av}}&\approx\frac{N}{4}+\frac{Ng^{2}\tau^{2}e^{-\Gamma\tau}}{8}\bigg{[}4(\bar{n}+1/2)-\frac{\sigma^{2}\tau^{2}}{3}(\bar{n}+1/2)+\frac{g^{2}\sigma^{2}\tau^{4}}{18}\bigg{]}\\\
\frac{\delta_{\eta}\braket{\hat{J}_{y}(T)}^{\text{av}}}{\eta}&\approx\frac{\sqrt{N}e^{-\Gamma\tau/2}g\tau(\tau-2T)}{2}\bigg{[}1-\frac{\sigma^{2}}{12}(2T^{2}-2\tau
T+\tau^{2})\bigg{]}.\end{split}$ (S51)

Hence the sensitivity is

$\displaystyle\begin{split}(\delta\eta)^{2}_{\text{classical}}&\approx\underbrace{\frac{e^{\Gamma\tau}}{g^{2}\tau^{2}(2T-\tau)^{2}}}_{\text{Ideal
result}}\hskip
28.45274pt+\underbrace{\frac{2\bar{n}+1}{(2T-\tau)^{2}}}_{\text{Ideal spin-
phonon correlations}}\\\\[10.0pt]
&+\underbrace{\frac{\sigma^{2}e^{\Gamma\tau}(2T^{2}-2\tau
T+\tau^{2})}{6g^{2}\tau^{2}(2T-\tau)^{2}}}_{\text{Signal reduction}}\hskip
14.22636pt+\underbrace{\frac{\sigma^{2}(\bar{n}+1/2)}{6}}_{\text{Spin-phonon
correlations}}+\underbrace{\frac{g^{2}\sigma^{2}\tau^{4}}{36(2T-\tau)^{2}}}_{\text{Spin-
spin correlations}}.\end{split}$ (S52)

Figure S3: Theoretical schematics of the electric field sensing protocols. The
classical protocol is shown in the left, with a drive $\eta$ applied for a
time $T$ and the ODF interaction used as a read-out and applied for a time
$\tau$ in the end. The quantum protocol is shown in the right and includes an
extra application of the ODF Hamiltonian for a time $\tau$ at the beginning.

As in the analysis of the displacement protocol, we have separated the
different contributions to the sensitivity according to their origin. The
first line includes effects present in the absence of COM frequency
fluctuations ($\sigma=0$): the first term is the increase in sensitivity due
to the acquisition of a bigger signal. However, as discussed in section S2.3,
this is accompanied by an increase in variance at the end of the protocol due
to spin-phonon entanglement, and this creates the second term in the first
line. The second line includes the effects due to COM frequency fluctuations
and are of a similar nature as those present in the displacement protocol.
Note that

$(\delta\eta)^{2}_{\text{classical}}>\frac{1}{4T^{2}}=(\delta\eta)^{2}_{\text{SQL}},$
(S53)

and so we never beat the SQL with this scheme. The achievable sensitivity is
even much worse:

$(\delta\eta)^{2}_{\text{classical}}>\frac{2\bar{n}+1}{4T^{2}}=(\delta\eta)^{2}_{\text{th}},$
(S54)

which includes thermal effects. We thus see that thermal occupation of the
bosons very negatively impacts the “Classical Protocol”.

Figure S4: We show the optimal ODF interaction time $\tau$ (left) and the
resulting optimal sensitivity (right) as a function of the protocol time $T$
for $g=2\pi\times 3.88$ kHz, $N=150$, $\bar{n}=5$, $\Gamma=520\text{ s}^{-1}$,
and $\sigma=2\pi\times 40$ Hz. For short $T$, $\tau$ is exactly half of the
protocol time, especially in the quantum protocol. For large $T$, $\tau$
saturates. The right panel shows that the classical protocol is generically
about 10 dB worse than the quantum one. We also include the SQL for
comparison, defined as $\delta\eta^{2}_{\text{SQL}}=1/(4T^{2})$, and note that
the quantum protocol does get below it.

####  S3.4.2 Quantum Protocol

The scheme is shown in Fig. S3 (right), with total protocol duration $T$ and
the ODF interaction is applied both at the beginning and at the end of the
protocol for a time $\tau$ in each section. As in the previous case,
$\eta(t)=\eta$, a constant throughout the evolution. Again, we calculate
$h(T)$, $p(T)$ and $q(T)$:

$\displaystyle\begin{split}|h(T)|^{2}&=\frac{16g^{2}}{\delta^{2}}\sin\bigg{(}\frac{\delta\tau}{2}\bigg{)}^{2}\sin\bigg{[}\frac{\delta(T-\tau)}{2}\bigg{]}^{2}\\\\[8.0pt]
p(T)&=-\frac{2g^{2}}{\delta^{2}}\bigg{\\{}\delta\tau-\sin(\delta\tau)-2\sin\bigg{(}\frac{\delta\tau}{2}\bigg{)}^{2}\sin[\delta(T-\tau)]\bigg{\\}}\\\\[8.0pt]
q(T)&=-\frac{4g\eta}{\delta^{2}}\sin\bigg{(}\frac{\delta\tau}{2}\bigg{)}\sin\bigg{[}\frac{\delta(T-\tau)}{2}\bigg{]}\cos\bigg{(}\frac{\delta
T}{2}\bigg{)}.\end{split}$ (S55)

As in the previous case, we numerically average the expectation values over
$\delta$ and calculate the sensitivity according to

$(\delta\eta)^{2}=\frac{\braket{\hat{J}_{y}^{2}(T)}^{\text{av}}}{\big{(}\partial\braket{\hat{J}_{y}(T)}^{\text{av}}/\partial\eta\big{)}^{2}}.$
(S56)

To second order in $\sigma$, we find after averaging that

$\displaystyle\begin{split}\braket{\hat{J}_{y}^{2}(T)}^{\text{av}}&\approx\frac{N}{4}+\frac{Ng^{2}\sigma^{2}\tau^{2}e^{-2\Gamma\tau}}{8}\bigg{[}4(\bar{n}+1/2)(T-\tau)^{2}+2g^{2}\tau^{2}\bigg{(}T-\frac{4\tau}{3}\bigg{)}^{2}\bigg{]}\\\
\frac{\delta_{\eta}\braket{\hat{J}_{y}(T)}^{\text{av}}}{\eta}&\approx-\sqrt{N}e^{-\Gamma\tau}g\tau(T-\tau)\bigg{[}1-\frac{\sigma^{2}}{12}(2T^{2}-\tau
T+\tau^{2})\bigg{]}.\end{split}$ (S57)

Hence the sensitivity is

$\displaystyle\begin{split}(\delta\eta)^{2}_{\text{quantum}}&\approx\underbrace{\frac{e^{2\Gamma\tau}}{4g^{2}\tau^{2}(T-\tau)^{2}}}_{\text{Ideal
result}}\\\\[10.0pt] &+\underbrace{\frac{\sigma^{2}e^{2\Gamma\tau}(2T^{2}-\tau
T+\tau^{2})}{24g^{2}\tau^{2}(T-\tau)^{2}}}_{\text{Signal
reduction}}+\underbrace{\frac{\sigma^{2}(2\bar{n}+1)}{4}}_{\text{Spin-phonon
correlations}}+\underbrace{\frac{g^{2}\sigma^{2}\tau^{2}(3T-4\tau)^{2}}{36(T-\tau)^{2}}}_{\text{Spin-
spin correlations}}.\end{split}$ (S58)

We have again included the effects of COM frequency fluctuations only in the
second line. Note that, in contrast to the “Classical Protocol”, the ideal
result is not bounded (when $\Gamma=0$) and thus can surpass the SQL [Eq.
(S53)] and the thermal bound [Eq. (S54)], irrespective of the thermal
occupation of the bosons. The COM frequency dependent contributions are of a
similar nature as those in the displacement sensing protocol.

In both the classical and quantum electric field sensing protocols we obtain
the sensitivity as a function of two variables, $T$ and $\tau$. For each $T$
we numerically optimize over $\tau$ for realistic values of $g$, $N$ and
$\sigma$. We show the optimal value of $\tau$ as a function of $T$ in Fig. S4
(left) and the optimal sensitivity in Fig. S4 (right). For the parameters used
in the main text, the perturbative expressions for
$(\delta\eta)_{\text{classical}}$ and $(\delta\eta)_{\text{quantum}}$ are
accurate only for $T\leq 0.6$ ms so the construction of the theory curves in
Fig. 4(b) of the main text required us to use the full numerical treatment.

##  S4 Experimental details

Here we outline some details of the Penning trap system used for the
measurements presented in the main text [8, 6]. Subsequent sections in this
supplementary material will address experimental calibrations and data
analysis.

The Penning trap consists of a stack of cylindrical electrodes (see cross-
sectional schematic in Fig. S5) which are used to generate an axially
confining electrostatic potential. The trap is inserted into the bore of a
$B=4.5$ T superconducting magnet oriented vertically such that the magnetic
field is aligned with the $\hat{z}-$axis defined by the trap. Radial
confinement of the ions results from the $E\times B$ induced rotation of the
ions at a frequency $\omega_{r}$. In a frame rotating at $\omega_{r}$, the
confining potential is

$q\Phi_{trap}=\frac{1}{2}M\omega_{z}^{2}(Z^{2}+\beta_{r}r^{2}),$ (S59)

where $Z$ is the axial distance and $r$ the cylindrical radius from the trap
center, $w_{z}/(2\pi)=1.59$ MHz is the COM mode frequency, which characterizes
the strength of the confining potential, and
$\beta_{r}=\frac{\omega_{r}(\Omega_{c}-\omega_{r})}{\omega_{z}^{2}}-\frac{1}{2}$
is the relative strength of the radial confinement for rotation frequency
$\omega_{r}$ and cyclotron frequency $\Omega_{c}/(2\pi)=7.6$ MHz. The rotation
frequency of the crystal is precisely with a rotating dipole electric field
[9] and, for the $N\approx 150$ ions typical of the work described in the main
text, set to $\omega_{r}/(2\pi)\approx 180$ kHz. At this rotation frequency,
the crystal is a 2-D single-plane array of ions.

Figure S6 shows the relevant energy levels for 9Be+ at high magnetic field
($4.5~{}T$). Though 9Be+ has nuclear spin $I=3/2$, we optically pump the
nuclear spin to the $m_{I}=+3/2$ level, where it remains throughout the
duration of an experiment and which allows for simplifying the level structure
as shown in Fig. S6. We define the qubit states as two levels in the $S_{1/2}$
manifold: $\ket{\uparrow}\equiv\ket{m_{J}=+1/2}$ and
$\ket{\downarrow}\equiv\ket{m_{J}=-1/2}$. The qubit is, essentially, the
${}^{9}\mathrm{Be}^{+}$ valence electron spin in the ground state. For a
magnetic field of $B=4.46$ T, the splitting between the qubit states
$\ket{\uparrow}$ and $\ket{\downarrow}$ is approximately 124 GHz. A low phase
noise microwave source at 124 GHz allows for global rotations of the spins.
The duration of a $\pi-$pulse (the time to flip $\ket{\downarrow}$ to
$\ket{\uparrow}$) is $45~{}\mu$s.

Doppler cooling and state-selective readout is performed on the cycling
transition $\ket{\uparrow}\rightarrow{}^{2}P_{3/2}(m_{J}=+3/2)$. State
initialization is performed with optical pumping into the $\ket{\uparrow}$
state with the Doppler cooling and repump lasers (see Fig. S6). For the work
presented here, all readout is made via projective measurements of the ion
spin states through state-dependent global resonance fluorescence. Photons are
collected on a photomultiplier tube (PMT) from the global fluorescence of ions
in the $\ket{\uparrow}$ state in a direction perpendicular to the magnetic
field (the ‘side view’ in Fig. S5). The PMT counts are integrated over $2.5$
ms to reduce the impact of photon shot noise.

Figure S5: Illustration of the Penning trap electrodes and imaging. The view
is cross-sectional and the electrodes are cylindrically symmetric. Doppler
cooling is performed along the trap axis and perpendicular to it. Photons are
imaged on to the ‘side view’ PMT for global state-dependent readout of the
spins. Figure S6: Relevant energy levels for 9Be+ at $B=4.5$ T. The splitting
between $\ket{\uparrow}$ and $\ket{\downarrow}$ is 124 GHz. All lasers are in
the UV with $\lambda\approx 313$ nm. Doppler cooling and state-dependent
projective measurements are performed on the cycling transition
$\ket{\uparrow}\rightarrow{}^{2}P_{3/2}(m_{J}=+3/2)$. A repump laser pumps
spin population out of the $\ket{\downarrow}$ state and prepares the state
$\ket{\uparrow}$ with high fidelity. A pair of off-resonant lasers couple the
qubit states to states in the ${}^{2}P_{3/2}$ manifold, producing a spin-
dependent optical-dipole force.

##  S5 Experimental calibrations

This section will address the various experimental calibrations that provide
input to the theoretical expressions.

###  S5.1 Applied electric field calibration

The COM mode is driven by applying a resonant rf voltage to a nearby trap
electrode. To calibrate the strength of this driving force at the ions, we
measure the axial displacement of the ion crystal when a large DC voltage is
applied. Applying hundreds of volts displaces the crystal by several microns,
which we measure on a sideview imaging system with a magnification of $11.6\pm
0.2$. Figure S7 shows these measured displacements over a range of voltages,
and the linear fit used to extract the voltage-to-displacement calibration of
$12.9\pm 0.46$ nm/V.

Figure S7: Voltage-to-displacement calibration of a nearby trap electrode. A
linear fit to the measured axial displacement of the crystal versus applied DC
voltage gives a calibration of $12.9\pm 0.46$ nm/V.

Since the ions are confined axially in a quadratic potential, this
displacement is linearly related to the force applied to the ions

$F_{D}=-m_{\rm Be}\omega_{z}^{2}z.$ (S60)

Through Eq. S60 and the voltage-to-displacement calibration, we are able to
calculate the driving force applied to the ions for a given rf voltage and
therefore the amplitude of the COM mode that is excited for a given drive
duration. The uncertainty in the voltage-to-displacement calibration is the
largest systematic uncertainty in these experiments.

###  S5.2 Spin depolarization due to photon scattering

Figure S8: Depolarization due to photon scattering from the off-resonance
laser beams that produce the ODF. The black points are data for the
decoherence rate due to ODF beams. The blue lines are fits to the data. In
this case, $\Gamma_{tot}=250\pm 3~{}s^{-1}$, respectively.

An important source of noise that must be calibrated is the spin
depolarization due to photon scattering from the off-resonance ODF laser
beams. For the conditions of the NIST Penning trap set-up, Rayleigh scattering
is the dominant source of decoherence over Raman scattering [10]. Typically,
$\Gamma_{el}\approx 4\Gamma_{ram}$. As a result, the dominant effect of photon
scattering here is spin depolarization rather than spin-flip errors. To
empirically measure this depolarization rate, a spin-echo Ramsey sequence is
used with the ODF turned on during the two free evolution periods each of time
$\tau$ (for a total interaction time of $2\tau$). The ODF detuning $\delta$ is
set $50$ kHz above the COM mode. In Fig. S8 the contrast is the measured
length of the composite Bloch vector of the spins. The contrast is given by
$\braket{\hat{\sigma}_{x}}=\exp(-2\Gamma_{\mathrm{tot}}\tau)$, where
$\Gamma_{\mathrm{tot}}=\frac{1}{2}(\Gamma_{ram}+\Gamma_{el})$. In the main
text and elsewhere in this supplementary material, the simplifying convention
$\Gamma\approx\Gamma_{el}$ is used. Figure S8 is an example of such a dataset.
As a result of the ACSS being nulled as well as the amplitudes stabilized for
each beam, decoherence due to laser power fluctuations is negligible.

###  S5.3 Spin-motion coupling strength

To couple the spin and motional degrees of freedom, a spin-dependent optical
dipole force (ODF) is generated by a pair of off-resonant laser beams. The two
beams intersect at the ions and create a one-dimensional traveling wave
potential at a frequency $\mu$ equal to the resulting beatnote frequency. The
relevant Hamiltonian for this interaction is:

$H_{\text{SDF}}=F_{0}\cos(\mu t)\sum_{i}\hat{z}_{i}\hat{\sigma}_{i}^{z},$
(S61)

where $\hat{\sigma}_{i}^{z}$ is the $z$-component of the Pauli spin matrix for
ion $i$, and $\hat{z}_{i}$ is the axial position of ion $i$. Details of the
generation of this spin-dependent force are given in the supplementary
information of [11, 12]. The axial motion of the crystal can be decomposed
into $N$ modes with frequencies $\omega_{m}$ (with $\omega_{z}$, the COM
frequency, being the highest) and normalized eigenfunctions $\vec{b}_{m}$
(with components $b_{m}^{i}$) that can be calculated as detailed in the
supplementary material of [11].

We measure the strength of $F_{0}$ by tuning the spin-dependent force to a
frequency close to the COM frequency $\omega_{z}$ such that the detuning
$\delta=\mu-\omega_{z}=2\pi/\tau$ (where $\tau$ is the interaction time) is
chosen to decouple the spin and COM motion degrees of freedom at the end of
the sequence [13]. This stroboscopically realizes a pure Ising spin
interaction. The effect of the remaining modes (other than the COM) is also
captured by an Ising interaction by virtue of their large detuning with
respect to the spin-dependent force frequency. The net effect of the ODF on
the spins is then well described by the Ising Hamiltonian

$\hat{H}_{\text{Ising}}=\frac{\hbar}{N}\sum_{i<j}J_{i,j}\hat{\sigma}_{i}^{z}\hat{\sigma}_{i}^{z}$
(S62)

where the amplitudes $J_{i,j}$ are related to the mode parameters through

$J_{i,j}=\frac{F_{0}^{2}}{4\hbar\,m_{Be}}\sum_{m}\frac{Nb_{m}^{i}b_{m}^{j}}{\omega_{m}(\mu-\omega_{m})}.$
(S63)

The only unknown quantity is $F_{0}$, the strength of the spin-dependent
force. We calibrate it through a measurement of the so-called mean-field spin
precession, the spin precession induced by the spin-motion coupling in a
regime where the spins decouple from the motion and the induced spin-spin
interaction is weak, resulting in minimal spin correlations. Under this
approximation, each spin precesses about an effective magnetic field created
by the rest of the spins:
$\hat{H}_{\text{mf}}=\sum_{i}B_{i}\hat{\sigma}_{z}^{i}/2$, with
$B_{i}=2\hbar\sum_{j\neq i}J_{i,j}\braket{\hat{\sigma}_{z}^{j}}/N$. By
applying the Hamiltonian in Eq. (S62) and measuring the mean-field precession
of the collective Bloch vector $\braket{\hat{J}}$ we can then extract $F_{0}$
[11, 12].

In the main text we concentrate only on the COM mode and so the spin-motion
Hamiltonian reduces to

$\hat{H}=\hbar\,\delta\hat{a}^{\dagger}\hat{a}+\frac{\hbar
g}{\sqrt{N}}(\hat{a}^{\dagger}+\hat{a})\hat{J}_{z}.$ (S64)

In this notation, the COM mode coupling strength is $g=F_{0}z_{0}/\hbar$ and
$z_{0}=\sqrt{\hbar/(2m_{\rm Be}\omega_{z})}$. For this work, typical values
range from $g/(2\pi)=3.7$ to $4.1$ kHz. Specifically, for the amplitude
sensing experiment in the main text $g/(2\pi)=3.91$ kHz, while $g/(2\pi)=3.88$
kHz for electric field sensing.

###  S5.4 Center-of-mass frequency fluctuations

Figure S9: a) Sequence used to measure the COM mode phonon occupation and
frequency fluctuations. b) Scan over COM mode lineshape, used to determine the
phonon occupation of the mode. By fitting to the data, a value of
$\bar{n}=4.4$ is extracted. c) Measurement of COM frequency fluctuations by
nominally setting $\delta=0$ and measuring the increase of $<P_{\uparrow}>$
above the expectation from spin depolarization as a function of the protocol
duration. The orange line is a fit to the data (black points), yielding rms
frequency fluctuations of $\sigma/(2\pi)=30$ Hz. The blue dashed line is
theory including only spin depolarization due to photon scattering.

The frequency fluctuations of the center-of-mass mode were determined by
measuring the inability of a time reversal to undo the spin-motion
entanglement in the sequence shown in Fig. S9(a). If no spin-dependent force
is applied the spins are rotated to the down state at the end of the sequence,
resulting in zero bright fraction. Figure S9(b) shows the signal line shape as
the frequency $\mu$ of the spin-dependent force is swept across the center-of-
mass mode frequency $\omega_{z}$. On resonance ($\delta=\mu-\omega_{z}$) the
spin-motion entanglement generated in the first pulse is taken out by the
second pulse, resulting in a sharp dip in the line shape.

However, if the COM mode frequency changes from one experimental trial to the
next, this cancellation will be imperfect. We characterize this effect by
measuring $\langle P_{\uparrow}\rangle$ with $\delta$ nominally set to zero
while varying the duration of the ODF (Fig. S9(c)). Frequency fluctuations of
the COM mode will result in the bright fraction $\langle P_{\uparrow}\rangle$
increasing above the expected background due to photon scattering. We obtain
the amplitude of frequency fluctuations by fitting this increase to
$\braket{P_{\uparrow}}=\frac{1}{2}[1-\braket{\hat{J}_{x}(2\tau)}]$, where
$\braket{\hat{J}_{x}(2\tau)}$ is obtained from Eqs. (S44) . Since the sequence
of steps used in the calibration is essentially the same as the displacement
sensing protocol without the displacement (see Fig. S2), $|h(2\tau)|^{2}$ and
$p(2\tau)$ are given by Eq. (S45), giving us the formula

$\braket{P_{\uparrow}}=\frac{1}{2}-\frac{1}{2}\frac{e^{-2\Gamma_{\mathrm{tot}}\tau}}{\sqrt{1+\frac{g^{2}\sigma^{2}\tau^{4}(2\bar{n}+1)}{N}+\frac{4g^{4}\sigma^{2}\tau^{2}}{9N}}},$
(S65)

which we find agreement with at both short and long times. The size of the
measured frequency fluctuations varied slowly over a time period of an hour
and from day to day. Typically, the rms amplitude of frequency fluctuations
obtained through this fitting is $\sigma/(2\pi)\approx 40\pm 20$ Hz.

###  S5.5 Center-of-mass mode damping & ring-down

A potential source of noise that could impact the amplitude and electric field
measurements is damping of the COM mode. This might occur by energy transfer
out of the COM mode and to the other axial or in-plane modes. To assess
whether this is an issue for this experiment, we perform a ring down
experiment wherein a large coherent motion is excited by a drive applied to an
endcap electrode and the excited COM mode amplitude is subsequently measured
following a variable wait duration. The amplitude of the motion is calibrated
to be $Z_{c}\approx 50$ nm, about an order of magnitude larger than the
thermal axial extent of the COM mode
$z_{\mathrm{cm}}=\sqrt{\frac{\hbar}{2Nm_{\mathrm{Be}}\omega_{z}}(2\bar{n}+1)}=5.5$
nm, with $N=120$, $\omega_{z}/(2\pi)=1.59$ MHz, and $\bar{n}=4.6$. The motion
is read out by mapping it onto the spin state with a spin-dependent ODF and
performing a projective measurement of the spins. The phase of the driven
motion relative to the ODF is left uncontrolled and varies randomly from one
trial of the experiment to the next [5]. The protocol used in [5] measures a
single quadrature of the motion, but here we measure the rms sum of both
quadratures with the sequence shown in Fig. S10 so that the sensitivity to the
relative phase $\phi$ is removed and the signal-to-noise ratio is not limited
by fluctuations in $\phi$. Note that the signal is sensitive to the rms sum of
the spin precession in each arm of the sequence. By applying a $\pi/2$
rotation about the $x$-axis in the middle of the sequence followed by a
$\pi/2$ phase shift of the ODF, both quadratures are mapped into rotations of
the spins: $\theta_{1}=\theta_{\mathrm{max}}\cos(\phi)$ and
$\theta_{2}=\theta_{\mathrm{max}}\sin(\phi)$, where
$\theta_{\mathrm{max}}=F_{0}Z_{c}\tau\hbar$ with $\tau$ the duration of each
ODF arm. The signal, then, is
$\left<P_{\uparrow}\right>=\frac{1}{2}\left[1-e^{-\Gamma_{\mathrm{tot}}2\tau}\left<\cos(\theta_{1})\cos(\theta_{2})\right>\right]\approx\frac{1}{2}\left[1-e^{-2\Gamma_{\mathrm{tot}}\tau}\cos(\theta_{\mathrm{max}})\right]$,
where this approximation is valid for small angles. Figure S10(b) shows the
results of an experiment where the previously described sequence is applied
and the wait duration following excitation of the COM mode is varied.
Typically, no reduction in the signal is observed for over $100$ ms. By
fitting to this data with an exponential decay, a value of $1/\kappa=300$ ms
is extracted, where $\kappa$ is the exponential decay rate. Calculations show
that this weak damping is sufficiently small that it does not impact the
displacement and electric field sensitivity studied in this manuscript.

Figure S10: Experimental data for the ring down sequence described in the
text. An exponential decay is fit to the data and yields a decay rate of
$\kappa=1/(300\>\mathrm{ms})$.

###  S5.6 Limits on stochastic electric fields

Ion traps, in particular small rf ion traps, are known to exhibit stochastic
electric fields that heat the motional degrees of freedom of a trapped ion
crystal. These stochastic electric fields could be due to noisy external
circuits connected to the trap electrodes. In general, this source of noise
can be filtered, but there remains an anomalous source of stochastic electric
fields due to fluctuating potentials on the surfaces of the trap electrodes
[14]. Measurements of the ion crystal heating rate can be used to place limits
on the size of a stochastic electric field sensed by the trapped-ion crystal.

The heating rate of the COM mode was determined by applying the sequence
discussed is Sec. S5.D to measure the COM mode phonon occupation, but with a
variable wait time after the state preparation. The mode will heat during this
wait, and by measuring the mode occupation for various wait times the heating
rate can be extracted. Figure 13 shows an example of this heating rate
measurement following EIT cooling to near the motional ground state. A fit to
this data yields a heating rate of $58\pm 5$ quanta/s. The measured heating
rate would vary from day to day, with an upper bound of $100$ quanta/s. Note
that the COM mode heating rate due to electric field noise increases linearly
with the number of ions. Therefore, the upper bound of 100 quanta/s performed
with $N\approx 100$ ions provides an upper bound on a single-ion heating rate
due to stochastic electric fields of $1$ quanta/s.

The $1$ quanta/s heating rate, if entirely due to stochastic electric fields,
corresponds to a power spectral density of electric field noise of
$S_{E}(2\pi\times 1.6\mathrm{MHz})\approx 0.25\times
10^{-14}~{}(\mathrm{V/m})^{2}~{}\mathrm{s}$ [14] or a $50$ nV
m${}^{-1}/\sqrt{\mathrm{Hz}}$ stochastic electric field spectral density. We
emphasize that this is an upper limit, as the measured heating rate could be
due to sources other than stochastic electric fields. Our demonstrated $240$
nV m${}^{-1}/\sqrt{\mathrm{Hz}}$ sensitivity for measuring a coherent electric
field is large compared to this stochastic electric field upper bound.

Figure S11: By measuring the COM mode occupation after a variable wait, the
heating rate can be extracted. Blue points are data following a fit to the COM
mode lineshape. The orange line is a fit to the data and gives a heating rate
of $58\pm 5$ quanta/s.

###  S5.7 Excess noise

To measure any excess experimental noise beyond projection noise, we perform
the displacement sensing experiment detuned $50$ kHz from the COM mode where
the effects of thermal and frequency fluctuations of this mode are negligible.
In this off-resonant experiment, an $18\%$ increase in the noise above
projection noise is measured. This excess noise is reduced as the strength of
the ODF is decreased, which points towards effects due to departure from the
Lamb-Dicke regime as the likely cause. For short ODF arm times, this noise
dominates over COM mode frequency fluctuations. When comparing theory and
experimental data, we have reduced the signal-to-noise of the theory by
including this excess noise in the theory curves.

##  S6 Experimental analysis

To compare with theory, we relate the unitless parameter $\beta$ to a
displacement $Z_{c}$ through $\beta=Z_{c}\sqrt{N}/(2z_{0})$, where
$z_{0}=\sqrt{\hbar/(2mw_{z})}$ is the ground state wavefunction size for a
single ion and $Z_{c}=F_{D}T/(2m_{Be}w_{z})$ is the zero-to-peak amplitude of
the COM motion after a drive interval $T$ with a force per ion $F_{D}$. The
force $F_{D}$ or equivalently the applied electric field at the ions is
calibrated as discussed in Sec. S5. A. The signal-to-noise ratio for a single
experiment $Z_{c}/\Delta Z_{c}\equiv\beta/\Delta\beta$ is extracted by
repeatedly measuring $\hat{J}_{y}$, which is related to the spin population in
the $\ket{\uparrow}$ state by
$\braket{P_{\uparrow}}=\frac{1}{2}[1-\braket{\hat{J}_{y}}]$. We take the
difference between pairs of subsequent measurements with a $\pi$ phase shift
in the amplitude $\beta$ such that $\hat{J}_{y}^{(1)}=-\hat{J}_{y}^{(2)}$. The
signal-to-noise data (Fig. 3(a) of the main text) is produced by taking the
average of the measurement pair differences
${\braket{\hat{J}_{y}^{(1)}}-\braket{\hat{J}_{y}^{(2)}}}$, dividing by the
standard deviation of the set of difference measurements, and further dividing
by $\sqrt{2}$ to produce the single-measurement signal-to-noise ratio. The
experimental single-measurement signal-to-noise is defined as

$\frac{Z_{c}}{\delta
Z_{c}}|_{\mathrm{single,exp}}\equiv\frac{\beta}{\Delta\beta}=\frac{\braket{J_{y}^{(1)}}-\braket{J_{y}^{(2)}}}{\sqrt{2}~{}\sigma(J_{y}^{(1)}-J_{y}^{(2)})}.$
(S66)

The variance $(\Delta\beta)^{2}$ of a single measurement of $\beta$, which is
plotted in Fig. 3(b) of the main text, is calculated from the experimentally
determined signal-to-noise ratio for a single measurement of an amplitude
$Z_{c}$ through the expression

$(\Delta\beta)^{2}=\frac{1}{4}\left(\frac{\sqrt{N}Z_{c}}{z_{0}\frac{Z_{c}}{\delta
Z_{c}}|_{\mathrm{single}}}\right)^{2}$ (S67)

where $N$ is the number of ions and $z_{0}=\sqrt{\hbar/2m_{\rm
Be}\omega_{z}}$.

Similarly, $(\Delta\eta)^{2}$ (see Fig. 4(b) of the main text) is calculated
through

$(\Delta\eta)^{2}=\frac{1}{4}\left(\frac{\sqrt{N}Z_{c}}{z_{0}T\frac{Z_{c}}{\delta
Z_{c}}|_{\mathrm{single}}}\right)^{2},$ (S68)

where $T$ is the duration of the spin-independent drive. The standard
deviation of a single measurement of the electric field, then, is

$\Delta\epsilon=\frac{2Z_{c}m_{\rm Be}\omega_{z}}{qT\frac{Z_{c}}{\delta
Z_{c}}|_{\mathrm{single}}}.$ (S69)

## References

  * Lewis-Swan  _et al._ [2020] R. J. Lewis-Swan, D. Barberena, J. A. Muniz, J. R. K. Cline, D. Young, J. K. Thompson, and A. M. Rey, Protocol for precise field sensing in the optical domain with cold atoms in a cavity, Phys. Rev. Lett. 124, 193602 (2020).
  * Zurek [2001] W. H. Zurek, Sub-planck structure in phase space and its relevance for quantum decoherence, Nature 412, 712 (2001).
  * Toscano  _et al._ [2006] F. Toscano, D. A. R. Dalvit, L. Davidovich, and W. H. Zurek, Sub-planck phase-space structures and heisenberg-limited measurements, Phys. Rev. A 73, 023803 (2006).
  * Walls and Milburn [2008] D. F. Walls and G. J. Milburn,  _Quantum Optics_ , 2nd ed. (Springer, 2008).
  * Gilmore  _et al._ [2017] K. A. Gilmore, J. G. Bohnet, B. C. Sawyer, J. W. Britton, and J. J. Bollinger, Amplitude sensing below the zero-point fluctuations with a two-dimensional trapped-ion mechanical oscillator, Phys. Rev. Lett. 118, 263602 (2017).
  * Affolter  _et al._ [2020] M. Affolter, K. A. Gilmore, J. E. Jordan, and J. J. Bollinger, Phase-coherent sensing of the center-of-mass motion of trapped-ion crystals, Phys. Rev. A 102, 052609 (2020).
  * Rey  _et al._ [2008] A. M. Rey, L. Jiang, M. Fleischhauer, E. Demler, and M. D. Lukin, Many-body protected entanglement generation in interacting spin systems, Phys. Rev. A 77, 052305 (2008).
  * Biercuk  _et al._ [2009] M. J. Biercuk, H. Uys, A. P. Vandevender, N. Shiga, W. M. Itano, and J. J. Bollinger, High-fidelity Quantum Control Using Ion Crystals in a Penning Trap, Quantum Info. Comput. 9, 920 (2009).
  * Huang  _et al._ [1998] X.-P. Huang, J. J. Bollinger, T. B. Mitchell, W. M. Itano, and D. H. E. Dubin, Precise control of the global rotation of strongly coupled ion plasmas in a Penning trap, Phys. Plasmas 5, 1656 (1998).
  * Uys  _et al._ [2010] H. Uys, M. J. Biercuk, A. P. VanDevender, C. Ospelkaus, D. Meiser, R. Ozeri, and J. J. Bollinger, Decoherence due to elastic Rayleigh scattering, Phys. Rev. Lett. 105, 200401 (2010).
  * Britton  _et al._ [2012] J. W. Britton, B. C. Sawyer, A. C. Keith, C.-C. J. Wang, J. K. Freericks, H. Uys, M. J. Biercuk, and J. J. Bollinger, Engineered two-dimensional Ising interactions in a trapped-ion quantum simulator with hundreds of spins, Nature 484, 489 (2012).
  * Bohnet  _et al._ [2016] J. G. Bohnet, B. C. Sawyer, J. W. Britton, M. L. Wall, A. M. Rey, M. Foss-Feig, and J. J. Bollinger, Quantum spin dynamics and entanglement generation with hundreds of trapped ions, Science 352, 1297 (2016).
  * Leibfried  _et al._ [2003] D. Leibfried, B. DeMarco, V. Meyer, D. Lucas, M. Barrett, J. Britton, W. M. Itano, B. Jelenković, C. Langer, T. Rosenband, and D. J. Wineland, Experimental demonstration of a robust, high-fidelity geometric two ion-qubit phase gate., Nature 422, 412 (2003).
  * Brownnutt  _et al._ [2015] M. Brownnutt, M. Kumph, P. Rabl, and R. Blatt, Ion-trap measurements of electric-field noise near surfaces, Rev. Mod. Phys. 87, 1419 (2015).