# Fast equilibrium reconstruction by deep learning on EAST tokamak

Jingjing Lu  University of Science and Technology of China, Hefei, 230026,
China  Institute of Plasma Physics, Chinese Academy of Sciences, Hefei 230031,
China  Youjun Hu<EMAIL_ADDRESS>Institute of Plasma Physics, Chinese Academy
of Sciences, Hefei 230031, China  Nong Xiang<EMAIL_ADDRESS>Institute of
Plasma Physics, Chinese Academy of Sciences, Hefei 230031, China  Youwen Sun
Institute of Plasma Physics, Chinese Academy of Sciences, Hefei 230031, China

###### Abstract

A deep neural network is developed and trained on magnetic measurements
(input) and EFIT poloidal magnetic flux (output) on the EAST tokamak. In
optimizing the network architecture, we use automatic optimization in
searching for the best hyperparameters, which helps the model generalize
better. We compare the inner magnetic surfaces and last-closed-flux surfaces
(LCFSs) with those from EFIT. We also calculated the normalized internal
inductance, which is completely determined by the poloidal magnetic flux and
can further reflect the accuracy of the prediction. The time evolution of the
internal inductance in full discharges is compared with that provided by EFIT.
All of the comparisons show good agreement, demonstrating the accuracy of the
machine learning model, which has the high spatial resolution as the off-line
EFIT while still meets the time constraint of real-time control.

##  I Introduction

Reconstructing magnetic configuration using magnetic measurements is a routine
task of tokamak operation. There are many equilibrium solvers, e.g.,
EFITLao_1985 (1, 2, 3, 4, 5, 6), that can do this kind of reconstruction by
solving the Grad-Shafranov equation under the constraint of magnetic
measurements. In recent years, accumulation of data resulting from these
reconstruction practices, along with the development of machine learning
algorithms, software frameworks and computing power, have made it possible to
train deep neural networks to provide reconstructions as accurate as those by
EFIT. This has been demonstrated on KSTARKSTAR_2020 (7) and DIII-DDIII-D_2022
(8).

On the EAST tokamakWan_2017 (9), EFIT has been routinely used in tokamak
operations for more than ten years and substantial equilibrium data have been
accumulatedQian2009 (10, 11, 12, 13). In this paper, we report the results of
magnetic reconstruction by a deep neural network trained on the magnetic
measurements and EFIT reconstructed 2D magnetic poloidal flux.

There are two versions of EFIT used on EAST, one is for real-time control and
one for off-line analysis. The former is often restricted to lower accuracy
due to the time constraint of real time control, while the latter is of higher
accuracy. In this work, we use the off-line EFIT data in training the neural
network. The model trained this way have the higher accuracy as the off-line
EFIT while still meets the time constraint of the real-time control.

There are many hyperparameters in a neural network that usually need to be set
manually, such as number of hidden layers, units per layer, mini-batch size,
learning rate, number of epochs of training. In recent years, there appear
optimization libraries that can automatically set the values of these
hyperparameters. In this work, we use the Optuna optimization
frameworkoptuna_2019 (14) in setting the hyperparameters. The hyperparameters
found this way turn out to be much better than our previously manually set
ones in terms of the model accuracy. The size of the network architecture
found by the automatic hyperparameter tuning turns out to be relatively small
(with less than 2 million parameters). This small size allows for very fast
equilibrium construction that can be easily deployed.

The input to the network is limited to only the magnetic measurements. The
output of the network is the 2D poloidal magnetic flux function,
$\Psi(R,Z)\equiv A_{\phi}R$, which is related to the poloidal magnetic field,
$B_{R}$ and $B_{Z}$, by

$B_{R}=-\frac{1}{R}\frac{\partial\Psi}{\partial Z},$ (1)
$B_{Z}=\frac{1}{R}\frac{\partial\Psi}{\partial R},$ (2)

where $(R,\phi,Z)$ are the cylindrical coordinates. The 2D contours of $\Psi$
in $(R,Z)$ plane correspond to the magnetic surfaces. We compare the inner
magnetic surfaces and the LCFSs with those given by EFIT, in order to evaluate
the accuracy of $\Psi$ predicted by the network. We also calculate the
normalized internal inductance, $l_{i}$, which is a quantity that is solely
determined by $\Psi$ and thus can reflect how accuracy the predicted $\Psi$
is. The time evolution of the internal inductance in full discharges is
compared with that provided by EFIT. All of the comparisons show good
agreement, demonstrating the accuracy of the machine learning model, which has
the high spatial resolution as the off-line EFIT while still meets the time
constraint of real-time control.

The rest of this paper is organized as follows. Section II presents how the
data are collected and normalized. Sec. III explains the structure of our
neural network and how the hyperparameters are chosen by automatic
optimization. In section IV, we test the predicting capability of the trained
network. Section V discusses a small network used to predict some volume-
averaged quantities, namely the plasma stored energy $W_{\operatorname{mhd}}$,
normalized toroidal beta $\beta_{N}$, and edge safety factor $q_{95}$. A brief
summary is given in section VI.

##  II Data collection and normalization

Figure 1a illustrates the poloidal locations of the magnetic measurements used
as inputs to our model. A typical time evolution of some of the magnetic
measurements from EAST discharge 113019 are plotted in Figure 1b-f.

Figure 1:  Left: location of the magnetic probes, flux loops, and poloidal
coils on EAST. A typical time evolution of some of the magnetic measurements
from EAST discharge 113019 are plotted in the right panel. There are total 38
magnetic probes measuring the equilibrium poloidal magnetic field, and only 34
of them are working and are used in this work.

The inputs (features) to the neutral network (NN) are 84 magnetic
measurements: 35 poloidal magnetic flux ($\Psi_{\operatorname{FL}}$) values
measured by flux loops, 34 equilibrium poloidal magnetic field (MP) values
measured by magnetic probes, 14 poloidal field (PF) coil currents and 1 plasma
current (Ip) measured by a Rogowski loop.

The outputs (targets) of the neutral network are the values of the poloidal
magnetic flux $\Psi$ at $R\times Z=129\times 129=16641$ spatial locations. The
$\Psi$ used in the training process is computed by off-line EFIT and are
downloaded from EAST MDSplus server (mds.ipp.ac.cn). The input and the output
signals are interpolated to the same time slices before they are fed to the
NN.

The inputs and outputs are summarized in Table 1.

Signal | Measure method | Signal meaning | Num. of values
---|---|---|---
Input |  |  | 84
$\Psi_{\operatorname{FL}}$ | Flux loop | Poloidal magnetic flux | 35
MP | Magnetic probe | Poloidal magnetic field | 34
PF | Rogowski loop | Poloidal field coil current | 14
Ip |  Rogowski loop | Plasma current | 1
Output |  |  | 16641
$\Psi$(R,Z) | EFIT | Poloidal magnetic flux | 16641

Table 1:  The inputs and outputs of the model.

The data used in training, validation and testing process were downloaded from
the EAST MDSplus server by using Python API, which scans a series of
discharges and automatically skips discharges where necessary signals are
missing. Specifically, we scan every 5 discharges among all the discharges
spanning from #114000 to #117000, resulting in total 45,544 equilibria (time
slices). These discharges are from experiments performed in one EAST campaign
from June to July in 2022. This range is casually chosen with no particular
criterion, except that we prefer recent discharges and avoid old discharges
because locations of some magnetic probes were changed in previous campaigns.
The auxiliary heating methods on EAST used in this campaign include neutral
beam injection (50-70keV Deuterium beam), lower-hybrid waves (2.45GHz and
4.6GHz), electron cyclotron waves (140GHz), and ion cyclotron waves
(25-70MHz). Typical values of total heating source power are between 4-10MW.

Figure 2 is the distribution of the EFIT equilibrium data in
$(l_{i},\beta_{N})$ plane, where $l_{i}$ is the normalized internal induction
and $\beta_{N}$ is the normalized plasma beta.

Figure 2:  Distribution of the EFIT equilibrium data (training set \+
validation set+ testing set ) in $(l_{i},\beta_{N})$ plane.

The collected data are split into three sets: training set (81%), validation
set (9%), and testing set (10%), where training set is used in training the
NN, validation set is used in monitoring potential overfitting and tuning
hyperparameters, and testing set is used in testing the predicting capability
of the trained model.

Figure 1b-f shows that there is a difference of up to six orders of magnitude
in the values of the input signals. In order to eliminate scale differences
among features, we use the min-max normalization method to normalize the input
data. The general formula is given by

$x^{\prime}=\frac{x-x_{\min}}{x_{\max}-x_{\min}}$ (3)

where $x$ is the original value of the feature, $x^{\prime}$ is the normalized
value, $x_{\min}$ and $x_{\max}$ are respectively the minimal and maximal
value of a feature in the data sets excluding the testing set. The $x_{\min}$
and $x_{\max}$ obtained here are then used to normalize the input data in the
testing set when doing prediction using the trained NN.

Figure 3 plots the time evolution of the normalized input signals
corresponding to those in Fig. 1b-f.

Figure 3: Normalized signals corresponding to those in Fig. 1b-f. The sudden
change in the poloidal magnetic flux (Fl Loop 3) after 10s is brought about by
the PF1 coil current, which is actively adjusted by the plasma control system
to provide Ohm field to maintain a constant plasma current.

The magnitude of the output ($\Psi$ in SI units) is near 1, so no
normalization is applied to it.

##  III Model architecture and automatic hyperparameter tuning

An artificial neural network is a kind of computational network, which usually
consists of multiple layers: input layer, one or more inner layers (known as
hidden layers), and a layer of outputs. Each layer is made of units. Each unit
in the computing layers (hidden and output layers) receives information and
process the information using some linear transform (matrix multiplication)
and some nonlinear transform (activation function).

A fully-connected feed-forward network showed in figure 4 is used here to
predict the poloidal magnetic flux $\Psi$ based on the magnetic measurements.
Here “fully-connected” means that each unit of a computing layer receives
information from all the units of the previous layer. “Feed-forward” means
that information move in only one direction (from the input layer to the
hidden layers, and to the output layers), no cycles or loops, and no intra-
layer connections.

Figure 4:  Fully connected feed-forward neural network used in this work.

Each unit (neuron or node) in the computing layers has trainable parameters,
often called weights and biases. Denote the output of the
$j^{\operatorname{th}}$ neuron in $l^{\operatorname{th}}$ layer by
$a_{j}^{l}$, then a neural network model assumes that $a_{j}^{l}$ is related
to the $a^{l-1}$ (output of the previous layer) via

$a^{l}_{j}=\sigma\left(\sum_{k}w_{jk}^{l}a^{l-1}_{k}+b^{l}_{j}\right),$ (4)

where $w^{l}_{jk}$ and $b_{j}^{l}$ are the weight and bias, the summation is
over all neurons in the $(l-1)^{\operatorname{th}}$ layer, and $\sigma$ is a
function called activation function. The weights and biases will be adjusted
in the training process by gradient descent methods to reduce the loss (cost
or error) function, which is defined in this work as

$L(\mathbf{w},\mathbf{b})\equiv\frac{1}{2n}\sum_{i=1}^{n}\|\mathbf{y}_{i}-\hat{\mathbf{y}}_{i}\|^{2},$
(5)

where $\mathbf{y}_{i}$ is the EFIT poloidal magnetic flux and
$\hat{\mathbf{y}}_{i}$ is the NN output, and the summation is over all the
samples in the training set. The loss function in Eq. (5) is the mean squared
error (MSE). The loss function measures the derivation of the approximate
solution away from the desired exact solution. So the goal of a learning
algorithm is to find weights and biases that minimize the loss function. To
minimize the loss function over $(\mathbf{w},\mathbf{b})$ using the gradient
descent method, we need to compute the partial derivatives $\partial
L/\partial w^{l}_{jk}$ and $\partial L/\partial b^{l}_{j}$, which can be
efficiently computed by the well known back-propagating methodRumelhart1986
(15, 16). The back-propagating algorithm and the corresponding gradient
descent method are the core algorithms in all deep learning software
frameworks.

Besides the trainable parameters, there are various hyperparameters in a NN
that usually need to be set manually, such as number of hidden layers, units
per layer, activation function, NN optimizers, learning rate, batch size,
number of epochs of training. In recent years, there appear automatic
optimization libraries that can search for the best combination of
hyperparameters. In this work, we use the Optuna optimization
frameworkoptuna_2019 (14) in setting the hyperparameters. Optuna automates the
hyperparameter optimization process by defining a search space of
hyperparameters and exploring the space using efficient searching algorithms.
The tree-structured Parzen estimator (TPE) algorithm is used in this work.
This algorithm models the relationship between hyperparameters and their
corresponding performance metrics and makes efficient decisions on which
hyperparameters to try next.

Optuna automate the selection of the best hyperparameter combination. After
multiple experiments, we have found that the model accuracy is not sensitive
to the number of hidden layers, activation functions, and optimizers (an
example showing the relative importance of these hyperparameters is given in
Fig. 5). Therefore, these hyperparameters are fixed in the fine tuning step,
in order to improve the speed of the model selection process, and explore more
hyperparameter regimes to which the model may be sensitive. For other
hyperparameters, we use Optuna framework to find the optimal combination of
hyperparameters. Relative importance of these hyperparameters are shown in
Fig. 6. The above results indicate that learning rate is the dominant factor
that determines the model performance.

Figure 5:  Hyperparameter importance distribution. Model accuracy is not
sensitive to activation functions for this case.

Figure 6:  Relative importance of the hyperparameters in determining the model
accuracy in the fine tuning.

The final values of the hyperparameters used in the model are shown in Table
2.

Hyperparameter | Meaning | Final values
---|---|---
n_layers* | Number of hidden layers | 4
n_units | Number of nodes per hidden layer | 86
Activation* | Activation function | tanh
Optimizer* | Optimizer type | Adam
$\eta$ | Learning rate | $2.26\times 10^{-5}$
Loss* | Loss function | MSE
batch_size | Number of samples used in a step | 16
Epochs | Number of epochs | 97

* Fixed hyperparameters during fine tuning

Table 2:  Final values of hyperparameters of the model. The hyperparameters
with an asterisk (*) are fixed during the fine tuning. The above activation
function refers to that used in the hidden layers. For the output layer, the
linear activation function is used.

The network are constructed and trained using Keras & TensorFlow2
chollet2015keras (17, 18), which is a broadly adopted open source deep
learning framework in industry and research community. Figure. 7 plots loss
function values as a function of the training epochs. The loss function is
also evaluated on the validation set, which serves as a monitor for the
possible overfitting. The validation loss follows the same trend as the
training loss, indicating no overfitting.

Figure 7:  Training history showing the training and validation loss function
values versus the training epochs. One epoch corresponds to go through all the
samples in the training set. In the mini-batch stochastic gradient descent
method used here, each gradient descent step uses randomly selected 16 samples
(a mini batch), and the loss function values shown here are obtained by
summation over a mini batch of samples in the training set or validation set.

##  IV Performance of the neural network

###  IV.1 Performance of the model on testing set

After the model is trained on the training set, we assess its prediction
capability on the data that are not seen by the training process. To evaluate
the reconstruction quality, we employ three widely adopted metrics: the
Pearson correlation coefficient $r$ (definition is given in appendix A), the
coefficient of determination $R^{2}$ (definition is given in appendix B), and
the peak signal-to-noise ratio (PSNR) (definition is given in appendix C).

Figure 8a plots the NN prediction of the poloidal flux $\Psi_{NN}$ vs. EFIT
results $\Psi_{\operatorname{EFIT}}$ for the testing set (total 4555
equilibria, each with 16641 values). The Pearson correlation coefficient $r$
and the coefficient of determination $R^{2}$ are also shown in the figure,
which are very close to 1, indicating a strong predictive capability. Figure
8b plots the distribution of the correlation coefficient $r$ between NN
predictions and EFIT results for each equilibrium of the 4555 equilibria in
the testing set. The results indicate that the majority of the values are
greater than 0.998, indicating good correlation between NN prediction and EFIT
result for each equilibrium.

Figure 8:  (a) Neural network prediction of the poloidal flux $\Psi_{NN}$ vs.
EFIT results $\Psi_{\operatorname{EFIT}}$ for the testing set (total 4555
equilibria). Red dashed line is the $y=x$ line. The total number of points
shown here is $4555\times 129\times 129$. Color represents density of data
points. (b) The distribution of the correlation coefficient for each
equilibrium of the 4555 equilibria in the testing set.

To test the accuracy of the model in predicting the plasma magnetic surface,
we compare the 2D contours of the poloidal magnetic flux predicted by the NN
with those given by EFIT. The results are shown in figures 9(a), (c), (e) and
(g), where the NN predictions of $\Psi$ contours are overlaid on the $\Psi$
contours of EFIT. It displays four randomly selected samples from the 4555
equilibria in the testing set (the four displayed samples may not necessarily
come from the same discharge). Since our reconstruction results take the form
of images with resolution determined by the spatial grid points, it is also
useful to use $\operatorname{PSNR}$ in evaluating the reconstruction quality
of the magnetic surface. The values of PSNR for the four equilibrium are shown
in the figure.

Figure 9:  Panels (a), (c), (e) and (g) compare the magnetic surfaces from
EFIT (solid black lines) with those from the NN (dashed red lines) for four
randomly selected time slices from the 4555 testing samples. Panels (b), (d,)
(f) and (h) show the corresponding normalized residua [$\Psi$(R,Z)
$\Psi^{\operatorname{NN}}$(R,Z)]/max(|$\Psi$ (R,Z)|).

To further assess the accuracy of the mode, we locate the LCFSs predicted by
the NN model and compare them with those given by EFIT. The LCFSs
corresponding to the four equilibrium of Fig. 9 are shown in Fig. 10, which
indicates that the NN and EFIT results are in good agreement. Minor
discrepancies appear near the X points.

Figure 10:  Comparison of LCFS between the NN prediction and EFIT result. The
four panel correspond to the four equilibria shown in fig. 9. In locating the
LCFS from the NN prediction, we analyze a series of $\Psi$ contours and
determine the outermost contour that are near the magnetic separatrix.

###  IV.2 Performance of the model on four complete discharges

In this section, we arbitrarily select 3 full discharges that are not in the
dataset used above to examine the time evolution of the magnetic configuration
during an entire discharge (from ramp-up to flat-top then to ramp-down).

Besides the plasma currents $I_{p}$, we also calculate the normalized internal
inductance $l_{i}$ (definition is given in Appendix D), which is a quantity
that is solely determined by $\Psi$ and thus can further reflect how accuracy
the predicted $\Psi$ is. We plot the time evolution of $l_{i}$ and compare it
with the EFIT results. By doing this, we can assess the accuracy of the NN in
predicting the time evolution of some key volume-integrated quantities
characterizing magnetic configuration.

Figure 11 compares the time evolution of $I_{p}$ and $l_{i}$ predicted by the
NN and that by EFIT for discharge #113388.

Figure 11:  (a) Comparison of time evolution of plasma current given by the NN
and EFIT. (b) Comparison of time evolution of $l_{i}$ given by the NN and
EFIT. Four time slices are indicated on the graph, which are time slices
selected for the magnetic configuration comparison shown in Figs. 12-13.

Figure 12 compares the contours of $\Psi$ given by the NN model and that given
by EFIT at 4 time slices (indicated in Fig. 11) in discharge #113388. The
results indicate the relative error between the NN and EFIT results is less
than $2\%$.

Figure 12:  Panels (a), (c), (e) and (g) compare the magnetic surfaces from
EFIT (solid black lines) with those from the NN (dashed red lines) for
discharge #113388. $\Psi$ at four different time slices during the discharge,
0.580s (early ramp-up), 2.680s,, 5.130s (flat top), 7.370s (ramp-down) are
shown. Panels (b), (d,) (f) and (h) show the corresponding relative error
$(\Psi_{\operatorname{EFIT}}-\Psi_{\operatorname{NN}})/\max(|\Psi_{\operatorname{EFIT}}|)$.

Figure 13 compares the LCFSs given by the NN model and that given by EFIT at 4
time slices in discharge #113388. The results show good agreement between the
two models. Minor differences usually appear in the ramp up/down phase, and
near the X-points.

Figure 13:  NN reconstructions of the LCFS (dotted-dashed red) for EAST shot
#113388 overlaid against the EFIT LCFS (solid black). LCFSs from four
different times in the discharge #113388 are shown.

Similar results for discharge #117016 are shown in Fig. 14-16.

Figure 14:  Same as figure 11, except that the discharge is #117016.

Figure 15:  Same as figure 12, except that the discharge is #117016.

Figure 16:  Same as figure 13, except that the discharge is #117016.

Similar results for discharge #113019 are shown in Fig. 17-19.

Figure 17:  Same as figure 11, except that the discharge is #113019.

Figure 18:  Same as figure 12, except that the discharge is #113019.

Figure 19:  Same as figure 13, except that the discharge is #113019.

##  V Neural network prediction of $W_{\operatorname{mhd}}$, $\beta_{N}$, and
$q_{95}$

Besides the $l_{i}$ discussed above, there are some other global parameters
that can be constructed from the magnetic measurements, namely the plasma
stored energy $W_{\operatorname{mhd}}$, normalized plasma beta $\beta_{N}$,
and edge safety factor $q_{95}$. These parameters depend on information beyond
the poloidal magnetic flux, namely the toroidal magnetic field and plasma
pressure. Therefore they can not be fully determined by using only the
poloidal magnetic flux predicted from the above network. Following Ref. DIII-
D_2022 (8), we construct a new NN for predicting these parameters (called NN2
in the following; the previous one will be called NN1), where the network has
only 3 output values, namely $W_{\operatorname{mhd}}$, $\beta_{N}$, $q_{95}$.
The input to NN2 includes a new signal, the current in the toroidal field (TF)
coils, which determines the toroidal field. (In the NN1, this signal is not
included because it has negligible effect on the prediction of the poloidal
magnetic flux.) The NN2 has only one hidden layer consisting of 16 units, and
uses the sigmoid as the activation function for both the hidden and output
layers. The input and output signals of NN2 are normalized by using the same
min-max scalar as used for NN1.

The training data consist of about $1/4$ randomly selected part of the data
used for NN1. We found that using larger dataset makes this small network
prone to overfitting. The testing set consists of 1000 time slices. Figure 20
plots the NN2 predictions against the EFIT values for the testing set. The
results indicate that the NN2 predictions are in reasonable agreement with the
EFIT values for all the 3 parameters. The NN2 predictions of $q_{95}$ are a
little worse than those of the other two parameters, judging from the values
of $r$ and $R^{2}$.

Figure 20: NN2 predictions of $W_{\operatorname{mhd}}$, $\beta_{N}$, and
$q_{95}$ against the EFIT values for the testing set.

To evaluate the accuracy of NN2 prediction for a full discharge, we
arbitrarily chosen a discharge and compare the time evolution of
$W_{\operatorname{mhd}}$, $\beta_{N}$, and $q_{95}$ between the NN2
predictions and EFIT values. The results are shown in Fig. 21, which shows
good agreement between the network predictions and EFIT values.

Figure 21: Time evolution of plasma stored energy $W_{\operatorname{mhd}}(J)$,
normalized beta $\beta_{N}$, and edge safety factor $q_{95}$ in EAST discharge
#100000. Solid lines are EFIT values and dash-dot lines are network
predictions.

##  VI Summary and discussion

In this work, we train a multiply-layer neural network on the magnetic
measurements (input) and EFIT poloidal magnetic flux (output) on EAST tokamak.
The prediction capability of the network is examined by comparing the
reconstructed magnetic surfaces, last closed flux surfaces, plasma current,
and normalized internal inductance with those of EFIT. The neural network
shows good agreement with EFIT for the data unseen in the training process.

In constructing the neural network, we use automatic optimization in searching
for the best hyperparameters of the model. The hyperparameters found this way
turn out to be better than our previously manually set hyperparameters in
terms of the model accuracy.

Based on the model’s good prediction capability and efficiency in terms of
computational time (about 0.5ms per equilibrium on a desktop computer, using
an 11th Gen Intel(R) Core(TM<EMAIL_ADDRESS>CPU with a single thread), it
looks promising to apply the neural network to real-time magnetic
configuration control. The above computational time does not include the time
used for tracing boundary/internal magnetic surfaces, and other related
calculations to obtain $l_{i}$. These computations (not optimized in this
work) seem too inefficient to be used in real-time control. The purpose of
computing $l_{i}$ for the NN1 model is to evaluate the accuracy of the
predicted $\Psi$. To predict these volume-integrated parameters, one usually
uses an additional small network, as we did in Sec. V, which is efficient
enough for real-time control because the network size is usually very small.

This work is limited to magnetic measurements. We plan to add more diagnostics
relating to the inner safety factor profiles and pressure profiles into the
model, in order to construct more realistic equilibria. This will rely on the
kinetic EFIT output. We are accumulating these kind of training data.

##  VII Data availability

The data that supports the findings of this study are available from the
corresponding author upon reasonable request.

##  VIII Acknowledgments

The authors thank Ting Lan, Tonghui Shi, Chengguang Wan, Zhengping Luo, Yao
Huang, Guoqiang Li, and Jingping Qian for useful discussions. This work was
supported by Comprehensive Research Facility for Fusion Technology Program of
China under Contract No. 2018-000052-73-01-001228, by Users with Excellence
Program of Hefei Science Center CAS under Grant No. 2021HSC-UE017, and by the
National Natural Science Foundation of China under Grant No. 11575251.

###  VIII.1 Conflict of interest

The authors have no conflicts to disclose.

##  Appendix A Pearson correlation coefficient $r$

The Pearson correlation coefficient $r$ is a statistical measure used to
assess the strength and direction of a linear relationship between the
predicted and true values of the data. It ranges from -1 to 1, where 1
indicates a perfect positive correlation, 0 indicates no correlation, and -1
indicates a perfect negative correlation. The formula for $r$ is

$r=\frac{\sum^{n}_{i=1}(y_{i}-\bar{y})(\hat{y}_{i}-\overline{\hat{y}})}{\sqrt{\sum^{n}_{i=1}(y_{i}-\bar{y})^{2}}\sqrt{\sum^{n}_{i=1}(\hat{y}_{i}-\overline{\hat{y}})^{2}}},$
(6)

where $n$ is the number of data in the testing set, $y_{i}$ is the value given
by EFIT, $\hat{y}_{i}$the prediction by the NN,
$\bar{y}=\frac{1}{n}\sum^{n}_{i=1}y_{i}$ is the mean value of the values given
by EFIT, and $\overline{\hat{y}}=\frac{1}{n}\sum_{i=1}^{n}\hat{y}_{i}$is the
mean value predicted by the NN.

##  Appendix B Coefficient of determination $R^{2}$

Another relevant metric used to assess how well a model fits the data is the
coefficient of determination $R^{2}$, which is defined by

$R^{2}=1-\frac{\sum_{i=1}^{n}(y_{i}-\hat{y}_{i})^{2}}{\sum_{i=1}^{n}(y_{i}-\bar{y})^{2}},$
(7)

where $n,y_{i}$, $\hat{y}_{i}$and $\bar{y}$ mean the same as in section 6. The
value of $R^{2}$ ranges from arbitrary negative values to 1, where 1
represents a perfect fit between the model predictions and the actual data
points. A higher value of $R^{2}$ suggests that the model is a better fit for
the data. The coefficient of determination $R^{2}$ is usually not equal to the
squared Pearson correlation coefficient except in some specific cases.

##  Appendix C Peak signal-to-noise ratio (PSNR)

The PSNR is a metric that measures the quality of an image by comparing the
original image to a reconstructed version. A higher PSNR value indicates a
higher quality reconstruction. It is defined by

$\displaystyle\operatorname{PSNR}$ $\displaystyle=$ $\displaystyle
10\times\log_{10}\left(\frac{\max(y_{i})^{2}}{\operatorname{MSE}}\right)$ (8)
$\displaystyle=$ $\displaystyle
10\times\log_{10}\left(\frac{\max(y_{i})^{2}}{\frac{1}{M}\sum^{M}_{i=1}(y_{i}-\hat{y}_{i})^{2}}\right)$

where $\max(y_{i})$ is the maximum value of $\Psi$ given by EFIT in $(R,Z)$
plane, and MSE is the mean squared error between the EFIT and NN.

##  Appendix D Normalized internal inductance

The normalized internal inductance $l_{i}$ is defined by

$l_{i}=\frac{\langle B^{2}_{\theta}\rangle_{P}}{\langle
B^{2}_{\theta}\rangle_{S}},$ (9)

where $P$ is the integration over the plasma volume, $\langle
B^{2}_{\theta}\rangle_{S}$ is the surface average of poloidal field over the
plasma boundary. $l_{i}$ reflects the peakness of the plasma current density
profile: a small value of $l_{i}$ corresponds to a broad current profile.

For circular cross section with minor radius $a$ and assuming $B_{\theta}$ is
independent of the poloidal angle, then, Ampere’s law gives
$B_{\theta}(a)=\mu_{0}I/(2\pi a)$. Then $\langle B^{2}_{\theta}\rangle_{S}$ is
approximated as

$\langle B^{2}_{\theta}\rangle_{S}\approx
B^{2}_{0}(a)=\frac{{\mu^{2}_{0}}I^{2}}{4\pi^{2}a^{2}},$ (10)

Using this and noting $V\approx\pi a^{2}2\pi R_{0}$, where $R_{0}$ is the
major radius of the device, Eq. (9) is written as

$l_{i}=\frac{4\pi^{2}a^{2}}{{\mu^{2}_{0}}I^{2}}\langle
B^{2}_{\theta}\rangle_{P}=\frac{4\pi^{2}a^{2}R_{0}}{{\mu^{2}_{0}}I^{2}R_{0}}\langle
B^{2}_{\theta}\rangle_{P}=\frac{2V}{{\mu^{2}_{0}}I^{2}R_{0}}\langle
B^{2}_{\theta}\rangle_{P}.$ (11)

Eq. (11) is used in this work to calculate $l_{i}$.

## References

  * (1) L. Lao, H. S. John, R. Stambaugh, A. Kellman, and W. Pfeiffer,  Nuclear Fusion 25, 1611 (1985).
  * (2) L. Lao, J. Ferron, R. Groebner, W. Howl, H. S. John, E. Strait, and T. Taylor,  Nuclear Fusion 30, 1035 (1990).
  * (3) L. L. Lao, H. E. St John, Q. Peng, J. R. Ferron, E. J. Strait, T. S. Taylor, W. H. Meyer, C. Zhang, and K. I. You,  Fusion Science and Technology 48 (2005).
  * (4) Q. Jinping, W. Baonian, L. L. Lao, S. Biao, S. A. Sabbagh, S. Youwen, L. Dongmei, X. Bingjia, R. Qilong, G. Xianzu, and L. Jiangang,  Plasma Science and Technology 11, 142 (2009).
  * (5) D. O’Brien, L. Lao, E. Solano, M. Garribba, T. Taylor, J. Cordey, and J. Ellis,  Nuclear Fusion 32, 1351 (1992).
  * (6) Y. Park, S. Sabbagh, J. Berkery, J. Bialek, Y. Jeon, S. Hahn, N. Eidietis, T. Evans, S. Yoon, J.-W. Ahn, et al.,  Nuclear Fusion 51, 053001 (2011).
  * (7) S. Joung, J. Kim, S. Kwak, J. Bak, S. Lee, H. Han, H. Kim, G. Lee, D. Kwon, and Y.-C. Ghim,  Nuclear Fusion 60, 016034 (2019).
  * (8) L. L. Lao, S. Kruger, C. Akcay, P. Balaprakash, T. A. Bechtel, E. Howell, J. Koo, J. Leddy, M. Leinhauser, Y. Q. Liu, S. Madireddy, J. McClenaghan, D. Orozco, A. Pankin, D. Schissel, S. Smith, X. Sun, and S. Williams,  Plasma Physics and Controlled Fusion 64, 074001 (2022).
  * (9) B. Wan, Y. Liang, X. Gong, J. Li, N. Xiang, G. Xu, Y. Sun, L. Wang, J. Qian, H. Liu, X. Zhang, L. Hu, J. Hu, F. Liu, C. Hu, Y. Zhao, L. Zeng, M. Wang, H. Xu, G. Luo, A. Garofalo, A. Ekedahl, L. Zhang, X. Zhang, J. Huang, B. Ding, Q. Zang, M. Li, F. Ding, S. Ding, B. Lyu, Y. Yu, T. Zhang, Y. Zhang, G. Li, T. Xia, the EAST team, and Collaborators,  Nuclear Fusion 57, 102019 (2017).
  * (10) Q. Jinping, W. Baonian, L. L. Lao, S. Biao, S. A. Sabbagh, S. Youwen, L. Dongmei, X. Bingjia, R. Qilong, G. Xianzu, and L. Jiangang,  Plasma Science and Technology 11, 142 (2009).
  * (11) G. Q. Li, Q. L. Ren, J. P. Qian, L. L. Lao, S. Y. Ding, Y. J. Chen, Z. X. Liu, B. Lu, and Q. Zang,  Plasma Physics and Controlled Fusion 55, 125008 (2013).
  * (12) L. Zhengping, X. Bingjia, Z. Yingfei, and Y. Fei,  Plasma Science and Technology 12, 412 (2010).
  * (13) C. Wan, Z. Yu, A. Pau, O. Sauter, X. Liu, Q. Yuan, and J. Li,  Nuclear Fusion 63, 056019 (2023).
  * (14) T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama,  CoRR abs/1907.10902 (2019).
  * (15) D. E. Rumelhart, G. E. Hinton, and R. J. Williams,  Nature 323, 533 (1986).
  * (16) M. A. Nielsen,  Neural networks and deep learning, volume 25,  Determination press San Francisco, CA, USA, 2015.
  * (17) F. Chollet et al.,  Keras,  https://keras.io, 2015.
  * (18) M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng,  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015,  Software available from tensorflow.org.


# Dynamic-based entanglement witnesses for harmonic oscillators

Pooja Jayachandran  Centre for Quantum Technologies, National University of
Singapore, 3 Science Drive 2, Singapore 117543  Lin Htoo Zaw  Centre for
Quantum Technologies, National University of Singapore, 3 Science Drive 2,
Singapore 117543  Valerio Scarani  Centre for Quantum Technologies, National
University of Singapore, 3 Science Drive 2, Singapore 117543  Department of
Physics, National University of Singapore, 2 Science Drive 3, Singapore 117542

###### Abstract

We introduce a family of entanglement witnesses for continuous variable
systems, which rely on the sole assumption that their dynamics is that of
coupled harmonic oscillators at the time of the test. Entanglement is inferred
from the Tsirelson nonclassicality test on one of the normal modes, without
any knowledge about the state of the other mode. In each round, the protocol
requires measuring only the sign of one coordinate (e.g. position) at one
among several times. This dynamic-based entanglement witness is more akin to a
Bell inequality than to an uncertainty relation: in particular, it does not
admit false positives from classical theory. Our criterion detects non-
Gaussian states, some of which are missed by other criteria.

##  I Introduction

Once troublesome to the founders of quantum mechanics [1, 2], entanglement is
now well-established as one of the defining features of quantum theory. While
entanglement in discrete systems has gone through much scrutiny [3, 4], the
field of continuous variable (CV) quantum entanglement has had its challenges
outside the Gaussian regime [5]. The obvious testbed for CV entanglement are
optical modes, but there has also been significant progress for massive
oscillators [6]. In this work, we present a criterion to certify the
entanglement of two CV degrees of freedom that exploits the knowledge of the
dynamics of the systems. We call it a dynamic-based entanglement witness (DEW)
111The name “dynamical entanglement witness” would have been more elegant, but
could have been confused with a “witness of dynamical entanglement”, which is
a notion of entanglement for channels [27]. Here we are concerned with
entanglement of states. Specifically, we consider degrees of freedom that
undergo types of harmonic dynamics, and build on a nonclassicality test for a
single oscillator [8, 9].

A priori, our criterion exhibits two elegant features: the absence of false
positives from classical theory, and the fact that only one observable needs
to be measured. These features, to be defined precisely below, can be
appreciated by contrast with existing entanglement witnesses, based on
generalized uncertainty relations [10, 11, 12, 13, 14]. Consider specifically
the criterion by Duan  _et al._ [10]. Generalizing the observables defined in
the Einstein–Podolsky–Rosen (EPR) argument [1], these authors defined the
commuting dimensionless variables
$u=\absolutevalue{c}\tilde{x}_{1}+\frac{1}{c}\tilde{x}_{2}$ (with
$\tilde{x_{j}}=x_{j}\sqrt{{m_{j}\omega_{j}}/{\hbar}}$) and
$v={\absolutevalue{c}}\tilde{p}_{1}-\frac{1}{c}\tilde{p}_{2}$ (with
$\tilde{p}_{j}=p_{j}/\sqrt{m_{j}\hbar\omega_{j}}$) for some $c\in\mathbb{R}$.
The two subsystems are then entangled if $\langle{\pqty{\Delta
u}^{2}}\rangle+\langle{\pqty{\Delta v}^{2}}\rangle<c^{2}+\frac{1}{c^{2}}$.
This requires measuring both positions and momenta, and with a precision set
by $\hbar$. At the precision of (say) human perception, two springs at
equilibrium are described by $\tilde{x}_{j}=\tilde{p}_{j}=0$, values which
would imply entanglement if plugged naively in the criterion above. Gross
though it is, this example shows the danger of false positives.

A posteriori, we find that our criterion detects states with negative Wigner
functions (thus, non-Gaussian), some of which are missed by all existing
criteria. Thus, besides being elegant, our DEW is also a useful addition to
the existing toolbox.

##  II The single-oscillator protocol

We review the Tsirelson nonclassicality test [8] following the generalization
given in [9]. The assumption is that the physical quantity $A_{1}$ is
undergoing a uniform precession at pulsation $\omega$, i.e.
$A_{1}(t)=A_{1}(0)\cos\omega t+A_{2}(0)\sin\omega t$, where $A_{2}$ is another
physical quantity. For classical systems, $A_{1}(t)$ is the value of $A_{1}$
at time $t$; for quantum systems, it is the corresponding observable in the
Heisenberg representation. The protocol for the test (which we call precession
protocol hereafter) goes as follows. In each round, the sign of $A_{1}$ is
measured at one of $K$ different times given by $t_{k}=(k/K)T$, where $K>0$,
$k=0,1,...,K-1$, and $T$ is the period of oscillation. After several rounds,
one estimates

$\displaystyle{P}_{K}$ $\displaystyle=$
$\displaystyle\frac{1}{K}\sum_{k=0}^{K-1}\Bqty{\Pr[A_{1}(t_{k})>0]+\frac{1}{2}\Pr[A_{1}(t_{k})=0]},\,\,\,$
(1)

where the second term in the bracket was introduced in [9] to avoid singular
behaviors for states with non-infinitesimal concentration on $A_{1}=0$. By
inspection [9], the upper bound $P_{K}\leq\mathbf{P}_{K}^{c}$ for a classical
theory is easily derived: $\mathbf{P}_{K}^{c}=\frac{1}{2}$ for $K$ even, and

$\displaystyle\begin{aligned}
\mathbf{P}_{K}^{c}&=\frac{1}{2}\left(1+\frac{1}{K}\right)\text{ for $K$
odd.}\end{aligned}$ (2)

Remarkably, in spite of the fact that the precessing dynamics is identical to
the classical one, there exist quantum states for which
$P_{K}>\mathbf{P}_{K}^{c}$ for any odd $K>1$.

For the remainder of the paper, we focus on the harmonic oscillator, i.e. a
material point, whose time evolution is governed by the Hamiltonian
$H=\frac{1}{2m}p^{2}+\frac{1}{2}m\omega^{2}x^{2}$. The pair
$(A_{1},A_{2})=(\tilde{x},\tilde{p})$ clearly precesses at pulsation $\omega$
and thus satisfies the assumption. On a given state $\rho$, quantum theory
predicts $P_{K}=\Tr(\rho Q_{K})$ where

$\displaystyle Q_{K}$ $\displaystyle=$
$\displaystyle\frac{1}{K}\sum_{k=0}^{K-1}\operatorname{pos}(X(t_{k}))\,,$ (3)

with $\operatorname{pos}(X)$ defined by
$\operatorname{pos}(X)\ket{x}=\frac{1}{2}(1+\operatorname{sgn}(x))\ket{x}$.

The maximum quantum score
$\mathbf{P}_{K}=\max_{\ket{\psi}}\expectationvalue{Q_{K}}{\psi}$ (denoted
$\mathbf{P}^{\infty}_{K}$ in [9]) is achieved by $\ket{\mathbf{P}_{K}}$, the
eigenstate of $Q_{K}$ with the largest eigenvalue. Tsirelson proved that
$\mathbf{P}_{3}\gtrsim 0.709>\mathbf{P}_{3}^{c}=2/3$ [8]; similar violations
are found for all odd $K$ [9]. The violation can be attributed to having
suitable patterns in the Wigner function, in particular, suitably distributed
negativities. A state with positive Wigner function cannot give any violation.

##  III Entanglement of two harmonic oscillators

Now, the key insight is that $x$ may be the position of an effective
oscillator, built out of two (or more) physical ones. We focus on the case of
two oscillators, with arbitrary masses and frequencies, and an $x$–$x$
coupling. The standard decomposition in normal modes yields

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{2}\pqty{\frac{1}{2m_{j}}p_{j}^{2}+\frac{1}{2}m_{j}\omega_{j}^{2}x_{j}^{2}}-\frac{1}{2}gx_{1}x_{2}$
(4) $\displaystyle=$
$\displaystyle\sum_{\sigma\in\\{+,-\\}}\frac{1}{2\mu}p_{\sigma}^{2}+\frac{1}{2}\mu\omega_{\sigma}^{2}x_{\sigma}^{2},$

where $\mu=\sqrt{m_{1}m_{2}}$,

$\displaystyle x_{+}(t)$
$\displaystyle=\pqty{\frac{m_{1}}{m_{2}}}^{1/4}\cos\theta\;x_{1}(t)+\pqty{\frac{m_{2}}{m_{1}}}^{1/4}\sin\theta\;x_{2}(t)$
(5) $\displaystyle x_{-}(t)$
$\displaystyle=\pqty{\frac{m_{2}}{m_{1}}}^{1/4}\cos\theta\;x_{2}(t)-\pqty{\frac{m_{1}}{m_{2}}}^{1/4}\sin\theta\;x_{1}(t),$

with mixing angle
$\theta=\operatorname{arctan2}\bqty{g,\mu(\omega_{1}^{2}-\omega_{2}^{2})}/2$,
and normal frequencies
$\omega^{2}_{\pm}=\frac{\omega_{1}^{2}+\omega_{2}^{2}}{2}\pm\sqrt{\pqty{\frac{\omega_{1}^{2}-\omega_{2}^{2}}{2}}^{2}+\frac{g^{2}}{4\mu^{2}}}$.
The time evolution of the $x_{\sigma}(t)$ is a uniform precession around phase
space with the period $T_{\sigma}=2\pi/\omega_{\sigma}$. Therefore, the
single-oscillator protocol can be performed for coupled oscillators with
different frequencies by measuring the normal modes $x_{\sigma}(t_{k}+t_{0})$
at times $t_{k}=(k/K)T_{\sigma}$ for $k=0,1,\dots,K-1$. Among the many ways to
estimate $P_{K}$ for $x_{\sigma}$, one consists in measuring separately
$x_{1}$ and $x_{2}$ in each round, forming the corresponding $x_{\sigma}(t)$,
and checking its sign. Notice also that, up to a multiplicative constant,
$x_{\sigma}(t)$ has the same form as $u$ defined in the entanglement criterion
by Duan and coworkers [10].

It is important to stress that $H$ describes the dynamics during the
certification protocol, not the interaction that prepared the state under
study. Thus, all values of the parameters are allowed. In particular, the
certification protocol can be performed when $g=0$ and
$\omega_{1}=\omega_{2}$. In this case, $\theta$ can take on any value: indeed,
for uncoupled oscillators precessing at the same frequency, all linear
combinations of $x_{1}$ and $x_{2}$ are normal modes at that same frequency.

If $P_{K}>\mathbf{P}_{K}^{c}$ for $x_{\sigma}$, the state of that mode has
certainly a negative Wigner function. We want to study when one can further
infer that the physical subsystems are entangled. This is not straightforward
because, by performing the protocol on one of the normal modes, we learn
nothing about the state of the other mode: the latter could be very mixed; or
the two normal modes may be even entangled. We are going to provide the
conditions under which entanglement can indeed be certified.

##  IV Results

For the quantum system, we denote the annihilation operators of the two
physical oscillators as $\\{a_{1},a_{2}\\}$: they are the subsystems whose
entanglement we want to certify. As hinted, $x_{\sigma}(t)$ is the position of
an effective oscillator denoted by the annihilation operator $a_{\sigma}$.
Specifically, let $\\{a_{+},a_{-}\\}$ be a new basis of modes, related to the
original by the passive transformation

$\pmqty{undef}=\pmqty{undef}\pmqty{undef}.$ (6)

Here, $\theta$ is the mixing angle as previously defined. We are going to
study the operator $Q_{K}$ given by Eq. (3) for the position operator
$X_{\sigma}(t)=\sqrt{\frac{\hbar}{2\mu\omega_{\sigma}}}\pqty{a_{\sigma}e^{-i\omega_{\sigma}t}+a_{\sigma}^{\dagger}e^{i\omega_{\sigma}t}}$.

When $\theta=\pi/4$, we have an analytical proof of existence of entanglement:

######  Result.

If the precession protocol is performed with $\theta=\pi/4$, all states that
violate the classical bound (2) for $a_{+}$ (or $a_{-}$) are entangled in
$\\{a_{1},a_{2}\\}$.

Proof. We show this for $a_{+}$ with a proof by contradiction. The proof for
$a_{-}$ proceeds in a similar way.

Take $\rho=\sum_{k}p_{k}\rho_{1}^{(k)}\otimes\rho_{2}^{(k)}$ separable in the
$\\{a_{1},a_{2}\\}$ subsystems. Its Wigner function is

$W_{\rho}(\alpha_{1},\alpha_{2})=\sum_{k}p_{k}W_{\rho_{1}}^{(k)}(\alpha_{1})W_{\rho_{2}}^{(k)}(\alpha_{2}),$
(7)

where $\\{\alpha_{1},\alpha_{2}\\}$ are the phase-space coordinates in the
$\\{a_{1},a_{2}\\}$ modes. For $\theta=\pi/4$ in (6), the Wigner function in
terms of $\\{\alpha_{+},\alpha_{-}\\}$, the phase-space coordinates in the
$\\{a_{+},a_{-}\\}$ modes, can be found with a straightforward coordinate
transformation:

$W_{\rho}(\alpha_{+},\alpha_{-})=\sum_{k}p_{k}W_{\rho_{1}}^{(k)}\\!\pqty{\tfrac{\alpha_{+}-\alpha_{-}}{\sqrt{2}}}W_{\rho_{2}}^{(k)}\\!\pqty{\tfrac{\alpha_{+}+\alpha_{-}}{\sqrt{2}}}.$
(8)

The measurement outcome in the $a_{+}$ basis is determined solely by the
reduced Wigner function
$W_{\tr_{-}(\rho)}(\alpha_{+})=\int\differential{\alpha_{-}}W_{\rho}(\alpha_{+},\alpha_{-})$.
By converting the passive coordinates in the arguments into active
transformations of the states, we find

$\displaystyle W_{\tr_{-}(\rho)}(\alpha_{+})$ $\displaystyle=$ $\displaystyle
2\sum_{k}p_{k}\int\differential[2]{\gamma}W_{\tilde{\rho}_{1}(\alpha_{+})}^{(k)}(\gamma)W_{\tilde{\rho}_{2}(\alpha_{+})}^{(k)}(\gamma)$
(9) $\displaystyle=$
$\displaystyle\frac{2}{\pi}\sum_{k}p_{k}\tr(\tilde{\rho}_{1}(\alpha_{+})\tilde{\rho}_{2}(\alpha_{+})).$

where

$\displaystyle\tilde{\rho}_{1}(\alpha_{+})=$ $\displaystyle
D\\!\pqty{\tfrac{1}{\sqrt{2}}\alpha_{+}}e^{-i\pi a^{\dagger}a}\rho_{1}e^{i\pi
a^{\dagger}a}D^{\dagger}\\!\pqty{\tfrac{1}{\sqrt{2}}\alpha_{+}}$
$\displaystyle\tilde{\rho}_{2}(\alpha_{+})=$ $\displaystyle
D\\!\pqty{-\tfrac{1}{\sqrt{2}}\alpha_{+}}\rho_{2}D^{\dagger}\\!\pqty{-\tfrac{1}{\sqrt{2}}\alpha_{+}}\,.$

and $D(\alpha)$ is the usual displacement operator. Thus
$W_{\tr_{-}(\rho)}\geq 0$, since it is a convex sum of inner products between
density operators. However, negativity in the Wigner function of $a_{+}$ is
necessary for a violation of the classical bound of the precession protocol
[9]. Therefore, when $\theta=\pi/4$, _any_ violation of the classical bound
witnesses entanglement of the $\\{a_{1},a_{2}\\}$ subsystems. ∎

Figure 1: (a) Heat map of logarithmic negativity as a function of the mixing
angle $\theta$ and the score $P_{3}$. Only the range $0\leq\theta\leq\pi/4$ is
shown: other values of $\theta$ corresponds to this range with a sign change
of $x_{1}$ and $x_{2}$, which can be effected with a local unitary on the
$\\{a_{1},a_{2}\\}$ basis and hence does not affect the amount of
entanglement. The dimension of the Hilbert space is truncated at
$\mathbf{n}=11$ for both modes $a_{+}$ and $a_{-}$. The dotted line separates
the states with $S_{N}>0$ (such that $S_{N}-\varepsilon>0$, where
$\varepsilon$ is the dual gap of the SDP) and those with $S_{N}=0$. (b)
Horizontal line cuts of the heat map. For fixed $P_{3}$, $S_{N}$ increases
with $\theta$. (c) Vertical line cuts of the heat map. For fixed $\theta$,
$S_{N}$ increases with $P_{3}$.

Next, we are going to study in detail the protocol with $K=3$. For any value
of $\theta$ and of $P_{3}>\mathbf{P}_{3}^{c}=\frac{2}{3}$, we are going to
compute a numerical lower bound on the amount of certifiable entanglement of
$\\{a_{1},a_{2}\\}$. We choose the logarithmic negativity
$S_{N}(\rho)=\log\tr\\!\absolutevalue{\rho^{\Gamma_{2}}}$ as the quantifier of
entanglement. Since
$\min_{\rho}S_{N}(\rho)=\log\min_{\rho}\tr\\!\absolutevalue{\rho^{\Gamma_{2}}}$,
and

$\displaystyle\min_{\rho}\tr\\!\absolutevalue{\rho^{\Gamma_{2}}}$ (10)
$\displaystyle\;\;\text{subject to}\;\tr(\rho Q_{3}(\theta))=P_{3}$

is a minimization of the trace norm under convex constraints, it can be cast
as a standard semidefinite program (SDP) when $\rho$ is truncated in the basis
of the excitations of $a_{+}$ and $a_{-}$ [15]. We run the SDP for truncation
$0\leq n_{+},n_{-}\leq\mathbf{n}=11$ (both, the form of the SDP and the choice
of the truncation are described in Appendix A). The results are plotted in
Fig. 1. We observe that the certifiable $S_{N}$ increases with $\theta$ and
$P_{3}$, as shown more explicitly by the line cuts. We already knew that every
$P_{3}$ certifies entanglement when $\theta={\pi}/{4}$, and the graph
indicates that this remains true down to $\theta\simeq{\pi}/{8}$. Below this
value, one needs a sufficiently large $P_{3}$, a low violation of the
classical bound being compatible with separable states. When $\theta=0$, the
precession protocol is performed on the first oscillator, and so no amount of
violation detects entanglement.

##  V Comparison with other witnesses

Now we put our DEW in the context of entanglement witnesses for continuous
variables (CVs), by comparing it to other criteria.

First of all, our DEW uses _quadrature measurements_ in the terminology of
quantum optics. Other measurements than quadratures can be used to witness CV
entanglement: for instance, one witness in Zhang  _et al._ [11] uses local
measurements of the generators of $\mathrm{SU}(N)$. In fact, any CV
entanglement can in principle be witnessed by projecting the state into a
finite dimensional subspace, then applying techniques to witness entanglement
of qudits [16]. Quadrature measurements have the appeal of having a classical
analog, are practical in many platforms, and in some setups may be even the
only available ones at this time (e.g. optomechanical systems with large
masses).

As already mentioned in the introduction, the other entanglement witnesses we
are aware of are open to false positives from classical theory [10, 11, 12,
13, 14]. By contrast, in the case of our DEW, poor precision or wrong
calibration may prevent the detection of entanglement, but will not lead to
false positives. This is similar to what happens with Bell inequalities, where
noise and lack of precision may decrease or cancel the violation, but not fake
it.

Having mentioned this, it is natural to move on to the comparison in terms of
characterization of the devices. Fully device-independent entanglement
witnesses (Bell inequalities) that use only quadrature measurements have been
hard to find: the few known examples are very specific [17, 18, 19]. Recently,
a measurement-device-independent criterion was introduced, under the
assumption that a trusted source of coherent states is available [20].
Meanwhile, our DEW is semi-device-independent: it works under the assumptions
that the dynamics is a uniform precession, and that the same quadrature is
measured whatever time is picked. Both assumptions are well-defined both in
classical and in quantum theory. In fact, the notions of “position” and
“uniform precession” are operational, and their meaning immediate by everyday
experience (what may not be immediate is the identification of a normal mode).

Lastly, let us compare the states that are detected. Many existing criteria
can witness entanglement of Gaussian states (like the two-mode squeezed state
and its limiting case, the EPR state). Our DEW misses these states, since
their Wigner function is positive. Conversely, there are states whose
entanglement is detected by our DEW and is missed by all the previous ones. A
concrete example is given in Appendix B, where we construct the class of
states $\\{\ket{\Psi_{\mathbf{n}}}\\}$. These states are always entangled, but
all are missed by [10, 11, 20], and some by [12]. In particular, the maximal
eigenstates of $Q_{K}$ (with or without truncation) are of this type. More
generally, our DEW is not a subset of any member of the family of uncertainty-
based entanglement witnesses defined by Shchukin and Vogel [13] and Nha and
Zubairy [14], which include [10, 12] as special cases. Indeed, those families
can only detect entangled states with negative partial transpose. By contrast,
for $\theta=\pi/4$ we can prove that our DEW is a nondecomposable entanglement
witness (proof in Appendix C): this implies that it can detect at least one
bound entangled state, although we have not been able to find an explicit
example.

##  VI Conclusion

We have introduced a dynamic-based entanglement witness for two harmonic
oscillators. It consists of certifying the quantumness of a normal mode using
the Tsirelon protocol [8, 9]: the entanglement of the physical oscillators can
then be inferred, without having any information about the other normal mode
(obviously, having also some information about it can only tighten the lower
bounds that we have obtained).

Our criterion detects a different set of states than those captured by
previous ones. Also, it does not rely on other features of quantum theory
(e.g. some form of uncertainty relations): it only assumes the form of the
dynamics. As only straightforward coordinate measurements are used, and false
positives from classical theory are excluded, our criterion is useful for
objects which are too massive to be fully tomographed. It is timely, as recent
advances in optomechanics have allowed for quantum control of objects with
masses in the mesoscopic and macroscopic scales [6]. The ability to generate
and detect quantum effects in such systems has proven useful in technologies
that aim to exploit the effects of quantum mechanics, in particular, quantum
sensing [21, 22] and metrology [23]. Moreover, it is also important in
experimental studies of fundamental physics, in tests of collapse theories
[24] and quantum-classical transitions [25].

In order to focus on the essentials of the idea, in this paper we have kept to
the simplest form of dynamics, that of two coupled harmonic oscillators. The
underlying Tsirelson protocol can be extended to a wide class of Hamiltonians,
possibly with an additional energy constraint [26]. Dynamic-based entanglement
certification can also be extended to dissipative dynamics: leaving
quantitative estimates for future work, it is clear that one will still be
able to certify entanglement provided the relaxation time is long enough.

## Acknowledgments

We acknowledge financial support from the National Research Foundation and the
Ministry of Education, Singapore, under the Research Centres of Excellence
programme. We also thank the National University of Singapore Information
Technology for the use of their high performance computing resources.

## References

  * Einstein  _et al._ [1935] A. Einstein, B. Podolsky, and N. Rosen, Can quantum-mechanical description of physical reality be considered complete?, Phys. Rev. 47, 777 (1935).
  * Schrödinger [1935] E. Schrödinger, Die gegenwärtige situation in der quantenmechanik, Naturwissenschaften 23, 823 (1935).
  * Horodecki  _et al._ [2009] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Rev. Mod. Phys. 81, 865 (2009).
  * Plenio and Virmani [2014] M. B. Plenio and S. S. Virmani, An introduction to entanglement theory, in  _Quantum Information and Coherence_, edited by E. Andersson and P. Öhberg (Springer International Publishing, Cham, 2014) pp. 173–209.
  * Adesso and Illuminati [2007] G. Adesso and F. Illuminati, Entanglement in continuous-variable systems: recent advances and current perspectives, Journal of Physics A: Mathematical and Theoretical 40, 7821 (2007).
  * Sarma  _et al._ [2021] A. K. Sarma, S. Chakraborty, and S. Kalita, Continuous variable quantum entanglement in optomechanical systems: A short review, AVS Quantum Science 3, 015901 (2021), https://doi.org/10.1116/5.0022349 .
  * Note [1] The name “dynamical entanglement witness” would have been more elegant, but could have been confused with a “witness of dynamical entanglement”, which is a notion of entanglement for channels [27]. Here we are concerned with entanglement of states.
  * Tsirelson [2006] B. Tsirelson, How often is the coordinate of a harmonic oscillator positive?, arXiv preprint quant-ph/0611147 (2006).
  * Zaw  _et al._ [2022] L. H. Zaw, C. C. Aw, Z. Lasmar, and V. Scarani, Detecting quantumness in uniform precessions, Phys. Rev. A 106, 032222 (2022).
  * Duan  _et al._ [2000] L.-M. Duan, G. Giedke, J. I. Cirac, and P. Zoller, Inseparability criterion for continuous variable systems, Physical Review Letters 84, 2722 (2000).
  * Zhang  _et al._ [2010] C.-J. Zhang, H. Nha, Y.-S. Zhang, and G.-C. Guo, Detection of bound entanglement in continuous-variable systems, Phys. Rev. A 82, 032323 (2010).
  * Hillery and Zubairy [2006] M. Hillery and M. S. Zubairy, Entanglement conditions for two-mode states, Physical review letters 96, 050503 (2006).
  * Shchukin and Vogel [2005] E. Shchukin and W. Vogel, Inseparability criteria for continuous bipartite quantum states, Phys. Rev. Lett. 95, 230502 (2005).
  * Nha and Zubairy [2008] H. Nha and M. S. Zubairy, Uncertainty inequalities as entanglement criteria for negative partial-transpose states, Phys. Rev. Lett. 101, 130402 (2008).
  * Vandenberghe and Boyd [1996] L. Vandenberghe and S. Boyd, Semidefinite programming (1996) pp. 49–95, https://doi.org/10.1137/1038003 .
  * Sperling and Vogel [2009] J. Sperling and W. Vogel, Verifying continuous-variable entanglement in finite spaces, Phys. Rev. A 79, 052313 (2009).
  * Banaszek and Wódkiewicz [1999] K. Banaszek and K. Wódkiewicz, Testing quantum nonlocality in phase space, Phys. Rev. Lett. 82, 2009 (1999).
  * Nha and Carmichael [2004] H. Nha and H. J. Carmichael, Proposed test of quantum nonlocality for continuous variables, Phys. Rev. Lett. 93, 020401 (2004).
  * García-Patrón  _et al._ [2004] R. García-Patrón, J. Fiurášek, N. J. Cerf, J. Wenger, R. Tualle-Brouri, and P. Grangier, Proposal for a loophole-free bell test using homodyne detection, Phys. Rev. Lett. 93, 130409 (2004).
  * Abiuso  _et al._ [2021] P. Abiuso, S. Bäuml, D. Cavalcanti, and A. Acín, Measurement-device-independent entanglement detection for continuous-variable systems, Phys. Rev. Lett. 126, 190502 (2021).
  * de Lépinay  _et al._ [2021] L. M. de Lépinay, C. F. Ockeloen-Korppi, M. J. Woolley, and M. A. Sillanpää, Quantum mechanics-free subsystem with mechanical oscillators, Science 372, 625 (2021).
  * Gut  _et al._ [2020] C. Gut, K. Winkler, J. Hoelscher-Obermaier, S. G. Hofer, R. M. Nia, N. Walk, A. Steffens, J. Eisert, W. Wieczorek, J. A. Slater, M. Aspelmeyer, and K. Hammerer, Stationary optomechanical entanglement between a mechanical oscillator and its measurement apparatus, Phys. Rev. Research 2, 033244 (2020).
  * Kwon  _et al._ [2022] H. Kwon, Y. Lim, L. Jiang, H. Jeong, and C. Oh, Quantum metrological power of continuous-variable quantum networks, Phys. Rev. Lett. 128, 180503 (2022).
  * Zhang  _et al._ [2017] J. Zhang, T. Zhang, and J. Li, Probing spontaneous wave-function collapse with entangled levitating nanospheres, Phys. Rev. A 95, 012141 (2017).
  * Schlosshauer [2007] M. Schlosshauer,  _Decoherence: And the Quantum-To-Classical Transition_ , The Frontiers Collection (Springer, 2007).
  * Zaw and Scarani [2022] L. H. Zaw and V. Scarani, Manuscript in preparation (2022), unpublished.
  * Gour and Scandolo [2020] G. Gour and C. M. Scandolo, Dynamical entanglement, Phys. Rev. Lett. 125, 180505 (2020).
  * Jiang  _et al._ [2013] Z. Jiang, M. D. Lang, and C. M. Caves, Mixing nonclassical pure states in a linear-optical network almost always generates modal entanglement, Physical Review A 88, 044301 (2013).
  * Lewenstein  _et al._ [2000] M. Lewenstein, B. Kraus, J. I. Cirac, and P. Horodecki, Optimization of entanglement witnesses, Phys. Rev. A 62, 052310 (2000).

##  Appendix A Semidefinite programming

In this appendix, we detail the numerical procedures performed to obtain the
results reported in the main text using semidefinite programming (SDP).

###  A.1 Minimizing entanglement for observed $P_{K}$ and given $\theta$

The entanglement of a state $\rho$ in the $\\{a_{1},a_{2}\\}$ basis can be
quantified with the logarithmic negativity
$S_{N}(\rho)=\log\tr\\!\absolutevalue{\rho^{\Gamma_{2}}}$. By bounding the
logarithmic negativity from below, we can deduce the minimum entanglement of a
state that violates the classical bound of the precession protocol.

To cast the minimization of $S_{N}(\rho)$ as an SDP, we introduce a Hermitian
basis of operators $\\{B_{j}\\}_{j=0}^{\pqty{\mathbf{n}+1}^{2}-1}$ for both
$a_{\sigma}$ systems, for some truncation $0\leq n_{\sigma}\leq\mathbf{n}$,
with the convention that $B_{0}\propto\mathbb{1}$ and
$\tr(B_{j}B_{k})=\delta_{j,k}$. In our numerical program, we used as $B_{j}$
the normalized form of the generalized Gell-Mann matrices. The operator basis
over the truncation of $\\{a_{+},a_{-}\\}$ is $\\{B_{j}\otimes
B_{k}\\}_{j,k=0}^{\pqty{\mathbf{n}+1}^{2}-1}$. Excluding $B_{0}\otimes B_{0}$,
we collect them into a vector as

$\displaystyle\overrightarrow{B\\!B}\equiv\Large($
$\displaystyle\phantom{B_{1}\otimes B_{0},}B_{0}\otimes
B_{1},\;\dots,\;B_{0}\otimes B_{\pqty{\mathbf{n}+1}^{2}-1},$ $\displaystyle
B_{1}\otimes B_{0},B_{1}\otimes B_{1},\;\dots,\;B_{1}\otimes
B_{\pqty{\mathbf{n}+1}^{2}-1},$ $\displaystyle\phantom{B_{1}\otimes
B_{0},B_{1}\otimes B_{0},}\;\dots,\;B_{\pqty{\mathbf{n}+1}^{2}-1}\otimes
B_{\pqty{\mathbf{n}+1}^{2}-1}\Large).$

Then $\rho={B_{0}\otimes
B_{0}}/{\tr(B_{0})^{2}}+\vec{x}\cdot\overrightarrow{B\\!B}$ and

$\displaystyle\rho^{\Gamma_{2}}=\frac{\pqty{B_{0}\otimes
B_{0}}^{\Gamma_{2}}}{\tr(B_{0})^{2}}+\vec{x}\cdot\overrightarrow{B\\!B}^{\Gamma_{2}},$
(11)

where $\vec{x}=\tr(\rho\overrightarrow{B\\!B})$, and
$\overrightarrow{B\\!B}^{\Gamma_{2}}$ is found by representing each operator
in $\overrightarrow{B\\!B}$ as a matrix in the $\\{a_{1},a_{2}\\}$ basis (here
is where $\theta$ is used) and taking the partial transpose in the usual way.
Denoting its positive eigenvalues as $\lambda_{+k}$ and its negative
eigenvalues as $-\lambda_{-k}$, so that $\lambda_{\pm k}\geq 0$, we have the
spectral decomposition

$\displaystyle\rho^{\Gamma_{2}}$
$\displaystyle=\underbrace{\sum_{\lambda_{+k}}\lambda_{+k}\outerproduct{\lambda_{+k}}{\lambda_{+k}}}_{\equiv\varrho_{+}}-\underbrace{\sum_{\lambda_{-k}}\lambda_{-k}\outerproduct{\lambda_{-k}}{\lambda_{-k}}}_{\equiv\varrho_{-}}$
(12)
$\displaystyle=\varrho_{+}-\varrho_{-},\qquad\text{where}\;\varrho_{\pm}\succcurlyeq
0.$

Considering that the restriction $0\leq n_{+},n_{-}\leq\mathbf{n}$ allows for
states in the $\\{a_{1},a_{2}\\}$ basis with $0\leq n_{1},n_{2}\leq
2\mathbf{n}$, the Hermitian operator basis for each $a_{i}$ is
$\\{A_{j}\\}_{j=0}^{\pqty{2\mathbf{n}+1}^{2}-1}$. The full space of
$\varrho_{\pm}$ is spanned by $\\{A_{j}\otimes
A_{k}\\}_{j,k=0}^{\pqty{2\mathbf{n}+1}^{2}-1}$, and we similarly define
$\overrightarrow{A\\!A}\equiv(A_{0}\otimes
A_{1},\dots,A_{\pqty{2\mathbf{n}+1}^{2}-1}\otimes
A_{\pqty{2\mathbf{n}+1}^{2}-1})$. Then,

$\displaystyle\varrho_{+}$ $\displaystyle=z\frac{A_{0}\otimes
A_{0}}{\tr(A_{0})^{2}}+\vec{y}_{+}\cdot\overrightarrow{A\\!A}$ (13)
$\displaystyle\varrho_{-}$ $\displaystyle=(z-1)\frac{A_{0}\otimes
A_{0}}{\tr(A_{0})^{2}}+\vec{y}_{-}\cdot\overrightarrow{A\\!A},$

where $\vec{y}_{\pm}=\tr(\varrho_{\pm}\overrightarrow{A\\!A})$,
$z=\tr(\varrho_{+})$, and $\tr(\varrho_{-})$ is fixed by
$\tr(\varrho_{+}-\varrho_{-})=1$. This means that the logarithmic negativity
of $\rho$ in the $\\{a_{1},a_{2}\\}$ basis is

$S_{N}(\rho)=\log\tr\\!\absolutevalue{\varrho_{+}-\varrho_{-}}=\log(2z-1).$
(14)

Finally, defining $\vec{q}\equiv\tr(Q_{K}\overrightarrow{B\\!B})$, we have

$P_{K}=\tr(\rho Q_{K})=\frac{1}{2}+\vec{x}\cdot\vec{q}.$ (15)

#### SDP Procedure

Define $A_{j}$ and $B_{j}$ according to the preceding discussion. Then,

$\displaystyle\min_{\begin{array}[]{rcl}\scriptstyle\vec{x}&\scriptstyle\in&\scriptstyle\mathbb{R}^{\pqty{\mathbf{n}+1}^{4}-1};\\\
\scriptstyle\vec{y}_{+},\vec{y}_{-}&\scriptstyle\in&\scriptstyle\mathbb{R}^{\pqty{2\mathbf{n}+1}^{4}-1};\\\
\scriptstyle z&\scriptstyle\in&\scriptstyle\mathbb{R}\end{array}}z$ (16)
subject to $\displaystyle\;\;\;\;P_{K}=\tfrac{1}{2}+\vec{x}\cdot\vec{q}$
$\displaystyle\;\;\;\;\tfrac{\pqty{B_{0}\otimes
B_{0}}^{\Gamma_{2}}}{\tr(B_{0})^{2}}+\vec{x}\cdot\overrightarrow{B\\!B}^{\Gamma_{2}}=\tfrac{A_{0}\otimes
A_{0}}{\tr(A_{0})^{2}}+\pqty{\vec{y}_{+}-\vec{y}_{-}}\cdot\overrightarrow{A\\!A},$
$\displaystyle\;\;\;\;\spmqty{undef}\succcurlyeq 0.$

With the minimum $z$, $S_{N}(\rho)=\log(2z-1)$. For the results in the main
text, we used $K=3$ and $\mathbf{n}=11$. The choice of $\mathbf{n}$ is
justified in the following appendix.

###  A.2 Choice of truncation $\mathbf{n}$

We consider how our choice of truncation $\mathbf{n}$ affects the entanglement
of states that maximally violate the classical bound. For each $\mathbf{n}$,
we perform an SDP optimization as laid out in the preceding appendix with
$P_{3}=\mathbf{P}_{3}$, where $\mathbf{P}_{3}$ is the largest eigenvalue of
$Q_{3}$ restricted to the truncated subspace $0\leq
n_{+},n_{-}\leq\mathbf{n}$.

Figure 2: Log negativity of maximally-violating states against truncation
$\mathbf{n}$ for various values of $\theta$, found with SDP optimization. In
all cases, the minimum $S_{N}$ occurs at $\mathbf{n}=11$.

In Fig. 2, we plot the logarithmic negativity against the truncation
$\mathbf{n}$ for various values of $\theta$. We find that the log negativity
increases with $\mathbf{n}$ in steps of $\Delta\mathbf{n}=6$. This reflects
how $\mathbf{P}_{K}$ increases in steps of $\Delta\mathbf{n}=2K$, which is due
to the symmetry of the protocol and the relationship between the odd and even
number states of $\ket{\mathbf{P}_{K}}$ [9]. Importantly, the minimum value of
$S_{N}$ is found for $\mathbf{n}=11$ for all values of $\theta$. As such,
choosing $\mathbf{n}=11$ for the numerical results in our main text provides
us with a conservative estimate of the amount of entanglement witnessed for a
given violation $P_{3}$.

##  Appendix B Class of states missed by uncertainty-based entanglement
witnesses

For $K\geq 3$ odd, consider the family of states

$\displaystyle\ket{\Psi_{\mathbf{n}}}$ $\displaystyle=$
$\displaystyle\sum_{j=0}^{\mathbf{n}K}\ket{\Psi^{j}}_{1}\otimes\ket{j}_{2}$
(17) $\displaystyle=$
$\displaystyle\big{(}\sum_{n=0}^{\mathbf{n}}\psi_{n}\ket{n}_{+}\big{)}\otimes\ket{0}_{-}\equiv\ket{\psi_{\mathbf{n}}}_{+}\otimes\ket{0}_{-},$

where $\mathbf{n}$ can go to infinity, $\absolutevalue{\psi_{0}}\neq 1$, and
$a_{+}=\cos\theta\,a_{1}+\sin\theta\,a_{2}$ with $\theta\in(0,\pi/2)$ defined
like in Section IV, and where

As these states are separable in the $\\{a_{+},a_{-}\\}$ basis, and do not
take the form of the states given in Jiang  _et al._ [28], they are entangled
in the $\\{a_{1},a_{2}\\}$ basis.

For these states, we have

$\displaystyle\expectationvalue{a_{1}}{\Psi_{\mathbf{n}}}$
$\displaystyle=\expectationvalue{a_{1}^{2}}{\Psi_{\mathbf{n}}}=0,$
$\displaystyle\expectationvalue{a_{2}}{\Psi_{\mathbf{n}}}$
$\displaystyle=\expectationvalue{a_{2}^{2}}{\Psi_{\mathbf{n}}}=0,$
$\displaystyle\expectationvalue{a_{1}a_{2}}{\Psi_{\mathbf{n}}}$
$\displaystyle=0,$
$\displaystyle\expectationvalue{a_{1}^{\dagger}a_{1}}{\Psi_{\mathbf{n}}}$
$\displaystyle=\pqty{\cos\theta}^{2}\underbrace{\expectationvalue{a_{+}^{\dagger}a_{+}}{\psi_{\mathbf{n}}}}_{\equiv\expectationvalue{n}_{\psi_{\mathbf{n}}}},$
$\displaystyle\expectationvalue{a_{2}^{\dagger}a_{2}}{\Psi_{\mathbf{n}}}$
$\displaystyle=\pqty{\sin\theta}^{2}\expectationvalue{n}_{\psi_{\mathbf{n}}},$
$\displaystyle\expectationvalue{a_{1}^{\dagger}a_{2}}{\Psi_{\mathbf{n}}}$
$\displaystyle=\sin\theta\cos\theta\expectationvalue{n}_{\psi_{\mathbf{n}}}.$

###  B.1 Application to criterion by Duan _et al._

In the criterion laid out by Duan  _et al._ [10], a state is entangled if for
some $c\in\mathbb{R}$,

$\expectationvalue{\pqty{\Delta u}^{2}}+\expectationvalue{\pqty{\Delta
v}^{2}}<c^{2}+\frac{1}{c^{2}},$ (19)

where

$\displaystyle u$
$\displaystyle=\absolutevalue{c}\sqrt{\frac{m_{1}\omega_{1}}{\hbar}}x_{1}+\frac{1}{c}\sqrt{\frac{m_{2}\omega_{2}}{\hbar}}x_{2}$
$\displaystyle=\frac{\absolutevalue{c}}{\sqrt{2}}\pqty{a_{1}+a_{1}^{\dagger}}+\frac{1}{c\sqrt{2}}\pqty{a_{2}+a_{2}^{\dagger}},$
$\displaystyle v$
$\displaystyle=\frac{\absolutevalue{c}}{\sqrt{m_{1}\omega_{1}\hbar}}p_{1}-\frac{1}{c\sqrt{m_{1}\omega_{1}\hbar}}p_{2}$
$\displaystyle=\frac{\absolutevalue{c}}{\sqrt{2}i}\pqty{a_{1}-a_{1}^{\dagger}}+\frac{1}{c\sqrt{2}i}\pqty{a_{2}-a_{2}^{\dagger}}.$

For the class of states in (17), we have

$\displaystyle\expectationvalue{\pqty{\Delta
u}^{2}}+\expectationvalue{\pqty{\Delta v}^{2}}$
$\displaystyle=\sin(2\theta)\expectationvalue{n}_{\psi_{\mathbf{n}}}\pqty{\frac{c^{2}}{\tan\theta}+\frac{\tan\theta}{c^{2}}+\operatorname{sgn}(c)}+c^{2}+\frac{1}{c^{2}}.$

Note that $\sin\theta,\cos\theta>0$ for the valid range of $\theta$. Using
$\operatorname{sgn}(x)\geq-1$ and $x+1/x\geq 2$ for $x\geq 0$,

$\displaystyle\expectationvalue{\pqty{\Delta
u}^{2}}+\expectationvalue{\pqty{\Delta v}^{2}}$
$\displaystyle\geq\sin(2\theta)\expectationvalue{n}_{\psi_{\mathbf{n}}}+c^{2}+\frac{1}{c^{2}}$
(20) $\displaystyle>c^{2}+\frac{1}{c^{2}},$

since $\absolutevalue{\psi_{0}}^{2}\neq
1\implies\expectationvalue{n}_{\psi_{\mathbf{n}}}>0$.

Therefore, all states in this family fails the Duan criterion for any choice
of $c$, while some—in particular, the maximally-violating state with or
without truncation—are detected by our criterion.

###  B.2 Application to criterion by Zhang et al.

For states with $\expectationvalue{x_{j}}=\expectationvalue{p_{j}}=0$, as is
the case for the family of states given in (17), the entanglement criteria
given in (6) of Zhang  _et al._ [11] becomes

$\displaystyle
4\expectationvalue{a_{1}^{\dagger}a_{1}}\expectationvalue{a_{2}^{\dagger}a_{2}}$
$\displaystyle<4\absolutevalue{\expectationvalue{a_{1}a_{2}}}^{2},$ (21)
$\displaystyle
4\expectationvalue{a_{1}^{\dagger}a_{1}}\expectationvalue{a_{2}^{\dagger}a_{2}}$
$\displaystyle<4\absolutevalue{\expectationvalue{a_{1}^{\dagger}a_{2}}}^{2},$

where the state is entangled if any of the two inequalities are met. For our
family of states, we have

$\begin{gathered}4\expectationvalue{a_{1}^{\dagger}a_{1}}\expectationvalue{a_{2}^{\dagger}a_{2}}=4\absolutevalue{\expectationvalue{a_{1}^{\dagger}a_{2}}}^{2}=\pqty{\sin(2\theta)}^{2}\expectationvalue{n}^{2}_{\psi_{\mathbf{n}}},\\\
4\absolutevalue{\expectationvalue{a_{1}a_{2}}}^{2}=0.\end{gathered}$ (22)

Hence, neither of the inequalities in (21) are satisified, so no member of our
family of states will be found to be entangled by Zhang  _et al._ ’s criteria.

###  B.3 Application to criterion by Hillery and Zubairy

In the criterion by Hillery and Zubairy [12], a state is entangled if

$\absolutevalue{\expectationvalue{a_{1}a_{2}^{\dagger}}}^{2}>\expectationvalue{a_{1}^{\dagger}a_{1}a_{2}^{\dagger}a_{2}},$
(23)

As before, for our class of states in (17),

$\displaystyle\expectationvalue{a_{1}^{\dagger}a_{1}a_{2}^{\dagger}a_{2}}$
$\displaystyle=\pqty{\sin\theta\cos\theta}^{2}\pqty{\expectationvalue{n^{2}}_{\psi_{\mathbf{n}}}-\expectationvalue{n}_{\psi_{\mathbf{n}}}}$
(24)
$\displaystyle\absolutevalue{\expectationvalue{a_{1}a_{2}^{\dagger}}}^{2}$
$\displaystyle=\pqty{\sin\theta\cos\theta}^{2}\expectationvalue{n}_{\psi_{\mathbf{n}}}^{2}.$
(25)

Therefore, any state with
$\expectationvalue{n^{2}}_{\psi_{\mathbf{n}}}-\expectationvalue{n}_{\psi_{\mathbf{n}}}^{2}\geq\expectationvalue{n}_{\psi_{\mathbf{n}}}$
will not be found to be entangled by Hillery and Zubairy’s criterion. In
particular, this includes maximally-violating states of $Q_{K}$ with
truncation $\mathbf{n}$.

###  B.4 Application to criterion by Abiuso et al.

Abiuso  _et al._ [20] introduces a trusted source of coherent states
$\ket{\alpha}_{\alpha}$ and $\ket{\beta}_{\beta}$, in the ancilla modes
$a_{\alpha}$ and $a_{\beta}$ respectively. This source produces coherent
states following the distribution

$P(\alpha)=\frac{1}{\pi\sigma^{2}}e^{-{\absolutevalue{\alpha}^{2}}/{\sigma^{2}}},\quad
P(\beta)=\frac{1}{\pi\sigma^{2}}e^{-{\absolutevalue{\beta}^{2}}/{\sigma^{2}}}.$
(26)

Then, a pure state of the total system is of the form
$\ket{\Psi_{\text{all}}}=\ket{\alpha}_{\alpha}\otimes\ket{\psi}_{1,2}\otimes\ket{\beta}_{\beta}$,
where the entanglement of $\ket{\psi}_{1,2}$ over the $\\{a_{1},a_{2}\\}$
basis is of interest.

This criterion can detect in principle any entangled state, if no restriction
is put on the measurements. When restricted to quadratures, the criterion
becomes similar to the one of Duan  _et al._ [10] and can detect Gaussian
entangled states. Let us now prove that it misses the family of states under
consideration here. Given the pair of observables

$\displaystyle U_{\kappa}$ $\displaystyle=$
$\displaystyle\frac{\kappa}{2}\pqty{a_{\alpha}+a_{\alpha}^{\dagger}+a_{1}+a_{1}^{\dagger}}-\frac{1}{2\kappa}\pqty{a_{\beta}+a_{\beta}^{\dagger}+a_{2}+a_{2}^{\dagger}}$
$\displaystyle{}-{}\frac{1}{\sqrt{2}}\pqty{\kappa\real(\alpha)-\frac{\imaginary(\beta)}{\kappa}},$
$\displaystyle V_{\kappa}$ $\displaystyle=$
$\displaystyle\frac{\kappa}{2i}\pqty{a_{\alpha}-a_{\alpha}^{\dagger}-a_{1}+a_{1}^{\dagger}}+\frac{1}{2i\kappa}\pqty{a_{\beta}-a_{\beta}^{\dagger}-a_{2}+a_{2}^{\dagger}}$
$\displaystyle{}-{}\frac{1}{\sqrt{2}}\pqty{\kappa\real(\alpha)+\frac{\imaginary(\beta)}{\kappa}},$

the criterion by Abiuso  _et al._ [20] states that $\ket{\psi}_{1,2}$ is
entangled if

$\expectationvalue{U_{\kappa}^{2}}+\expectationvalue{V_{\kappa}^{2}}<\frac{\kappa^{2}+\kappa^{-2}}{2}\frac{\sigma^{2}}{1+\sigma^{2}}.$
(27)

Here, the expectation value of $U_{\kappa}^{2}$ is over the distribution of
the coherent states

$\expectationvalue{U_{\kappa}^{2}}=\int\differential[2]{\alpha}\int\differential[2]{\beta}P(\alpha)P(\beta)\expectationvalue{U_{\kappa}^{2}}{\Psi_{\text{all}}},$
(28)

with a similar expression for $V_{\kappa}^{2}$. Replacing $\ket{\psi}_{1,2}$
with the states $\ket{\Psi_{\mathbf{n}}}$ as defined in Eq. (17), we have

$\displaystyle\expectationvalue{U_{\kappa}^{2}}{\Psi_{\text{all}}}$
$\displaystyle=\frac{1}{2}\pqty{3-2\sqrt{2}}\pqty{\kappa\real(\alpha)-\frac{\real(\beta)}{\kappa}}^{2}+\frac{1}{2}\pqty{\kappa^{2}+\kappa^{-2}}+\frac{1}{2}\pqty{\cos\theta-\sin\theta}^{2}\expectationvalue{n}_{\psi_{\mathbf{n}}}$
(29) $\displaystyle\expectationvalue{V_{\kappa}^{2}}{\Psi_{\text{all}}}$
$\displaystyle=\frac{1}{2}\pqty{3-2\sqrt{2}}\pqty{\kappa\imaginary(\alpha)+\frac{\imaginary(\beta)}{\kappa}}^{2}+\frac{1}{2}\pqty{\kappa^{2}+\kappa^{-2}}+\frac{1}{2}\pqty{\cos\theta+\sin\theta}^{2}\expectationvalue{n}_{\psi_{\mathbf{n}}}.$
(30)

Therefore,

$\displaystyle\expectationvalue{U^{2}_{\kappa}}+\expectationvalue{V^{2}_{\kappa}}$
$\displaystyle=\frac{\kappa^{2}+\kappa^{-2}}{2}\pqty{\pqty{3-2\sqrt{2}}\sigma^{2}+2}+\expectationvalue{n}_{\psi_{\mathbf{n}}}$
(31)
$\displaystyle\geq\frac{\kappa^{2}+\kappa^{-2}}{2}\frac{\sigma^{2}}{1+\sigma^{2}},$

where we have used that $\expectationvalue{n}_{\psi_{\mathbf{n}}}\geq 0$ and
$\pqty{3-2\sqrt{2}}\sigma^{2}+2\geq\frac{\sigma^{2}}{1+\sigma^{2}}$ for all
$\sigma$. As such, none of these states will be found to be entangled by this
criterion.

##  Appendix C The dynamic-based entanglement witness is an optimal,
nondecomposable entanglement witness

The canonical form of an entanglement witness is an operator $W=W^{\dagger}$
such that $\tr(W\rho)\geq 0$ for all $\rho$ separable and $\tr(W\rho)<0$ for
at least one $\rho$ entangled. We can rewrite our DEW for $\theta=\pi/4$ into
the canonical form of an entanglement witness with

$W=\frac{1}{2}\pqty{1+\frac{1}{K}}\mathbbm{1}-Q_{K},$ (32)

as in this case, any violation of the classical bound $\pqty{1+1/K}/2$
witnesses entanglement. Following this, we shall prove a few properties of
$W$.

###  C.1 $W$ is an optimal entanglement witness

From Theorem 1 of Lewenstein  _et al._ [29], we call $W$ optimal iff for all
$P\succcurlyeq 0$ and $\epsilon>0$, $W^{\prime}=(1+\epsilon)W-\epsilon P$ is
not an entanglement witness. Let us consider the state

$\displaystyle\ket{r}$
$\displaystyle\equiv\ket{\alpha=-r}_{1}\otimes\ket{\alpha=-r}_{2}$ (33)
$\displaystyle={|\alpha=-\sqrt{2}r\rangle}_{+}\otimes\ket{\alpha=0}_{-}.$

The Wigner function of this state in $a_{+}$ is a Gaussian centered around
$x_{+}=-2r/\sqrt{m\hbar\omega}$, $p_{+}=0$. A direct calculation yields

$\expectationvalue{W}{r}=\frac{1}{6}\pqty{1-2\operatorname{erf}(r)+\operatorname{erf}(2r)}\leq\frac{1}{6}\pqty{1-\operatorname{erf}(r)}.$
(34)

Notice that $\expectationvalue{W}{r}$ can be chosen to be arbitrarily small by
choosing $r$ sufficiently large. In particular, by choosing
$r>\operatorname{erf}^{-1}\pqty{1-\frac{6\epsilon}{1+\epsilon}\expectationvalue{P}{r}}$,
we have

$\expectationvalue{W}{r}<\frac{\epsilon}{1+\epsilon}\expectationvalue{P}{r}\implies\expectationvalue{W^{\prime}}{r}<0.$
(35)

Since $\ket{r}$ is separable, $W^{\prime}$ is not an entanglement witness. As
an appropriate $r$ can be chosen for any $P\succcurlyeq 0$ and $\epsilon>0$,
$W$ is an optimal entanglement witness.

###  C.2 $W$ is a nondecomposable entanglement witness

We can now proceed to show that $W$ is not decomposable—that is, it _cannot_
be written as $W=cP+(1-c)Q^{\Gamma_{2}}$ for some $c\in[0,1]$ and
$P,Q\succcurlyeq 0$. Theorem 2 of Lewenstein  _et al._ [29] states that if a
decomposable entanglement witness is optimal, then it can be written as
$W=Q^{\Gamma_{2}}$, where $Q\succcurlyeq 0$ contains no product vector in its
range.

As we have shown that $W$ is optimal, this means that
$W^{\Gamma_{2}}\succcurlyeq 0$ if $W$ is decomposable. Here, we consider the
smallest eigenvalue of
$W^{\prime\prime}=\Pi_{\lceil\frac{K}{2}\rceil}W^{\Gamma_{2}}\Pi_{\lceil\frac{K}{2}\rceil}$,
where

$\Pi_{\lceil\frac{K}{2}\rceil}=\pqty{\sum_{n=0}^{\lceil\frac{K}{2}\rceil}\outerproduct{n}{n}_{1}}\otimes\pqty{\sum_{m=0}^{\lceil\frac{K}{2}\rceil}\outerproduct{m}{m}_{2}}$
(36)

is a projection of $W^{\Gamma_{2}}$ onto a truncated space. The matrix
elements of $Q_{K}$ are known analytically, so $W^{\prime\prime}$ can be
constructed explicitly. By diagonalizing $W^{\prime\prime}$ with standard
techniques, the smallest eigenvalue of $W^{\prime\prime}$ is found to be
negative. Hence, $W^{\prime\prime}{\not{\\!\succcurlyeq\;}}0\implies
W^{\Gamma_{2}}{\not{\\!\succcurlyeq\;}}0$.

Therefore, since $W$ is an optimal entanglement witness and $W\neq
Q^{\Gamma_{2}}$ for some $Q\succcurlyeq 0$, $W$ is not decomposable.

###  C.3 $W$ can detect positive partial transposed entangled states

Theorem 3 of Lewenstein  _et al._ [29] states that an entanglement witness is
nondecomposable iff it detects positive partial transposed entangled states.
Having established in the preceeding discussion that $W$ is both optimal and
nondecomposable, $W$ detects entangled positive partial transpose states.