11institutetext: Centre for Medical Image Computing, Wellcome/EPSRC Centre for
Interventional and Surgical Sciences, Department of Medical Physics and
Biomedical Engineering, University College London, London, U.K.
22institutetext: State Key Laboratory of Robotics and System, Harbin Institute
of Technology, Harbin, China 33institutetext: School of Biomedical Engineering
& Imaging Sciences, King’s College London, London, U.K.

# Privileged Anatomical and Protocol Discrimination in Trackerless 3D
Ultrasound Reconstruction

Qi Li 11 Ziyi Shen 11 Qian Li 1122 Dean C. Barratt 11 Thomas Dowrick 11
Matthew J. Clarkson 11 Tom Vercauteren 33 Yipeng Hu 11

###### Abstract

Three-dimensional (3D) freehand ultrasound (US) reconstruction without using
any additional external tracking device has seen recent advances with deep
neural networks (DNNs). In this paper, we first investigated two identified
contributing factors of the learned inter-frame correlation that enable the
DNN-based reconstruction: anatomy and protocol. We propose to incorporate the
ability to represent these two factors - readily available during training -
as the privileged information to improve existing DNN-based methods. This is
implemented in a new multi-task method, where the anatomical and protocol
discrimination are used as auxiliary tasks. We further develop a
differentiable network architecture to optimise the branching location of
these auxiliary tasks, which controls the ratio between shared and task-
specific network parameters, for maximising the benefits from the two
auxiliary tasks. Experimental results, on a dataset with 38 forearms of 19
volunteers acquired with 6 different scanning protocols, show that 1) both
anatomical and protocol variances are enabling factors for DNN-based US
reconstruction; 2) learning how to discriminate different subjects (anatomical
variance) and predefined types of scanning paths (protocol variance) both
significantly improve frame prediction accuracy, volume reconstruction
overlap, accumulated tracking error and final drift, using the proposed
algorithm.

###### Keywords:

Freehand ultrasound Privileged information Multi-task learning.

##  1 Introduction

3D ultrasound (US) reconstruction is a promising technique both for
diagnostics and image guidance, such as image fusion with other image
modalities [5], registration with preoperative data during surgery [3], volume
visualisation and measurement [2]. Currently, most clinically applied 3D
reconstruction of 2D freehand US imaging use spatial tracking devices, such as
optical, electromagnetic or mechanical positioning [12]. Research for reducing
dependency on external devices in such 3D US reconstruction has been motivated
for its portability, accessibility and low-cost. Previous non-learning-based
methods utilised the speckle correlation between US frames [1, 3]. In recent
learning-based approaches, convolutional neural networks and their variants
have been proposed to use two adjacent frames as network input [10, 15, 16],
for predicting spatial transformation between them. More generally, sequential
modelling methods, e.g. using recurrent neural network [4, 7, 8, 11] and
transformer [14], have also been tested with the same goal of localising
relative positions between two or more US frames.

With the promising results from deep learning, one might question what was
learned in these data-driven approaches for predicting inter-frame
transformations. Specifically, in addition to speckle patterns (which holds
within a limited spatial scale), what are other factors that generated
correlation between US frames, such that their relative locations can be
inferred?

We first hypothesized that two factors, common anatomical characteristics
between subjects and predefined scanning protocols, are responsible for such
predictability in this application. A variance-reduction study is presented in
Sec. 3.2 to demonstrate that sufficient anatomical and protocol variance in
training data is indeed required for the reconstruction.

In this work, we propose to encode the anatomical and protocol patterns using
two classification tasks, discriminating between subjects and between types of
scanning paths, respectively. We then investigate methods to train these tasks
together with the main reconstruction task to improve the performance of the
main task.

In addition to US images as network input, previous studies have also
integrated additional information or signals, such as optical flow [21]
between US frames and acceleration, orientation, angular rates from inertial
measurement units (IMU) [8, 9]. Optical flow was derived from image sequence
itself, while the IMU-measured signals are required at both training and
inference. Different to these applications, the proposed discrimination task
labels, subject and protocol indices, are in general available during
training, but not required at inference, thus known as privileged information
[19], further discussed in Sec. 2.1.

In summary, our contributions include 1) a multi-task learning approach to
formulate two factors in freehand US as privileged information, for improving
reconstruction accuracy; 2) a mixture model formulation for optimisable
branching locations for auxiliary tasks, implemented with a differentiable
network and a gradient-based bi-level optimisation; 3) extensive experimental
results to quantify the improved performance due to the privileged tasks; and
4) open-source code and data 111https://github.com/ucl-candi/freehand for
public access and reproducibility.

##  2 Method

###  2.1 Preliminaries: Privileged Auxiliary Tasks and Shared Parameters

Assume a main task $f_{\theta}(y|x)$ that predicts ${y}$ with input data $x$
(here, using a $\theta$-parameterised neural network), which is optimised
alongside $\mathcal{J}$ auxiliary tasks
$f_{\theta^{s}_{j},\theta^{a}_{j}}(y^{a}_{j}|x)$, $j=1,...,\mathcal{J}$, where
$\theta^{s}_{j}$ are shared parameters (with the main task), and
$\theta^{a}_{j}$ are task-specific parameters, both for predicting
$y^{a}_{j}$. Therefore, $\theta^{s}_{j}\subseteq\theta$ and
$\theta^{a}_{j}\cap\theta=\emptyset$. When the goal is to predict $y$ from
$x$, the supervision and prediction of auxiliary task are only required during
training. The main task benefits from this privileged information, which is
not required at inference.

Such multi-task learning incorporates the privileged information but may
suffer from absolute negative transfer [20] \- here, when the additional
auxiliary tasks negatively impact the main task performance, and/or relative
negative transfer - one auxiliary task reduces the main-task-improving
potentials from the other auxiliary task(s). One approach for reducing
negative transfer, or optimising the transferability, is adjusting the ratio
between the shared $\theta^{s}_{j}$ and task-specific $\theta^{a}_{j}$
parameters [13]. As no parameters are shared (i.e.
$\theta^{s}_{j}=\emptyset$), no negative transfer is possible (although any
benefit from this auxiliary task also diminishes).

Assume a binary task descriptor $z_{j}$ (extended to a one-hot vector in Sec.
2.3) indicating the use of the $j^{th}$ auxiliary task. The main task is thus
conditioned on the use of the auxiliary task
$f_{\theta,\theta^{a}_{j}}(y|x,z_{j})$. Importantly, where to place the task
descriptor $z_{j}$ determines which and how many network parameters
$\theta^{s}_{j}$ are shared. As in Fig. 1 (a), the closer the task descriptor
is to the input, i.e. early branching, the less shared parameters.

(a)

(b)

Figure 1: Example of network architectures with (a) various possible branching
locations with different shared parameters and (b) two auxiliary tasks, each
modelled as a mixture of $\mathcal{I}$ candidate tasks.

###  2.2 Example Auxiliary Tasks: Anatomical and Protocol Discrimination

Each auxiliary task is trained with a cross-entropy loss, between predicted
class probability vectors $y^{a}_{j}$ and ground-truth targets $t^{a}_{j}$ in
one-hot vectors.

$\mathcal{L}^{CE}_{j}=-\sum_{N_{j}}(t^{a}_{j}\cdot
log[f_{\theta^{s}_{j},\theta^{a}_{j}}(y^{a}_{j}|x)])$ (1)

The two tasks, minimising $\mathcal{L}^{CE}_{j=1}$ and
$\mathcal{L}^{CE}_{j=2}$, classify different training subjects and types of
scanning protocols, with $N_{j=1}$ and $N_{j=2}$ number of classes,
respectively. An underlying assumption is that the anatomical and protocol
variance can impact the 3D US reconstruction performance. This assumption is
tested by quantifying the reconstruction accuracy changes, as number of
subjects and/or protocol types are reduced. Results are reported in Sec. 3.2.

###  2.3 Parameterised Task Descriptor Locations

Assume $\mathcal{I}$ locations in the main task network, at which a task
descriptor $z_{i,j}$ can be conditioned, where $i=1,...,\mathcal{I}$ for
$j^{th}$ auxiliary task. For a single main task in this case, each auxiliary
task can branch out from these $\mathcal{I}$ locations, as illustrated in Fig.
1 (b). The task descriptor thus represents the probability of branching
location, with additional constraints $z_{i,j}\in[0,1]$ and
$\sum_{i}z_{i,j}=1$.

With the task descriptor $z_{i,j}$, the $j^{th}$ auxiliary task is
parameterised by a mixture model of $\mathcal{I}$ candidate tasks
$f_{\theta^{s}_{i,j},\theta^{a}_{i,j}}(y^{a}_{i,j}|x)$, each performed by one
branch, where the additional subscript $i$ is the candidate task index.

$f_{\theta,\theta^{a}_{1,j},...,\theta^{a}_{\mathcal{I},j}}(y^{a}_{j}|x,z_{j})=\sum_{i=1}^{\mathcal{I}}[z_{i,j}\cdot
f_{\theta^{s}_{i,j},\theta^{a}_{i,j}}(y^{a}_{i,j}|x)]$ (2)

where $z_{j}=[z_{i,j}]^{\top}_{i=1,...,\mathcal{I}}$ is the task descriptors
for all the auxiliary tasks, with shared parameters
$\theta^{s}_{j}=[\theta^{s}_{i,j}]^{\top}_{i=1,...,\mathcal{I}}$ and task-
specific parameters
$\theta^{a}_{j}=[\theta^{a}_{i,j}]^{\top}_{i=1,...,\mathcal{I}}$. For $j^{th}$
auxiliary task, the final prediction is therefore the matrix product of all
candidate (branch) predictions
$y^{a}_{j}=[y^{a}_{i,j}]^{\top}_{i=1,...,\mathcal{I}}$ (a
$N_{j}\times\mathcal{I}$ matrix) and the location weights $z_{j}$ (a
$\mathcal{I}\times 1$ vector, generated by a softmax function). Multiplicative
task conditioning is also known as gating. The loss defined in Eq. 1 remains.

This formulation allows all candidate branches trained together without a
predefined architecture decision. We show in Sec. 2.5 that the task descriptor
may be considered as a hyperparameter and optimised using a gradient-based bi-
level optimisation algorithm, for efficient inference using a single branch.

###  2.4 The Main Task and Evaluation

We adopt the approach in our previous work [4] for the main task,
reconstructing a 3D US scan by sequentially predicting inter-frame
transformations with inputs of frame sequences. The same reconstruction loss
is used, but is now conditioned on the two auxiliary tasks, denoted as
$\mathcal{L}^{Rec}(t,f_{\theta,\theta^{a}_{1},\theta^{a}_{2}}(y|x,z_{1},z_{2}))$
where the network output is the six parameters for rigid spatial
transformation. It is important to clarify that this proposed loss function
also utilises a multi-task learning for predicting transformations between
nearby frames, but different to and independent of that used in this work.
Unless specified otherwise, the methodology and the implementation of the main
task network, based on EfficientNet (b1) [18], remain the same, with further
details in the original publication [4].

To test the generalisation and reconstruction performance of the proposed
method, four evaluation metrics are used. For each frame, frame prediction
accuracy ($\epsilon_{frame}$) is used to evaluate the generalisation of the
method, denoting the Euclidean distance between the ground-truth and
prediction- transformed four corner points of each frame. The scan
reconstruction performance is quantified by an accumulated tracking error
($\epsilon_{acc.}$), indicating the averaged point distance on all frame
pixels, a volume reconstruction overlap ($\epsilon_{dice}$), denoting the
overlap of all pixels between the ground-truth and prediction volume, and a
final drift ($\epsilon_{drift}$), denoting the frame prediction accuracy of
the last frame in a scan.

###  2.5 Bi-level Optimisation of Task Descriptor

Given a loss for the main task
$\mathcal{L}^{Rec}(t,f_{\theta,\theta^{a}_{1},\theta^{a}_{2}}(y|x,z_{1},z_{2}))$
and those for the two auxiliary tasks defined in Eq. 1, the overall loss
function thus is:

$\mathcal{L}(\theta,\theta^{a}_{1},\theta^{a}_{2},z_{1},z_{2}|D)=\mathcal{L}^{Rec}(\theta|x,t)+\sum_{j=1}^{2}\mathcal{L}^{CE}_{j}(\theta,\theta^{a}_{j},z_{j}|x,t^{a}_{1},t^{a}_{2})$
(3)

where the loss is rearranged as a function of relevant network parameters
$\\{\theta,\theta^{a}_{1},\theta^{a}_{2}\\}$ and task descriptors
$\\{z_{1},z_{2}\\}$, with observed data $D=\\{x,t,t^{a}_{1},t^{a}_{2}\\}$.
Given a training data set $\mathcal{D}_{train}$ and a validation data set
$\mathcal{D}_{val}$, the empirical losses are
$\mathcal{L}_{train}(\cdot)=\mathbf{E}_{D\in\mathcal{D}_{train}}[\mathcal{L}(\cdot|D)]$
and
$\mathcal{L}_{val}(\cdot)=\mathbf{E}_{D\in\mathcal{D}_{val}}[\mathcal{L}(\cdot|D)]$,
respectively. Optimising the task descriptors subject to optimised network
parameters leads to the following meta-learning task with a bi-level
optimisation problem: ^z_1, ^z_2 = argmin_z_1, z_2 L_val (^θ,^θ^a_1,^θ^a_2;
z_1, z_2)
s.t. ^θ,^θ^a_1,^θ^a_2 = argmin_θ,θ^a_1,θ^a_2 L_train (θ,θ^a_1,θ^a_2; z_1, z_2)

Since the $\frac{\partial\mathcal{L}_{val}}{\partial z_{j}}$ can be estimated,
the task descriptors can be optimised by gradient-based updates alternating
between minimising the two empirical loss functions [17]. The numerical
algorithm is summarised in Algorithm 1. The three tasks are weighed equally in
Eq. 3 as the task descriptor hyperparameters make explicit optimising or
predefining these weights redundant.

Algorithm 1: The bi-level optimisation algorithm.
---
1\. Initialise task descriptors $\\{z_{1},z_{2}\\}$ for two auxiliary tasks.
2\. Optimise the task descriptors using meta-learning algorithm in a
differentiable network:
while _not converged_ do
(1) Update network parameters {$\theta,\theta^{a}_{1},\theta^{a}_{2}$} for the
main task and all auxiliary tasks
by descending
$\nabla_{\theta,\theta^{a}_{1},\theta^{a}_{2}}\mathcal{L}_{train}(\theta,\theta^{a}_{1},\theta^{a}_{2};z_{1},z_{2})$
(2) Update task descriptors {$z_{1},z_{2}$}
by descending
$\nabla_{z_{1},z_{2}}\mathcal{L}_{val}((\theta,\theta^{a}_{1},\theta^{a}_{2})-\xi\nabla_{\theta,\theta^{a}_{1},\theta^{a}_{2}}\mathcal{L}_{train}(\theta,\theta^{a}_{1},\theta^{a}_{2};z_{1},z_{2});z_{1},z_{2})$.
3\. Finalise the network architecture using the optimised task descriptors
{$z_{1},z_{2}$}.
11footnotetext: In this work, the first-order approximation was used when
optimising task descriptors with $\xi=0$, i.e.,
$\nabla_{z_{1},z_{2}}\mathcal{L}_{val}(\theta,\theta^{a}_{1},\theta^{a}_{2};z_{1},z_{2})$
[6].

##  3 Experiments and Results

###  3.1 Dataset and Network Settings

The US data used in this study was from our previously published data set in
[4], acquired at 20 frame per second (fps) by an Ultrasonix machine (BK,
Europe) with a curvilinear probe (4DC7-3/40, 6 MHz, 9 cm depth and a median
level of speckle reduction) and an NDI Polaris Vicra (Northern Digital Inc.,
Canada) tracker. All US images with a size of 480×640 pixels, after spatial
and temporal calibration, acquired from 38 forearms of 19 volunteers were used
in this study, more than 40,000 US frames in total. For each forearm, three
predefined scanning protocols (straight line shape, ‘C’ shape and ‘S’ shape,
as in Fig. 2 (a)), in a distal-to-proximal direction, with the US probe
perpendicular of and parallel to the forearm, were acquired, resulting in 6
different protocols and 228 scans in total. The US scans have a various number
of frames, from 36 to 430 frames, equivalent to a probe travel distance of
between 100 and 200 mm. The data was split into train, validation and test
sets by a ratio of 3:1:1 on a scan level.

With the EfficientNet (b1) for the main reconstruction task [4], nine
locations were used for candidate branches, with each being a single fully-
connected classification layer, denoted as Branches 1-9 and Branch 1 being
closest to network input. The nine candidate predictions are weighted by the
softmax-generated task predictor, as described in Sec. 2.3, to form each final
auxiliary task prediction. The two auxiliary tasks resulted in total of 18
branches.

A minibatch of 32, an Adam optimizer were used for model training, with a
learning rate of $10^{-4}$ (tested among $\\{10^{-3},10^{-4},10^{-5}\\}$) and
the input sequence length being $M=\\{100,140\\}$ (tested among
$M=\\{49,75,100,140\\}$), selected based on the validation set performance.
Each model was trained for at least 15,000 epochs until convergence, for up to
4 days, on Ubuntu 18.04.6 LTS with a single NVIDIA Quadro P5000 GPU card. The
model with the best validation set performance was selected to evaluate test
set performance. Other hyperparameters were found relatively insensitive to
model performance and configured empirically based on validation performance.

###  3.2 The Effect of Anatomical and Protocol Variances

The impact on performance of the main task was tested by training models with
various anatomical and protocol variance in the training data: 1) all training
data (All); 2) straight only (Straight); 3) C-shape and S-shape (C-S); 4) 25%
subjects in training data (Sub 25%); 5) 50% training subjects (Sub 50%); 6)
75% training subjects (Sub 75%); 7) 50% of frames in a scan (Frm 50%) and 8)
75% of frames (Frm 75%), and tested on the same test set. Fig. 2 (b) plotted
$\epsilon_{acc.}$ with $M=20$ as an example of the reconstruction performance
using different training sets. The other models and metrics yielded a
consistent trend, which saw $\epsilon_{acc.}$ increased with both reduced
anatomical and protocol variances. It indicates both are factors impacting the
main reconstruction task performance.

Figure 2: Illustration of different US acquisition protocols and their imapct
to reconstruction accuracy. (a) Protocols of straight line shape (p1, p4), ‘C’
shape (p2, p5) and ‘S’ shape (p3, p6) with US probe perpendicular of and
parallel to the forearm, (b) $\epsilon_{acc.}$ changes due to reduced
anatomical and protocol variance.

###  3.3 Ablation Studies and Comparison

In our experiments, two types protocol discrimination tasks are considered,
six-class classification (as in Fig. 2) and three-class classification
(combining the perpendicular and parallel scans for the same scan shapes).
Each is combined with the 38-class classification task (38 training subjects)
as the anatomical discrimination task. Results from different $M={100,140}$
are also included, together with those from the main task network without any
branches (no-branch).

The optimised $z_{i,j}$ values versus the training epochs are plotted in Fig.
3 (a) and (b), where the optimum branch was selected by the maximum $z_{i,j}$
value, at the epoch. Table. 1 summarised the reconstruction performance of the
proposed method, using the optimised branches (their indices are denoted with
asterisks) for the two auxiliary tasks. Comparing with the no-branch models,
the improved performances from the proposed methods can be seen, for both
$M=100$ and $M=140$, regardless of the three- or six-class were used as
protocol discrimination tasks. For example, at $M=100$, $\epsilon_{drift}$ was
lowered from 14.52 to 6.56 mm, using the proposed methods with optimised
Branches 4 and 4, for protocol and anatomical discrimination tasks, with
$p$-value = 0.013 (unpaired t-test, at a significance level $\alpha=0.05$).
Statistical significance was found in performance improvement using all four
evaluation metrics.

Although no previous work utilise these discrimination tasks as privileged
information for assisting this application, two baseline models were
implemented as alternatives to the main task network. We have re-implemented
the proposed method using one of the first proposed approaches for this
application [16] as the main task network, with two adjacent frames as input
and outputting the transformation between them. In addition, the same network
architecture using more (10) input frames were also tested, denoted as
‘[16]-10’. These results are also summarised in Table. 1, with and without
optimised branches. However, we would like to emphasize that the reported
inferior results from these compared baselines need to be interpreted with
caution, as they were neither designed for nor tuned with incorporating these
auxiliary tasks used in this study. They are included for completeness and
reference values.

Table 1: Mean and standard deviation of four metrics of the proposed method and no-branch method. $M$ |  Num. of protocols |  Branch index (protocol/anatomy) | $\epsilon_{frame}$ | $\epsilon_{acc.}$ | $\epsilon_{dice}$ | $\epsilon_{drift}$
---|---|---|---|---|---|---
100 | n/a | No-branch | $0.18\pm 0.05$ | $7.03\pm 3.97$ | $0.84\pm 0.08$ | $14.52\pm 10.51$
100 | 3 |  Branches $5^{*}/4^{*}$ | $0.19\pm 0.06$ | $3.92\pm 3.50$ | $0.73\pm 0.21$ | $7.06\pm 7.30$
100 | 6 |  Branches $4^{*}/4^{*}$ | $0.17\pm 0.08$ | $3.80\pm 3.97$ | $0.76\pm 0.24$ | $6.56\pm 7.53$
140 | n/a | No-branch | $0.14\pm 0.05$ | $3.68\pm 3.10$ | $0.62\pm 0.28$ | $7.30\pm 7.40$
140 | 3 |  Branches $9^{*}/4^{*}$ | $0.15\pm 0.08$ | $3.36\pm 3.26$ | $\bm{0.94\pm 0.00}$ | $\bm{6.20\pm 6.31}$
140 | 6 |  Branches $5^{*}/7^{*}$ | $\bm{0.13\pm 0.05}$ | $\bm{2.90\pm 2.10}$ | $0.89\pm 0.00$ | $6.53\pm 5.98$
[16] | n/a | No-branch | $0.59\pm 0.28$ | $29.03\pm 9.15$ | $0.43\pm 0.32$ | $35.67\pm 11.20$
[16] | 3 |  Branches $1^{*}/9^{*}$ | $0.70\pm 0.50$ | $32.71\pm 18.10$ | $0.60\pm 0.22$ | $59.53\pm 36.87$
[16] | 6 |  Branches $4^{*}/9^{*}$ | $0.68\pm 0.46$ | $30.29\pm 17.44$ | $0.67\pm 0.16$ | $55.05\pm 32.69$
[16]-10 | n/a | No-branch | $0.38\pm 0.21$ | $17.60\pm 9.77$ | $0.67\pm 0.22$ | $22.64\pm 12.47$
[16]-10 | 3 |  Branches $4^{*}/4^{*}$ | $0.42\pm 0.30$ | $19.37\pm 10.63$ | $0.50\pm 0.30$ | $26.32\pm 13.33$
[16]-10 | 6 |  Branches $9^{*}/4^{*}$ | $0.43\pm 0.38$ | $21.72\pm 12.74$ | $0.50\pm 0.25$ | $29.64\pm 16.29$
Figure 3: The trend of task descriptor for anatomy (a) and protocol (b),
$M=100$, and the reconstruction performance (c). The epoch indicating best
performance of the model on validation set is denoted by a gray dotted line.
Scans with various scanning path are reconstructed using no-branch, optimised
branches, and random branches strategies, $M=100$.

##  4 Conclusion and Discussion

This work demonstrated the impact of anatomical and protocol variance towards
the 3D reconstruction of trackerless freehand US and formulated two respective
discrimination tasks for taking advantage these privileged information during
training. Using the proposed algorithm, substantially improved reconstruction
performance was achieved, which may indicate a promising new direction for
improving the potentials of this application for clinical adoption. Future
work includes testing clinical applications with specific challenges, such as
those without predefined protocol classes (where a clustering task may be used
instead), and comparison with approaches such as gradient surgery [22], which
may need adaptation for a single main task.

####  4.0.1 Declarations.

This work was supported by the EPSRC [EP/T029404/1], a Royal Academy of
Engineering / Medtronic Research Chair [RCSRF1819\7\734] (TV), Wellcome/EPSRC
Centre for Interventional and Surgical Sciences
[203145Z/16/Z], and the International Alliance for Cancer Early Detection, an
alliance between Cancer Research UK [C28070/A30912; C73666/A31378], Canary
Center at Stanford University, the University of Cambridge, OHSU Knight Cancer
Institute, University College London and the University of Manchester. TV is
co-founder and shareholder of Hypervision Surgical. Qi Li was supported by the
University College London Overseas and Graduate Research Scholarships. For the
purpose of open access, the authors have applied a CC BY public copyright
licence to any Author Accepted Manuscript version arising from this
submission. This study was performed in accordance with the ethical standards
in the 1964 Declaration of Helsinki and its later amendments or comparable
ethical standards. Approval was granted by the Ethics Committee of local
institution (UCL Department of Medical Physics and Biomedical Engineering) on
$20^{th}$ Jan. 2023 [24055/001].

## References

  * [1] Chang, R.F., Wu, W.J., Chen, D.R., Chen, W.M., Shu, W., Lee, J.H., Jeng, L.B.: 3-d us frame positioning using speckle decorrelation and image registration. Ultrasound in medicine & biology 29(6), 801–812 (2003)
  * [2] Guo, H., Chao, H., Xu, S., Wood, B.J., Wang, J., Yan, P.: Ultrasound volume reconstruction from freehand scans without tracking. IEEE Transactions on Biomedical Engineering 70(3), 970–979 (2022)
  * [3] Lang, A., Mousavi, P., Gill, S., Fichtinger, G., Abolmaesumi, P.: Multi-modal registration of speckle-tracked freehand 3d ultrasound to ct in the lumbar spine. Medical image analysis 16(3), 675–686 (2012)
  * [4] Li, Q., Shen, Z., Li, Q., Barratt, D.C., Dowrick, T., Clarkson, M.J., Vercauteren, T., Hu, Y.: Trackerless freehand ultrasound with sequence modelling and auxiliary transformation over past and future frames. arXiv preprint arXiv:2211.04867 (2022)
  * [5] Lindseth, F., Kaspersen, J.H., Ommedal, S., Langø, T., Bang, J., Hokland, J., Unsgaard, G., Nagelhus Hemes, T.A.: Multimodal image fusion in ultrasound-based neuronavigation: improving overview and interpretation by integrating preoperative mri with intraoperative 3d ultrasound. Computer Aided Surgery 8(2), 49–69 (2003)
  * [6] Liu, H., Simonyan, K., Yang, Y.: Darts: Differentiable architecture search. arXiv preprint arXiv:1806.09055 (2018)
  * [7] Luo, M., Yang, X., Huang, X., Huang, Y., Zou, Y., Hu, X., Ravikumar, N., Frangi, A.F., Ni, D.: Self context and shape prior for sensorless freehand 3d ultrasound reconstruction. In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2021: 24th International Conference, Strasbourg, France, September 27–October 1, 2021, Proceedings, Part VI 24. pp. 201–210. Springer (2021)
  * [8] Luo, M., Yang, X., Wang, H., Du, L., Ni, D.: Deep motion network for freehand 3d ultrasound reconstruction. In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2022: 25th International Conference, Singapore, September 18–22, 2022, Proceedings, Part IV. pp. 290–299. Springer (2022)
  * [9] Mikaeili, M., Bilge, H.Ş.: Trajectory estimation of ultrasound images based on convolutional neural network. Biomedical Signal Processing and Control 78, 103965 (2022)
  * [10] Miura, K., Ito, K., Aoki, T., Ohmiya, J., Kondo, S.: Localizing 2d ultrasound probe from ultrasound image sequences using deep learning for volume reconstruction. In: Medical Ultrasound, and Preterm, Perinatal and Paediatric Image Analysis: First International Workshop, ASMUS 2020, and 5th International Workshop, PIPPI 2020, Held in Conjunction with MICCAI 2020, Lima, Peru, October 4-8, 2020, Proceedings 1. pp. 97–105. Springer (2020)
  * [11] Miura, K., Ito, K., Aoki, T., Ohmiya, J., Kondo, S.: Pose estimation of 2d ultrasound probe from ultrasound image sequences using cnn and rnn. In: Simplifying Medical Ultrasound: Second International Workshop, ASMUS 2021, Held in Conjunction with MICCAI 2021, Strasbourg, France, September 27, 2021, Proceedings 2. pp. 96–105. Springer (2021)
  * [12] Mozaffari, M.H., Lee, W.S.: Freehand 3-d ultrasound imaging: a systematic review. Ultrasound in medicine & biology 43(10), 2099–2124 (2017)
  * [13] Newell, A., Jiang, L., Wang, C., Li, L.J., Deng, J.: Feature partitioning for efficient multi-task architectures. arXiv preprint arXiv:1908.04339 (2019)
  * [14] Ning, G., Liang, H., Zhou, L., Zhang, X., Liao, H.: Spatial position estimation method for 3d ultrasound reconstruction based on hybrid transfomers. In: 2022 IEEE 19th International Symposium on Biomedical Imaging (ISBI). pp. 1–5. IEEE (2022)
  * [15] Prevost, R., Salehi, M., Jagoda, S., Kumar, N., Sprung, J., Ladikos, A., Bauer, R., Zettinig, O., Wein, W.: 3d freehand ultrasound without external tracking using deep learning. Medical image analysis 48, 187–202 (2018)
  * [16] Prevost, R., Salehi, M., Sprung, J., Ladikos, A., Bauer, R., Wein, W.: Deep learning for sensorless 3d freehand ultrasound imaging. In: Medical Image Computing and Computer-Assisted Intervention- MICCAI 2017: 20th International Conference, Quebec City, QC, Canada, September 11-13, 2017, Proceedings, Part II. pp. 628–636. Springer (2017)
  * [17] Rajeswaran, A., Finn, C., Kakade, S.M., Levine, S.: Meta-learning with implicit gradients. Advances in neural information processing systems 32 (2019)
  * [18] Tan, M., Le, Q.: Efficientnet: Rethinking model scaling for convolutional neural networks. In: International conference on machine learning. pp. 6105–6114. PMLR (2019)
  * [19] Vapnik, V., Vashist, A.: A new learning paradigm: Learning using privileged information. Neural networks 22(5-6), 544–557 (2009)
  * [20] Wu, S., Zhang, H.R., Ré, C.: Understanding and improving information transfer in multi-task learning. arXiv preprint arXiv:2005.00944 (2020)
  * [21] Xie, Y., Liao, H., Zhang, D., Zhou, L., Chen, F.: Image-based 3d ultrasound reconstruction with optical flow via pyramid warping network. In: 2021 43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC). pp. 3539–3542. IEEE (2021)
  * [22] Yu, T., Kumar, S., Gupta, A., Levine, S., Hausman, K., Finn, C.: Gradient surgery for multi-task learning. Advances in Neural Information Processing Systems 33, 5824–5836 (2020)


# Harnessing the center-of-mass excitations in quantum metrology

Karol Gietka<EMAIL_ADDRESS>Quantum Systems Unit, Okinawa Institute of
Science and Technology Graduate University, Onna, Okinawa 904-0495, Japan

###### Abstract

In quantum metrology, one typically creates correlations among atoms or
photons to enhance measurement precision. Here, we show how one can use other
excitations to perform quantum-enhanced measurements on the example of center-
of-mass excitations of a spin-orbit coupled Bose-Einstein condensate and a
Coulomb crystal. We also present a method to simulate a homodyne detection of
center-of-mass excitations in these systems, which is required for optimal
estimation.

##  I Introduction

Quantum metrology relies on exploiting non-classical correlations between
system constituents to overcome the standard quantum limit of precision when
estimating an unknown parameter [1, 2]. In typical situations, these
correlations are being created among atoms [3] or photons [4, 5] which are
subsequently used to imprint the information about the unknown parameter.
Then, the final state of the system is measured, and the resulting measurement
outcomes are used to estimate the parameter. Unfortunately, in some cases,
obtaining a suitable correlated state of atoms or photons might be
challenging, hindering the possibility of performing quantum-enhanced
measurements. However, atoms and photons do not exhaust the set of possible
excitations among which the correlations could be generated and potentially
used [6, 7, 8, 9, 10].

Let us, for example, consider a single atom in a harmonic trap for which it
seems impossible to create quantum correlations for metrological purposes.
Nevertheless, if the harmonic trap has a very high frequency, the wavefunction
of the atom will be very localized, which means that the center-of-mass
excitations must be correlated (squeezed). Such squeezing is known to enhance
the precision of measurement in the context of photons [11]; unfortunately,
most of the forces that one would want to measure do not couple to the center-
of-mass motion (except gravity), but to the atomic internal degrees of
freedom, thus such correlations might seem indeed metrologically useless.

To overcome this apparent limitation, one could couple the center-of-mass
motion of an atom with its internal degree of freedom (for example, via spin-
orbit coupling) and potentially transfer the information from an atom whose
internal structure is affected by some force to the center-of-mass motion and
subsequently measure the position or momentum distribution to estimate the
unknown parameter. Then the precision of such estimation should scale with the
number of center-of-mass excitations which for a very squeezed state might be
very large even for a single atom. A similar approach is typically exploited
in light-matter systems. There, some force affects an ensemble of atoms, and
due to the light-matter interactions, the state of the light is also affected
by that force [12, 13, 14]. Performing then quadrature measurements on the
light usually allows one to estimate the strength of the force. In this case,
the precision scales typically with the number of field excitations (photons).

In this manuscript, on the example of center-of-mass excitations, we show how
one can harness less obvious excitations in quantum-enhanced measurements in
analogy to light-matter systems. To this end, we consider an adiabatic
protocol involving a spin-orbit coupled Bose-Einstein condensate and a
(sudden) quench protocol involving a Coulomb crystal. We also show how one can
simulate a homodyne detection of center-of-mass excitations in these systems
which is required for optimal estimation.

##  II Adiabatic protocol with a spin-orbit coupled Bose-Einstein condensate

A spin-orbit coupled Bose-Einstein condensate is a system in which atoms are
coupled to different momenta depending on their spin state. Although it is
well understood [15] and has been realized experimentally [16, 17], its
metrological potential remains unexplored. In the simplest possible case of a
two-mode non-interacting condensate, the behavior of each atom can be
described by the following single-particle Hamiltonian (we set $\hbar=1$
throughout the entire manuscript)

$\displaystyle\hat{H}=\frac{\left(\hat{p}+k\hat{\sigma}_{z}\right)^{2}}{2m}+\frac{m\omega^{2}\hat{x}^{2}}{2}+\frac{\Omega}{2}\hat{\sigma}_{x}$
(1)

where $\hat{x}$ and $\hat{p}\equiv-i\partial_{x}$ are position and momentum
operators of an atom, $\omega$ is the frequency of the harmonic trap, $m$ is
the atomic mass, $\Omega$ is the atomic frequency, $\hat{\sigma}_{i}$ are the
standard Pauli matrices, and $k$ is the spin-orbit interaction strength.
Depending on the strength of the coupling between the two modes of the
condensate $\Omega$, one can distinguish two phases. The single minimum phase
is characterized by $\Omega>2k^{2}/m$, and a double minimum phase (stripe
phase) is characterized by $\Omega<2k^{2}/m$.

For the sake of comparison with light-matter systems, we can now rotate
Hamiltonian (1) and rewrite it in the form of a quantum Rabi model using the
annihilation and creation operators of the harmonic oscillator

$\displaystyle\hat{H}=\omega\hat{a}^{\dagger}\hat{a}+\frac{\Omega}{2}\hat{\sigma}_{z}+k\sqrt{\frac{\omega}{2m}}(\hat{a}+\hat{a}^{\dagger})\hat{\sigma}_{x}.$
(2)

In other words, a non-interacting spin-orbit coupled Bose-Einstein condensate
of $N$ atoms trapped in a harmonic potential is equivalent to $N$
indistinguishable copies of a quantum Rabi model with critical coupling
strength $k_{c}=\sqrt{\Omega m/2}$. The latter has been extensively studied in
the context of light-matter interactions [18, 19] and (critical) quantum
metrology [20, 21, 22, 23].

Let us now discuss an equilibrium protocol which relies on adiabatic
preparation of the ground state of the system in the single minimum phase
[20]. For the sake of brevity, we assume now that $\Omega$ is an unknown
parameter, and we can control all other parameters to an arbitrary extent. In
order to find the ground state (in a suitable form), we can apply the
Schrieffer-Wolff transformation [24]

$\displaystyle\hat{U}=\exp\left(i\frac{k}{k_{c}}\frac{\sqrt{\omega}}{2\sqrt{\Omega}}(\hat{a}+\hat{a}^{\dagger}){\hat{\sigma}_{y}}\right),$
(3)

and for $1-k^{2}/k_{c}^{2}\gg({\omega}/{\Omega})^{2/3}$ rewrite the
Hamiltonian as (see Appendix A for details)

$\displaystyle\hat{H}\approx\omega\hat{a}^{\dagger}\hat{a}+\frac{\Omega}{2}\hat{\sigma}_{z}+k^{2}{\frac{\omega}{2m\Omega}}(\hat{a}+\hat{a}^{\dagger})^{2}\hat{\sigma}_{z}.$
(4)

whose ground state is a squeezed vacuum coupled to the spin down state
$|{\psi}\rangle=\hat{S}(\xi)|0\rangle\otimes|\\!\\!\downarrow\,\rangle$ with
$\hat{S}(\xi)\equiv\exp\\{(\xi/2)(\hat{a}^{\dagger})^{2}-(\xi^{*}/2)\hat{a}^{2}\\}$
being the squeeze operator with $\xi=-\frac{1}{4}\ln\\{1-(k/k_{c})^{2}\\}$. We
can now use this ground state to calculate the quantum Fisher information
$\mathcal{I}_{\Omega}=\langle\partial_{\Omega}\psi|\partial_{\Omega}\psi\rangle-\langle\psi|\partial_{\Omega}\psi\rangle^{2}$
which is related to the sensitivity of estimating an unknown parameter through
the Cramér-Rao bound $\Delta^{2}\Omega\geq 1/N\mathcal{I}_{\Omega}$ [25] with
$N$ being the number of system’s copies. The quantum Fisher information
becomes

$\displaystyle\mathcal{I}_{\Omega}=\frac{1}{8\Omega^{2}\left(1-\frac{k^{2}}{k_{c}^{2}}\right)^{2}}\frac{k^{4}}{k_{c}^{4}}.$
(5)

Since the ground state is a squeezed vacuum state, which is entirely
characterized by the width of the distribution (a consequence of Gaussianity),
the measurement that leads to such Fisher information is the measurement of
the second moment of position (or momentum) [26, 27]. The above expression
tells us, in fact, how distinguishable neighboring ground states (as a
function of $k$) are [28]. In order to make a metrological sense out of this
expression, we now assume that the coupling was adiabatically changed from
$k=0$ to $k=k_{f}<k_{c}$ which took time [20]
$T_{k_{f}}\approx\left(2\gamma\omega\sqrt{1-{k_{f}^{2}}/{k_{c}^{2}}}\right)^{-1}$
where $\gamma\ll 1$. Here we see the effect of critical slowing down;
increasing $k$ towards $k_{c}$ closes the energy gap which means that
preparation time also has to increase. Plugging this time into the expression
for the quantum Fisher information yields

$\displaystyle\mathcal{I}_{\Omega}=8\frac{\omega^{2}}{\Omega^{2}}\gamma^{2}\langle\hat{n}\rangle^{2}T_{k_{f}}^{2},$
(6)

where
$\langle\hat{n}\rangle=\sinh^{2}\xi\approx\left(4\sqrt{1-{k_{f}^{2}}/{k_{c}^{2}}}\right)^{-1}$
is the average number of center-of-mass excitations [29]. Thus, the Cramér-Rao
bound becomes

$\displaystyle\Delta^{2}\Omega\geq\left(8N\frac{\omega^{2}}{\Omega^{2}}\gamma^{2}\langle\hat{n}\rangle^{2}T_{k_{f}}^{2}\right)^{-1}$
(7)

and exhibits a Heisenberg scaling with respect to the center-of-mass
excitations (similar scalings can be shown for other parameters of the
system).

The above expression for the quantum Fisher information (and Cramér-Rao bound)
is valid as long as $1-k^{2}/k_{c}^{2}\gg({\omega}/{\Omega})^{2/3}$; in order
to check what happens once this condition is no longer satisfied, we have to
perform numerical calculations. These are presented in Fig. 1 where we
calculate the number of excitations, quantum Fisher information, and the
normalized energy gap $\Delta_{\epsilon}/\omega$ as a function of
$1-k^{2}/k_{c}^{2}$ for a fixed $\omega/\Omega$ using the Hamiltonian (1) (see
also Appendix A).

Figure 1: Comparison between the effective model (dashed black line) and
Hamiltonian (1) (solid red line) for $\Omega/\omega=10^{6}$. The vertical
dotted line indicates $1-k^{2}/k_{c}^{2}=({\omega}/{\Omega})^{2/3}$. (a), (b),
and (c) shows the number of excitations, quantum Fisher information, and
energy gap, respectively.

As can be seen in Fig. 1, once the condition
$1-k^{2}/k_{c}^{2}\gg({\omega}/{\Omega})^{2/3}$ is no longer satisfied, the
state does not change anymore and the number of excitations (Fig. 1a), quantum
Fisher information (Fig. 1b), and the energy gap (Fig. 1c) does not alter.

The above results prove that adiabatic preparation of the ground-state of a
spin-orbit coupled Bose-Einstein condensate can, in principle, lead to
sensitivities exhibiting Heisenberg scaling with the number of center-of-mass
excitations without exploiting any correlations between the atoms. In
practice, however, overcoming the standard quantum limit with an adiabatic
approach might be substantially limited due to the critical slowing down
effect. This is caused by the fact that during adiabatic dynamics, the
(squeezed) excitations are created very slowly—and thus ineffectively—as we
are closing the energy gap while increasing $k$ towards $k_{c}$. In order to
remove the problem of critical slowing down, one can perform a fast quench
beyond the transition point $k_{c}$ to increase the rate at which the
excitations are being created and thus make an advantage out of closing energy
gap. Although one can apply the sudden quench protocol to the spin-orbit
coupled Bose-Einstein condensate, in order to show the generality of using
center-of-mass excitations, we will now demonstrate this mechanism on the
example of a Coulomb crystal.

##  III Sudden quench protocol with a Coulomb crystal

A Coulomb crystal [30] is a crystal-like system composed of cold ions which
can serve as a quantum simulator. For example, in a recent experiment [31] a
Coulomb crystal in the presence of an external transverse field was used to
engineer a simulator of the paradigmatic Dicke Hamiltonian [32]. This was
achieved by coupling vibrational modes of the crystal with the spin degree of
freedom via a spin-dependent optical dipole force generated by the
interference of a pair of lasers. In the presence of an additional transverse
field, the system is described by the Dicke Hamiltonian

$\hat{H}=-\delta\hat{a}^{\dagger}\hat{a}+B\hat{S}_{x}-\frac{g}{\sqrt{N}}(\hat{a}+\hat{a}^{\dagger})\hat{S}_{z},$
(8)

where $\hat{a}$ ($\hat{a}^{\dagger}$) is the annihilation (creation) operator
of the center-of-mass excitation and $\hat{S}_{i}$ are the collective spin
operators. In the above Hamiltonian $\delta\equiv\omega_{R}-\omega$ is
detuning between the laser frequency $\omega_{R}$ and the center-of-mass mode
frequency $\omega$, $g$ is the homogeneous coupling between each ion and the
center-of-mass mode, and $B$ is the strength of the applied transverse field.
Depending on the transverse field strength, one can distinguish two phases.
The paramagnetic (normal) phase for $B>g^{2}/|\delta|\equiv B_{c}$ and the
ferromagnetic (superradiant) phase for $B<B_{c}$.

Dynamic metrological protocol relies on performing a sudden quench from the
single minimum (normal) phase to the double minimum (superradiant) phase and
[33] exploits the dynamics of an inverted harmonic oscillator [34]. This
provides an exponential speed up with respect to the adiabatic protocol [35].
In the following, we treat $\delta$ as the unknown parameter and change $B$ to
perform a quench between the two phases. In order to see the effect of an
inverted harmonic oscillator, we first apply the Schrieffer-Wolff
transformation

$\hat{U}=\exp\left(i\frac{B_{c}}{B}\frac{\delta}{g\sqrt{N}}(\hat{a}+\hat{a}^{\dagger}){\hat{S}_{y}}\right),$
(9)

to the Hamiltonian (8) and rewrite it with the help of quadrature operators
(see Appendix B)

$\hat{H}=\frac{\delta}{2}\hat{P}^{2}+\frac{\delta}{2}\left(1-\frac{B_{c}}{B}\right)\hat{X}^{2}.$
(10)

The above Hamiltonian describes a harmonic oscillator with a tunable frequency
which can become negative. It gives an accurate ground state description only
in the paramagnetic phase ($B>B_{c}$), as for $B<B_{c}$ the ground state does
not exist. Nevertheless, this Hamiltonian can be still used to describe
quenches from paramagnetic to the ferromagnetic phase until a certain number
of center-of-mass excitations is reached ($n_{\mathrm{max}}\propto
NB/\delta$), and we assume now that this number is not exceeded (see Appendix
B).

Figure 2: Comparison between the effective model (10) (dashed-black line) and
the Hamiltonian (8) (solid red line) for $B_{f}=B_{c}/2$. (a), (b), and (c)
shows $\langle\hat{X}^{2}\rangle$, number of excitations, and normalized
quantum Fisher information, respectively. For the numerical calculations, we
have set $N=1$, $B=B_{c}/2$, $g/\delta=80$, and $B/\delta=3200$.

In order to calculate the quantum Fisher information, we assume that the
system is quenched from the ground state at $B\gg B_{c}$, where the center-of-
mass excitations are decoupled from the spin, to some $B_{f}<B_{c}$.
Unfortunately, apart from dependence on hyperbolic functions (see Appendix C)
which reflects the fact that inverted harmonic oscillator dynamics stimulates
squeezing, the full expression for quantum Fisher information is complicated
and one does not clearly see its advantage. For a special case of
$B_{f}=B_{c}/2$, one can show that the quantum Fisher information is equal to

$\mathcal{I}_{\delta}=\frac{B_{c}^{2}}{2\delta^{2}B_{f}^{2}}\sinh^{4}\delta
t=\frac{2}{\delta^{2}}\langle\hat{n}\rangle^{2},$ (11)

where $\langle\hat{n}\rangle=\sinh^{2}\delta t$ is the instantaneous number of
center-of-mass excitations.

In order to inspect what happens once the inverted harmonic oscillator
dynamics is no longer valid, we have to resort to numerical calculations.
These are presented in Fig. 2, where we can see $\langle\hat{X}^{2}\rangle$,
$\langle\hat{n}\rangle^{2}$, and $\mathcal{I}_{\delta}$ as a function of time
after the quench. As can be seen in Fig. 2(b), once the number of excitations
does not grow exponentially, the same happens with the quantum Fisher
information [Fig. 2(c)] and squeezing [Fig. 2(a)].

In the case of a quench protocol, the final state is Gaussian, and the optimal
measurement that can lead to saturation of the Cramér-Rao bound is also a
homodyne detection but at an angle that is dependent on time and parameters of
the system [35]. Therefore, one needs a method to simulate a homodyne
detection at an arbitrary angle for the center-of-mass excitations.

##  IV Simulating homodyne detection

In the case of light-matter systems and Gaussian states of light, it is well
known that in order to extract the maximum amount of information from the
system one has to perform a homodyne measurement in the optimal direction
[36], which can be controlled with the phase of the local oscillator. For a
Bose-Einstein condensate and a Coulomb crystal, it seems that one is limited
to the quadrature measurement in only two directions, namely $\hat{x}$
($\hat{X}$) and $\hat{p}$ ($\hat{P}$). In the following, we will explain how
an equivalent of a generalized quadrature operator
$\hat{Q}(\theta)=\hat{x}\cos\theta+\hat{p}\sin\theta$ ($\theta$ is the local
oscillator phase) for the center-of-mass excitations could be measured in a
spin-orbit coupled Bose-Einstein condensate and a Coulomb crystal. To this
end, let us consider Hamiltonian (1) for $k=0$. Such a Hamiltonian describes
two subsystems that are decoupled. This means that for $k=0$, the dynamics of
center-of-mass excitations is governed solely by
$\hat{H}=\omega\hat{a}^{\dagger}\hat{a}$. Such a Hamiltonian generates
rotations in the quadrature phase space around its origin. Therefore, by
performing a sudden quench back to $k=0$ after the final state was reached, it
is possible to rotate the wavefunction by an arbitrary angle $\theta$ in the
phase space and perform thus a measurement of $\hat{Q}(\theta)$ by measuring
either momentum or position distribution.

In the case of a Coulomb crystal [Hamiltonian (8)], decoupling of the
collective spin from the center-of-mass motion can be done by turning off the
lasers or, alternatively, quenching the magnetic field back to $B\gg B_{c}$
after the final state is obtained. Therefore an arbitrary rotation can be
performed in the phase space of center-of-mass excitations by exploiting an
appropriate parameter quench. Then $\hat{Q}(\theta)$ could be measured by
inspecting either momentum or position distribution (see Appendix D for a
visualization) similarly as in the case of a spin-orbit-coupled Bose-Einstein
condensate.

##  V Conclusions

In this work, we have shown how to exploit excitations other than photonic or
atomic to enhance measurement precision. To this end, we have proposed how to
exploit center-of-mass excitations of a Bose-Einstein condensate and a Coulomb
crystal. We presented two protocols based on an adiabatic and a sudden quench,
originally developed for light-matter systems, that could be used to estimate
unknown parameters of the system with sensitivities exhibiting Heisenberg
scaling with respect to the center-of-mass excitations. We have also presented
how one can simulate the measurement of a generalized quadrature operator,
which is required for extracting the maximal amount of information about the
unknown parameter from these systems. This could be achieved by decoupling the
center-of-mass mode from the spin degree of freedom, waiting until the state
_aligns_ with the position or momentum axis, and finally measuring its
distribution. Then, as the resulting states are Gaussian, the second moment of
position and momentum can serve as the optimal estimator.

In the case of a spin-orbit coupled Bose-Einstien condensate, we did not
consider any interactions. In general, small spin-spin repulsive interactions
will shift the critical spin-orbit coupling strength to a larger value and
_change the shape_ of the wavefunction (it will no longer be Gaussian). This,
however, should not prevent using center-of-mass excitations in quantum
metrology. The major difference will be the necessity to use the entire
distribution of position (or momentum) to construct an optimal estimator.

In both cases, for finite temperature systems, where the system occupies
incoherently excited states, it can be shown that while the effective
Hamiltonian description is valid, the quantum Fisher information increases
with growing temperature (see Appendix E) which is caused by the increased
number of excitations in a thermal state. Once the effective harmonic
oscillator model is not valid anymore, the quantum Fisher information
decreases with increasing temperature. In this case, the entire distribution
of position (or momentum) is also required for optimal estimation as the state
is no longer Gaussian.

The presented protocols could be used in any system that can be mapped to a
harmonic oscillator coupled to a (collective) spin; one intriguing example is
electrons trapped on the surface of liquid Helium, where the role of the
harmonic oscillator is played by Landau levels [37]. We hope that this work
will stimulate the search for new metrological protocols (in known and well-
understood systems) involving various kinds of excitations that might be more
robust to decoherence and losses than atoms or photons. In particular,
protocols involving spin-orbit coupled fermionic atoms [38] (or a Tonks-
Girardeau gas [39, 40, 41]) might lead to protocols robust to atom losses and
benefiting from both correlated center-of-mass excitations and naturally
existing fermionic correlations [42].

##  VI Acknowledgements

K.G. is pleased to acknowledge Mohamed Hatifi, Thomás Fogarty, Ayaka Usui,
Yongping Zhang, and Thomas Busch for fruitful discussions. Simulations were
performed using the open-source QuantumOptics.jl framework in Julia [43]. This
work was supported by the Okinawa Institute of Science and Technology Graduate
University.

##  Appendix A Effective Hamiltonian for the spin-orbit coupled Bose-Einstein
condensate

Here we derive the effective Hamiltonian (4). The Schrieffer-Wolff
transformation

$\hat{U}=\exp\left(i\frac{k}{k_{c}}\frac{\sqrt{\omega}}{2\sqrt{\Omega}}(\hat{a}+\hat{a}^{\dagger}){\hat{\sigma}_{y}}\right),$
(12)

is a simultaneous rotation of the spin by an angle
$\phi=\frac{k}{k_{c}}\frac{\sqrt{\omega}}{\sqrt{\Omega}}\left(\hat{a}^{\dagger}+\hat{a}\right)$
and displacement of the harmonic oscillator by
$\alpha=i\frac{k}{k_{c}}\frac{\sqrt{\omega}}{2\sqrt{\Omega}}\hat{\sigma}_{y}$.
Applying this transformation to the Hamiltonian (2) yields

$\displaystyle\hat{H}=\omega\left(\hat{a}^{\dagger}-\alpha\right)\left(\hat{a}+\alpha\right)+\frac{\Omega}{2}\left(\cos\phi\hat{\sigma}_{z}-\sin\phi\hat{\sigma}_{x}\right)+k\sqrt{\frac{\omega}{2m}}\left(\hat{a}+\hat{a}^{\dagger}\right)\left(\cos\phi\hat{\sigma}_{x}+\sin\phi\hat{\sigma}_{z}\right).$
(13)

Assuming $\omega/\Omega\rightarrow 0$, we obtain Hamiltonian (4)

$\displaystyle\hat{H}\approx\omega\hat{a}^{\dagger}\hat{a}+\frac{\Omega}{2}\hat{\sigma}_{z}+k^{2}{\frac{\omega}{2m\Omega}}(\hat{a}+\hat{a}^{\dagger})^{2}\hat{\sigma}_{z}.$
(14)

Figure 3: Comparison between the Hamiltonian (2) and (4). For
$(\omega/\Omega)^{2/3}\ll 1-k^{2}/k_{c}^{2}$, the two Hamiltonians have the
same ground state (the difference in negligible). Once this condition is not
satisfied, the ground state of Hamiltonian (2) does not change anymore by
further increasing $k$ towards $k_{c}$. (a) and (b) shows logarithm of
$2m\omega\langle\hat{x}^{2}\rangle$ and $m\omega/2\langle\hat{p}^{2}\rangle$.
(c) and (d) shows second moments of $\hat{x}$ and $\hat{p}$ for Hamiltonian
(2) normalized to second moments of $\hat{x}$ and $\hat{p}$ for Hamiltonian
(4).

By projecting the above Hamiltonian onto the spin down sector we finally get

$\displaystyle\hat{H}=\omega\hat{a}^{\dagger}\hat{a}-k^{2}{\frac{\omega}{2m\Omega}}(\hat{a}+\hat{a}^{\dagger})^{2},$
(15)

which can be rewritten using the position and momentum operators in the form
of a harmonic oscillator with a tunable frequency.

$\hat{H}=\frac{p^{2}}{2m}+\frac{m\omega^{2}}{2}\left(1-\frac{k^{2}}{k^{2}_{c}}\right)\hat{x}^{2}.$
(16)

We now expect that there exists a weaker condition than
$\omega/\Omega\rightarrow 0$ that relates $\omega/\Omega$ and $k/k_{c}$ such
that the effective description is still valid. In order to find this condition
and check the limits of applicability of the effective Hamiltonian we now
compute its ground state as a function of $1-k^{2}/k_{c}^{2}$ and
$\Omega/\omega$ and compare with the ground state of the original Hamiltonian.
We choose to compare the second moments of $\hat{x}$ and $\hat{p}$ as they
characterize a Gaussian state centered at the origin of the phase space. The
results of the simulations are presented in Fig. 3 where we also indicate the
$1-k^{2}/k_{c}^{2}=({\omega}/{\Omega})^{2/3}$ condition with a dashed-black
line. Similar calculations can be performed for the Coloumb crystal.

##  Appendix B Effective Hamiltonian for the Coulomb crystal

Here we derive the effective Hamiltonian (10). Similarly as in the case of a
spin-orbit-coupled Bose-Einstein condensate, the Schrieffer-Wolff
transformation

$\hat{U}=\exp\left(i\frac{B_{c}}{B}\frac{\delta}{g\sqrt{N}}(\hat{a}+\hat{a}^{\dagger}){\hat{S}_{y}}\right)$
(17)

is a rotation of the collective spin around the $y$ axis by an angle
$\varphi=\frac{B_{c}}{B}\frac{\delta}{g\sqrt{N}}(\hat{a}+\hat{a}^{\dagger})$
and displacement of the harmonic oscillator by
$\beta=i\frac{B_{c}}{B}\frac{\delta}{g\sqrt{N}}{\hat{S}_{y}}$. Applying the
transformation to Hamiltonian (8) yields

$\hat{H}=-\delta\left(\hat{a}^{\dagger}-\beta\right)\left(\hat{a}+\beta\right)+B\left(\cos\varphi\hat{S}_{x}+\sin\varphi\hat{S}_{z}\right)-\frac{g}{\sqrt{N}}\left(\hat{a}+\hat{a}^{\dagger}\right)\left(\cos\varphi\hat{S}_{z}-\sin\varphi\hat{S}_{x}\right).$
(18)

Assuming $\delta/B\rightarrow 0$ we get

$\hat{H}=-\delta\hat{a}^{\dagger}\hat{a}+B\hat{S}_{x}-\frac{g^{2}}{2NB}(\hat{a}+\hat{a}^{\dagger})^{2}\hat{S}_{x},$
(19)

which upon projecting on the low energy sector and using quadrature operators
gives

$\hat{H}=\frac{\delta}{2}\hat{P}^{2}+\frac{\delta}{2}\left(1-\frac{B_{c}}{B}\right)\hat{X}^{2},$
(20)

and describes a harmonic oscillator with a tunable frequency which can become
negative. By quenching from paramagnetic to ferromagnetic phase this
Hamiltonian will give a proper description until the number of excitations
approaches the number of excitations characterizing the ground state of the
quenched Hamiltonian which can be calculated to be (assuming $B<B_{c}$)

$\langle\hat{a}^{\dagger}\hat{a}\rangle_{\mathrm{GS}}=\frac{NB}{4\delta}\left(\frac{B_{c}}{B}-\frac{B}{B_{c}}\right),$
(21)

where the subscript GS indicates the ground state. Therefore, for a fixed $B$
the time for which the effective description is valid can be extended by
increasing $\frac{NB}{\delta}$. This is illustrated in Fig. 4 where we plot
the number of excitations as a function of time for different parameters.

Figure 4: The dashed black line is the inverted harmonic oscillator dynamics.
In (a) the number of ions is fixed to $N=1$; yellow, blue, and red line is a
simulation for $NB/\delta=0.1875$, $NB/\delta=18.75$, and $NB/\delta=1875$,
respectively. In (b) the $B/\delta$ ratio is fixed to $B/\delta=18.75$;
yellow, blue, and red line is a simulation for $NB/\delta=18.75$ ($N=1$),
$NB/\delta=187.5$ ($N=10$), and $NB/\delta=1875$ ($N=100$), respectively.

##  Appendix C Quantum Fisher information for a sudden quench in a Coulomb
crystal

Here we calculate the quantum Fisher information for a sudden quench from the
paramagnetic phase to the ferromagnetic phase. The quantum Fisher information
can be calculated as a variance of the local generator $\hat{h}_{\lambda}$,
where $\lambda$ is the unknown parameter. The local generator
$\hat{h}_{\lambda}$ for a Hamiltonian
$\hat{H}=\hat{H}_{0}+\lambda\hat{H}_{\lambda}$ can be expressed as

$\begin{split}\hat{h}_{\lambda}=ie^{i\hat{H}t}\left(\partial_{\lambda}e^{-i\hat{H}t}\right)=\int_{0}^{t}e^{i\hat{H}s}\hat{H}_{\lambda}e^{-i\hat{H}s}\mathrm{d}s=\int_{0}^{t}\sum_{n=0}^{\infty}\frac{(is)^{n}}{n!}\left[\hat{H},\hat{H}_{\lambda}\right]_{n}\mathrm{d}s=-i\sum_{n=0}^{\infty}\frac{(it)^{n+1}}{(n+1)!}\left[\hat{H},\hat{H}_{\lambda}\right]_{n},\end{split}$
(22)

where
$[\hat{H},\hat{H}_{\lambda}]_{n+1}=[\hat{H},[\hat{H},\hat{H}_{\lambda}]]_{n}$
and $[\hat{H},\hat{H}_{\lambda}]_{0}=\hat{H}_{\lambda}$. Now it can be shown
that if $[\hat{H},\hat{H}_{\lambda}]=i\hat{C}$, $[\hat{H},i\hat{C}]=-\hat{D}$,
and $[\hat{H},-\hat{D}]=i\epsilon\hat{C}$, then

$[\hat{H},\hat{H}_{\lambda}]_{2n+1}=i\epsilon^{n}\hat{C}\quad\mathrm{and}\quad[\hat{H},\hat{H}_{\lambda}]_{2n+2}=-\epsilon^{n}\hat{D}.$
(23)

Therefore, one can express the local generator $\hat{h}_{\lambda}$ as

$\hat{h}_{\lambda}=\hat{H}_{\lambda}t-i\sum_{n=0}^{\infty}\frac{(it)^{2n+2}}{(2n+2)!}i\epsilon^{n}\hat{C}+i\sum_{n=0}^{\infty}\frac{(it)^{2n+3}}{(2n+3)!}\epsilon^{n}\hat{D}=\hat{H}_{\lambda}t+\frac{\cos(\sqrt{\epsilon}t)-1}{\epsilon}\hat{C}-\frac{\sin(\sqrt{\epsilon}t)-\sqrt{\epsilon}t}{\sqrt{\epsilon}\epsilon}\hat{D}.$
(24)

Taking the effective Hamiltonian for the Coulomb crystal

$\hat{H}=\frac{\delta}{2}\hat{P}^{2}+\frac{\delta}{2}\left(1-\frac{B_{c}}{B}\right)\hat{X}^{2},$
(25)

assuming that $\delta$ is unknown, choosing
$\hat{H}_{\delta}=\frac{1}{2}(\hat{X}^{2}+\hat{P}^{2})$ and
$\hat{H}_{0}=-\frac{g^{2}}{2B}\hat{X}^{2}$, one can show that

$\begin{split}\hat{h}_{\delta}=&\frac{t}{2}\left(\hat{X}^{2}+\hat{P}^{2}\right)-\frac{g^{2}}{2B}\frac{\cos\left(\sqrt{\epsilon}t\right)-1}{\epsilon}\left(\hat{X}\hat{P}+\hat{P}\hat{X}\right)+\frac{g^{2}\delta}{B}\frac{\sin\left(\sqrt{\epsilon}t\right)-\sqrt{\epsilon}t}{\epsilon^{3/2}}\left(\hat{X}^{2}-\hat{P}^{2}\right)-\frac{g^{4}}{B^{2}}\frac{\sin\left(\sqrt{\epsilon}t\right)-\sqrt{\epsilon}t}{\epsilon^{3/2}}\hat{X}^{2},\end{split}$
(26)

where $\epsilon=4\delta^{2}({1-B_{c}/B})$. Once $\epsilon$ becomes imaginary
($B<B_{c}$), the sine and cosine functions become hyperbolic which reflects
the exponential generation of excitations following a quench and leads to
exponential dependence of the quantum Fisher information
$\mathcal{I}_{\delta}=4\Delta^{2}\hat{h}_{\delta}$ on time. A similar
calculation can be performed treating other parameters as unknown as well as
for the spin-orbit-coupled Bose-Einstein condensate.

##  Appendix D simulating homodyne detection

Here we present simulations explaining how one can simulate homodyne detection
of center-of-mass excitations. In order to simulate the homodyne detection, we
decouple the center-of-mass mode from the spin degree of freedom after the
final state is obtained. Waiting subsequently for a proper amount of time
makes it possible to rotate the state in the phase space and align it with the
real (position) or imaginary (momentum) axis. Therefore by measuring the
distribution of position or momentum, one can measure an arbitrary quadrature.
This mechanism is presented in Fig. 3, where we plot the Husimi representation
$Q(\alpha)=\frac{1}{\pi}\langle\alpha|\hat{\rho}_{\psi}(t)|\alpha\rangle$
(where $|\alpha\rangle$ is a coherent state with amplitude $\alpha$ and
$\hat{\rho}_{\psi}(t)\equiv|\psi(t)\rangle\langle\psi(t)|$ is the state of the
system) for the (a) initial ($t=0$) and (b) final ($t=T$) state, as well as
the rotated final state aligned with the (c) $\hat{x}$ and (d) $\hat{p}$
quadratures.

Figure 5: Simulating homodyne detection of center-of-mass excitations. (a) is
the initial state (vacuum), (b) is the final state (squeezed state), (c) is
the final state rotated to align with the $\hat{x}$ quadrature, (d) is the
final state rotated to align with the $\hat{y}$ quadrature.

##  Appendix E Finite temperature calculations

In the main text, we have focused on the zero temperature case where the
system resides in its ground state and is described by a pure state. From a
practical point of view, however, it is relevant to consider a finite
temperature case where the system occupies (incoherently) excited levels as
well. In such a case, the state of the system can be described using a density
matrix operator (statistical mixture)

$\hat{\rho}=\sum_{n=0}^{\infty}p_{n}|n\rangle\langle n|,$ (27)

where $n$ labels the excited states and $p_{n}$ is given by

$p_{n}=\frac{1}{Z}\exp(-\beta E_{n}),$ (28)

where $Z=\sum_{n}\exp(-\beta E_{n})$ is the partition function (statistical
sum) with $E_{n}$ being the energy of the $n$th excited state, and
$\beta=1/k_{b}T$ where $k_{b}$ is the Boltzmann constant and $T$ is the
temperature. In order to calculate the quantum Fisher information for a
statistical mixture, we have to use the following formula

$\mathcal{I_{\lambda}}=2\sum_{n,m}\frac{(p_{n}-p_{m})^{2}}{p_{n}+p_{m}}|\langle
n|\hat{h}_{\lambda}|m\rangle|^{2}.$ (29)

The results of numerical simulations for the spin-orbit coupled Bose-Einstein
condensate are presented in Fig 6. In panel (a), we plot the quantum Fisher
information for the limiting case where $\omega/\Omega\rightarrow 0$ (solid
lines) and for finite $\omega/\Omega$ (dashed lines). Interestingly, the
quantum Fisher information increases with growing temperature until the
maximum amount of squeezing is obtained. Once the plateau is reached,
increasing the temperature decreases the quantum Fisher information. This is
further illustrated in panel (b), where we plot the quantum Fisher information
as a function of temperature for the case where the maximum amount of
squeezing has been reached ($1-k^{2}/k_{c}^{2}=10^{-5}$). The (initial) growth
of the quantum Fisher information is caused by the increased number of center-
of-mass excitations for a thermal state.

Figure 6: Quantum Fisher information for various temperatures of a spin-orbit
coupled Bose-Einstein condensate. In panel (a) the solid lines represent
$\Omega/\omega\rightarrow\infty$ case and the dashed lines represent
$\Omega/\omega=10^{4}$ case. Black, blue, and red colors correspond to
$\omega\beta=2,0.3,0.1$, respectively. In panel (b), we plot the final value
of the quantum Fisher information as a function of $\omega\beta$. Once
$\omega\lesssim 2$ the plateau value of the quantum Fisher information starts
to decrease.

The results of numerical simulations for the Coulomb crystal are presented in
Fig 7. Similarly as in Fig. 6, in panel (a) we plot the quantum Fisher
information for the limiting case where $NB/\delta\rightarrow 0$ (solid lines)
and for finite $NB/\delta$ (dashed lines). The quantum Fisher information
increases (again) with growing temperature until the inverted harmonic
oscillator picture is no longer applicable. Once quantum Fisher information
does not grow exponentially anymore, the effect of increased temperature
decreases the quantum Fisher information. The growth of the quantum Fisher
information is further illustrated in panel (b), where we plot the quantum
Fisher information as a function of temperature for the case when the inverted
harmonic oscillator picture is still valid (without the loss of generality, we
choose $\delta t=4$ point).

Figure 7: Quantum Fisher information for various temperatures of a Coulomb
crystal. In panel (a) the solid lines represent $NB/\delta\rightarrow\infty$
case and the dashed lines represent $NB/\delta=200$ case. Black, blue, and red
colors correspond to $\delta\beta=2,1,0.5$, respectively. In panel (b), we
plot the value of the quantum Fisher information for the infinite case
($NB/\delta\rightarrow\infty$) as a function of $\delta\beta$ for $\delta
t=4$. Once $\delta\beta\lesssim 4$, the final value of the quantum Fisher
information starts to increase.

## References

  * Giovannetti  _et al._ [2006] V. Giovannetti, S. Lloyd, and L. Maccone, Quantum metrology, Phys. Rev. Lett. 96, 010401 (2006).
  * Giovannetti  _et al._ [2011] V. Giovannetti, S. Lloyd, and L. Maccone, Advances in quantum metrology, Nat. Photonics 5, 222 (2011).
  * 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).
  * Demkowicz-Dobrzański  _et al._ [2015] R. Demkowicz-Dobrzański, M. Jarzyna, and J. Kołodyński, Chapter four - quantum limits in optical interferometry (Elsevier, 2015) pp. 345–435.
  * Polino  _et al._ [2020] E. Polino, M. Valeri, N. Spagnolo, and F. Sciarrino, Photonic quantum metrology, AVS Quantum Science 2, 024703 (2020).
  * Pirkkalainen  _et al._ [2015] J.-M. Pirkkalainen, E. Damskägg, M. Brandt, F. Massel, and M. A. Sillanpää, Squeezing of quantum noise of motion in a micromechanical resonator, Phys. Rev. Lett. 115, 243601 (2015).
  * Zhang  _et al._ [2018] J. Zhang, M. Um, D. Lv, J.-N. Zhang, L.-M. Duan, and K. Kim, Noon states of nine quantized vibrations in two radial modes of a trapped ion, Phys. Rev. Lett. 121, 160502 (2018).
  * Xin  _et al._ [2021] M. Xin, W. S. Leong, Z. Chen, Y. Wang, and S.-Y. Lan, Rapid quantum squeezing by jumping the harmonic oscillator frequency, Phys. Rev. Lett. 127, 183602 (2021).
  * Kustura  _et al._ [2022] K. Kustura, C. Gonzalez-Ballestero, A. d. l. R. Sommer, N. Meyer, R. Quidant, and O. Romero-Isart, Mechanical squeezing via unstable dynamics in a microcavity, Phys. Rev. Lett. 128, 143601 (2022).
  * Wollack  _et al._ [2022] E. A. Wollack, A. Y. Cleland, R. G. Gruenke, Z. Wang, P. Arrangoiz-Arriola, and A. H. Safavi-Naeini, Quantum state preparation and tomography of entangled mechanical resonators, Nature 604, 463 (2022).
  * Lawrie  _et al._ [2019] B. J. Lawrie, P. D. Lett, A. M. Marino, and R. C. Pooser, Quantum sensing with squeezed light, ACS Photonics 6, 1307 (2019).
  * Gietka  _et al._ [2019] K. Gietka, F. Mivehvar, and H. Ritsch, Supersolid-based gravimeter in a ring cavity, Phys. Rev. Lett. 122, 190801 (2019).
  * Gietka  _et al._ [2021] K. Gietka, F. Mivehvar, and T. Busch, Cavity-enhanced magnetometer with a spinor bose–einstein condensate, New J. Phys. 23, 043020 (2021).
  * Amorós-Binefa and Kołodyński [2021] J. Amorós-Binefa and J. Kołodyński, Noisy atomic magnetometry in real time, New J. Phys. 23, 123030 (2021).
  * Zhang  _et al._ [2016] Y. Zhang, M. E. Mossman, T. Busch, P. Engels, and C. Zhang, Properties of spin–orbit-coupled bose–einstein condensates, Front. Phys. 11, 118103 (2016).
  * Lin  _et al._ [2011] Y. J. Lin, K. Jiménez-García, and I. B. Spielman, Spin–orbit-coupled bose–einstein condensates, Nature 471, 83 (2011).
  * Li  _et al._ [2016] J. Li, W. Huang, B. Shteynas, S. Burchesky, F. i. m. c. b. u. i. e. i. f. Top, E. Su, J. Lee, A. O. Jamison, and W. Ketterle, Spin-orbit coupling and spin textures in optical superlattices, Phys. Rev. Lett. 117, 185301 (2016).
  * Crespi  _et al._ [2012] A. Crespi, S. Longhi, and R. Osellame, Photonic realization of the quantum rabi model, Phys. Rev. Lett. 108, 163601 (2012).
  * Braak  _et al._ [2016] D. Braak, Q.-H. Chen, M. T. Batchelor, and E. Solano, Semi-classical and quantum rabi models: in celebration of 80 years, J. Phys. A: Math. Theor. 49, 300301 (2016).
  * Garbe  _et al._ [2020] L. Garbe, M. Bina, A. Keller, M. G. A. Paris, and S. Felicetti, Critical quantum metrology with a finite-component quantum phase transition, Phys. Rev. Lett. 124, 120504 (2020).
  * Chu  _et al._ [2021] Y. Chu, S. Zhang, B. Yu, and J. Cai, Dynamic framework for criticality-enhanced quantum sensing, Phys. Rev. Lett. 126, 010502 (2021).
  * Ilias  _et al._ [2022] T. Ilias, D. Yang, S. F. Huelga, and M. B. Plenio, Criticality-enhanced quantum sensing via continuous measurement, PRX Quantum 3, 010354 (2022).
  * Gietka [2022] K. Gietka, Squeezing by critical speeding up: Applications in quantum metrology, Phys. Rev. A 105, 042620 (2022).
  * Schrieffer and Wolff [1966] J. R. Schrieffer and P. A. Wolff, Relation between the anderson and kondo hamiltonians, Phys. Rev. 149, 491 (1966).
  * Pezzé and Smerzi [2009] L. Pezzé and A. Smerzi, Entanglement, nonlinear dynamics, and the heisenberg limit, Phys. Rev. Lett. 102, 100401 (2009).
  * Lvovsky and Raymer [2009] A. I. Lvovsky and M. G. Raymer, Continuous-variable optical quantum-state tomography, Rev. Mod. Phys. 81, 299 (2009).
  * Weedbrook  _et al._ [2012] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Gaussian quantum information, Rev. Mod. Phys. 84, 621 (2012).
  * Braunstein and Caves [1994] S. L. Braunstein and C. M. Caves, Statistical distance and the geometry of quantum states, Phys. Rev. Lett. 72, 3439 (1994).
  * Walls [1983] D. F. Walls, Squeezed states of light, Nature 306, 141 (1983).
  * Thompson [2015] R. C. Thompson, Ion coulomb crystals, Contemp. Phys. 56, 63 (2015).
  * Safavi-Naini  _et al._ [2018] A. Safavi-Naini, R. J. Lewis-Swan, J. G. Bohnet, M. Gärttner, K. A. Gilmore, J. E. Jordan, J. Cohn, J. K. Freericks, A. M. Rey, and J. J. Bollinger, Verification of a many-ion simulator of the dicke model through slow quenches across a phase transition, Phys. Rev. Lett. 121, 040503 (2018).
  * Garraway [2011] B. M. Garraway, The dicke model in quantum optics: Dicke model revisited, Philos. Trans. Royal Soc. A 369, 1137 (2011).
  * Hayes  _et al._ [2018] A. J. Hayes, S. Dooley, W. J. Munro, K. Nemoto, and J. Dunningham, Making the most of time in quantum metrology: concurrent state preparation and sensing, Quantum Sci. Technol. 3, 035007 (2018).
  * Gietka and Busch [2021] K. Gietka and T. Busch, Inverted harmonic oscillator dynamics of the nonequilibrium phase transition in the dicke model, Phys. Rev. E 104, 034132 (2021).
  * Gietka  _et al._ [2022] K. Gietka, L. Ruks, and T. Busch, Understanding and Improving Critical Metrology. Quenching Superradiant Light-Matter Systems Beyond the Critical Point, Quantum 6, 700 (2022).
  * Monras [2006] A. Monras, Optimal phase measurements with pure gaussian states, Phys. Rev. A 73, 033821 (2006).
  * Yunusova  _et al._ [2019] K. M. Yunusova, D. Konstantinov, H. Bouchiat, and A. D. Chepelianskii, Coupling between rydberg states and landau levels of electrons trapped on liquid helium, Phys. Rev. Lett. 122, 176802 (2019).
  * [38] K. Gietka, (unpublished).
  * Tonks [1936] L. Tonks, The complete equation of state of one, two and three-dimensional gases of hard elastic spheres, Phys. Rev. 50, 955 (1936).
  * Lieb and Liniger [1963] E. H. Lieb and W. Liniger, Exact analysis of an interacting bose gas. i. the general solution and the ground state, Phys. Rev. 130, 1605 (1963).
  * Kinoshita  _et al._ [2004] T. Kinoshita, T. Wenger, and D. S. Weiss, Observation of a one-dimensional tonks-girardeau gas, Science 305, 1125 (2004).
  * Benatti  _et al._ [2014] F. Benatti, R. Floreanini, and U. Marzolino, Entanglement in fermion systems and quantum metrology, Phys. Rev. A 89, 032326 (2014).
  * Krämer  _et al._ [2018] S. Krämer, D. Plankensteiner, L. Ostermann, and H. Ritsch, QuantumOptics.jl: A Julia framework for simulating open quantum systems, Comput. Phys. Commun 227, 109 (2018).