# Massive Trajectory Matching and Construction from Aerial Videos based on
Frame-by-Frame Vehicle Detections

Ruyi Fenga, Zhibin Lia and Changyan Fana CONTACT Zhibin Li Email:
<EMAIL_ADDRESS>aSchool of Transportation, Southeast University, Nanjing,
China

###### Abstract

Vehicle trajectory data provides critical information for traffic flow
modeling and analysis. Unmanned aerial vehicles (UAV) is an emerging
technology for traffic data collection because of its flexibility and
diversity on spatial and temporal coverage. Vehicle trajectories are
constructed from frame-by-frame detections. The increase of vehicle counts
makes multiple-target matching more challenging. Errors are caused by pixel
jitter, vehicle shadows, road marks as well as some missing detections. This
research proposes a novel framework for construction of massive vehicle
trajectories from aerial videos by matching vehicle detections based on
traffic flow dynamic features. The You Look Only Once (YOLO) v4 is used for
vehicle detection in UAV videos based on Convolution Neural Network (CNN).
Trajectory construction is proposed in detected bounding boxes with trajectory
identification, integrity enhancement, and coordinate transformation from
image coordinates to the Frenet coordinates. The raw trajectory obtained is
then denoised by the ensemble empirical mode decomposition (EEMD). Our
framework is tested on two aerial videos taken by a UAV on city expressway
covering congested and free-flow traffic conditions. The results show that the
proposed framework achieves a Recall of 93.00% and 86.69%, and a Precision of
98.86% and 98.83% for vehicle trajectories in the free-flow and congested
traffic conditions.The trajectory processing speed is about 30s per track.

###### keywords:

Vehicle trajectory; UAV; Vehicle detection; YOLOv4; Trajectory construction

##  1 Introduction

Vehicle trajectories provide vital data support for traffic flow studies. It
contains both free-flow and congested traffic states. From the trajectory map,
we can calculate macroscopic traffic flow parameters such as average speed,
density, and volume, which are the outputs of the widely used inductive loop
detectors. The trajectory map also provides an intuitive view of traffic
phenomena such as kinematic wave and traffic breakdown. Moreover, it includes
traffic parameters at the microscopic scope such as speed; time headway and
space headway of individual vehicles. Previously, numerous studies have
utilized vehicle trajectory data for traffic flow model calibration ( e.g.
Tordeux et al., 2010; Zheng et al., 2013; Talebpour et al., 2015), traffic
feature/phenomena exploration, (e.g. Chen et al., 2012; Coifman, 2015), car-
following and lane-changing behavior analysis, (e.g. Siqueira et al., 2016;
Herrera and Bayen, 2010), and driving strategy development, (e.g. Li et al.,
2012; Rhoades et al., 2016; Aghabayk et al., 2013).

The most famous vehicle trajectory dataset is the Next Generation Simulation
(NGSIM) launched by the Federal Highway Administration. The project installs
multiple fixed cameras on the top of a nearby building to collect traffic
videos. It provides vehicle trajectory data of four road segments, containing
information such as instantaneous vehicle velocity, acceleration, position
coordinates, vehicle length, and vehicle type. This dataset has been widely
used since its release (He, 2017). However, the dataset has some limitations
such as fixed road segment, insufficient coverage range, limited traffic flow
condition, limited vehicle component, as well as erroneous speed and
acceleration information (Montanino and Punzo, 2013). Such limitations are
associated with the fixed locations of cameras at inclined shooting angles as
well as the trajectory extraction methods applied on those videos (e.g.
Montanino and Punzo 2013; Punzo et al. 2011; Zheng and Washington 2012). A
multi-step reconstruction algorithm was provided for verifying abnormal
trajectories in this dataset, e.g.(Montanino and Punzo, 2015).

In recent years, unmanned aerial vehicles (UAV) technology has emerged in
traffic data collection because of its high flexibility, low cost and broad
covering range. Previous studies have proposed approaches to extract vehicle
trajectories from UAV videos. The critical challenges include small-sized
objects, indistinguishable features, the influence of shadows and simultaneous
operation for multiple vehicles. e.g.(Barmpounakis et al., 2016; Feng et al.,
2020). Some studies apply feature-based tracking algorithms to track vehicles
in consecutive frames based on the maximum similarity criterion. They usually
need the initial positions of the to-be-tracked targets in the region of
interest (ROI) which locates in upstream. For example, Xiao et al. (2007)
apply a Sniff object-tracking algorithm with the optical flow for vehicle
position initialization. Saleemi and Shah (2013) propose a multi-frame many-
many point correspondence-tracking algorithm using the Harris corner detector.
Chen et al. (2019) combines different feature descriptors to form a multi-view
learning and sparse representation-based target tracker. However, the
performance of the tracking-based trajectory extraction algorithms is very
sensitive to the initialization of vehicle position. In other words, the
trajectory will be entirely lost if it is not detected in the ROI.

To avoid such issue, studies construct the vehicle trajectories by detecting
vehicles in the entire spatial-temporal region of the videos. For example, Cao
et al. (2011) use a support vector machine (SVM) classifier with the bLPS-HOG
features for vehicle detections. They then applied a spatiotemporal
appearance-related similarity method (STARS) for target correlation. Apeltauer
et al. (2015) use the AdaBoost classifier to learn the multi-scale block local
binary patterns features for detection and tracked vehicles with the single-
particle filter. Azevedo et al. (2014) adopt the background subtraction for
vehicle detection, and the k-shortest disjoints paths algorithm for data
correlation. The performance is highly affected by the vehicle detection
accuracy. Specifically, if the loss rate of detections exceeds some
thresholds, vehicles could be mismatched in different frames and erroneous
trajectories are produced. Thus, how to accurately match vehicles from massive
frame-by-frame detection records to construct true trajectories is rather
challenging.

This research aims at proposing an accurate vehicle matching and trajectory
reconstruction framework using frame-by-frame vehicle detection data. The
traffic dynamic features in massive vehicle flow are applied to generate
constraints for trajectories matching and construction. The performance of our
proposed method is tested on two aerial videos with different target
quantities and traffic conditions. We also estimate the processing speed of
our trajectory construction method.

##  2 METHODS

###  2.1 Overall Framework

Our framework contains four modules, as shown in Figure 1, which are image
stabilization, vehicle detection, trajectory construction, and trajectory
denoising. In the first module, shaking in UAV videos is eliminated by Scale
Invariant Feature Transform (SIFT) operator matching. YOLOv4 algorithm for
vehicle detection is applied for position obtaining as the foundation of
trajectory construction. Then a data correlation algorithm based on traffic
flow theory combined with regional traffic dynamics is proposed to construct
vehicle trajectories. In the end, noise in rough vehicle trajectories is
separated and removed by the EEMD-based denoising algorithm. The construction
details of each section are described below.

Figure 1: Framework flow chart of our algorithm.

###  2.2 Image Stabilization

UAV may shake due to the interference of airflow and wind (e.g.Li et al.,
2019).The external disturbances could cause viewing angle changes in the
image. Such movement will overlap with the target movement in the video and
cause an error. In image stabilization, a feature point selection method using
the scale-invariant feature as decision method, SIFT (Lowe, 2004) is applied.
SIFT is an image local feature description operator based on scale-space that
maintains invariance to image scaling, rotation, and even affine
transformation. The key points that SIFT looks for are some very prominent
”stable” feature points that will not change due to factors such as
illumination, affine transformation and noise, such as corner points, edge
points, bright spots in dark areas, and dark points in bright areas.

The selection of feature point pairs uses the joint restriction of Euclidean
distance and location range. The Euclidean distance is to ensure the
similarity of the two feature points to the internal scale-invariant vector.
And the location range is to exclude the mismatch caused by different objects
had similar feature vectors based on no violent motion of UAV. Based on these
two constraints, all video frames will be affine transformed by the good-key-
point pairs found by the SIFT operator and aligned with the first frame. For
stronger fault tolerance, the multiple good-key-point pairs found will all be
used for matching operations. The calculated affine transformation matrix
searched by the minimum sum of errors for all good-key-point.

###  2.3 CNN-based Vehicle Detection

YOLO series (Redmon et al., 2016; Bochkovskiy et al., 2020) and R-CNN (Fan et
al., 2016) series are two well-known masterpieces of detection in the field of
CNN (Benjdira et al., 2019).Some scholars have conducted an experimental
comparison between respective outstanding works of the two series, YOLOv3 and
Faster-RCNN. This research shows more suitability for vehicle detection in
aerial videos of YOLOv3 and achieves higher detection accuracy and faster
processing speed than R-CNN (Benjdira et al., 2019). The current YOLO has been
upgraded to version 4 (Bochkovskiy et al., 2020), which fully integrates
computing skills and has a significant improvement in speed and edge accuracy.
Note that the latest updated version, YOLOv5 (Jocher et al., 2020), emphasizes
its lightweight and training speed, and its accuracy can almost reach the same
as YOLOv4. As the detection accuracy is concerned most, YOLOv4 is used in this
research.

In YOLOv4, the convolutional layer and sampling layer segment operate the
image alternating continuously (Figure 2). The input image first goes to
convolutional layers to get samples of features. And then, 1×1 or 3×3 layers
resetting is applied for quicker calculating speed. The shortcut layers are
arranged among convolutional layers at specific intervals to divide the CNN
into dozens of pieces. The dividing operation is necessary for controlling the
propagation of gradients as well as avoiding the problems of gradient
diffusion and gradient exploding. At last, the yolo layers perform
classification and position prediction of the targets. The entire calculation
process and target characteristics will be continuously summarized and updated
as a weights file. This weights file is the database of target detection.

Figure 2: YOLOv4 natural work schematic diagram.

The model performs instable in dark vehicle and large vehicle detection
because of their distinguished features and small training samples. To solve
this problem, we propose a merged training strategy to improve detection
performance (see Figure 3). We first train the YOLOv4 model with all vehicle
samples in the training database. In the second step, we further enhance our
method by separately training two additional neural network models for dark
vehicles and large vehicles. After that, the vehicle detections from the three
models are merged to form the raw detection pool.

Figure 3: The merged basic and enhanced training strategy.

The raw data may contain duplicated detection boxes of the same vehicle, which
need to be reduced. Intersection-over-union (IoU) is used in the duplication
reduction method, the calculation of IoU is shown in Figure 4, and the
recognition mechanism is expressed as (Eq. 1, 2,3):

Figure 4: The calculation of IoU.

$max(IoU_{obj1},IoU_{obj2})\leq area_{overlap}\;$ (1) $area_{overlap}\leq
min(IoU_{obj_{1}},IoU_{obj_{2}})\;$ (2) $IoU_{obj1}\leq
area_{overlap}<IoU_{obj2}\;$ (3)

where $IoU_{obj1},IoU_{obj2}$ indicate the IoU value of two overlapped object
obj1 and obj2, $area_{overlap}$ is the maximum acceptable percentage of target
bounding box overlapping. In Eq. 1-3, only one item will be established for
each judgment. When it suits Eq. 1 or Eq. 2, the higher confidence object will
be reserved.And if Eq. 3 holds, the lower overlapped obj1 will be choose.

After completing the repeated screening, the location of the vehicle at each
moment is obtained. At this time, we only have a scattered distribution of
points without any links. In the trajectory construction part, these points
will be connected in series to construct the trajectory.

###  2.4 Trajectory Construction

A large number of vehicle targets are contained in every frame. Errors are
also contained which could be caused by pixel jitter, vehicle shadows, road
marks as well as some missing detections. Such errors make it more difficult
to matching construction in massive vehicle candidates. Directly linking these
detection frames with errors result in irreparable effects on traffic data.
Existing association methods are processed in a mixture of these right and
wrong candidate boxes and result in a high error rate (Azevedo et al., 2014;
Luetteke et al., 2012).

It is crucial to choose the only true position of one vehicle among the
candidates mixed with interference target. And in an individual vehicle, the
position during travel space must be determined by vehicle dynamics and
surroundings. In other words, by combining vehicle dynamics and surroundings,
we can predict a more reliable position in multiply candidates. In this
section, traffic flow theory is used to deduce the correct way to match
vehicles. In addition, under the constraints of the time section, there will
still be some unreasonable fragments in the trajectory. The unreasonable
fragments should be cut, and the vacant parts should be compensated. In the
following section, trajectory identification integrity enhancement and
coordinate transformation are introduced.

####  2.4.1 Step 1: Trajectory Identification

At the beginning of this step, it is necessary to apply an explanation of
traffic flow theory in differential conditions proposed by Newell (1993) The
original equation is as Eq. 4. Then we do a conventional deformation and it
looks like Eq. 5, where $\Delta Q$ is the inflow (or increase) flow in the
differential area, $\Delta speed$ is the speed increment in this area and
$\Delta k$ indicates the density increment of this area.

$\frac{\partial k(x,t)}{\partial t}+\frac{\partial q(x,t)}{\partial x}=0\;$
(4) $\Delta Q=\Delta speed\times(-\Delta k)$ (5)

In the proposed explanation, we will redefine Eq. 5 from the perspective of
tendency. That is, $\Delta Q$ is defined as the inflow tendency of vehicles,
$-\Delta k$ is redefined as the regional gravity. When either the $\Delta
speed$ or the $-\Delta k$ becomes larger, the stronger tendency of the vehicle
flows into. Besides, we can also redefine regional gravity as Eq. 6. It shows
that the definition of density. For a wider use of this equation, we assume
that there are N (N is a constant number) invisible vehicles in every
differential region. As a result, in the gap area, the larger of space, the
larger of $-\Delta k$. Then, it causes a positive value of regional gravity
($-\Delta k$) and a stronger vehicle inflow tendency.

$-k=-\frac{n}{space}$ (6)

Such an inflow tendency explanation can also be confirmed from the car
following and lane changing models. Whether in the IDM, Gipps model of car-
following(Treiber et al., 2000; Gipps, 1981) or the MOBIL model, CA
model(Treiber and Helbing, 2016; Nagel et al., 1998) of lane-changing, the
driver will continuously determine the vehicle’s running speed and direction
according to the gap size during driving.

In the proposed inflow tendency explanation, two critical factors affect the
inflow trend. One is regional speed, and the other is the density of the
differential area. In this step of trajectory construction, the two factors
will be both used as the limitation for an accurate candidate selection.

We first use speed as a broad constrain. The vehicle position distribution of
the former and the later frames are mainly applied, and the procedure is shown
in Fig. 5.For the detection boxes in the former frame, we conduct a kinematic
prediction and find its possible position in the next frame. In the later
frame, the single object of position will be directly categorized as the same
vehicle. But sometimes there are multiple choices, maybe a shadow detection
box. Thus scoring methods are proposed to determine the most likely one.

Figure 5: Schematic diagram of identification in speed searching step.

Within a certain speed range, it needs to select the driving direction. In the
proposed method, we assume the vehicle will move in the direction of greater
clearance. So, the vehicle distribution around the target in the former frame
is collected and used to calculate the probability of movement in each
direction. As shown in Fig. 6, in the direction decision model, the largest
gap will be standardized as 1, and the link to other vehicles will be
expressed in a ratio. This ratio will be regarded as the potential of each
direction, and the candidate box with a high potential in the later frame will
be selected. Only the angle is considered, and the candidate frame at the
angle between the two inflection points will be scored by interpolation. Due
to the different positions and different influences of surrounding vehicles,
even if the candidate shadow box is very close to the target vehicle, its
direction of arrival should be different from the candidate of the target
vehicle.

Figure 6: Regional vehicle cooperative direction decision model.

It needs to notice that once a frame is not detected, the trajectory will end.
Therefore, to be better compatible with the identification problem, we will
continue to look for the future trajectory of the target vehicle. This search
cannot be unlimited, so $TH$ is proposed to restrict the maximum number of
future frames in the search process. In our study, we conduct multiple tests
and finally set $TH$ to be 5, which has the best performance in our
experiments. If a detection box is not detected in the next five frames, it is
considered a false one and removed from the trajectory data. If the box can be
continuously identified, its position is recorded to form the trajectory track
for the vehicle in the following frames. The procedure repeats until there is
no detection box to be correlated, and the trajectory is considered ended. The
above process is conducted for every detection box in all frames. In this way,
the raw trajectory data is produced.

####  2.4.2 Step 2: Integrity Enhancement

The trajectory obtained in the above steps may have false positives. These
false positive parts are called ghost track hereinafter. The following two
steps are applied for processing these unreasonable fragments.

1) Ghost track fragments are taken out first. The movement of ghost track
fragment does not match the characteristics of driving. It may have an unfixed
moving direction and abnormal changes of speed. Non-linear stimulus car-
following model (Eq. 7 was applied in each trajectory to check whether the
position and speed are valid.$X_{n}(t)$ indicates the position of vehicle n in
time $t$, $\dot{X}_{n}(t)$ is speed and $\ddot{X}_{n+1}$ is acceleration. The
first second of the trajectory is used to initialize the parameter
$\lambda_{1}$ and reaction time $T$. During the verification, the leader
vehicle will update as its location changing. When the leader is not found or
changes, it will be treated as no stimulation until the vehicle in front is
stable. Then the entire trajectory with anomalies and the abnormal trajectory
parts at the beginning and end will be removed.

$\ddot{X}_{n+1}(t+T)=\frac{\lambda_{1}}{[X_{n}(t)-X_{n+1}(t)]^{2}}[\dot{X}_{n}(t)-\dot{X}_{n+1}(t)]\;$
(7)

2) Reasonable trajectories will be connected accordingly. After the above
processing, the trajectory fragments that come to this step are reasonable but
may be broken. These broken pieces will be compensated using quantic
polynomial fitting based on the assumption that rate of acceleration change is
in a uniform change speed.

####  2.4.3 Step 3: Coordinate Transformation

The purpose of the coordinate transformation is to convert the vehicle
position information from the video image coordinates (where the axes are the
border of the video) to the so-called Frenet coordinates (e.g.Mao et al.,
2015; Malik et al., 2003), whose axes are along and vertical to the road lane.
The details of the coordinate transformation algorithm are presented as
follows.

Firstly, road-lane-curve information is needed for establishing axes of Frenet
coordinates. Considering the UAV camera is stationary during data collection,
we manually collect sample points of lane boundaries and use the polyfit
function to fit the lane curves. Each lane center curve is obtained by taking
the average of its upper and lower lane boundaries. The UAV may move slightly
due to the interference of wind or signal instability. As a result, we perform
the lane fitting procedure multiple times and update the coordinates to ensure
their accuracy.

Since we have known the lane information, the current vehicle position under
Frenet coordinates can be calculated as follows (Fig. 7). Where
$\\{x(t),f(x(t))\\},(t=0,1,2...)$ represents vehicle position at time t, and
$\\{s(t),g(s(t))\\},(t=0,1,2,..)$ is the point closest to the current vehicle
position on the lane center curve. $q(t)$ is the vehicle movement along the
road lane, m(q(t)) is the vehicle movement vertical to the road lane.
$\triangle q(t)$ is the movement of q(t) in each frame. Note that if taken
each frame as a microelement, $\triangle q(t)$ can be approximated to the
distance between current and last mapping point $\\{s(t),g(s(t))\\}$. After
the above steps, the trajectory data in the Frenet coordinates is established.

Figure 7: Sketch map of transforming image coordinates to Frenet coordinates.

###  2.5 Trajectory Denoising

The trajectories contain noises with features of high frequency and small
magnitude due to the limitations of integer pixels and background
interference. Although these deviations are composed of white noise, they will
still have an impact on the calculation of instantaneous velocity and
acceleration. The ensemble empirical mode decomposition (EEMD) is applied for
our research (Wu and Huang, 2009). EEMD have a good performance for both
linear and non-linear signal decomposition. Besides, the decomposed multiple
intrinsic mode functions (IMF) in different frequencies can robustly avoid
modal mixing, which indicates either an IMF consisting of signals of widely
disparate scales, or a signal of a similar scale residing in different IMF
components.

In the trajectory decomposition, in order to better determine the frequency
band of each IMF, the additional white noise of finite-amplitude is added to
the IMF. Accurately, after adding white noise, the current observation can be
described as follows (Eq. 8):

$X_{i}(t)=x(t)+w_{i}(t)\;$ (8)

where $i$ is the iteration number, $x(t)$ is the raw trajectory data,
$w_{i}(t)$ is the added white noise, $X_{i}(t)$ is the noise-aided trajectory
data.

We decompose the trajectory through the following rules. We first line maximum
points and minimum points with cubic splines to form the upper and lower
envelopes. Then the local means of the upper and lower envelope are taken as
the next observation $x_{i}(t)$. The IMF is found if $x_{i}(t)$ meets the
following conditions (Eq. 9, 10):

$N_{z}-1\leq N_{e}\leq N_{z}+1\;$ (9) $\frac{f_{max}(t)+f_{min}(t)}{2}=0\;$
(10)

Define the extrema of $x_{i}(t)$ as $N_{z}$, the zero-crossing point as
$N_{e}$, the upper and lower envelope as $f_{max}(t)$ and $f_{min}(t)$.

When IMFs can no longer be found, this searching will end and original
$X_{i}(t)$ can be expressed as Eq. 11:

$X_{i}(t)=\sum_{j=1}^{n}imf_{j}(t)+r_{n}(t)\;$ (11)

where $imf_{j}(t)$ is the IMFs separated from $X_{i}(t)$, andv$r_{n}(t)$ is
the residual.

Appropriate IMFs are determined with an energy-based IMF selection method Chen
et al. (2018). In this research, it shows that the energy of the signal
represents the amount of information stored in it. The energy of signal can be
calculated as Eq. 12. The energy of noise signal only accounts for a small
proportion while the majority of energy is concentrated in the meaningful
signals. A log function is used to distinguish the noise signal energy which
processes a negative result, and the other meaningful signals need to meet the
Eq. 13.

$E_{j}=\frac{1}{num}\sum_{k=1}^{num}[c_{j}(k)]^{2}$ (12) $log_{2}E_{j}>0$ (13)

where $E_{j}$ is the energy of the jth IMF, $c_{j}(k)$ is the point collection
of jth IMF, num is the total number of collection points. The sum of the
selected IMFs is the denoised trajectory of the vehicle. After the above
steps, accurate vehicle trajectories are extracted via our framework.

##  3 EXPERIMENT DESIGN

The proposed framework is evaluated on two aerial videos captured by high
definition camera mounted on a UAV (DJI Mavic professional). Videos are
captured at 24 frames-per-second (fps) and with the resolution of 4096 pixel
×2160 pixel. The two test videos were shot at different times of the urban
expressway in Nanjing, China. Test video #1 is taken at a 280m altitude in a
free flow traffic condition, while test video #2 is taken at a 310m altitude
in a congested condition. In the course of data acquisition, the UAV is set to
be stabilized in the air. The video information is shown in Table 3.

Test Video Information  Video Information Test Video #1 Test Video #2 Traffic
condition Free flow Congested Road geometry Curve Curve Frame rate 24fps 24fps
Resolution $4096\times 2160$ $4096\times 2160$ Length 386 427 Duration 333s
300s Frame 7992 7200 Focal length 23mm 23mm Flying height 280m 310m

The experiment is launched on a workstation with a 2.5 GHz E5-2678v3 dual
processor and a NVIDIA 2080 Ti graphics card. The calculation and processing
are developed on the platform of Visual Studio 2015 and Matlab 2016a.

In terms of model evaluation, we use counting methods to evaluate the
completeness of large-scale trajectories. The confusion matrix and two
goodness-of-fit measures are quoted to evaluate the accuracy of the trajectory
counting results (Table 3 and Eq. 14, 15). And also, for the accuracy of
trajectory details, the research group manually mark the ground truth and
calculate the IoU with the extracted trajectory.

Confusion Matrix  Ground Truth Model Results True False True True positive
(TP) False Positive (FP) False False Negative (FN) True Negative (TN)

$Recall=\frac{TP}{GT}$ (14) $Precesion=\frac{TP}{TP+FP}$ (15)

##  4 RESULTS

###  4.1 Results of Image Stabilization

All frames of UAV videos will be matched with the first frame. 200 strongest
initial feature points are selected and matched in every frame by the SIFT
calculator. Not all matched couples are effectively the same point, so the
Euclidean distance of the feature invariant vectors of the two matching points
is used to test and screen the effective feature point pairs. In order to
better obtain effective feature points, we limit the search range of feature
points to outside the road area. Most of the images in these areas are static
buildings, which can provide effective feature points.

As shown in Figure 8(a), the colored circles represent the acquired feature
points, and the colored lines connect the effective feature point pairs. The
proportion of retained results after the screening is not large, but they are
all exact matches. Registration requires at least three pairs of points to
achieve affine transformation. Now, these dozens of matching points have been
able to ensure the registration of the picture, which will match based on the
principle of minimum error. The stabilized image is shown in Figure 8(b). It
can be seen that the characteristic objects are all in the same position after
the image changes, and the non-pixel portions in border areas are filled with
black.

Figure 8: After image stabilization, each pixel matched and the non-pixel
portions in borders area are filled with black.

###  4.2 Model Training and Parameter Calibration

The basic network and an enhanced network of YOLOv4 are trained for aerial
vehicle detection. According to the creator and some researchers who have
applied YOLOv4, the training settings and training samples almost play a
decisive role in the performance of the network. In our study, 20000 vehicle
samples are applied for the basic training. Besides, 4000 dark vehicle samples
and 1000 large vehicle samples are applied for enhanced training. These
datasets are manually calibrated to ensure that the external edge error is
within three pixels.

Critical parameters of YOLOv4 and the values we choose are listed in Table4.2.
The Batch indicates the number of samples that participated in each iteration.
A higher batch will bring a better accuracy but a lower training speed. The
subdivision is set to relieve the usage of GPU by separating training datasets
into groups. Width and Height are the resized image specifications for
training. Large Width and Height will retain more image information but
consume more memory. Max-batches is the upper limit of training times. The
best value of the Max-batch is to get less loss, but overlarge Max-batch can
easily overfit. The learning rate is to changes the model in each based on
training results. A higher learning rate is required at beginning to get a
quick decay in loss and low down in the back to ensure learning accuracy and
prevent overfitting.

YOLOv4 Training Parameters  Configuration parameters Basic detection Enhanced
detection Dark vehicle detection Large vehicle detection Number of training
samples 20000 4000 1000 Batch 128 128 128 Subdivisions 32 32 32 Width 672 672
672 Height 672 672 672 Max batches 60000 20000 20000 Learning rate
0.001/0.0001 0.001/0.0001 0.001/0.0001

The training parameters are set for a balance of training performance and
capability of the hardware. A high accuracy requires a high batch as 128, and
the limitation of hardware sets the subdivisions and training size as 32 and
672×672, respectively. Max-batches are set to 60000 in the initial training
and 20000 in the enhanced training. The minimum loss is found near 48000
batches in the initial training and 12000 batches in the enhanced training.
The learning rate is set to 0.001 at the beginning of training, and a step
decay policy is applied to decrease the learning rate from 0.001 to 0.0001
after 40000 batches in initial training and 8000 batches in enhanced training.

###  4.3 Results of Vehicle Detections

Vehicles are detected on the test videos using the calibrated model. When
evaluating the integrity of the trajectory, we sample of 500 frames randomly,
for counting and manually labeling the bounding boxes. In order to evaluate
the performance of our training strategy, we calculate the results of each
step in the training procedure, i.e. basic training, enhanced training, and
duplication reduction. The vehicle detection performance is shown in Table4.3.

YOLOv4 Detection Results  Performance Test video #1 Test video #2 Bsc train
Enh train Dup reduc Bsc train Enh training Dup reduc Ground truth 50998 125993
True positive 46520 47698 47535 105091 120046 119681 False negative 6385 2954
3100 20902 5947 6312 False positive 1132 6119 551 7257 24934 1424 Recall
91.22% 93.53% 93.21% 83.41% 95.28% 94.99% (sd=0.0017) (sd=0.0006) (sd=0.0012)
(sd=0.0027) (sd=0.0011) (sd=0.0012) Precision 97.63% 88.63% 98.86% 93.54%
82.80% 98.83% Average IoU 0.8377 0.8135

It shows that the YOLOv4 detection model after parameter selection and
enhancement performs stable and shows the comprehensive result in UAV videos.
Notably, the Recall in test video #1 is 91.22%, which is higher than that in
test video #2 (83.41%). The main reason for the difference is that compared
with video #1 which is under the free-flow condition, video #2 is under the
congested condition and possesses larger vehicle amounts and smaller spacing
which lowers the confidence score in the detection section. And the average
IoU is only calculated by the true positive detection boxes. The two test
videos are at a similar and good level although a vehicle only covers about
1500 pixels. We find that with basic training models, some dark vehicles and
large vehicles are not accurately detected (Fig. 3).

Table 4.3 shows that the Recall in both videos increases apparently to 93.53%
and 95.28% after the enhanced training. The main reason is that the separate
neural network models successfully detect more dark and large vehicles.
However, the Precision shows a significant decrease resulting from the overlap
of vehicle detection boxes from different models (Fig. 3).

After the duplication reduction, as shown in Table 4.3, the Precision
successfully increases from 88.63% to 98.86% in test video #1 and 82.80% to
98.83% in test video #2, respectively. It means the false vehicle detections
are successfully removed from the data. The Recall maintains at a high level,
indicating that the duplication reduction does not wrongly delete true
detections. The final Recall is 93.21% in test video #1 and 94.99% in test
video #2, respectively. It indicates that the detection method of our
framework works well in both congested and free-flow traffic conditions.

In the confusion matrix(Table 3), the two types of errors are false negative
and false positive. False negative indicates the missed detections of true
vehicles. It usually happens when the lost vehicle has a special shape, color
or vehicle type, vehicle type, or the vehicle is obscured by road decorations
such as signs. This problem may be reduced when more training samples are
equipped, and the number of samples of each type is the same. False positive
indicates fake vehicle detections, are caused by the mix-up detection of the
road surface or road markings with vehicles because of confusing features.

Two typical examples are shown in Figure 9. In Figure 9(a1), the road marking
is very similar to a white vehicle so that the detector makes a false positive
detection. The pap similarity of the color histograms is 0.6690 in Figure 9
(a2), which shows that the two pictures are very similar. In Figure 9 (b1),
the vehicle contour is not clear so that the detector did not recognize it and
made a false negative mistake. Their color histograms possess a pap similarity
of 0.3846 in Figure 9 (b2).

Figure 9: Two typical false examples detection errors.

###  4.4 Results of Trajectory Construction

Through trajectory construction, raw vehicle trajectories are extracted from
two test videos. The results are shown in Table 4.4. The research team
manually counts the true count of trajectories, which is 500 and 541 in the
two videos. After the trajectory construction, 484 and 497 trajectories are
correctly obtained, while 4 and 43 trajectories are false positive ones. The
improvement by integrity enhancement on trajectory count is evaluated in Table
4.4 (Trajectory Identification is simplified to T Identification and Integrity
enhancement is simplified to I Enhancement). Results show that Recall
increases from 92.20% to 93.00% and 83.92% to 86.69% in test video #1 and test
video #2, respectively. Mainly, fifteen trajectories are put into integrity
enhancement in test video #2 while only four trajectories are enhanced in test
video #1.

This result shows a good performance in trajectory identification steps based
on traffic flow theory. The majority of trajectories have been matched and
constructed correctly. Besides, a large number of false positive and false
negative in detection steps have been eliminated in the trajectory
construction step.

Trajectory Construction Performances  Performances Test Video #1 Test Video #2
T Identification I Enhancement T Identification I Enhancement Ground Truth 500
541 True positive 461 465 454 469 False negative 39 35 87 72 False positive 4
4 43 43 Recall 92.20% 93.00% 83.92% 86.69% Precision 99.13% 99.15% 91.35%
91.60%

In the two videos, 107 trajectories were missing or incorrectly constructed.
We carefully checked these errors and found that most of the missing
trajectories are caused by a long period of detection failure. The continuum
of detection boxes in the trajectories is insufficient, so our methods are not
able to form the complete trajectories. Only a few missing trajectories are
caused by the detection box in one trajectory identifying wrongly to the
trajectory of another vehicle. The mismatch occurs when the to-be-identified
detection box is lost in the next several frames and too much position
derivation is needed. At this time, the interference of adjacent vehicles will
increase, making it easier to misidentify.

###  4.5 Results of Trajectory Denoising

It is very difficult, or impossible, to obtain the ground truth of the vehicle
trajectories for the validation of the accuracy of constructed trajectories.
As a result, in our study, we estimate the validity of the denoising procedure
by checking the reasonableness of the distribution range of parameters such as
vehicle speed, acceleration, space headway, time headway, and gap.

A typical vehicle trajectory denoising result using EEMD is shown in Figure
10. The contribution of smoothing out the noise can be seen in Figure 10. The
raw trajectory data has large fluctuation, while the denoised data is much
smoother. The acceleration change resulting from the denoising process is
shown in Figure 10 (b). According to Montanino and Punzo (2013), the physical
limitations of vehicle acceleration and deceleration rate are $9m/s^{2}$ and
$-9m/s^{2}$, respectively. And as our know, the comfortable acceleration are
within $\pm 3m/s^{2}$. Our denoised acceleration data is all within a
reasonable range (between the pink dotted lines), and most of them are
distribute in the range of comfortable accelerations (between the purple
dotted lines). The cumulative percentage curve of the acceleration of all
vehicles from our data and from Montanio’s method is shown in Figure 10 (c).
It can be clearly seen that the distributions of the acceleration of two
trajectory data after denoising are consistent.

Figure 10: The vehicle trajectory denoising result.

To further validate the reasonableness of our data, we make a comparison of
the value of key parameters (speed, space headway, time headway, and gap) with
the NGSIM dataset from the southbound traffic on US-101 in Los Angeles, CA,
recorded from 14:50-15:05, on June 15, 2005. The visualization of the
comparison is shown in Figure 11. As seen in Figure 11 (c-d), the space
headway with the highest frequency in our data (the red bar) is around 12m,
while that in NGSIM data (the blue bar) is around 40m. Similarly, in Figure
11(e-h), the time headway with the highest frequency is around 1.6s in both
data, and the gap with the highest frequency is 8m and 40m in our data and
NGSIM data, respectively. The difference in the value of the gap and the space
headway maybe because our test videos partially suffer from congestions. For
the same reason, the speed of our data (Figure 11 (c)) has two peaks around
$30km/h$ and $55km/h$. The speed of NGSIM data has one peak around 40km
because the data were collected from the free-flow traffic state.

Figure 11: Comparison of key traffic parameters with NGSIM.

###  4.6 Evaluation of Processing Speed

The processing speed of our framework for trajectory extraction is calculated.
The majority of time is spent on the offline training of the YOLOv4 algorithm.
The training time is about 70h, which is similar to the previous study
(Benjdira et al., 2019). After the training, the online processing speed is
very fast. The processing time in each step of our framework is shown in Table
4.6. The detection step only costs 0.1 h, suggesting the remarkable advantage
of applying YOLOv4. The trajectory construction costs 3.8 h (26s/track) and
4.0 h (27s/track), including the trajectory identification, integrity
enhancement and coordinate transformation. The denoising step is fast, which
costs only 0.3 h (2s/track).

In summary, the total processing time is 4.2 h and 4.3 h on the test video #1
and #2. The slightly longer processing time in video #2 is caused by the fact
that there are more vehicles in congested traffic conditions. The average
online processing speed is 29s/track and 30s/track, or in other words 0.6
s/frame and 0.7 s/frame, on the two test videos.

Processing Time of Our Framework  Section Processing time Processing speed
Test video #1 Test video #2 Test video #1 Test video #2 Training 70h
13s/vehicle Image stabilization 40 minutes 0.3s/frame Detection 0.1h 0.1h
$6*10^{-3}s/vehicle$ $3*10^{-3}s/vehicle$ T Identification 1.3h 1.5h 9s/track
10s/track I Enhancement 0.5h 0.5h 3s/track 3s/track C transformation 2h 2h
14s/track 14s/track Denosing 0.3h 0.3h 2s/track 2s/track whole framework 4.2h
4.4h 29s/track 30s/track (without training)

##  5 CONCLUSIONS AND DISCUSSION

A novel framework for massive vehicle trajectory matching and construction
from aerial videos is proposed in this research. This framework is based on
the vehicle dynamic constraints derived from Newell’s traffic flow theory. In
the framework, an integrated high-precision vehicle detection by YOLOv4 is
arranged and strengthened. Unique trajectory decision-making and
reconstruction method is applied to detected bounding boxes frame-by-frame.
This step is equipped with trajectory identification, integrity enhancement
and coordinate transformation from image coordinates to the Frenet
coordinates. Then the raw trajectory obtained from the construction section
will be denoised by EEMD. Our framework is tested on two aerial videos taken
by a UAV on a city expressway covering congested and free-flow traffic
conditions. The extracted vehicle trajectories are compared with manual counts
and labeled bounding boxes.

The experimental results showed that the average Recall of vehicle detection
count in the two videos are 93.21% and 94.99%, respectively. And the average
Precision of vehicle detection are 98.86% and 98.83%. respectively. IoU of the
vehicle detection is 0.8377 and 0.8135 respectively. The Recall of trajectory
construction count was 93.00% and 86.69% while the Precision of trajectory
construction is 99.15% and 91.60% in the two videos. The denoising algorithm
performed effectively in eliminating outliners and keeping the traffic
parameters within a reasonable range. The total processing time of our
framework was 4.2 h and 4.4 h on test videos #1 and #2.

Our methodology can help enrich the trajectory data for traffic flow studies.
Based on our video data, we construct the space-time trajectory map of each
lane, as shown in Figure 12, 13. It can be observed that the traffic is free-
flowing in video #1 with a relatively constant speed on each lane. While in
video #2, the formation of congestion as well as kinematic waves is evident.
The critical car-following and lane-changing parameters can be estimated from
the trajectory map for more comprehensive traffic flow modeling and analysis.

Figure 12: Space-time trajectory maps (test video #1). Figure 13: Space-time
trajectory maps (test video #2).

## References

  * Aghabayk et al. (2013) Kayvan Aghabayk, Majid Sarvi, William Young, and Lukas Kautzsch.  A novel methodology for evolutionary calibration of vissim by multi-threading.  In _Australasian Transport Research Forum_ , volume 36, pages 1–15, 2013.
  * Apeltauer et al. (2015) Jiřĺ Apeltauer, Adam Babinec, David Herman, and Tomáš Apeltauer.  Automatic vehicle trajectory extraction for traffic analysis from aerial video data.  _The International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences_ , 40(3):9, 2015.
  * Azevedo et al. (2014) Carlos Lima Azevedo, João L Cardoso, Moshe Ben-Akiva, João P Costeira, and Manuel Marques.  Automatic vehicle trajectory extraction by aerial remote sensing.  _Procedia-Social and Behavioral Sciences_ , 111:849–858, 2014.
  * Barmpounakis et al. (2016) Emmanouil N Barmpounakis, Eleni I Vlahogianni, and John C Golias.  Unmanned aerial aircraft systems for transportation engineering: Current practice and future challenges.  _International Journal of Transportation Science and Technology_ , 5(3):111–122, 2016.
  * Benjdira et al. (2019) Bilel Benjdira, Taha Khursheed, Anis Koubaa, Adel Ammar, and Kais Ouni.  Car detection using unmanned aerial vehicles: Comparison between faster r-cnn and yolov3.  In _2019 1st International Conference on Unmanned Vehicle Systems-Oman (UVS)_ , pages 1–6. IEEE, 2019.
  * Bochkovskiy et al. (2020) Alexey Bochkovskiy, Chien-Yao Wang, and Hong-Yuan Mark Liao.  Yolov4: Optimal speed and accuracy of object detection.  _arXiv preprint arXiv:2004.10934_ , 2020.
  * Cao et al. (2011) Xianbin Cao, Changxia Wu, Jinhe Lan, Pingkun Yan, and Xuelong Li.  Vehicle detection and motion analysis in low-altitude airborne video under urban environment.  _IEEE Transactions on Circuits and Systems for Video Technology_ , 21(10):1522–1533, 2011.
  * Chen et al. (2012) Danjue Chen, Jorge Laval, Zuduo Zheng, and Soyoung Ahn.  A behavioral car-following model that captures traffic oscillations.  _Transportation research part B: methodological_ , 46(6):744–761, 2012.
  * Chen et al. (2018) Xinqiang Chen, Zhibin Li, Yinhai Wang, Jinjun Tang, Wenbo Zhu, Chaojian Shi, and Huafeng Wu.  Anomaly detection and cleaning of highway elevation data from google earth using ensemble empirical mode decomposition.  _Journal of Transportation Engineering, Part A: Systems_ , 144(5):04018015, 2018.
  * Chen et al. (2019) Xinqiang Chen, Shengzheng Wang, Chaojian Shi, Huafeng Wu, Jiansen Zhao, and Junjie Fu.  Robust ship tracking via multi-view learning and sparse representation.  _The Journal of Navigation_ , 72(1):176–192, 2019\.
  * Coifman (2015) Benjamin Coifman.  Empirical flow-density and speed-spacing relationships: Evidence of vehicle length dependency.  _Transportation Research Part B: Methodological_ , 78:54–65, 2015.
  * Fan et al. (2016) Quanfu Fan, Lisa Brown, and John Smith.  A closer look at faster r-cnn for vehicle detection.  In _2016 IEEE intelligent vehicles symposium (IV)_ , pages 124–129. IEEE, 2016.
  * Feng et al. (2020) Ruyi Feng, Changyan Fan, Zhibin Li, and Xinqiang Chen.  Mixed road user trajectory extraction from moving aerial videos based on convolution neural network detection.  _IEEE Access_ , 8:43508–43519, 2020.
  * Gipps (1981) Peter G Gipps.  A behavioural car-following model for computer simulation.  _Transportation Research Part B: Methodological_ , 15(2):105–111, 1981.
  * He (2017) Zhengbing He.  Research based on high-fidelity ngsim vehicle trajectory datasets: A review.  _Research Gate_ , pages 1–33, 2017.
  * Herrera and Bayen (2010) Juan C Herrera and Alexandre M Bayen.  Incorporation of lagrangian measurements in freeway traffic state estimation.  _Transportation Research Part B: Methodological_ , 44(4):460–481, 2010.
  * Jocher et al. (2020) Glenn Jocher et al.  Yolov5.  _Code repository https://github. com/ultralytics/yolov5_ , 2020.
  * Li et al. (2012) Xiaopeng Li, Xin Wang, and Yanfeng Ouyang.  Prediction and field validation of traffic oscillation propagation under nonlinear car-following laws.  _Transportation research part B: methodological_ , 46(3):409–423, 2012.
  * Li et al. (2019) Zhibin Li, Xinqiang Chen, Lei Ling, Huafeng Wu, Wenzhu Zhou, and Cong Qi.  Accurate traffic parameter extraction from aerial videos with multi-dimensional camera movements.  Technical report, 2019.
  * Lowe (2004) David G Lowe.  Distinctive image features from scale-invariant keypoints.  _International journal of computer vision_ , 60(2):91–110, 2004.
  * Luetteke et al. (2012) Felix Luetteke, Xu Zhang, and Jörg Franke.  Implementation of the hungarian method for object tracking on a camera monitored transportation system.  In _ROBOTIK 2012; 7th German Conference on Robotics_ , pages 1–6. VDE, 2012.
  * Malik et al. (2003) J Malik et al.  High-quality vehicle trajectory generation from video data based on vehicle detection and description.  In _Proceedings of the 2003 IEEE International Conference on Intelligent Transportation Systems_ , volume 1, pages 176–182. IEEE, 2003.
  * Mao et al. (2015) Tianlu Mao, Hua Wang, Zhigang Deng, and Zhaoqi Wang.  An efficient lane model for complex traffic simulation.  _Computer Animation and Virtual Worlds_ , 26(3-4):397–403, 2015.
  * Montanino and Punzo (2013) Marcello Montanino and Vincenzo Punzo.  Making ngsim data usable for studies on traffic flow theory: Multistep method for vehicle trajectory reconstruction.  _Transportation Research Record_ , 2390(1):99–111, 2013.
  * Montanino and Punzo (2015) Marcello Montanino and Vincenzo Punzo.  Trajectory data reconstruction and simulation-based validation against macroscopic traffic patterns.  _Transportation Research Part B: Methodological_ , 80:82–106, 2015.
  * Nagel et al. (1998) Kai Nagel, Dietrich E Wolf, Peter Wagner, and Patrice Simon.  Two-lane traffic rules for cellular automata: A systematic approach.  _Physical Review E_ , 58(2):1425, 1998.
  * Newell (1993) Gordon F Newell.  A simplified theory of kinematic waves in highway traffic, part i: General theory.  _Transportation Research Part B: Methodological_ , 27(4):281–287, 1993.
  * Punzo et al. (2011) Vincenzo Punzo, Maria Teresa Borzacchiello, and Biagio Ciuffo.  On the assessment of vehicle trajectory data accuracy and application to the next generation simulation (ngsim) program data.  _Transportation Research Part C: Emerging Technologies_ , 19(6):1243–1262, 2011.
  * Redmon et al. (2016) Joseph Redmon, Santosh Divvala, Ross Girshick, and Ali Farhadi.  You only look once: Unified, real-time object detection.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 779–788, 2016.
  * Rhoades et al. (2016) Christine Rhoades, Xin Wang, and Yanfeng Ouyang.  Calibration of nonlinear car-following laws for traffic oscillation prediction.  _Transportation research part C: emerging technologies_ , 69:328–342, 2016.
  * Saleemi and Shah (2013) Imran Saleemi and Mubarak Shah.  Multiframe many–many point correspondence for vehicle tracking in high density wide area aerial videos.  _International journal of computer vision_ , 104(2):198–219, 2013.
  * Siqueira et al. (2016) Adriano F Siqueira, Carlos JT Peixoto, Chen Wu, and Wei-Liang Qian.  Effect of stochastic transition in the fundamental diagram of traffic flow.  _Transportation Research Part B: Methodological_ , 87:1–13, 2016.
  * Talebpour et al. (2015) Alireza Talebpour, Hani S Mahmassani, and Samer H Hamdar.  Modeling lane-changing behavior in a connected environment: A game theory approach.  _Transportation Research Part C: Emerging Technologies_ , 59:216–232, 2015.
  * Tordeux et al. (2010) Antoine Tordeux, Sylvain Lassarre, and Michel Roussignol.  An adaptive time gap car-following model.  _Transportation research part B: methodological_ , 44(8-9):1115–1131, 2010.
  * Treiber and Helbing (2016) Martin Treiber and Dirk Helbing.  Mobil: General lane-changing model for car-following models.  _Disponıvel Acesso Dezembro_ , 2016.
  * Treiber et al. (2000) Martin Treiber, Ansgar Hennecke, and Dirk Helbing.  Congested traffic states in empirical observations and microscopic simulations.  _Physical review E_ , 62(2):1805, 2000.
  * Wu and Huang (2009) Zhaohua Wu and Norden E Huang.  Ensemble empirical mode decomposition: a noise-assisted data analysis method.  _Advances in adaptive data analysis_ , 1(01):1–41, 2009.
  * Xiao et al. (2007) Jiangjian Xiao, Changjiang Yang, Feng Han, and Hui Cheng.  Vehicle and person tracking in aerial videos.  In _Multimodal Technologies for Perception of Humans_ , pages 203–214. Springer, 2007.
  * Zheng et al. (2013) Jian Zheng, Koji Suzuki, and Motohiro Fujita.  Car-following behavior with instantaneous driver–vehicle reaction delay: A neural-network-based methodology.  _Tramsportation Research Part C Emerging Technologies_ , 36(nov.):339–351, 2013.
  * Zheng and Washington (2012) Zuduo Zheng and Simon Washington.  On selecting an optimal wavelet for detecting singularities in traffic and vehicular data.  _Transportation Research Part C: Emerging Technologies_ , 25:18–33, 2012.


# Non-Gaussian entanglement via splitting of a few thermal quanta

Pradip Laha<EMAIL_ADDRESS>Department of Optics, Palacký University, 17.
listopadu 1192/12, 771 46 Olomouc, Czech Republic  Darren W. Moore
<EMAIL_ADDRESS>Department of Optics, Palacký University, 17. listopadu
1192/12, 771 46 Olomouc, Czech Republic  Radim Filip<EMAIL_ADDRESS>Department of Optics, Palacký University, 17. listopadu 1192/12, 771 46
Olomouc, Czech Republic

###### Abstract

Quanta splitting is an essential generator of Gaussian entanglement,
apparently the most commonly occurring form of entanglement in linearised
systems. In general, such Gaussian entanglement results from the strong
pumping of a nonlinear process with the highly coherent and low-noise external
drive necessary for the linearisation. In contrast, recent experiments
involving efficient trilinear processes in trapped ions and superconducting
circuits have opened the complementary possibility to test the splitting of a
few thermal quanta without a linearisation process. Stimulated by such small
thermal energy, the strong degenerate trilinear coupling generates large
amounts of nonclassicality, detectable by more than 3 dB of distillable
quadrature squeezing, in contrast with the nondegenerate trilinear coupling.
Substantial entanglement can be generated in tandem with the nonclassicality
generation via linear coupling to a third mode, achieving the entanglement
potential of the nonclassical states. Such entanglement, although not captured
by the covariance matrix, grows with the mean number of split thermal quanta;
something which does not happen with Gaussian entanglement. We extend our
analysis to higher-order degenerate quanta splitting processes to demonstrate
that this non-Gaussian entanglement is already common for strongly nonlinear
systems at the low-energy thermal level. We thus shed light on the generic
entanglement mechanism for nonlinear bosonic systems operating in such low-
energy regimes.

Introduction— Entangled quantum harmonic oscillators form the backbone of many
fields of quantum physics and applications in quantum technology. The most
frequently occurring entanglement seems to be the Gaussian entanglement of
weakly nonlinear systems, where a strong, coherent, and low-noise external
drive is used to linearise nonlinear systems Braunstein and van Loock (2005).
However, strongly nonlinear quantum interactions exhibit a much richer and
exciting variety of phenomena which can potentially produce entanglement from
low thermal energy, although typically at the cost of ease of detection,
simple explanations and analytical descriptions. Such unexplored phenomena may
stimulate new nonlinear optics experiments using second- and high-harmonic
generation Grosse  _et al._ (2008); Olsen (2018), three and four wave mixing
Boyer  _et al._ (2008); Coelho  _et al._ (2009); González  _et al._ (2018) or
multiphoton Kerr processes Silberhorn  _et al._ (2001); Gabriel  _et al._
(2011). In the microwave domain, superconducting circuits may take advantage
of the nonlinearity of Josephson junctions to generate parametric
amplification Perelshtein  _et al._ (2021), three-photon downconversion Chang
_et al._ (2020) and induce trilinear Hamiltonians (SNAILs) Frattini  _et al._
(2017); Quijandría  _et al._ (2020) outside the rotating wave approximation.
Among the growing number of such processes, entanglement might be much more
common than anticipated and even appear more autonomously without the need for
a strong, coherent, and low-noise external drive.

The mechanical motion of trapped atoms is an excellent candidate for a proof-
of-principle investigation of such entanglement due to several properties: (i)
the strong nonlinear coupling of highly linear mechanical oscillators, (ii)
low damping and substantial shielding from ambient background noises and (iii)
the potential scalability of the trapped atoms Leibfried  _et al._ (2003);
Tiecke  _et al._ (2014). For the mechanical motion of electrically trapped
ions, the intrinsic nonlinearity of the Coulomb coupling can be leveraged into
several multimode versions of a trilinear Hamiltonian Nie  _et al._ (2009) of
the type historically studied to investigate generalisations of squeezing to
higher order moments Braunstein and McLachlan (1987). Recently, such
Hamiltonians, including fully degenerate, partially degenerate and
nondegenerate interactions, were also achieved in superconducting circuits
Chang  _et al._ (2020); Vrajitoarea  _et al._ (2020). Scalability and the
ultimate removal of small parasitic nonlinearities in the individual
oscillators is under successful development. The trapped ion versions of such
Hamiltonians, both partially degenerate and nondegenerate, comprising two-mode
and three-mode nonlinearities respectively, have been experimentally achieved
Ding  _et al._ (2017) and shown to be effective for tasks such as quantum
simulation Ding  _et al._ (2018) and quantum refrigeration Maslennikov  _et
al._ (2019). These systems are particularly suitable for proof-of-principle
studies of new instances of non-Gaussian entanglement and further
exploitation.

In this Letter, we analyse the deterministic creation of generic non-Gaussian
entanglement arising from the splitting of a few thermal quanta in trilinear
interactions and increasing with an increase in thermal occupation. In order
to produce sufficiently large entanglement of two oscillators modes, induced
by few thermal quanta, there are necessarily three modes, with one thermal
mode implementing the internal thermal driving while the remaining two,
receiving quanta from the thermal energy via the nonlinear interaction, become
entangled. We begin our discussion from the lowest order nonlinear
interactions between linear oscillators, usually the most feasible in
experiments. Nondegenerate trilinear couplings, while producing Gaussian
entanglement under a classical pump Ou  _et al._ (1992), do not produce
bipartite entanglement from a thermal pump among any combination of the modes
González  _et al._ (2018); Agustí  _et al._ (2020) (see Supplementary Material
(SM) for details). Hence, we focus on a degenerate trilinear interaction
generating significant nonclassicality from just a few thermal quanta.
Nonclassicality is qualitatively visible both in the discrete energy basis via
Klyshko’s criteria Klyshko (1996); Innocenti  _et al._ (2022) and
quantitatively in the continuous variable basis, where it manifests via
distillable squeezing surpassing the 3dB limit Filip (2013). It is remarkably
robust, remaining detectable both by Klyshko’s criteria and distillable
squeezing despite an initial thermal occupation greater than one phonon in the
degenerate mode. When the degenerate interaction is continuously combined with
a linear coupling to an auxiliary mode the resulting non-Gaussian enhancement
of logarithmic negativity can unconditionally reach $\sim 0.77$, saturating
the entanglement potential of the nonclassical states. Non-Gaussian
entanglement thus materialises straightforwardly via splitting of a few
thermal quanta in a regime very different to the Gaussian regime involving a
strong coherent pump. To make this attractive for proof-of-principle tests, we
navigate towards possible experimental tests with trapped ions, or other
alternative platforms, leaving the details for further experimental
evaluation.

Trilinear Hamiltonians with Trapped Ions— A pair of ions with equal mass $m$
and charge $e$, contained in a harmonic trap, are coupled via the Coulomb
interaction. If the radial trapping frequencies $\omega_{x}$ and $\omega_{y}$
are assumed to be much greater than the axial trapping frequency $\omega_{z}$,
then the ions distribute themselves along the $z$-axis with equilibrium
positions $(0,0,z_{i}^{0})$ Marquet  _et al._ (2003), where
$z_{1}^{0}=-\sqrt[3]{\frac{e^{2}}{4m\omega_{z}^{2}}}=-z_{2}^{0}$. The Coulomb
interaction can be expanded to second order, inducing a natural transformation
to the normal modes of the motion, wherein the collective vibrational modes
are decoupled from each other. The first step beyond this harmonic
approximation, involving third order terms in the expansion, removes the
decoupling of the spatial motion leading to interactions involving both the
radial and axial modes. Components of these nonlinear interactions are
resonant properties of the collective motion of the ions.

Recent experiments have achieved regimes in which these higher order
interactions become relevant. The third order term in the Coulomb potential
expansion is expressed, with $Q_{\alpha}=\frac{1}{2}(r_{1\alpha}-r_{2\alpha})$
the relative position for axes $\alpha={x,y,z}$, as

$V^{(3)}=\frac{e^{2}}{4|Q_{0}|^{4}}Q_{z}\left(3\left(Q_{x}^{2}+Q_{y}^{2}\right)-2Q_{z}^{2}\right)\,,$
(1)

where $Q_{0}=\frac{1}{2}(z_{1}^{0}-z_{2}^{0})$. Evidently, the relative motion
of the ions becomes anharmonic and, importantly, the axial motion becomes
coupled to the radial motion. From here, we will assume that the most
significant components of the potential are those contributing to the
interaction between the $x$ and $z$ spatial components. In this approximation,
the potential energy can be expressed as

$V\approx\sum_{\alpha=x,z}m\Omega_{\alpha}^{2}Q_{\alpha}^{2}+\frac{3e^{2}}{4|Q_{0}|^{4}}Q_{z}Q_{x}^{2}\,,$
(2)

where we have defined the normal mode frequencies
$\Omega_{x}^{2}=\omega_{x}^{2}-\frac{e^{2}}{4m|Q_{z}^{0}|^{3}}$ and
$\Omega_{z}^{2}=\omega_{z}^{2}+\frac{e^{2}}{2m|Q_{z}^{0}|^{3}}$. Alternatively
for superconductng circuits trilinear interactions using SNAILs Frattini  _et
al._ (2017) can be considered to reach Eq. (1).

The motional variables can be quantized by promoting them to operators
satisfying the canonical commutation relations
$[\hat{Q}_{\alpha},\hat{P}_{\alpha}]=i\hbar$, where $\hat{P}_{\alpha}$ is the
momentum operator corresponding to the normal mode $\hat{Q}_{\alpha}$. The
canonical variables of the relative motion can be linked to annihilation
operators by the relations
$\hat{Q}_{z}=\sqrt{\frac{\hbar}{2m\Omega_{z}}}(\hat{a}+\hat{a}^{\dagger})$ and
$\hat{P}_{z}=i\sqrt{\frac{\hbar m\Omega_{z}}{2}}(\hat{a}^{\dagger}-\hat{a})$,
and $\hat{Q}_{x}=\sqrt{\frac{\hbar}{2m\Omega_{x}}}(\hat{b}+\hat{b}^{\dagger})$
and $\hat{P}_{x}=i\sqrt{\frac{\hbar
m\Omega_{x}}{2}}(\hat{b}^{\dagger}-\hat{b})$, with
$[a,\,a^{\dagger}]=[b,\,b^{\dagger}]=\mathbb{I}$. Expanding $V$ in terms of
the annihilation operators and moving to a frame rotating with the respect to
the free motion generates an interaction Hamiltonian $H$ of the form

$H=\frac{3e^{2}}{4|Q_{0}|^{4}}\left(ae^{-i\Omega_{z}t}+a^{\dagger}e^{i\Omega_{z}t}\right)\\\
\times\left(b^{2}e^{-i2\Omega_{x}t}+b^{\dagger}{}^{2}e^{i2\Omega_{x}t}+b^{\dagger}b+1\right).$
(3)

When $\Omega_{z}=2\Omega_{x}$, several terms become resonant and in the
rotating wave approximation the interaction Hamiltonian takes the trilinear
form

$H=\Omega_{T}\left(ab^{\dagger}{}^{2}+a^{\dagger}b^{2}\right),$ (4)

where $\Omega_{T}=\frac{3e^{2}}{4|Q_{0}|^{4}}$. For a more detailed
Hamiltonian analysis see, for instance, Refs. Marquet  _et al._ (2003); Lemmer
(2018). This is a partially degenerate trilinear interaction in which the
creation (annihilation) of a phonon in the axial direction results in the
annihilation (creation) of a pair of phonons in the radial direction.

Thermally Induced Nonclassicality and Non-Gaussian Entanglement— The quanta
splitting and ensuing nonclassicality from the nonlinear dynamics may proceed
from a small amount of classical thermal noise present in one of the
subsystems. By initialising the axial ($a$) or radial ($b$) components of the
motion in thermal states while the remainder is initialised in the ground
state, the degree to which the nonlinearity converts the initial resource of
thermal noise into nonclassical behaviour can be examined (Fig. 1). In
general, if the radial mode is used as a thermal pump of the axial mode then
the nonclassicality and corresponding entanglement potential is small. In
contrast, our starting point is that $a$, using the small incoherent energy of
a thermal pump, can drive nonclassicality in $b$.

Figure 1: The upper panels demonstrate the production of nonclassicality due
to thermal quanta splitting. Mode $a$ is coupled to mode $b$ through the
trilinear coupling given by Eq. (4). Initially, $a$ is in a thermal state
(upper left) with mean thermal energy $\bar{n}$, while $b$ is in the ground
state (for thermal occupation of $b$, see SM). Thermal quanta in mode $a$ are
dynamically split into pairs in mode $b$, which generates nonclassical phase-
insensitive phonon statistics (upper right), verified by Klyshko’s criteria
Klyshko (1996); Innocenti  _et al._ (2022). Furthermore, increasing $\bar{n}$
makes the nonclassical effects more pronounced. Purple and orange correspond
to $\bar{n}=1$ and 3, respectively, while green corresponds to a squeezed
state with 3 dB of squeezing. The nonclassicality of mode $b$ after the
trilinear interaction can also be quantified in a phase sensitive manner via
distillable squeezing. The probability densities $P(q_{x})$, of the
dimensionless quadrature $\hat{q}_{x}=\frac{b+b^{\dagger}}{\sqrt{2}}$, are
compared on the lower left for the ground state (dashed black) and with
thermally induced nonclassicality, $\bar{n}=3$, (red) demonstrating the sub-
Planck structure around the origin generated by the nonlinear interaction. The
yellow curve indicates the probability distribution of the distillable
squeezed state close to the asymptotic limit, corresponding to the yellow
curve in the lower right panel. Universal distillation of squeezing from
multiple copies of the nonclassical state in mode $b$ is shown on the lower
right, with the asymptotic limit in grey Filip (2013) and the 3 dB squeezing
limit in dashed green. Consuming progressively more copies produces a variance
moving towards and then below the shot-noise level (dashed black). The
distillation procedure saturates (blue), limited by the nonclassicality
present in mode $b$ after the trilinear interaction. Note, we have selected
the time to coincide with the first peak of the entanglement dynamics.

In the ideal case, the nonclassicality of mode $b$ can be immediately seen
from the phonon distribution $P_{k}=\braket{k}{\rho}{k}$; the reduced state of
mode $b$ is nonclassical due to its population of only even Fock states
irrespective of the strength of the thermal pump. This is verified directly by
Klyshko’s criteria Klyshko (1996); Innocenti  _et al._ (2022), which are
automatically satisfied for odd number occupations and satisfied strongly for
even occupations close to the peaks of the distribution. The nonclassicality
is phase insensitive, reflecting the incoherent phase of the thermal pump (see
SM for comparison with a coherent pump). Simultaneously, the non-Gaussian
position probability density, despite having a variance greater than the
ground state, displays sub-Planck structure Zurek (2001) around the global
maximum verified by the presence of universally distillable squeezing Filip
(2013), with the degree of squeezing available increasing with the thermal
energy. Notably, despite the presence of such sub-Planck structures the Wigner
function is positive at the level of the nonclassical state of $b$. The lower
right panel of Fig. 1 shows that the degree of distillable squeezing due to
quanta splitting is stimulated by $\bar{n}$ beyond 3 dB. Remarkably, these
nonclassical aspects are observable also for an initial thermal occupation of
$b$ larger than a single quantum (see SM).

The amount of nonclassicality capable of being converted to entanglement in a
particular mode can be further characterised through the entanglement
potential (EP), denoted $\mathcal{E}$ Asbóth  _et al._ (2005). The EP of mode
$a$ or $b$ is evaluated through its capacity to become entangled via passive
interactions with the ground state of an oscillator. This can be quantified by
determining the logarithmic negativity (LN) Plenio (2005), of the following
state:

$\displaystyle\rho_{\mathcal{E}}=U_{\text{BS}}^{\dagger}\,\left(\rho\otimes\ket{0}\bra{0}\right)U_{\text{BS}},$
(5)

where $U_{\text{BS}}=e^{\frac{\pi}{4}(c^{\dagger}A-cA^{\dagger})}$, where $c$
is an auxiliary mode prepared in the ground state, $A=a,b$ and $\rho$ is the
reduced state corresponding to the mode $A$.

The nonclassicality converted from the thermal noise resource can be leveraged
to create entanglement with another mode $c$ which interacts simultaneously
with $b$ via the simple linear interaction

$\displaystyle H_{\text{aux}}$
$\displaystyle=H+g(b\,c^{\dagger}+b^{\dagger}c)\,.$ (6)

In our analysis, we have selected $g=\Omega_{T}$, which appears to give the
best performance. Deviation from equality simply results in reduced LN,
although the qualitative effects remain the same. Note that the beamsplitter
coupling in Eqs. (5,6) in our case is incapable of generating any amount of
entanglement from classical states. The entanglement is not a short-term
effect, in the sense that it arises after an initial quiescent period in the
dynamics, and optimisation of the interaction time necessaerily occurs over
anharmonically oscillating dynamics (see SM). To compromise we select the
first peak of EP in order to balance finding significant entanglement while
remaining closer to the short times available to experimental settings, where
decoherence will be relevant.

Figure 2: For the trilinear coupling given by Eq. (4), the entanglement
potential $\mathcal{E}_{b}$ (blue) of mode $b$ generated by the thermal noise
in $a$ increases as the noise is increased. In the presence of a linear
coupling Eq. (6) to an auxiliary mode $c$, we show the logarithmic negativity
$\mathcal{L}$ between $b$ and $c$ (red), $a$ and $b$ (yellow), and $a$ and $c$
(green). $a$ is initially in a thermal state while $b$ and $c$ are initially
in their respective ground states. The values plotted here correspond to the
first peak in the temporal evolution of their respective quantities
$\mathcal{E}$ and $\mathcal{L}$ which otherwise undergo rather complex
evolution. $\mathcal{L}_{bc}$ closely mimics the behavior of $\mathcal{E}_{b}$
indicating that the entanglement potential of $b$ can be realised dynamically,
despite residual entanglement among the remaining modes. $\mathcal{E}_{b}$ is
further compared against the EP of a squeezed vacuum state (black), whose
degree of squeezing is selected in accord with the asymptotic distillation of
squeezing available for the corresponding trilinear state.

The EP generated in $b$ by a thermal pump in $a$ is shown in Fig. 2 as the
blue line. As a signature of nonclassicality, this is connected to the
previous discussion of distillable squeezing by noting that if the non-
Gaussian state with large EP is used to distill a Gaussian squeezed state (in
the asymptotic limit), then this squeezed state generates similar EP to the
non-Gaussian state from the trilinear interaction (black line). By introducing
an auxiliary mode $c$, as in Eq. (6), the resulting LN, denoted
$\mathcal{L}_{b:c}$, dynamically fulfills the EP of $b$ as the red line in
Fig. 2 achieves at least the value of $\mathcal{E}_{b}$. This persists even
though residual LN is generated with the thermal pump mode
($\mathcal{L}_{a:c}$ and $\mathcal{L}_{a:b}$) and is subsequently traced out.
The Gaussian entanglement, extracted from the covariance matrix of the state,
is always zero during the evolution, indicating that the thermally induced
entanglement involves correlations beyond the covariance matrix.

To illustrate the robustness of these phenomena we examined the deviation of
mode $b$ from the ideal ground state. Thermal occupation of mode $b$, while
reducing the absolute value of the EP, maintains the enhancement of
nonclassicality due to increasing the thermal energy in $a$ (see SM).
Remarkably, this occurs even when the thermal occupation of $b$ is greater
than 1 and can occur when greater than $\bar{n}$. This allows experimental
tests of these phenomena without being bound to the ground state of mode $b$.
Additionally, we examined the effect on the EP and LN when the trilinear
interaction occurs outside the RWA, so that the free evolution of the
oscillators contributes (see SM). While the EP is robust to these effects, the
LN decreases greatly, while still retaining the property that it increases
with $\bar{n}$. Finally we also examined the effect of coupling the modes to a
thermal environment. While both EP and LN are reduced by the interaction with
the environment, significant values of both quantities remain (see SM).

Higher Order Nonlinearities— It is possible that higher nonlinearities will
exhibit different phenomena and so we examine highly nonlinear Hamiltonians of
the form $H_{\text{NL}}\propto a(b^{\dagger})^{N}+a^{\dagger}b^{N}$, where
$N>2$. Continuing the logic, the internal thermal driving of $a$ will drive
higher order quanta splitting, resulting in generalised squeezing Braunstein
and McLachlan (1987) in $b$ and producing highly nonclassical and non-Gaussian
states with large entanglement potential.

Figure 3: (a) The probability densities $P(q_{x})$ from
$H_{\text{NL}}=\Omega_{T}(a^{\dagger}b^{N}+ab^{\dagger\,N})$ with $N=2$ (red),
3 (blue), 4 (green) and 5 (yellow) are compared at the optimal time (see
text), compared with the ground state distribution (dashed black). (b)
Universal (squares) and non-universal (crosses) squeezing distillation as a
function of number of copies of the state. For odd $N$ the universal method
does not find distillable squeezing around the global maximum of the
distribution, unlike the case of even $N$. Instead, the non-universal method
searches the whole distribution for nonclassical structures and distills them
to find distillable squeezing. In both panels mode $a$ is initially in a
thermal state with $\bar{n}=1$ and mode $b$ is in the ground state. Figure 4:
Entanglement potential $\mathcal{E}_{b}$ (solid curves) and
$\mathcal{L}_{b:c}$ (dashed curves) generated by incoherent thermal noise in
$a$ for Hamiltonians $H_{\text{NL}}$. The evolution times and coupling
strengths are as previously described in the main text. The entanglement
extracted by the auxiliary linear coupling ceases to closely follow the EP for
$N>2$.

In Fig. 3 we show the probability densities in $q_{x}$ for the nonclassical
states produced by these higher order nonlinearities. They are clearly non-
Gaussian, with embedded sub-Planck structures. However the distillable
squeezing shows features distinct from the $N=2$ case considered above. The
universal method of distillation Filip (2013), which extracts nonclassical
features from the curvature around the global maximum of the distribution,
does not produce squeezing for odd $N$. Indeed this is reflected in the
distributions which have low curvature around their maxima. Instead the non-
universal method Filip (2013); Filip and Zapletal (2014), which optimises the
distillation around arbitrary features of the distribution, can extract
squeezing even from the odd parity Hamiltonians. By measuring the copies in
the distillation around sub-Planck features in the distribution, optimising
the measurement at each distillation step, distillable squeezing, and thus
nonclassicality, can still be detected. However the distillable squeezing
found through the non-universal method does not increase with $\bar{n}$, in
contrast with the universal method analysed in detail in Fig 1.

This nonclassicality can again be used to generate entanglement. Fig. 4 shows
the entanglement potential and logarithmic negativity via coupling to an
auxiliary state for various $N$. Increasing $N$ tends to increase the
nonclassicality generated by the same thermal energy, however, similar to the
failure of universal squeezing distillation, the anticipated generation of
entanglement via linear coupling in parallel with the generation of
nonclassicality eventually fails to reach the EP for larger values of
$\bar{n}$. For $N=3$, replacing the linear coupling with a trilinear coupling
does not improve $\mathcal{L}_{b:c}$. This suggests that among these the
trilinear Hamiltonian of Eq. (4), combined with a linear coupling to another
mode, is the most effective in terms of generating non-Gaussian entanglement
from small thermal energy.

Discussion— We have singled out the simplest nonlinearity coupling linear
oscillators, the degenerate trilinear interaction, combined with a linear
coupling in order to transmute the thermal energy of a few quanta into greater
nonclassicality and non-Gaussian entanglement. Other possibilities among three
mode interactions are surprisingly limited in their capacity to generate
bipartite entanglement González  _et al._ (2018); Agustí  _et al._ (2020) (see
SM). That larger thermal fluctuations are converted into increasing coherent
quantum phenomena, such as entanglement Bose  _et al._ (2001); Kim  _et al._
(2002); Marek  _et al._ (2016); Laha  _et al._ (2022), detected by logarithmic
negativity, functions as a direct proof of genuine nonlinear quantum effects,
as it cannot happen for linearised dynamics within the Gaussian approximation
provided by the covariance matrix. While such dynamics intrinsically involves
higher than quadratic interaction terms, and shows quantum non-Gaussian
effects Ding  _et al._ (2017) such a direct experimental proof is yet to be
demonstrated. In the past, such features were predicted in systems taking
advantage of discrete nonlinearities in the level structure of atoms or
superconducting circuits with very strong saturation Marek  _et al._ (2016);
Slodička  _et al._ (2016); Laha  _et al._ (2022) while here entanglement
materialises without the direct use of saturation in such discrete
nonlinearities. It is also conceptually different to thermally induced
entanglement using controlled-SWAP gates Filip  _et al._ (2002); Gao  _et al._
(2019).

This proposal, inspired by recent experimental achievements regarding
trilinear interactions in trapped ion Ding  _et al._ (2017, 2018); Maslennikov
_et al._ (2019) and superconducting circuit platforms Chang  _et al._ (2020);
Vrajitoarea  _et al._ (2020), provides a clear example of the hidden power of
multimode nonlinear systems to achieve, autonomously and unconditionally,
increasing nonclassical quantum behaviour leveraging only small thermal
energy. Furthermore, the linear coupling required to directly test the
parallel generation of entanglement has already been proposed for optical
systems Krastanov  _et al._ (2021) and is likely to be easily adapted to
phononic systems and superconduting circuits.

##  I Acknowledgements

The authors acknowledge funding from Project 21-13265X of the Czech Science
Foundation. This project has received funding from the European Union’s 2020
research and innovation programme (CSA—Coordination and support action,
H2020-WIDESPREAD-2020-5) under grant agreement No. 951737 (NONGAUSS).

## References

  * Braunstein and van Loock (2005) S. L. Braunstein and P. van Loock, Rev. Mod. Phys. 77, 513 (2005).
  * Grosse  _et al._ (2008) N. B. Grosse, S. Assad, M. Mehmet, R. Schnabel, T. Symul, and P. K. Lam, Phys. Rev. Lett. 100, 243601 (2008).
  * Olsen (2018) M. K. Olsen, Phys. Rev. A 97, 033820 (2018).
  * Boyer  _et al._ (2008) V. Boyer, A. M. Marino, R. C. Pooser, and P. D. Lett, Science 321, 544 (2008).
  * Coelho  _et al._ (2009) A. S. Coelho, F. A. S. Barbosa, K. N. Cassemiro, A. S. Villar, M. Martinelli, and P. Nussenzveig, Science 326, 823 (2009).
  * González  _et al._ (2018) E.-A. R. González, A. Borne, B. Boulanger, J.-A. Levenson, and K. Bencheikh, Phys. Rev. Lett. 120, 043601 (2018).
  * Silberhorn  _et al._ (2001) C. Silberhorn, P. K. Lam, O. Weiß, F. König, N. Korolkova, and G. Leuchs, Phys. Rev. Lett. 86, 4267 (2001).
  * Gabriel  _et al._ (2011) C. Gabriel, A. Aiello, W. Zhong, T. G. Euser, N. Y. Joly, P. Banzer, M. Förtsch, D. Elser, U. L. Andersen, C. Marquardt, P. S. J. Russell, and G. Leuchs, Phys. Rev. Lett. 106, 060502 (2011).
  * Perelshtein  _et al._ (2021) M. Perelshtein, K. Petrovnin, V. Vesterinen, S. H. Raja, I. Lilja, M. Will, A. Savin, S. Simbierowicz, R. Jabdaraghi, J. Lehtinen, L. Grönberg, J. Hassel, M. Prunnila, J. Govenius, S. Paraoanu, and P. Hakonen, arXiv:2111.06145 [quant-ph] (2021).
  * Chang  _et al._ (2020) C. W. S. Chang, C. Sabín, P. Forn-Díaz, F. Quijandría, A. M. Vadiraj, I. Nsanzineza, G. Johansson, and C. M. Wilson, Phys. Rev. X 10, 011011 (2020).
  * Frattini  _et al._ (2017) N. E. Frattini, U. Vool, S. Shankar, A. Narla, K. M. Sliwa, and M. H. Devoret, Appl. Phys. Lett. 110, 222603 (2017).
  * Quijandría  _et al._ (2020) F. Quijandría, G. Johansson, A. Ferraro, S. Gasparinetti, G. Ferrini, and T. Hillmann, Phys. Rev. Lett. 125, 160501 (2020).
  * Leibfried  _et al._ (2003) D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, Rev. Mod. Phys. 75, 281 (2003).
  * Tiecke  _et al._ (2014) T. G. Tiecke, J. D. Thompson, N. P. de Leon, L. R. Liu, V. Vuletić, and M. D. Lukin, Nature 508, 241 (2014).
  * Nie  _et al._ (2009) X. R. Nie, C. F. Roos, and D. F. V. James, Phys. Lett. A 373, 422 (2009).
  * Braunstein and McLachlan (1987) S. L. Braunstein and R. I. McLachlan, Phys. Rev. A 35, 1659 (1987).
  * Vrajitoarea  _et al._ (2020) A. Vrajitoarea, Z. Huang, P. Groszkowski, J. Koch, and A. A. Houck, Nat. Phys. 16, 211 (2020).
  * Ding  _et al._ (2017) S. Ding, G. Maslennikov, R. Hablützel, H. Loh, and D. Matsukevich, Phys. Rev. Lett. 119, 150404 (2017).
  * Ding  _et al._ (2018) S. Ding, G. Maslennikov, R. Hablützel, and D. Matsukevich, Phys. Rev. Lett. 121, 130502 (2018).
  * Maslennikov  _et al._ (2019) G. Maslennikov, S. Ding, R. Hablützel, J. Gan, A. Roulet, S. Nimmrichter, J. Dai, V. Scarani, and D. Matsukevich, Nat. Commun. 10, 202 (2019).
  * Ou  _et al._ (1992) Z. Y. Ou, S. F. Pereira, H. J. Kimble, and K. C. Peng, Phys. Rev. Lett. 68, 3663 (1992).
  * Agustí  _et al._ (2020) A. Agustí, C. W. S. Chang, F. Quijandría, G. Johansson, C. M. Wilson, and C. Sabín, Phys. Rev. Lett. 125, 020502 (2020).
  * Klyshko (1996) D. N. Klyshko, Phys. Lett. A 213, 7 (1996).
  * Innocenti  _et al._ (2022) L. Innocenti, L. Lachman, and R. Filip, npj Quantum Inf. 8, 1 (2022).
  * Filip (2013) R. Filip, Phys. Rev. A 88, 063837 (2013).
  * Marquet  _et al._ (2003) C. Marquet, F. Schmidt-Kaler, and D. James, Appl. Phys. B 76, 199 (2003).
  * Lemmer (2018) A. Lemmer,  _Quantum dynamics with trapped ions_ , Dissertation, Universität Ulm (2018).
  * Zurek (2001) W. H. Zurek, Nature 412, 712 (2001).
  * Asbóth  _et al._ (2005) J. K. Asbóth, J. Calsamiglia, and H. Ritsch, Phys. Rev. Lett. 94, 173602 (2005).
  * Plenio (2005) M. B. Plenio, Phys. Rev. Lett. 95, 090503 (2005).
  * Filip and Zapletal (2014) R. Filip and P. Zapletal, Phys. Rev. A 90, 043854 (2014).
  * Bose  _et al._ (2001) S. Bose, I. Fuentes-Guridi, P. L. Knight, and V. Vedral, Phys. Rev. Lett. 87, 050401 (2001).
  * Kim  _et al._ (2002) M. S. Kim, J. Lee, D. Ahn, and P. L. Knight, Phys. Rev. A 65, 040101 (2002).
  * Marek  _et al._ (2016) P. Marek, L. Lachman, L. Slodička, and R. Filip, Phys. Rev. A 94, 013850 (2016).
  * Laha  _et al._ (2022) P. Laha, L. Slodička, D. W. Moore, and R. Filip, Opt. Express 30, 8814 (2022).
  * Slodička  _et al._ (2016) L. Slodička, P. Marek, and R. Filip, Opt. Express 24, 7858 (2016).
  * Filip  _et al._ (2002) R. Filip, M. Dušek, J. Fiurášek, and L. Mišta, Phys. Rev. A 65, 043802 (2002).
  * Gao  _et al._ (2019) Y. Y. Gao, B. J. Lester, K. S. Chou, L. Frunzio, M. H. Devoret, L. Jiang, S. M. Girvin, and R. J. Schoelkopf, Nature 566, 509 (2019).
  * Krastanov  _et al._ (2021) S. Krastanov, M. Heuck, J. H. Shapiro, P. Narang, D. R. Englund, and K. Jacobs, Nat. Commun. 12, 191 (2021).

##  II Supplementary Material

##  III Non-Gaussian entanglement via splitting of a few thermal quanta

###  .1 Nondegenerate and Multiple Trilinear Hamiltonians

In the main analysis, we focused on degenerate trilinear coupling as the best
candidate for the generic entanglement generation from a few thermal quanta.
Alternative configurations of trilinear Hamiltonians are briefly analysed here
and demonstrated to be inefficient in producing entanglement. In the
nondegenerate case, the Hamiltonian is given by

$H^{\mathrm{ND}}=\Omega(a^{\dagger}\,b\,c+a\,b^{\dagger}\,c^{\dagger}),$ (S1)

where $\Omega$ is the coupling strength. This interaction Hamiltonian has
already been realised using three ${}^{171}\mathrm{Yb}^{+}$ ions in a Paul
trap, where $a$ corresponds to the axial mode, and $b$ and $c$ to the two
radial modes Ding  _et al._ (2018). Thermally pumping any one of the modes
$a,\,b,\,c$, does not produce entanglement among the remaining pair. However,
we can check the EP of these modes individually and extract the entanglement
via continuous linear coupling to an auxiliary mode $d$. Setting $\tau=\Omega
t$ as before, and studying the induced nonclassicality in mode $b$ (equivalent
to $c$ by symmetry) by thermally driving mode $a$, while both modes $b$ and
$c$ are initially in their respective ground states we find (Fig. S1) that EP
can be generated in mode $b$ ($c$), in negligible quantities compared to the
degenerate case analysed in the main text. As before, we also see that the
nonclassicality converted from the thermal noise resource can be leveraged via
coupling $b$ to an auxiliary mode $d$ (i.e.
$H^{\mathrm{ND}}_{\textrm{aux}}=H^{\mathrm{ND}}+g(bd^{\dagger}+b^{\dagger}d)$,
or similarly with $c$) to create entanglement which follows the behavior of
$\mathcal{E}_{b}$ ($\mathcal{E}_{c}$). Again we set $g=\Omega$. We note that
$\mathcal{E}_{a}$ is greater than $\mathcal{E}_{b}$ (or equivalently,
$\mathcal{E}_{c}$) which is not the case for the degenerate case investigated
in the main text.

Figure S1: EP of mode $a$ (black), and identical modes $b$ and $c$ (blue) are
plotted as a function of $\bar{n}$ for the nondegenerate trilinear system
given in Eq.(S1). $a$ is prepared in a thermal state characterised by
$\bar{n}$ while $b$ and $c$ are in the ground state. Adding a linear coupling
between mode $b$ and an auxiliary mode $d$, as in
$H^{\mathrm{ND}}_{\textrm{aux}}$, the LN $\mathcal{L}_{b:d}$ increases with
$\bar{n}$ but is negligible compared to the degenerate case.

Alternatively, multiple trilinear interactions with one common mode are
described by the following Hamiltonians,

$\displaystyle H_{1}$
$\displaystyle=\Omega_{1}(ab^{\dagger}{}^{2}+a^{\dagger}b^{2})+\Omega_{2}(cb^{\dagger}{}^{2}+c^{\dagger}b^{2})$
(S2) $\displaystyle H_{2}$
$\displaystyle=\Omega_{1}(ab^{\dagger}{}^{2}+a^{\dagger}b^{2})+\Omega_{2}(ac^{\dagger}{}^{2}+a^{\dagger}c^{2}).$
(S3)

In Fig. S2 we show the behavior of $\mathcal{L}_{b:c}$ for $H_{1}$ (olive) and
$H_{2}$ (green). In both the cases the induced entanglement is much smaller
compared to $\mathcal{L}_{b:c}$ of the individual interaction given by Eq.
(4).

Figure S2: $\mathcal{L}_{b:c}$ for $H_{1}$ (olive), $H_{2}$ (green) is
compared with $\mathcal{L}_{b:c}$ for $H$ (red) in Fig. 2. The LN for multiple
trilinear interactions is much smaller than for the single degenerate version
preferred in the main text.

###  .2 Coherent and Incoherent Internal Driving

In the main text we focus on the incoherent resource of thermal states to
generate coherent quantum phenomena. With coherent resources the amount of
nonclassicality generated can be greatly increased, interestingly however the
qualitative behaviour can already be observed for the thermal states.

Figure S3: The role of different initial states on the entanglement
production. $\mathcal{E}_{b}$ (blue) and $\mathcal{L}_{b:c}$ (red) are
calculated using different resources to drive $a$. Represented are the thermal
states (solid lines), coherent states (solid circles), Poissonian noise
(dashed) and Fock states (solid triangles). $b$ is always initialised in the
ground state. Note that for coherent states and Poissonian noise we have
$\bar{n}=|\alpha|^{2}$.

In Fig. S3 we show the EP and LN with $a$ prepared in a coherent state (solid
circles) $\ket{\alpha}$ with $|\alpha|^{2}=\bar{n}$. This is much greater than
that achievable with thermal states (solid lines). If the phase is uncertain
then the state has Poissonian noise in the Fock basis, and we obtain the phase
randomised coherent states (PRCS) defined by

$\frac{1}{2\pi}\int_{0}^{2\pi}d\theta\,\ket{\alpha e^{i\theta}}\bra{\alpha
e^{i\theta}}=e^{-|\alpha|^{2}}\sum_{n=0}^{\infty}\frac{|\alpha|^{2n}}{n!}\ket{n}\bra{n}.$
(S4)

Here, it is instructive to compare with less noisy and more coherent drives
under the same conditions. Poissonian noise (dashed lines) produces less
nonclassicality and entanglement than the coherent state, but more than the
thermal states. Furthermore, the Fock states, which have zero noise around the
phonon number and no quantum coherence, but are extremely nonclassical, induce
larger nonclassicality than the other examples provided.

Figure S4: The variance of the generalised quadrature $X_{\theta}$ of mode $b$
as a function of scaled time $\tau=\Omega_{T}t$ and phase $\theta$, when $a$
is driven coherently (with strength $\bar{n}=|\alpha|^{2}=3$).

If mode $a$ is prepared in a coherent state then the nonclassical state
prepared in $b$ is phase sensitive, in contrast with that generated by thermal
energy, with the phase determined by the phase of the coherent state. For
large $\alpha$ the mode $a$ in the trilinear Hamiltonian can be approximated
with the c-number $\alpha=|\alpha|e^{i\phi}$ as

$H\sim\alpha(b^{\dagger})^{2}+\alpha^{*}b^{2}=|\alpha|\left(b^{\dagger\,2}e^{i\phi}+b^{2}e^{-i\phi}\right)\,.$
(S5)

The coherent phase $\phi$ can be connected with the quadrature angle
$\frac{\theta}{2}$ via $\phi=\theta-\frac{\pi}{2}$, so that the Hamiltonian
takes the standard squeeze operator form

$H=-i|\alpha|\left(b^{\dagger\,2}e^{i\theta}-b^{2}e^{-i\theta}\right)\,.$ (S6)

Fig. S4 shows the squeezing in the generalised quadrature
$X_{\theta}=\frac{1}{\sqrt{2}}\left(be^{-i\frac{\theta}{2}}+b^{\dagger}e^{i\frac{\theta}{2}}\right)$
for the coherent state with $|\alpha|^{2}=\bar{n}=3$ and $\phi=0$. The
squeezing is most manifest along the quadrature axis defined by
$\theta=\frac{\pi}{2}$, in accordance with the prediction of the saturated
trilinear system. In contrast, for a phase random initial state such as a
thermal state or PRCS, the nonclassicality is phase insensitive.

###  .3  Influence of the Environment

In our treatment of the trilinear dynamics we have assumed unitary evolution
due to the high quality of the dynamics and sufficient decoupling from the
environment for both trapped ions and superconducting circuits over time
required for such experiments. In this Section, we detail our numerical
observations on the robustness of the thermally induced entanglement
generation against such effects. If we ignore the contributions of the Lamb
shift leading to a small renormalization of the system energy levels, the
Lindblad master equation describing such environmental effects is given by

$\displaystyle\displaystyle\frac{d\rho_{S}}{dt}=-i\left[H,\rho_{S}\right]+\sum_{k}\left(L_{k}\rho_{S}L_{k}^{\dagger}-\frac{1}{2}\left[\rho_{S},\,L_{k}^{\dagger}L_{k}\right]\right).$
(S7)

Here, $\rho_{S}(t)$ is the reduced density matrix of the system. The operators
$L_{k}=\sqrt{\lambda}_{k}A_{k}$ are the Lindblad jump operators, and the
environment couples to the system through the operators $A_{k}$ with coupling
rates $\lambda_{k}$. For the two mechanical modes coupled to the environment
with same temperature $n_{\text{th}}$, we have the Lindblad operators
($L_{k}$’s) given by $\sqrt{\lambda_{a}(1+\bar{n}_{\text{th}})}\,a$,
$\sqrt{\lambda_{a}\,\bar{n}_{\text{th}}}\,a^{\dagger}$, and
$\sqrt{\lambda_{b}(1+\bar{n}_{\text{th}})}\,b$,
$\sqrt{\lambda_{b}\,\bar{n}_{\text{th}}}\,b^{\dagger}$. When considering the
LN there are additional terms for the mode $c$ involved in the linear
coupling.

In Fig. S5 the effect of decoherence on both the EP and LN is examined in
terms of the strength of the coupling to the environment. The system evolves
with mode $a$ in a thermal state with $\bar{n}=2$ and the remaining modes in
the ground state, until a fixed time set by the unitary evolution in the main
text. Both quantities experience a sharp decay to zero with similar behaviour,
with the rate increasing with the temperature of the environment. However in
Fig. S6 the direct evolution against time shows that while the nonclassicality
is largely reduced, significant quantities are still visible, albeit at
slightly shorter times.

Importantly, as explained in the following section, the enhancement due to
initial thermal noise remains visible despite the reduction in logarithmic
negativity.

Figure S5: Effect of environmental induced decoherence on $\mathcal{E}_{b}$
and $\mathcal{L}_{b:c}$ when all modes are connected to the same thermal bath
with different temperatures $\bar{n}_{\text{th}}=0$ (blue), 0.1 (red), 0.5
(yellow), and 1.0 (green). Initially, $a$ is in a thermal state with
$\bar{n}=2$, and $b$ and $c$ are in the ground states. The black dashed lines
represents the unitary evolution ($\lambda=0$). We fix the time as the
optimised time of the unitary case, as described in the main text.
Figure S6: Effect of environmental induced decoherence on $\mathcal{E}_{b}$
and $\mathcal{L}_{b:c}$ when all modes are connected to the same bath with
$\bar{n}_{\text{th}}=0$. Initially, $a$ is in a thermal state with
$\bar{n}=2$, and $b$ and $c$ are in the ground state. We fix the time as the
optimised time of the unitary case, as described in the main text.

###  .4 Thermal Occupation of the Radial Mode

In our analysis of the nonclassicality and EP for the trilinear system we
assumed the idealisation that mode $b$ is prepared in the ground state with
the only thermal energy in mode $a$. This is because both trapped ions
experiments and superconducting circuits can achieve good ground states of
these modes. However, to check the stability outside this limit we now briefly
consider the effect of thermal occupation of modes $b$ and $c$.
Advantangeously, the thermally induced effects persist even if the thermal
energy in mode $b$, characterised by an average occupation $\bar{n}_{b}$, is
greater than 1. Fig. S7 shows a recreation of Fig. 1 from the main text with
$\bar{n}_{b}=1.5$. The clear demonstration of nonclassicality from the phonon
statistics in Fig. 1 is lost however Klyshko’s criteria still verify
nonclassicality. Similarly, the distillable squeezing in $q_{x}$ can still be
reduced below the shot noise by very large thermal energy in $a$.

Figure S7: Analysis of the nonclassicality in mode $b$ as in Fig. 1 with
$\bar{n}_{b}=1.5$. The clear visual generation of nonclassicality by the
occupation of only even Fock states is lost, however Klyshko’s criteria still
indicate nonclassicality. The distillable squeezing can also still achieve
values below the shot noise by preparing $a$ in thermal states with
$\bar{n}_{a}>3$.

In Fig. S8 we show how the distillable squeezing varies as a function of
$\bar{n}_{b}$. The amount of thermal energy in $a$ required to see reduction
in the position variance below the shot noise increases with $\bar{n}_{b}$,
however it is still possible to see nonclassicality even for large values of
$\bar{n}_{b}$.

Figure S8: Universal distillable squeezing when mode $b$ deviates from the
ideal ground state configuration to a thermal state with average phonon number
$\bar{n}_{b}=0$ (blue), 0.5 (red), 1 (green) and 1.5 (yellow) as a function of
$\bar{n}_{a}$. The distillable squeezing indicates nonclassicality for large
enough thermal energy in $a$.
Figure S9: $\mathcal{E}_{b}$ as a function of the thermal occupation
$\bar{n}_{b}$ of the radial mode $b$ for several values of $\bar{n}_{a}$.
While the overall effect is to reduce the amount of EP produced, the EP is
still visibly increased by an increase in the thermal occupation of mode $a$.

Finally in Fig. S9 we also show the EP as a function of $\bar{n}_{b}$.
Increasing $\bar{n}_{b}$ (Fig S9) drastically decreases the absolute value of
the EP. However our main result, that the EP is increased by increasing the
thermal energy in mode $a$, is still visible up to saturation at
$\bar{n}\approx 3$.

###  .5 Effect of Free Motion on Nonclassicality and Entanglement

Our analysis takes place in a frame of reference that eliminates the free
motion of the component oscillators, which is suited to the platforms of
trapped ions that we consider as well as similar systems such as
superconducting circuits. If such trilinear Hamiltonians occur in a context in
which such transformations cannot be performed then the effect of the free
motion must be taken into account. The Hamiltonian, including the linear
coupling required to study the LN, is modified to

$H_{\text{FM}}=\sum_{k=a,b,c}\omega_{k}k^{\dagger}k+\Omega_{T}(a^{\dagger}b^{2}+ab^{\dagger}{}^{2})+g(b^{\dagger}c+bc^{\dagger})\,.$
(S8)

Figure S10: The entanglement potential and logarithmic negativity in the
presence of free motion. Blue and red correspond to the EP and LN analysed in
Fig. 2 of the main text. The green and yellow correspond to the
nonclassicality and entanglement dynamics involving the free motion.

We therefore briefly verify if the observed effects qualitatively persist if
the free evolution is present. As in the main text, we select $\Omega_{T}=g$.
Due to energy conservation, there is only one parameter for the frequency, as
$\omega_{a}=2\omega_{b}=2\omega_{c}$. While the effects of free motion terms
are typically destructive to quantum phenomena, here we find that the EP is in
fact quite resilient to these effects but the LN suffers as shown in Fig S10.
We note however that LN is still thermally stimulated and reaches observable
values for experimental tests within this regime.