# Analyzing and Improving Fault Tolerance of Learning-Based Navigation Systems

Zishen Wan1∗, Aqeel Anwar1, Yu-Shun Hsiao2, Tianyu Jia2 and  Vijay Janapa
Reddi2, Arijit Raychowdhury1 1Georgia Institute of Technology, Atlanta, GA
2Harvard University, Cambridge, MA *Email<EMAIL_ADDRESS>
###### Abstract.

Learning-based navigation systems are widely used in autonomous applications,
such as robotics, unmanned vehicles and drones. Specialized hardware
accelerators have been proposed for high-performance and energy-efficiency for
such navigational tasks. However, transient and permanent faults are
increasing in hardware systems and can catastrophically violate tasks safety.
Meanwhile, traditional redundancy-based protection methods are challenging to
deploy on resource-constrained edge applications. In this paper, we
experimentally evaluate the resilience of navigation systems with respect to
algorithms, fault models and data types from both RL training and inference.
We further propose two efficient fault mitigation techniques that achieve
2$\times$ success rate and 39% quality-of-flight improvement in learning-based
navigation systems.

##  1\. Introduction

Autonomous navigation continues to be deployed and attains significant
traction in the field of robotics, unmanned drones and autonomous vehicles at
all computing scales (liu2021robotic, ; wan2021survey, ; suleiman2019navion, ;
gao2021ielas, ; krishnan2021autosoc, ; wan2021energy, ). It helps an agent
avoid navigating in unknown environments and situations; and take safe actions
as needed. Recently, end-to-end learning-based techniques (anwar2018navren, ;
krishnan2019air, ; duisterhof2019learning, ) have demonstrated considerable
potential in navigation along with specialized hardware accelerators
(amaravati201855, ; krishnan2021machine, ), where the agent processes raw
sensor information and uses policy trained by reinforcement learning (RL) to
directly produce output actions.

While various approaches to achieve energy-efficient and accurate navigation
have been extensively studied, the inherent resilience of learning-based
navigation systems to hardware faults is not well understood. Hardware
accelerators are susceptible to permanent and transient faults, as a
consequence of manufacturing defects or single event upset. The increasing
rate of faults in technology nodes can cause application failures and safety
violations.

Faults may greatly impact the safety of navigation systems in both training
and inference. In the context of training, different from conventional offline
training and deployment method (i.e., train the policy on cloud and conduct
inference on edge), the learning-based system usually requires real-time
training and fine-tuning on edge nodes because the discrepancy between offline
simulated and new unknown environments will degrade the policy’s performance
(anwar2020autonomous, ). Faults occurred during training might impact the
learning process and even negatively affect the policy convergence due to the
complexity of learning algorithms and optimization procedures. In the context
of inference, different from DNN models in supervised learning, RL policy’s
quality depends on how effective it is in long-term sequential decision
making. Faults at one stage might propagate to subsequent stages, suggesting
that policies might be more vulnerable to faults than traditional DNN
applications.

Traditional methods to protect compute systems from hardware faults typically
adopt ECC (tawada2015bit, ) or replicate hardware components (matsuo2018dual,
; bannon2019computer, ). While these approaches are effective, they bring
large overhead in the hardware cost and energy which are challenging to deploy
on edge devices. In addition, this level of protection imposes a high design
effort at both micro-architecture and device level, which may be an overkill
for learning-based navigation policies.

In this paper, we perform an in-depth fault characterization study on the
learning-based navigation system from small (grid-based navigation task) to
large (drone navigation task) computing scales. The resilience against
permanent and transient faults is explored in both training and inference
stages. Our results demonstrate that different learning algorithms exhibit
different sensitivities to faults. We disclose the relationship between fault
type, bit error rate, policy convergence and exploration rate in the learning-
based system. Moreover, the system resilience depends on the model topology,
data type and faulty bit location. Based on our observations, we devise two
cost-effective solutions to mitigate faults impact: adjustment of the ratio of
exploration-to-exploitation during RL training and range-based anomaly
detection during inference.

This paper, therefore, makes the following contributions:

  * •

We develop a fault injection tool-chain that emulates various types of
hardware faults and enables rapid fault analysis of learning-based navigation
systems.

  * •

We conduct a large-scale fault injection study in both training and inference
stages of learning-based systems against permanent and transient faults. We
show how faults can impact policy convergence and exploration with different
training algorithms and evaluate the system resilience dependent on model
structure, data type and faulty bit position. To the best of our knowledge,
this is the first such comprehensive study of hardware faults in RL.

  * •

We present two low-overhead techniques to mitigate transient and permanent
fault impacts in both RL training and inference. By dynamically adjusting
exploration rate in training and detecting anomaly values in inference, the
learning-based system achieves up to 2$\times$ resilience improvement.

##  2\. Related Work

Reliability of autonomous system. The resilience of navigation system has been
recently studied in the autonomous vehicle field. (toschi2019characterizing, )
assesses the noise in sensor, while we focus on the hardware faults during
compute. (jha2019ml, ; hsiao2021mavfi, ) evaluate hardware faults on the
traditional kernel-based pipeline consisting of several function modules such
as perception, localization and planning. However, the reliability of the
emerging end-to-end learning-based navigation system has not been adequately
explored, which is the focus of our work.

Fault characterization. The impact of faults in neural network, a core
component in the learning-based navigation system, have been studied in a wide
body of recent work. (reagen2018ares, ; mahmoud2020pytorchfi, ;
chen2020tensorfi, ) build fault injection frameworks to quantify the error
resilience of DNN applications. (li2017understanding, ) evaluated the fault
propagation of DNN focused on the vulnerability of different layers in an ASIC
model. However, these works exclusively focus on the resiliency of neural
networks to hardware faults without considering the inherent resilience for
end-to-end learning-based safety-critical systems. Unlike a single neural
network model in supervised learning, the resiliency of models (i.e.,
policies) in the learning-based system depends on the sequential decision-
making process. Therefore, we take the first step to explore how hardware
faults affect the training process and long-term decision-making capability of
learning-based systems.

Fault mitigation techniques. Several general techniques have been proposed to
mitigate faults impacts, e.g., DMR (matsuo2018dual, ), TMR (hudson2018fault, )
and ECC (tawada2015bit, ). Some NN-specified techniques, such as explicit
redundancy and retraining, are proposed in (torres2017fault, ). These
techniques can be effective for fault detection and recovery, but incur area,
power or timing overhead which is challenging to deploy in resource-limited
edge nodes. Instead of costly fault tolerance operations, we present
application-aware fault mitigation techniques while considering the
characteristics of the learning-based system and bit significance. The
proposed techniques don’t require redundant bits and can be applied to both RL
training and inference stages.

##  3\. Methodology and Fault Model

###  3.1. Learning-Based Navigational System

We consider the task of reinforcement learning based navigational system
(anwar2018navren, ). The goal of such system is to learn a navigational policy
by interacting with the underlying environment achieving user-defined goals.
We learn these navigational policies with the help of RL where the agent
constantly interacts with the environment $\mathcal{E}$ characterized by a
Markov Decision Process (MDP). The MDP of the environment can be described by
the tuple
$\mathcal{M}=(\mathcal{S},\mathcal{A},\mathcal{P},\mathcal{R},\gamma)$ where
$\mathcal{S}$ is the state space, $\mathcal{A}$ is the action space,
$\mathcal{P}:\mathcal{S}\times\mathcal{A}\rightarrow\mathcal{S}$ is the MDP
transition probabilities,
$\mathcal{R}:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$ is the reward
function, and $\gamma\in(0,1)$ is the discount factor. With each interaction
$i$ with the environment $\mathcal{E}$, the agent observes the data tuple
$\mathcal{D}_{i}=(s_{i},a_{i},s_{i+1},r_{i})$ where
$s_{i},s_{i+1}\in\mathcal{S}$ is the current and next state,
$a_{i}\in\mathcal{A}$ is the action taken at state $s_{i}$, and
$r_{i}=\mathcal{R}(s_{i},a_{i})$ is the reward obtained. Given a policy $\pi$,
the value function (for a given state) and the Q-function (for a given state-
action pair) is given by

(1) $\displaystyle
V^{\pi}\left(s\right)=\mathbb{E}\left[\sum_{k=0}^{\infty}\gamma^{k}\mathcal{R}\left(s_{k},a_{k}\right)\mid
s_{0}=s\right]\quad,a_{k}\sim\pi(\cdot|s_{k})$ (2) $\displaystyle
Q^{\pi}\left(s,a\right)=\mathbb{E}\left[\sum_{k=0}^{\infty}\gamma^{k}\mathcal{R}\left(s_{k},a_{k}\right)\mid
s_{0}=s,a_{0}=a\right]$

We denote by $\rho$ the initial state distribution of state-space
$\mathcal{S}$. The long-term discounted return under $\rho$ can be written as
$V^{\pi}\left(\rho\right)=\mathbb{E}_{s\sim\rho}V^{\pi}\left(s\right).$ The
goal of RL is to learn an optimal policy $\pi^{*}$ given the observed data
tuples $\mathcal{D}_{i}$ that maximizes this return i.e.

(3) $\pi^{*}=\mathop{\arg\\!\max}_{\pi}V^{\pi}(\rho)$

We use Q-learning based RL where bellman backup operator is used to update the
Q-function until it converges.

(4) $\displaystyle
Q\left(s_{i},a_{i}\right)\leftarrow\mathbb{E}_{s^{\prime}\sim\mathcal{P}(\cdot\mid
s_{i},a_{i})}\left[\mathcal{R}(s_{i},a_{i})+\gamma\max_{a^{\prime}}Q\left(s_{i+1},a^{\prime}\right)\right]$

and the underlying policy at any iteration is given by

(5) $\displaystyle a_{k}=\mathop{\arg\\!\max}_{a^{\prime}}Q(s_{k},a^{\prime})$

The above policy converges to an optimal policy $\pi^{*}$ (deterministic) when
the Q-function converges under the application of the Bellman backup operator.
In order to estimate the Q-function, we consider both the tabular method (for
simpler problems) and neural network-based function approximation (for complex
problems). In tabular methods, we maintain a Q table of size
$|\mathcal{S}|\times|\mathcal{A}|$. For more complex problems with a larger
state/action space, we use neural network-based function approximation
$f_{\theta}:\mathcal{S}\rightarrow\mathcal{A}$ to estimate the Q-function
parameterized by $\theta$, where $\theta$ are the trainable weights of the
pre-defined neural network structure.

###  3.2. Fault Model

(a) Low obstacle density

(b) Middle obstacle density

(c) High obstacle density

Figure 1. Grid World problems with various obstacle densities. Red, green and
dark blue cells are obstacles, goal and agent respectively. Light blue path is
the mostly followed route by agent after training.

Fault type. We consider both permanent and transient faults in this paper.
Permanent faults usually result from manufacturing defects at the transistor
and interconnected structures in the integrated circuit and are continuous and
stable with time. They can cause bits to be held exclusively low (i.e., stuck-
at-0) or high (i.e., stuck-at-1), and introduce bit-flips if the logical value
of faulty bits differ from the value where the stuck happens. Transient faults
(i.e., soft error) mainly arise from external perturbations, such as particle
strikes or voltage violations, and may only exist for a short period. They
manifest as bit-flips in the computing and storage elements and result in
incorrect random values. We use widely adopted stuck-at and random bit-flip
fault models as abstractions of permanent and transient physical defect
mechanisms (salami2018resilience, ).

Fault location. We consider faults in the memory, which are pertinent to the
learning-based system due to its large intermediate states and values storage
requirement. For tabular-based policy, we consider data buffer storing tabular
values. For neural network-based policy, we consider input buffer, filter
buffer and output buffer which are used to store feature maps, weights and
activations, respectively. Faults can also happen in the datapath, and we
assume that faults in MAC units end up corrupting the operation outputs and
are manifested as random values in the output buffer. We do not consider
faults in control logic because the scheduling is mainly done by the host
(i.e., CPU). The fault model is derived from (reagen2018ares, ;
li2017understanding, ) and interested readers are pointed to the references
for details.

###  3.3. Fault Injection Method

Fault injection is performed in two modes: static and dynamic. Static faults
are injected offline before training and before inference is executed. Dynamic
injection is conducted during the execution, and the overheads are minimized
by implementing the fault model as tensor operations. This method is in line
with (reagen2018ares, ) and its accuracy has been validated on silicon
(whatmough201714, ). Permanent faults are injected statically since they are
independent of execution. Transient faults in weights during inference are
static injection since they are known after training and can be manipulated
before inference begins. In contrast, transient faults in activations are
injected dynamically since they are input-dependent.

##  4\. Fault Characterization

In this section, we evaluate the impact of faults on two learning-based
template problems: Grid World (Sec. 4.1) and drone autonomous navigation (Sec.
4.2).

###  4.1. Grid World Navigation Problem

We begin our experimentation with a simpler but popular navigation problem of
Grid World. The environment consists of a $n\times n$ grid. Each cell in the
grid is marked as either of the four categories: {source, goal, hell, free}.
The agent is placed at the source and is required to reach the goal and avoid
getting trapped in hell. Each cell corresponds to a unique state
($|\mathcal{S}|=n^{2}$). At each occupied cell, the agent samples from one of
the four actions from its action space $\mathcal{A}$: {move-up, move-down,
move-left, move-right} ($|\mathcal{A}|=4$) and moves to a new cell. The reward
is 1, -1 or 0 if the new cell is a goal, hell or free respectively. The
performance of policy is quantified by the agent’s success-rate which is
defined as the ratio of the number of times the agent reached the goal to the
total number of trials and cumulative reward (sum of rewards in an episode).
We consider three Grid World settings with various obstacle densities, as
shown in Fig. 1. Widely used tabular-based and neural network (NN)-based
methods (Sec. 3.1) are adopted and compared. The policies are quantized to
8-bit without loss of performance during the training and inference given the
memory and power constraints of edge devices. We repeat each fault injection
campaign 1000 times, which can lead to 95% confidence-level within 1% error
margin. Due to space limit, we only show results for Fig. 1b. We however
confirm that similar observations apply to both Fig. 1a and 1c.

####  4.1.1. Training in Grid World

Transient faults in training. Fig. 2a and Fig. 2c show the impact of transient
faults on both tabular and NN-based methods during training. Bit-flips are
injected in a single random episode with different bit error rates (BERs). It
is well observed that transient faults occurred in earlier episodes with low
BERs have a negligible impact on the agent’s final success rate (color-coded).
However, when the BER and fault injection episode are introduced beyond a
certain threshold, the performance and accuracy degrades quickly.

We find that with low BER, the agent can learn from these faults during
training and generate a good policy. In a few cases, policy trained with bit-
flips even yields a higher success rate than runs without faults. We believe
that these bit-flips can be seen as injecting an amount of noise that was
small enough to maintain a good policy and large enough to improve exploration
which benefits the convergence. This is akin to regularization in RL and
previous work in quantization report similar effect (krishnan2019quantized, ).
However, with high BER, the agent fails to learn from these faults which
decreases the efficacy of the policy and can even impact the convergence.

Another takeaway from Fig. 2a and Fig. 2c is that NN-based approach exhibits
higher resilience to transient faults than the tabular approach. We
hypothesize the following two reasons: NN-based training process is not
significantly impacted by transient faults and it can quickly recover from
faults. Fig. 3 shows that the latter case contributes more to the policy
resiliency. After bit-flips get injected (blue and orange curves), the
cumulative rewards of both tabular and NN-based training drop quickly,
indicating that faults destroy existing learned experience. However, the NN-
based approach recovers from faults quickly after training several episodes
due to its advanced optimization techniques, whereas the tabular approach
takes a longer time to recover or even fails to converge at the end of the
training.

(a) Tabular-based approach

(b) Tabular value histogram

(c) NN-based approach

(d) NN weight histogram

Figure 2. (a)(c) The impact of transient and permanent faults on Grid World
training for both tabular and NN-based approaches. The darker, the worse of
performance. (b)(d) The histogram of tabular values and NN parameters.

Permanent faults in training. Permanent fault is stuck at 0 or 1 at the same
memory location throughout training, and its impacts are shown in Fig. 2a and
Fig. 2c. We find that despite its permanent nature, these fault have
comparable impact as transient faults on the policies. This indicates that the
agent can learn and adapt itself to the permanent fault patterns in training.
For example, if bits are stuck at 0 and affect a cell which will be passed
with the uncontaminated policy, the agent will gradually learn to avoid this
cell and find a new path. Due to system’s intrinsic long-term learning
capability, permanent faults do not corrupt the learned policies as much as
transient faults that occur in the latter training episodes.

Interestingly, stuck-at-1 has a slightly higher impact than stuck-at-0 on
policy performance for tabular-based approach, as shown in Fig. 2a. This can
be explained from Fig. 2b where tabular values contain more 0 bits than 1
bits, therefore stuck-at-1 bring more faulty bits in the system. But due to
the tabular values’ characteristics (not centered at 0 and a wide range) and
cell’s functionality (most cells don’t contribute to the final path and one
cell only affects its adjacents), the distinction between the impact of stuck-
at-0 and stuck-at-1 faults is small.

Different from tabular method, stuck-at-0 and stuck-at-1 exhibit significantly
divergent impact on NN-based methods. As shown in Fig. 2c, stuck-at-0 has a
negligible impact within 0.7% BER, while stuck-at-1 noticeably fails the
training even with 0.2% BER. This observation is due to the sparsity and
limited range of NN data as visualized in Fig. 2d. At a bit level, there are
more 0s than 1 bits in the weights by 7.17$\times$, whereas this ratio is only
3.18$\times$ for tabular values. Further, most NN model data is centered
around 0 and have a much narrower range than tabular values. From NN
characteristics, weights with higher magnitude play more critical roles in
final output. Therefore, the greater magnitude faulty values arose from stuck-
at-1 have more impacts on the system than stuck-at-0.

(a) Tabular-based approach

(b) NN-based approach

Figure 3. Exampled cumulative return during training under transient and
permanent faults. BER represents bit error rate (%). EI represents faults
injected episode index.

The training processes under permanent faults are demonstrated in Fig. 3
(green and red curves). The agent can learn across stuck-at-faults and finally
converge itself although it takes a longer time in the beginning phase of
training, except for stuck-at-1 in the NN-based method which destroys the
convergence of system.

(a) Transient faults in Tabular-based training

(b) Permanent faults in Tabular-based training

(c) Transient faults in NN-based training

(d) Permanent faults in NN-based training

Figure 4. (a)(c) Episodes taken to converge ($>$95% success rate) after
transient faults injected in 900th episode. (b)(d) The policy’s success rate
after training 1000 more episodes under permanent faults. EI represents faults
injected episode index.

Convergence. We now ask the following question: whether the policy can finally
converge after faults are injected? Answering this question will enable us to
further understand the system’s long-term learning capability on different
types of faults. Fig. 4a and Fig. 4c demonstrate that both tabular and NN-
based approach can finally achieve convergence after transient faults
injected, although tabular method needs 2$\times$ episodes than NN method.
This further implicates the self-healing ability and higher resiliency to
transient faults of NN-based learning system. In contrast, Fig. 4b and Fig. 4d
indicate that extra training doesn’t bring obvious improvements for both
stuck-at-0 and stuck-at-1 when BER surpasses a threshold, and may even degrade
the performance due to over-training. This is because a steady low exploration
rate impedes the agent learn the permanent fault, which we will discuss more
in Sec. 5.

####  4.1.2. Inference in Grid World

Transient faults in inference. Inference of learning-based system is a
sequential decision-making process. The trained policy will be utilized in
each action step based on the agent’s current state. Hence transient faults
have two modes: occur in the read register which only impacts the current
single step (Transient-1) or occur in memory which impacts the whole
executions (Transient-M). As shown in Fig. 5, Transient-1 has a negligible
effect on navigation success rate compared to transient-M fault, this is
because a single wrong decision doesn’t necessarily results in task failure
and can be remedied in the following actions of learning-based system.

Permanent faults in inference. Permanent faults are stuck at 0 or 1 during
whole inference steps. For tabular-based inference, stuck-at-0 and stuck-at-1
have similar impacts. For NN-based inference, however, stuck-at-1 has a much
severe impact on policy than stuck-at-0 due to the sparsity nature of NN
parameters.

(a) Tabular-based inference

(b) NN-based inference

Figure 5. The impact of transient and permanents faults on Grid World
inference for both tabular and NN-based approaches.

(a) PEDRA environments

(b) C3F2 neural network policy

Figure 6. Training details for the drone autonomous navigation task

(a) Drone navigation training

(b) Environment

(c) Data location

(d) NN layer

(e) Data type

Figure 7. Fault characterization on drone navigation tasks. (a) is training
result, (b)-(e) are inference results.

###  4.2. Drone Autonomous Navigation

In this section, we experiment on a more complex drone navigation problem in
3D realistic environments. Drone is a complex system (krishnan2020sky, ), and
we use PEDRA (anwar2020autonomous, ) as the simulation platform which is
powered by Unreal Engine. During the task, the drone is initialized at a start
location and required to navigate across the environment avoiding obstacles.
There is no destination position, and the drone is required to fly as long as
it can. The monocular image captured from the front-facing camera is taken as
the state of the RL problem. We use a perception-based probabilistic action
space $\mathcal{A}$ with 25 actions. The reward is designed to encourage the
drone to stay away from the obstacles. We use the indoor-long and indoor-
vanleer environments provided by PEDRA, which can be seen in Fig. 6a. To
quantify the policy performance, we use Mean Safe Flight (MSF) which is the
average distance traveled by the drone before collision.

####  4.2.1. Training in Drone Navigation Problem

Transient faults in training. The model (Fig. 6b) is first trained offline
using Double DQN with experience replay and then fined-tuned last two layers
online using transfer learning (anwar2020autonomous, ). Transient fault is
injected in a single random step with various BER during online training, and
its impact is shown in Fig. 7a. We can observe that the flight quality
degrades with higher BER and latter fault occurrence, which is the same trend
as Sec. 4.1.1.

Permanent faults in training. Permanent faults exist throughout online
training. Fig. 7a demonstrates that stuck-at-1 faults (BER=$10^{-3}$) severely
affect the drone’s safe flight distance, while stuck-at-0’s impact is
moderate. These are the same trends shown in Sec. 4.1.1, so we omit the
explanation here due to space restriction.

####  4.2.2. Inference in Drone Navigation Problem

Different environments. We first compare the policy resilience between two
environments, shown in Fig. 7b. Transient faults are injected in weights once
the policy is trained. Experiments are repeated 100 times for each case. We
observe that faults have similar impact trends on both indoor-long and indoor-
vanleer, which brings more flight quality degradation with higher BER. Due to
this little difference, we focus on indoor-long for the rest of Fig. 7.

Different fault locations. The fault sensitivity of different data locations
(i.e., feature map, weights and activations) is evaluated in Fig. 7c. Weight
buffer is read-only once the policy is trained, while input and activation
buffer will be rewritten during execution. Fig. 7c shows that input buffer is
relatively immune to transient faults since bit-flips only affect values in a
single frame and image information has redundancy. Transient faults occurred
in ReLU activation (Act (T)) bring more flight distance drop than inputs
because an activation faulty value will be propagated to all neurons in the
next level without any masking. Permanent faults occurred in activations (Act
(P)) destroy the task since all action steps will be greatly affected
throughout the flight. Weights are more sensitive to transient faults than
activations (Act (T)) due to its massive reuse in convolutional layers. A
fault in filter will be multiplied across whole feature map, therefore impacts
larger value subsets.

Different layer types. The sensitivity of different positions and types of
layers is evaluated in Fig. 7d by injecting bit-flips into each layer
individually. We find that the first two convolutional layers exhibit higher
resiliency to transient faults. The reason is max-pooling and ReLU activation
are implemented at the end of the first two layers, which can discard the
faulty values and mask errors. However, there’s no pooling layer after the
third convolutional layer. Therefore, latter layers are generally more
vulnerable to faults since faults have a relatively lower probability of being
masked. The second fully-connected layer is the most vulnerable to bit-flips.
This is because faults can spread to all neurons right after occurrence, and
FC2 is the last layer that directly dictates the drone actions.

Different data types. Fig. 7e demonstrates the fault tolerance of three fixed-
point data types: Q (1,4,11), Q (1,7,8) and Q (1,10,5) (Q (sign, integer,
fraction)) with the same flight quality baseline. We observe that Q (1,4,11)
consistently exhibit higher resilience than Q (1,7,8) and Q (1,10,5). This is
because Q (1,4,11) has a narrower range compared to others, and are already
able to capture the range of weights without performance loss. Hence, a bit-
flip in MSB part of Q (1,10,5) is likely to result in a larger deviation from
zero and higher impacts. It is noted that to achieve higher resiliency, data
types should be able to optimally capture the value range rather than pursuing
an unnecessarily large range. This is in line with some recent works
(tambe2019adaptivfloat, ; langroudi2020adaptive, ; tambe2020algorithm, ).

##  5\. Fault Mitigation Techniques

In consideration of the significant impact of faults, we propose and evaluate
two fault mitigation schemes for training (Sec. 5.1) and inference (Sec. 5.1),
relying on our characterization observations.

###  5.1. Training: Adaptive Exploration Rate Adjustment

Detection. Based on our observations in Sec. 4, transient faults can bring a
sudden drop in reward and permanent faults can result in a continuous low
reward. Therefore, we use the change in cumulative reward as an indicator of
fault occurrence during training. For transient fault detection, if the
cumulative reward drop exceeds $x$% within $y$ consecutive episodes (where $x$
and $y$ are parameters), then it is supposed that fault occurs in the system.
For permanent fault detection, after the exploration rate drops to a constant
value and the agent goes into a stable greedy exploitation state, if the
reward is still lower than 50% of maximum value, then the permanent fault is
supposed to be detected.

Recovery. One of fundamental trade-offs impacting the RL training process is
the exploration and exploitation ratio, which determines how agents explore
the action space and exploit the knowledge at any given time step. Decaying-
epsilon-greedy strategy is commonly used where the agent performs more
explorations in the beginning and then explores less over time until a steady
rate. This can make the policy converge fast and facilitate the learning
process.

However, once faults occur, the agent needs more exploration actions to avoid
being stuck in the wrong states and adapt itself to fault pattern. Therefore,
we propose a dynamic exploration rate adjustment scheme to mitigate the
negative fault impact during training. If transient faults are detected, the
system will automatically increase its exploration rate (ER) followed by

(6) $\displaystyle ER_{new}=ER_{old}+\delta(ER)=ER_{old}+\alpha\times
min(f{(r)},f{(r)}f(t))$

where $f{(r)}$ is normalized reward drop $f{(r)}=\Delta{r}/r_{max}$. $f{(t)}$
characterizes fault occurrence time $t$ defined as $f{(t)}=t/T$ ($T$ is
episodes taken to reach steady exploitation in normal training). $ER_{new}$ is
the adjusted exploration rate, and $\alpha$ is adjustment coefficient.

The agent needs more episodes to learn permanent fault patterns. Hence when
permanent faults are detected, the agent will autonomously revert the
exploration rate to initial and slow down its decreasing speed by
$2^{n}\times$. This detection and recovery process will be triggered whenever
the agent goes into the steady exploitation state, and $n$ is the number of
detection times.

Evaluation. We measure the effectiveness of the proposed dynamic exploration
rate adjustment technique by repeating experiments in Fig. 2 where bit-flips
are randomly injected with different BERs during training. We choose $x$=25,
$y$=50 for detection, and $\alpha$=0.8 and $\alpha$=0.4 for tabular and NN-
based methods respectively due to their self-healing ability difference shown
in Sec. 4.1. The agent takes $T$=100 episodes to reach the steady exploitation
state.

(a) Tabular-based approach

(b) NN-based approach

Figure 8. The fault mitigation effect of dynamic exploration adjustment scheme
on training.

(a) Tabular-based approach

(b) NN-based approach

(c) Recovery Speed

Figure 9. The correlation between bit error rate, exploration ratio and
transient fault recovery speed during training.

The evaluation results are shown in Fig. 8. Compared with Fig. 2a and Fig. 2c,
we find that almost all transient faults happened before 800 episodes become
benign to the system, and the impact of faults happened in the last 200
episodes has also been greatly mitigated. The negative impact of permanent
faults is also relieved by 10%.

We also disclose the correlation between bit error rate, adjusted exploration
rate and fault recovery speed. As shown in Fig. 9a and Fig. 9b, the agent
automatically adjusts to a higher exploration rate when more transient faults
detected (red), and adopts lower exploration decreasing speed when more
permanent faults exist (blue). Note that if stuck-at-1 exists in NN-based
system, the agent tends to conduct much more exploration to prevent itself
stuck in the wrong states. Fig. 9c demonstrates that agent usually needs more
episodes to converge back if adjusting to a higher exploration rate. This is
reasonable because if agent chooses low exploration rates under high BER
scenarios, it may fail to recover from faults. The trade-off between recovery
speed and general recovery success rate can be dynamically determined by agent
in our proposed scheme.

(a) Grid World problem

(b) Drone navigation problem

Figure 10. The effectiveness of anomaly detection on inference.

###  5.2. Inference: Range-Based Anomaly Detection

Detection. To detect unusual values, we adopt statistically anomaly detection
technique that is widely used in safety-critical systems
(rajbahadur2018survey, ). Note that we perform value-level fault detection
rather than conventional bit-level mismatch detection due to the observation
that many bit-flips may eventually get masked or only result in small value
deviations which won’t impact final performance.

Based on the finding that faulty values with higher magnitude have a stronger
destructive impact, we use the value range of parameters as an indicator to
detect outliers. Once the policy is trained, the value ranges of each layer
will be instrumented and derived (i.e. ($a_{i},b_{i}$)), and 10% detection
margin will be applied on bound (i.e. ($1.1a_{i},1.1b_{i}$)) during inference.
We only compare values composed of sign and integer bits based on the
observation that fractional part has little impact, which can reduce the
hardware cost further.

Recovery. Once any anomaly data is detected, an alarm signal will be raised,
and the operations around this data will be skipped. The intuition is the
sparsity nature of NN parameters, and the data with small magnitude has a
higher possibility to become an outlier after bit-flip due to two’s complement
encoding method.

Evaluation. The effectiveness of our proposed anomaly detection scheme is
evaluated in both Grid World and drone navigation problems, shown in Fig. 10.
Transient faults are injected in NN weights. Compared to no mitigation cases,
the agent’s success rate increases by 2x in Grid World and drone’s flight
quality increases by 39% after the fault mitigation is applied. Moreover, its
runtime overhead is less than 3% without requiring redundant bits as opposed
to ECC.

##  6\. Conclusion

This paper empirically investigates the end-to-end learning-based navigation
system from the fault tolerance perspective, i.e., fault characterization and
mitigation. We characterize the impact of transient and permanent faults on
both training and inference stages through large-scale fault injection
experiments on various training algorithms, data types and model structures.
Based on the characterizations, we further propose two cost-effective fault
mitigation techniques by dynamically adjusting the exploration rate during
training and detecting anomaly ranges during inference. Our proposed method
can improve system resilience by up to 2$\times$.

## Acknowledgements

We thank Muya Chang (from Georgia Tech) for his technical support. This work
was supported in part by C-BRIC and ADA, two of six centers in JUMP, a
Semiconductor Research Corporation (SRC) program sponsored by DARPA.

## References

  * (1) S. Liu, Z. Wan, B. Yu, and Y. Wang, “Robotic computing on fpgas,” Synthesis Lectures on Computer Architecture, vol. 16, no. 1, pp. 1–218, 2021\.
  * (2) Z. Wan, B. Yu, T. Y. Li, J. Tang, Y. Zhu, Y. Wang, A. Raychowdhury, and S. Liu, “A survey of fpga-based robotic computing,” IEEE Circuits and Systems Magazine, vol. 21, no. 2, pp. 48–74, 2021.
  * (3) A. Suleiman, Z. Zhang, L. Carlone, S. Karaman, and V. Sze, “Navion: A 2-mw fully integrated real-time visual-inertial odometry accelerator for autonomous navigation of nano drones,” IEEE Journal of Solid-State Circuits, vol. 54, no. 4, pp. 1106–1119, 2019.
  * (4) T. Gao, Z. Wan, Y. Zhang, B. Yu, Y. Zhang, S. Liu, and A. Raychowdhury, “ielas: An elas-based energy-efficient accelerator for real-time stereo matching on fpga platform,” in 2021 IEEE 3rd International Conference on Artificial Intelligence Circuits and Systems (AICAS), pp. 1–4, IEEE, 2021\.
  * (5) S. Krishnan, T. Tambe, Z. Wan, and V. J. Reddi, “Autosoc: Automating algorithm-soc co-design for aerial robots,” arXiv preprint arXiv:2109.05683, 2021.
  * (6) Z. Wan, Y. Zhang, A. Raychowdhury, B. Yu, Y. Zhang, and S. Liu, “An energy-efficient quad-camera visual system for autonomous machines on fpga platform,” in 2021 IEEE 3rd International Conference on Artificial Intelligence Circuits and Systems (AICAS), pp. 1–4, IEEE, 2021.
  * (7) M. A. Anwar and A. Raychowdhury, “Navren-rl: Learning to fly in real environment via end-to-end deep reinforcement learning using monocular images,” in 2018 25th International Conference on Mechatronics and Machine Vision in Practice (M2VIP), pp. 1–6, IEEE, 2018.
  * (8) S. Krishnan, B. Borojerdian, W. Fu, A. Faust, and V. J. Reddi, “Air learning: An ai research platform for algorithm-hardware benchmarking of autonomous aerial robots,” arXiv preprint arXiv:1906.00421, 2019.
  * (9) B. P. Duisterhof, S. Krishnan, J. J. Cruz, C. R. Banbury, W. Fu, A. Faust, G. C. de Croon, and V. J. Reddi, “Learning to seek: Autonomous source seeking with deep reinforcement learning onboard a nano drone microcontroller,” arXiv preprint arXiv:1909.11236, 2019.
  * (10) A. Amaravati, S. B. Nasir, J. Ting, I. Yoon, and A. Raychowdhury, “A 55-nm, 1.0–0.4 v, 1.25-pj/mac time-domain mixed-signal neuromorphic accelerator with stochastic synapses for reinforcement learning in autonomous mobile robots,” IEEE Journal of Solid-State Circuits, vol. 54, no. 1, pp. 75–87, 2018.
  * (11) S. Krishnan, Z. Wan, K. Bharadwaj, P. Whatmough, A. Faust, S. Neuman, G.-Y. Wei, D. Brooks, and V. J. Reddi, “Machine learning-based automated design space exploration for autonomous aerial robots,” arXiv preprint arXiv:2102.02988, 2021.
  * (12) A. Anwar and A. Raychowdhury, “Autonomous navigation via deep reinforcement learning for resource constraint edge nodes using transfer learning,” IEEE Access, vol. 8, pp. 26549–26560, 2020.
  * (13) M. Tawada, S. Kimura, M. Yanagisawa, and N. Togawa, “A bit-write reduction method based on error-correcting codes for non-volatile memories,” in The 20th Asia and South Pacific Design Automation Conference, pp. 496–501, IEEE, 2015.
  * (14) I. B. M. Matsuo, L. Zhao, and W.-J. Lee, “A dual modular redundancy scheme for cpu–fpga platform-based systems,” IEEE Transactions on Industry Applications, vol. 54, no. 6, pp. 5621–5629, 2018.
  * (15) P. Bannon, G. Venkataramanan, D. D. Sarma, and E. Talpes, “Computer and redundancy solution for the full self-driving computer,” in 2019 IEEE Hot Chips 31 Symposium (HCS), pp. 1–22, IEEE Computer Society, 2019.
  * (16) A. Toschi, M. Sanic, J. Leng, Q. Chen, C. Wang, and M. Guo, “Characterizing perception module performance and robustness in production-scale autonomous driving system,” in IFIP International Conference on Network and Parallel Computing, pp. 235–247, Springer, 2019.
  * (17) S. Jha, S. Banerjee, T. Tsai, S. K. Hari, M. B. Sullivan, Z. T. Kalbarczyk, S. W. Keckler, and R. K. Iyer, “Ml-based fault injection for autonomous vehicles: A case for bayesian fault injection,” in 2019 49th annual IEEE/IFIP international conference on dependable systems and networks (DSN), pp. 112–124, IEEE, 2019.
  * (18) Y.-S. Hsiao, Z. Wan, T. Jia, R. Ghosal, A. Raychowdhury, D. Brooks, G.-Y. Wei, and V. J. Reddi, “Mavfi: An end-to-end fault analysis framework with anomaly detection and recovery for micro aerial vehicles,” arXiv preprint arXiv:2105.12882, 2021.
  * (19) B. Reagen, U. Gupta, L. Pentecost, P. Whatmough, S. K. Lee, N. Mulholland, D. Brooks, and G.-Y. Wei, “Ares: A framework for quantifying the resilience of deep neural networks,” in 2018 55th ACM/ESDA/IEEE Design Automation Conference (DAC), pp. 1–6, IEEE, 2018.
  * (20) A. Mahmoud, N. Aggarwal, A. Nobbe, J. R. S. Vicarte, S. V. Adve, C. W. Fletcher, I. Frosio, and S. K. S. Hari, “Pytorchfi: A runtime perturbation tool for dnns,” in 2020 50th Annual IEEE/IFIP International Conference on Dependable Systems and Networks Workshops (DSN-W), pp. 25–31, IEEE, 2020\.
  * (21) Z. Chen, N. Narayanan, B. Fang, G. Li, K. Pattabiraman, and N. DeBardeleben, “Tensorfi: A flexible fault injection framework for tensorflow applications,” in 2020 IEEE 31st International Symposium on Software Reliability Engineering (ISSRE), pp. 426–435, IEEE, 2020.
  * (22) G. Li, S. K. S. Hari, M. Sullivan, T. Tsai, K. Pattabiraman, J. Emer, and S. W. Keckler, “Understanding error propagation in deep learning neural network (dnn) accelerators and applications,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, pp. 1–12, 2017.
  * (23) S. Hudson, R. S. Sundar, and S. Koppu, “Fault control using triple modular redundancy (tmr),” in Progress in Computing, Analytics and Networking, pp. 471–480, Springer, 2018.
  * (24) C. Torres-Huitzil and B. Girau, “Fault and error tolerance in neural networks: A review,” IEEE Access, vol. 5, pp. 17322–17341, 2017.
  * (25) B. Salami, O. S. Unsal, and A. C. Kestelman, “On the resilience of rtl nn accelerators: Fault characterization and mitigation,” in 2018 30th International Symposium on Computer Architecture and High Performance Computing (SBAC-PAD), pp. 322–329, IEEE, 2018.
  * (26) P. N. Whatmough, S. K. Lee, H. Lee, S. Rama, D. Brooks, and G.-Y. Wei, “A 28nm soc with a 1.2 ghz 568nj/prediction sparse deep-neural-network engine with¿ 0.1 timing error rate tolerance for iot applications,” in 2017 IEEE International Solid-State Circuits Conference (ISSCC), pp. 242–243, IEEE, 2017\.
  * (27) S. Krishnan, S. Chitlangia, M. Lam, Z. Wan, A. Faust, and V. J. Reddi, “Quantized reinforcement learning (quarl),” arXiv preprint arXiv:1910.01055, 2019.
  * (28) S. Krishnan, Z. Wan, K. Bhardwaj, P. Whatmough, A. Faust, G.-Y. Wei, D. Brooks, and V. J. Reddi, “The sky is not the limit: A visual performance model for cyber-physical co-design in autonomous machines,” IEEE Computer Architecture Letters, vol. 19, no. 1, pp. 38–42, 2020.
  * (29) T. Tambe, E.-Y. Yang, Z. Wan, Y. Deng, V. J. Reddi, A. Rush, D. Brooks, and G.-Y. Wei, “Adaptivfloat: A floating-point based data type for resilient deep learning inference,” arXiv preprint arXiv:1909.13271, 2019.
  * (30) H. F. Langroudi, V. Karia, J. L. Gustafson, and D. Kudithipudi, “Adaptive posit: Parameter aware numerical format for deep learning inference on the edge,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops, pp. 726–727, 2020.
  * (31) T. Tambe, E.-Y. Yang, Z. Wan, Y. Deng, V. J. Reddi, A. Rush, D. Brooks, and G.-Y. Wei, “Algorithm-hardware co-design of adaptive floating-point encodings for resilient deep learning inference,” in 2020 57th ACM/IEEE Design Automation Conference (DAC), pp. 1–6, IEEE, 2020.
  * (32) G. K. Rajbahadur, A. J. Malton, A. Walenstein, and A. E. Hassan, “A survey of anomaly detection for connected vehicle cybersecurity and safety,” in 2018 IEEE Intelligent Vehicles Symposium (IV), pp. 421–426, IEEE, 2018.


# Virtual excitations and entanglement dynamics and polygamy
in three ultra-strongly coupled systems

Radouan Hab-arrih  Laboratory of Theoretical Physics, Faculty of Sciences,
Chouaïb Doukkali University, PO Box 20, 24000 El Jadida, Morocco  Ahmed Jellal
<EMAIL_ADDRESS>Laboratory of Theoretical Physics, Faculty of Sciences,
Chouaïb Doukkali University, PO Box 20, 24000 El Jadida, Morocco  Canadian
Quantum Research Center, 204-3002 32 Ave Vernon, BC V1T 2L7, Canada

###### Abstract

The Milburn dynamics of three ultra-strongly coupled oscillators are
investigated. We look at the Milburn dynamics of virtual excitations and how
they interact with pairwise entanglement. It is found that the dynamics of
excitations and entanglement experience similar profiles against time,
physical parameters, and decoherence rate. As well, we show that the
extinction of excitations entails separability. Additionally, we analyze the
effects of physical parameters on the redistribution of virtual excitations
among the three bi-partitions. As a result, we show the violation of the
monogamy of excitations as in quantum discord. This implies that excitations
can be considered as signatures of quantum correlations beyond entanglement.

###### pacs:

03.65.Fd, 03.65.Ge, 03.65.Ud, 03.67.Hk
Keywords: Two harmonic oscillators, magnetic coupling, Schmidt decomposition,
entanglement, dynamics.

###### pacs:

03.65.Fd, 03.65.Ge, 03.65.Ud, 03.67.Hk
Keywords:Ultra-strong coupling, oscillators, Milburn dynamics, virtual
excitations, entanglement, polygamy.

##  I Introduction

Ultra strong coupling (USC) physics, in which the coupling becomes comparable
to the system’s frequencies, has emerged as a new and increasingly popular
domain [1, 2, 3, 4, 5, 6, 7]. The importance of the USC lies in its potential
effects on several physical phenomena, including photon blockade [8], Purcell
effect [9], Zeno effect [10], and maintaining entanglement by virtual
excitations of a vacuum state [11, 12]. The fluctuations of the vacuum are one
of the most striking phenomena of modern physics and have colossal
applications in quantum technologies [13, 14]. The vacuum in USC systems
becomes populated by short-lived particles called virtual (not being absorbed
by detectors) excitations [1]. These excitations are responsible for mediating
forces (i.e., the Casimir effect) [15] and the transfer of energy between
oscillators and atoms [16]. Moreover, we point out that the known counter-
rotating terms (CRT) are responsible for simultaneous creation and
annihilation, as well as the amplification of excitations in a vacuum state.
As a result, the rotating wave approximation (RWA), where CRT is neglected,
breaks down [17, 2, 7].

Multipartite systems made of coupled oscillators are central in physics. In
fact, they lead to the modeling of coupled ions in ion traps [18], arrayed
coupled nano-sized electromechanical devices [19], light propagation in
inhomogeneous media [20], and a nitrogen vacancy ensemble embedded in a
diamond nanobeam [21]. Besides, the separability of a tripartite gaussian
state was studied in [22]. The stationary and non-stationary entanglement of
three oscillators were exhaustively studied in [23, 24], respectively.
Moreover, avoiding external decoherence and disentanglement in oscillator
systems was studied in [25, 26]. The interplay between coherence and
entanglement and their redistribution was discussed in [24, 27]. Additionally,
it was shown that with three coupled parametric oscillators, it is possible to
generate entanglement even at high temperatures [28]. Using the optimal
control theory, maximizing entanglement in coupled oscillators is discussed by
Stefanatos in [29].

Motivated by our last achievement with the Milburn dynamics of two ultra-
strongly coupled oscillators [30], showing the interplay between steering and
entanglement with virtual excitations. The current paper surveys the Milburn
dynamics of three ultra-strongly coupled oscillators. We harness the
covariance matrix formalism beyond RWA. Thus, the interconnection of
excitations with pairwise entanglement will be legitimately discussed in USC.
Additionally, this work sheds light on whether excitations are monogamous or
polygamous. The monogamy constraints are fundamental in quantum information
science [31] because they capture the quantumness of correlation. The
entanglement, in particular, is monogamous, i.e., $E(a|b)+E(a|c)\leq E(a|bc)$,
with $a,b$ and $c$ are representing three distinct parties. Nevertheless, it
is not the case for quantum discord [32]. By the way, we ask whether
excitations are monogamous or not. Addressing such questions of monogamy is
important to understanding the role of virtual excitations as signatures of
genuine non-classical correlations beyond entanglement.

The presented paper is formulated as follows. In Sec. II, we define the system
made of three oscillators as well as present the diagonalization scheme.
While, in Sec. III, we write out the unitary transformations in their
symplectic form and resolve the Milburn dynamics beyond the RWA. Sec. IV
investigates quantum entanglement and virtual excitations. In Sec. V, we
present our numerical results together with some discussions. Finally, we
conclude our work.

##  II Model and Diagonalization

We consider three interacting harmonic oscillators $(a,b,c)$, in which each
oscillator has its own angular frequency $\omega_{i}$ ($i=1,2,3$), as shown in
Fig. 1. The interaction between two oscillators $k$ and $l$ is carried out by
a time-independent constant $J_{kl}$, with $1\leq k<l\leq 3$.

Figure 1: (color online) Schematic diagram of three coupled harmonic
oscillators initially prepared in $\left|{000}\right\rangle$. Each of the two
oscillators $j$ and $k$ are coupled via ”position-position” $x_{j}x_{k}$ type
interaction with coupling strength $g_{jk}$, $(j<k)$. The dynamics of the
system is affected by the intrinsic decoherence quantified by the Milburn
decoherence parameter $\Gamma$.

The total Hamiltonian describing the isolated system is

$\displaystyle
H=\frac{p_{1}^{2}}{2}+\frac{p_{2}^{2}}{2}+\frac{p_{3}^{2}}{2}+\frac{1}{2}\omega_{1}^{2}x_{1}^{2}+\frac{1}{2}\omega_{2}^{2}x_{2}^{2}+\frac{1}{2}\omega_{3}^{2}x_{3}^{2}+J_{12}x_{1}x_{2}+J_{13}x_{1}x_{3}+J_{23}x_{2}x_{3}$
(1)

where the position $x_{i}$ and momentum $p_{i}$ operators are conjugate,
satisfying the commutation relations $\left[x_{k},p_{l}\right]=i\delta_{kl}$
and $\left[x_{k},x_{l}\right]=\left[p_{k},p_{l}\right]=0$. The masses will be
treated as a unit, and a simple transformation validates the assumption [33],
and setting $\hbar=1$. In matrix form, we have

$\displaystyle H$ $\displaystyle=$
$\displaystyle\frac{1}{2}\mathbb{P}^{t}\mathbb{P}+\frac{1}{2}\mathbb{X}^{t}\mathbb{V}\mathbb{X}.$
(2)

where $u^{t}$ stands for the transpose of the vector $u$,
$\mathbb{P}=(p_{1},p_{2},p_{3})^{t}$, $\mathbb{X}=(x_{1},x_{2},x_{3})^{t}$ and
$\mathbb{V}$ is the potential matrix given by

$\displaystyle\mathbb{V}=\left(\begin{array}[]{ccc}\omega_{1}^{2}&J_{12}&J_{13}\\\
J_{12}&\omega_{2}^{2}&J_{23}\\\
J_{13}&J_{23}&\omega_{3}^{2}\end{array}\right).$ (6)

Using the Sylvester criterion of the positivity of $\mathbb{V}$, we obtain a
manifold of physical parameters that gives a bound state. This is

$\displaystyle\mathcal{B}=\\{(\omega_{j},J_{jk})_{1\leq j<k\leq
3}\in\mathbb{C}^{6},\quad\max(\omega_{1}^{2},0)\times\max(\omega_{2}^{2}\omega_{3}^{2}-J_{23}^{2},0)\times\max(\det(\mathbb{V}),0)>0\\}.$
(7)

Now, to diagonalize the Hamiltonian using a unitary transformation, we address
the time-independent Euler unitary transformation
$\mathbb{T}(\alpha,\beta,\gamma)$ defined as [23, 24, 34]

$\displaystyle\mathbb{T}(\alpha,\beta,\gamma):=e^{-i\alpha\mathcal{J}_{3}}\times
e^{-i\beta\mathcal{J}_{2}}\times e^{-i\gamma\mathcal{J}_{3}}$ (8)

where $\mathcal{J}_{2,3}$ are the $y$ and $z$ components of the angular
momentum operator $\mathcal{J}$, respectively. By making the following choices
of the Euler angles, $\alpha,\beta$ and $\gamma$

$\displaystyle\cos(2\alpha)=2\left[\frac{\left(\Omega_{2}^{4}-(\omega_{1}^{2}+\omega_{2}^{2})\Omega_{2}^{2}+\omega_{1}^{2}\omega_{2}^{2}-J_{12}^{2}\right)\frac{\left(\Omega_{1}^{2}-\Omega_{3}^{2}\right)}{\left(\Omega_{3}^{2}-\Omega_{2}^{2}\right)}}{\left(\Omega_{1}^{4}-(\omega_{1}^{2}+\omega_{2}^{2})\Omega_{1}^{2}+\omega_{1}^{2}\omega_{2}^{2}-J_{12}^{2}\right)}+1\right]^{-1}-1$
(9)
$\displaystyle\cos(2\beta)=2\frac{\Omega_{3}^{4}-(\omega_{1}^{2}+\omega_{2}^{2})\Omega_{3}^{2}+\omega_{2}^{2}\omega_{1}^{2}-J_{12}^{2}}{(\Omega_{3}^{2}-\Omega_{1}^{2})(\Omega_{3}^{2}-\Omega_{2}^{2})}-1$
(10)
$\displaystyle\cos(2\gamma)=2\left[\frac{\Omega_{3}^{4}-(\omega_{1}^{2}+\omega_{3}^{2})\Omega_{3}^{2}+\omega_{1}^{2}\omega_{3}^{2}-J_{13}^{2}}{\Omega_{3}^{4}-(\omega_{2}^{2}+\omega_{3}^{2})\Omega_{3}^{2}+\omega_{2}^{2}\omega_{3}^{2}-J_{23}^{2}}+1\right]^{-1}-1$
(11)

we end up with a transformed Hamiltonian

$\displaystyle
H_{D}=\frac{1}{2}\mathbb{P}^{t}\mathbb{P}+\frac{1}{2}{\mathbb{X}}^{t}\mathbb{V}_{D}\mathbb{X}$
(12)

including the diagonalized potential matrix

$\displaystyle\mathbb{V}_{D}=\text{diag}\left(\Omega_{1}^{2},\Omega_{2}^{2},\Omega_{3}^{2}\right)$
(13)

and the normal frequencies

$\displaystyle\Omega_{1}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\left[\varpi+2\sqrt{p}\cos(\Phi)\right]$ (14)
$\displaystyle\Omega_{2}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\left[\varpi+2\sqrt{p}\cos\left(\Phi+\frac{2\pi}{3}\right)\right]$
(15) $\displaystyle\Omega_{3}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{3}\left[\varpi+2\sqrt{p}\cos\left(\Phi-\frac{2\pi}{3}\right)\right]$
(16)

where the parameters $\varpi$, $p$ and $\Phi$ read as

$\displaystyle
c_{0}=\sum\limits_{(i,i)\neq(j,k),j<k}^{3}\omega^{2}_{i}J_{jk}^{2}-\prod\limits_{i=1}^{3}\omega^{2}_{i}+2\prod\limits_{i<j}^{3}J_{ij},\quad\varpi=\sum\limits_{i=1}^{3}\omega^{2}_{i},\quad
c_{1}=\sum\limits_{i<j}^{3}\omega_{i}^{2}\omega_{j}^{2}-\sum\limits_{i<j}^{3}J_{ij}^{2}$
(17)
$\displaystyle\Phi=\frac{1}{3}\arctan\left(\frac{\sqrt{p^{3}-q^{2}}}{q}\right),\quad
p=\varpi^{2}-3c_{1},\quad
q=-\frac{27}{2}c_{0}+\varpi^{3}-\frac{9}{2}c_{1}\varpi.$ (18)

As a result, we can write (12) in the following diagonal form

$\displaystyle H_{D}$ $\displaystyle=$
$\displaystyle\frac{p_{1}^{2}}{2}+\frac{p_{2}^{2}}{2}+\frac{p_{3}^{2}}{2}+\frac{1}{2}\Omega_{1}^{2}x_{1}^{2}+\frac{1}{2}\Omega_{2}^{2}x_{2}^{2}+\frac{1}{2}\Omega_{3}^{2}x_{3}^{2}.$
(19)

To discuss entanglement together with virtual excitations, we employ the
annihilation $\mathfrak{a}_{j}$ and creation $\mathfrak{a}^{+}_{j}$ operators
via the mapping

$\displaystyle\mathfrak{a}_{j}=\left(\mathfrak{a}_{j}^{+}\right)^{+}=\frac{1}{\sqrt{2\omega_{j}}}\left(\omega_{j}x_{j}+ip_{j}\right),\quad
j=1,2,3.$ (20)

These enable us to cast the non-diagonalized Hamiltonian (1) as

$\displaystyle
H=\sum\limits_{j=1}^{3}\omega_{j}\left(\mathfrak{a}^{+}_{j}\mathfrak{a}_{j}+\frac{1}{2}\right)+\sum\limits_{j<k}^{3}g_{jk}(\mathfrak{a}_{j}^{+}+\mathfrak{a}_{j})(\mathfrak{a}^{+}_{k}+\mathfrak{a}_{k})$
(21)

where the coupling strengths have been set

$\displaystyle g_{jk}:=\frac{J_{jk}}{2\sqrt{\omega_{j}\omega_{k}}},\qquad
1\leq j<k\leq 3.$ (22)

Also, the previous Euler transformation (8) leads to a diagonalized
Hamiltonian

$\displaystyle H_{D}$ $\displaystyle=$
$\displaystyle\sum\limits_{j=1}^{3}\Omega_{j}\left(\mathfrak{A}^{+}_{j}\mathfrak{A}_{j}+\frac{1}{2}\right)$
(23)

in terms of the operators

$\displaystyle\mathfrak{A}_{j}=\left(\mathfrak{A}_{j}^{+}\right)^{+}=\frac{1}{\sqrt{2\Omega_{j}}}\left(\Omega_{j}x_{j}+ip_{j}\right)$
(24)

which satisfy the commutation relation
$\left[\mathfrak{A}_{k},\mathfrak{A}^{+}_{l}\right]=\delta_{kl}$, and the rest
are null. As a result, the Hamiltonian in terms of the diagonal operator reads
as

$\displaystyle
H=\mathbb{T}^{-1}(\alpha,\beta,\gamma)H_{D}\mathbb{T}(\alpha,\beta,\gamma)$
(25)

or equivalently

$\displaystyle H=e^{i\gamma\mathcal{J}_{3}}\ e^{i\beta\mathcal{J}_{2}}\
e^{i\alpha\mathcal{J}_{3}}\ H_{D}\ e^{-i\alpha\mathcal{J}_{3}}\
e^{-i\beta\mathcal{J}_{2}}\ e^{-i\gamma\mathcal{J}_{3}}.$ (26)

Now, it is worthwhile to note that $H_{D}$ is diagonal in the representation
$\left\\{\mathfrak{A}_{j},\mathfrak{A}_{j}^{+}\right\\}$. To diagonalize the
Hamiltonian in terms of original operators $\mathfrak{a}_{j}$ and
$\mathfrak{a}_{j}^{+}$, we rewrite (23) as

$\displaystyle
H_{D}=\sum\limits_{j=1}^{3}\tilde{\Omega}_{j}\left(\mathfrak{a}_{j}^{+}\mathfrak{a}_{j}+\frac{1}{2}\right)+\mathfrak{g}_{j}\left(\left(\mathfrak{a}_{j}^{+}\right)^{2}+\mathfrak{a}_{j}^{2}\right)$
(27)

where we have defined

$\displaystyle\tilde{\Omega}_{j}:=\frac{1}{2}\left(\frac{\Omega_{j}^{2}}{\omega_{j}}+\omega_{j}\right),\quad\mathfrak{g}_{j}:=\frac{1}{4}\left(\frac{\Omega_{j}^{2}}{\omega_{j}}-\omega_{j}\right)$
(28)

and perform the following squeezing operators

$\displaystyle
S_{j}(r_{j}):=e^{\frac{r_{j}}{2}\left[\mathfrak{a}_{j}^{2}-(\mathfrak{a}_{j}^{+})^{2}\right]},\quad
r_{j}:=\frac{1}{2}\ln\left(\frac{\Omega_{j}}{\omega_{j}}\right)$ (29)

giving rise to the following diagonal Hamiltonian in the original basis

$\displaystyle H_{d}=\prod\limits_{j=1}^{3}S_{j}(-r_{j})\ H_{D}\
\prod\limits_{j=1}^{3}S_{j}(r_{j})=\sum\limits_{j=1}^{3}\Omega_{j}\left(\mathfrak{a}_{j}^{+}\mathfrak{a}_{j}+\frac{1}{2}\right).$
(30)

Combining all to map the Hamiltonian (1) as

$\displaystyle H$ $\displaystyle=$ $\displaystyle e^{i\gamma\mathcal{J}_{3}}\
e^{i\beta\mathcal{J}_{2}}\ e^{i\alpha\mathcal{J}_{3}}\ S_{1}(-r_{1})\
S_{2}(-r_{2})\ S_{3}(-r_{3})\ H_{d}\ S_{1}(r_{1})S_{2}(r_{2})\ S_{3}(r_{3})\
e^{-i\alpha\mathcal{J}_{3}}\
e^{-i\beta\mathcal{J}_{2}}e^{-i\gamma\mathcal{J}_{3}}.$ (31)

More analysis will be conducted to glean further information from these
results and in accordance with the core components of the current system.

##  III Milburn dynamics and covariance matrix

To study the intrinsic decoherence effect on our quantum system, we follow the
Milburn Model (MM) [35]. In MM, the evolution of the density $\rho(t)$ on a
sufficiently small time scale is uncertain. As a result, the system evolves
from the state $\rho(t)$ to $\rho(t+\tau)$ via

$\displaystyle\rho(t+\tau)=e^{-i\mathfrak{P}(\tau)H}\rho(t)e^{i\mathfrak{P}(\tau)H}$
(32)

with the probability $p(\tau)$. It is worth emphasizing that if $t=0$ and
$\rho(0)$ is Gaussian, the state $\rho(t)$ is Gaussian for all $\tau$. This is
due to the stochastic unitarity of Milburn evolution and the quadratic form of
the Hamiltonian [35, 36]. In addition, we have
$\mathfrak{P}(\tau)\rightarrow\tau$ and $p(\tau)\rightarrow 1$ in standard
quantum mechanics. Assuming a Poisson model for the stochastic dynamics of
time steps, the state will be governed by the following Milburn master
equation

$\dot{\rho}(t)=\Gamma\left[e^{-\frac{iH}{\Gamma}}\rho(0)e^{\frac{iH}{\Gamma}}-\rho(t)\right]$
(33)

where $\Gamma^{-1}$ is the decoherence rate defined as

$\displaystyle\Gamma^{-1}=\lim\limits_{\tau\rightarrow 0}\mathfrak{P(\tau)}.$
(34)

The formal solution to Milburn’s equation is given by [36, 37]

$\displaystyle\rho(t,\Gamma)$ $\displaystyle=$ $\displaystyle e^{-\Gamma
t}\sum_{k=0}^{\infty}\frac{(\Gamma
t)^{k}}{k!}e^{-\frac{ikH}{\Gamma}}\rho(0)e^{\frac{iHk}{\Gamma}}.$ (35)

Additionally, by using (31), one can show the result

$\displaystyle e^{-i\frac{k}{\Gamma}H}=e^{i\gamma\mathcal{J}_{3}}\
e^{i\beta\mathcal{J}_{2}}\ e^{i\alpha\mathcal{J}_{3}}\ S_{1}(-r_{1})\
S_{2}(-r_{2})\ S_{3}(-r_{3})\ e^{-i\frac{k}{\Gamma}H_{d}}\ S_{1}(r_{1})\
S_{2}(r_{2})\ S_{3}(r_{3})\ e^{-i\alpha\mathcal{J}_{3}}\
e^{-i\beta\mathcal{J}_{2}}\ e^{-i\gamma\mathcal{J}_{3}}.$ (36)

Since the coupling becomes ultra-strong, the ground state becomes populated
with virtual excitations. This is due to the counter-rotating terms [11, 12],
which simultaneously create the excitations in two modes. Then we will
consider the ground state

$\displaystyle\rho(0)=|000\rangle\langle 000|.$ (37)

To go further, we collect the annihilation and creation operators of the three
modes in the vector

$\displaystyle\mathbb{A}=(\mathbb{A}_{n})=(\mathfrak{a}_{1},\mathfrak{a}_{1}^{+},\mathfrak{a}_{2},\mathfrak{a}_{2}^{+},\mathfrak{a}_{3},\mathfrak{a}_{3}^{+})^{t}.$
(38)

The commutation relations
$[\mathfrak{a}_{i},\mathfrak{a}^{+}_{j}]=\delta_{ij}$, reduce to the following
compact form

$\displaystyle[\mathbb{A}_{n},\mathbb{A}_{m}]=i\mathbb{J}_{n,m},\qquad
i\mathbb{J}_{6\times 6}=\bigoplus\tilde{\mathbb{I}}_{2\times
2},\qquad\tilde{\mathbb{{I}}}_{2\times 2}=\left(\begin{array}[]{cc}0&1\\\
-1&0\end{array}\right).$ (41)

Now, to put the transformations used above in their symplectic form, we recall
that

$\displaystyle\mathcal{U}\mathbb{A}\mathcal{U}^{-1}=\mathbb{S}\mathbb{A}$ (42)

where the matrix $\mathbb{S}$ is called the symplectic representation of
$\mathcal{U}$ and verifies $\det(\mathbb{S})=1$ with
$\mathbb{S}^{+}\mathbb{J}\mathbb{S}=\mathbb{J}$. In our formalism, the
covariance matrix reduces to

$\displaystyle\sigma_{n,m}$ $\displaystyle=$
$\displaystyle\langle\\{\mathbb{A}_{n},\mathbb{A}^{+}_{m}\\}\rangle_{\rho}-2\langle\mathbb{A}_{n}\rangle_{\rho}\langle\mathbb{A}^{+}_{m}\rangle_{\rho}.$
(43)

We mention that our state is centered, thus,
$\langle\mathbb{A}_{n}\rangle_{\rho}=\langle\mathbb{A}^{+}_{m}\rangle_{\rho}=0$.
On the other hand, we rewrite (35) in its symplectic form

$\displaystyle\sigma(t,\Gamma)=e^{-\Gamma
t}\sum\limits_{k=0}^{\infty}\frac{(\Gamma
t)^{k}}{k!}\mathbb{H}\sigma(0)\mathbb{H}^{+}$ (44)

by showing that

$\displaystyle\sigma(0)=\text{diag}(1,1,1,1,1,1)$ (45)

where $(\cdot)^{+}$ is the Hermitian conjugate of $(\cdot)$ and $\mathbb{H}$
is the symplectic representation of the operator $e^{-\frac{ikH}{\Gamma}}$. We
use (36) to demonstrate that the matrix $\mathbb{H}$ is

$\displaystyle\mathbb{H}=\mathbb{S}_{12}(\gamma)\mathbb{S}_{13}(\beta)\mathbb{S}_{12}(\alpha)\mathbb{S}_{123}(-r_{1},-r_{2},-r_{3})\mathbb{D}(\Gamma)\mathbb{S}_{123}(r_{1},r_{2},r_{3})\mathbb{S}_{12}(-\alpha)\mathbb{S}_{13}(-\beta)\mathbb{S}_{12}(-\gamma)$
(46)

where the symplectic matrices of $e^{iz\mathcal{J}_{3}}$ and
$e^{i\beta\mathcal{J}_{2}}$ are given by

$\displaystyle\mathbb{S}_{12}(z)$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\cos(z)&0&\frac{R+R^{-1}}{2}\sin(z)&\frac{R-R^{-1}}{2}\sin(z)&0&0\\\
0&\cos(z)&\frac{R-R^{-1}}{2}\sin(z)&\frac{R+R^{-1}}{2}\sin(z)&0&0\\\
-\frac{R+R^{-1}}{2}\sin(z)&\frac{R-R^{-1}}{2}\sin(z)&\cos(z)&0&0&0\\\
\frac{R-R^{-1}}{2}\sin(z)&-\frac{R+R^{-1}}{2}\sin(z)&0&\cos(z)&0&0\\\
0&0&0&0&1&0\\\ 0&0&0&0&0&1\end{pmatrix}$ (47)
$\displaystyle\mathbb{S}_{13}(\beta)$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\cos(\beta)&0&0&0&-\frac{\tilde{R}+\tilde{R}^{-1}}{2}\sin(\beta)&\frac{\tilde{R}-\tilde{R}^{-1}}{2}\sin(\beta)\\\
0&\cos(\beta)&0&0&\frac{\tilde{R}-\tilde{R}^{-1}}{2}\sin(\beta)&-\frac{\tilde{R}+\tilde{R}^{-1}}{2}\sin(\beta)\\\
0&0&1&0&0&0\\\ 0&0&0&1&0&0\\\
\frac{\tilde{R}+\tilde{R}^{-1}}{2}\sin(\beta)&\frac{\tilde{R}-\tilde{R}^{-1}}{2}\sin(\beta)&0&0&\cos(\beta)&0\\\
\frac{\tilde{R}-\tilde{R}^{-1}}{2}\sin(\beta)&\frac{\tilde{R}+\tilde{R}^{-1}}{2}\sin(\beta)&0&0&1&\cos(\beta)\\\
\end{pmatrix}$ (48)

with $z=\alpha,\gamma$, $R=\sqrt{\frac{\omega_{1}}{\omega_{2}}}$ and
$\tilde{R}=\sqrt{\frac{\omega_{3}}{\omega_{1}}}$ have been set. The symplectic
representation of the three squeezors is reduced to

$\displaystyle\mathbb{S}_{123}(r_{1},r_{2},r_{3})=\bigoplus\limits_{j=1}^{3}\begin{pmatrix}\cosh(r_{j})&\sinh(r_{j})\\\
\sinh(r_{j})&\cosh(r_{j})\\\ \end{pmatrix},\quad j=1,2,3.$ (49)

Finally, the symplectic form of the transformation
$e^{-i\frac{k}{\Gamma}H_{d}}$, which is a diagonal matrix $\mathbb{D}$ [38],
is straightforwardly demonstrated. This is

$\displaystyle\mathbb{D}(\Gamma)=\bigoplus\limits_{j=1}^{3}\begin{pmatrix}e^{-ik\frac{\Omega_{j}}{\Gamma}}&0\\\
0&e^{ik\frac{\Omega_{j}}{\Gamma}}\\\ \end{pmatrix}.$ (50)

We mention that after performing the matrix product (44), the infinite sum can
be exactly computed by using the identity

$\displaystyle e^{-\Gamma t}\sum\limits_{k=0}^{+\infty}\frac{(\Gamma
t)^{k}}{k!}e^{if(\Omega_{1},\Omega_{2},\Omega_{3})\frac{k}{\Gamma}}=e^{\Gamma
t\left(e^{i\frac{f(\Omega_{1},\Omega_{2},\Omega_{3})}{\Gamma}}-1\right)}.$
(51)

where $f$ is a well defined function of normal frequencies $\Omega_{j}$.

##  IV Quantum entanglement and virtual excitations

We analyze the dynamics and redistribution of virtual excitations between the
three oscillators, together with the dynamics of bipartite entanglement. Both
quantities are derived from the covariance matrix (44), that is

$\displaystyle\sigma(t,\Gamma)=\begin{pmatrix}\sigma_{a}&\sigma_{ab}&\sigma_{ac}\\\
\sigma_{ab}^{t}&\sigma_{b}&\sigma_{bc}\\\
\sigma_{ac}^{t}&\sigma_{bc}^{t}&\sigma_{c}\end{pmatrix}$ (52)

where the single mode states are given by

$\displaystyle\sigma_{a}=\left(\begin{array}[]{cc}\sigma_{11}&\sigma_{12}\\\
\sigma_{12}^{\ast}&\sigma_{11}\end{array}\right),\quad\sigma_{b}=\left(\begin{array}[]{cc}\sigma_{33}&\sigma_{34}\\\
\sigma_{34}^{\ast}&\sigma_{33}\end{array}\right),\quad\sigma_{c}=\left(\begin{array}[]{cc}\sigma_{55}&\sigma_{56}\\\
\sigma_{56}^{\ast}&\sigma_{55}\end{array}\right)$ (59)

and the correlation matrices between two modes are

$\displaystyle\sigma_{ab}=\left(\begin{array}[]{cc}\sigma_{13}&\sigma_{14}\\\
\sigma_{23}^{\ast}&\sigma_{24}\end{array}\right),\ \
\sigma_{ac}=\left(\begin{array}[]{cc}\sigma_{15}&\sigma_{16}\\\
\sigma_{25}^{\ast}&\sigma_{26}\end{array}\right),\ \
\sigma_{bc}=\left(\begin{array}[]{cc}\sigma_{35}&\sigma_{36}\\\
\sigma_{45}^{\ast}&\sigma_{46}\end{array}\right)$ (66)

where $t$ stands for transpose. By making use of the tracing-out prescription,
we end up with the bipartite states [39]

$\displaystyle\sigma^{ab}=\begin{pmatrix}\sigma_{a}&\sigma_{ab}\\\
\sigma_{ab}^{t}&\sigma_{b}\end{pmatrix},\quad\sigma^{ac}=\begin{pmatrix}\sigma_{a}&\sigma_{ac}\\\
\sigma_{ac}^{t}&\sigma_{c}\end{pmatrix},\quad\sigma^{bc}=\begin{pmatrix}\sigma_{b}&\sigma_{bc}\\\
\sigma_{bc}^{t}&\sigma_{c}\end{pmatrix}.$ (67)

Due to the length of the explicit expressions for the matrix elements
$\sigma_{nm}$ ($n,m=1,\cdots 6$), we only report their formal relations and
numerical results.

We investigate the impact of virtual excitations on quantum entanglement in a
harmonic oscillator system with ultra-strongly coupled. We show that the
average number of excitations in the ground state is

$\displaystyle\langle N_{j}\rangle$ $\displaystyle=$
$\displaystyle\langle\mathfrak{a}_{j}^{+}\mathfrak{a}_{j}\rangle=\frac{1}{2}(\sigma_{2j-1,2j-1}-1),\
\ j=1,2,3.$ (68)

Additionally, to quantify the bipartite entanglement $(k|l)$ while
$k,l=a,b,c$. The bipartite covariance matrix $\sigma^{kl}$ (67) are partially
transposed to obtain $\tilde{\sigma}$. The positive minimal symplectic
eigenvalue is given by [40]

$\displaystyle\tilde{\nu}_{min}^{kl}=\sqrt{\frac{1}{2}\left(\Delta_{kl}-\sqrt{\Delta_{kl}^{2}-4\det\sigma^{kl}}\right)},\quad
k,l=a,b,c$ (69)

where the seralian $\Delta_{kl}$ is defined as

$\displaystyle\Delta_{kl}=\det\sigma_{k}+\det\sigma_{l}-2\det\sigma_{kl}.$
(70)

Consequently, the logarithmic negativity reduces to

$\displaystyle E_{kl}=\max\left(0,-\ln(\tilde{\nu}^{kl}_{min})\right).$ (71)

In contrast to non-classical correlations that go beyond inseparability,
entanglement is monogamous among the parties of the system. This entails that
entanglement is not freely shared between the parties of the whole system. In
particular, for a tripartite system, the monogamy of entanglement was
formulated in [41] and shown in [42], that is

$\displaystyle E_{k|lm}\geq E_{kl}+E_{km},\quad k,l,m\in\\{a,b,c\\}$ (72)

where $E_{kl}$ represents the bipartite entanglement of $k$ and $l$, the
vertical bar denotes the bipartite split. Motivated by the results of [12,
11], which demonstrate the interaction between entanglement and virtual
excitations, we ask a fundamental question: whether virtual excitations are
distributed in the same way that entanglement is.

##  V Numerical results and discussions

###  V.1 Isotropic case

Assume you have a resonant system with
$\omega_{1}=\omega_{2}=\omega_{3}=\omega_{r}$ and $J_{12}=J_{13}=J_{23}=J$. As
a result, the standard frequencies become

$\displaystyle\Omega_{1}=\sqrt{\omega_{r}^{2}+2J},\quad\Omega_{2}=\Omega_{3}=\sqrt{\omega_{r}^{2}-J}.$
(73)

The matrix form of the rotation operator that leads to the diagonalized
Hamiltonian is

$\displaystyle\mathcal{R}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{ccc}\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}\\\
0&\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{3}}\\\
-\sqrt{\frac{2}{3}}&\frac{1}{\sqrt{6}}&\frac{1}{\sqrt{6}}\end{array}\right)$
(77)

and the rotation angles are reduced to

$\displaystyle\gamma$ $\displaystyle\to$
$\displaystyle-\arctan\left(\frac{1}{3}\right)$ (78) $\displaystyle\alpha$
$\displaystyle\to$ $\displaystyle\arctan\left(\sqrt{\frac{3}{2}}\right)$ (79)
$\displaystyle\beta$ $\displaystyle\to$
$\displaystyle\arccos\left(\frac{1}{\sqrt{6}}\right).$ (80)

Then, we show the following virtual excitations

$\displaystyle\langle N_{1}\rangle(t)$ $\displaystyle=$
$\displaystyle\frac{24}{50}{\cosh^{2}(s_{1})}{\sinh^{2}(s_{1})}\left[1-\exp[-\Gamma
t(1-\cos(2\Omega_{1}/\Gamma))]\cos(\Gamma t\sin(2\Omega_{1}/\Gamma))\right]$
$\displaystyle+$
$\displaystyle\frac{76}{50}{\cosh^{2}(s_{2})}{\sinh^{2}(s_{2})}\left[1-\exp[-\Gamma
t(1-\cos(2\Omega_{2}/\Gamma))]\cos(\Gamma t\sin(2\Omega_{2}/\Gamma))\right]$
$\displaystyle\langle N_{2}\rangle(t)$ $\displaystyle=$
$\displaystyle\frac{1}{50}{\cosh^{2}(s_{1})}{\sinh^{2}(s_{1})}\left[1-\exp[-\Gamma
t(1-\cos(2\Omega_{1}/\Gamma))]\cos(\Gamma t\sin(2\Omega_{1}/\Gamma))\right]$
$\displaystyle+$
$\displaystyle\frac{99}{50}{\cosh^{2}(s_{2})}{\sinh^{2}(s_{2})}\left[1-\exp[-\Gamma
t(1-\cos(2\Omega_{2}/\Gamma))]\cos(\Gamma t\sin(2\Omega_{2}/\Gamma))\right]$
$\displaystyle\langle N_{3}\rangle(t)$ $\displaystyle=$
$\displaystyle\frac{3}{2}{\cosh^{2}(s_{1})}{\sinh^{2}(s_{1})}\left[1-\exp[-\Gamma
t(1-\cos(2\Omega_{1}/\Gamma))]\cos(\Gamma t\sin(2\Omega_{1}/\Gamma))\right]$
$\displaystyle+$
$\displaystyle\frac{1}{2}{\cosh^{2}(s_{2})}{\sinh^{2}(s_{2})}\left[1-\exp[-\Gamma
t(1-\cos(2\Omega_{2}/\Gamma))]\cos(\Gamma t\sin(2\Omega_{2}/\Gamma))\right].$

Additionally, by virtue of the above expressions, we show the steady values of
excitations, such that

$\displaystyle\langle N_{1}\rangle_{\infty}$ $\displaystyle=$
$\displaystyle\frac{24}{50}{\cosh^{2}(s_{1})}{\sinh^{2}(s_{1})}+\frac{76}{50}{\cosh^{2}(s_{2})}{\sinh^{2}(s_{2})}=\frac{J^{2}(43\omega_{r}^{2}+14J)}{200\omega_{r}^{2}(\omega_{r}^{2}+2J)(\omega_{r}^{2}-J)}$
(84) $\displaystyle\langle N_{2}\rangle_{\infty}$ $\displaystyle=$
$\displaystyle\frac{1}{50}{\cosh^{2}(s_{1})}{\sinh^{2}(s_{1})}+\frac{99}{50}{\cosh^{2}(s_{2})}{\sinh^{2}(s_{2})}=\frac{J^{2}(103\omega_{r}^{2}+149J)}{800\omega_{r}^{2}(\omega_{r}^{2}+2J)(\omega_{r}^{2}-J)}$
$\displaystyle\langle N_{3}\rangle_{\infty}$ $\displaystyle=$
$\displaystyle\frac{3}{2}{\cosh^{2}(s_{1})}{\sinh^{2}(s_{1})}+\frac{1}{2}{\cosh^{2}(s_{2})}{\sinh^{2}(s_{2})}=\frac{J^{2}(13\omega_{r}^{2}-10J)}{32\omega_{r}^{2}(\omega_{r}^{2}+2J)(\omega_{r}^{2}-J)}$

which are independent of $\Gamma$. Note also that the three steady values will
be equal when the coupling takes $J=\omega_{r}^{2}/2$. Furthermore, it is
simple to demonstrate that $\langle N_{j}\rangle$ are increasing functions
with respect to $J$. Furthermore, we show $\langle
N_{3}\rangle_{\infty}>\langle N_{1}\rangle_{\infty}>\langle
N_{2}\rangle_{\infty}$ for all $J<\omega_{r}^{2}/2$, and the hierarchy will be
inverted, i.e., $\langle N_{3}\rangle_{\infty}<\langle
N_{1}\rangle_{\infty}<\langle N_{2}\rangle_{\infty}$ for $J>\omega_{r}^{2}/2$.
We also mention that the dynamics reaches a steady state during times of order

$\displaystyle
t_{steady}\sim\max\left[\frac{1}{\Gamma(1-\cos(2\Omega_{1}/\Gamma))},\frac{1}{\Gamma(1-\cos(2\Omega_{2}/\Gamma))}\right].$
(85)

As shown in Fig. 2, as the coupling increases, $t_{steady}$ increases.
Additionally, we observe that excitations increase as the coupling becomes
ultra-strong. Furthermore, as previously discussed, the steady values of
excitations will be equal as long as $J=\omega_{r}^{2}/2$ and will become
hierarchical if $J\neq\omega_{r}^{2}/2$.

Figure 2: (color online) The effects of coupling $J$ on the time of the
transient state. The parameters $\omega_{r}=1$ and $\Gamma=100$ were set.

###  V.2 Effects of Milburn decoherence

To show the effects of the intrinsic decoherence on the dynamics of
excitations, we plot in Fig. 3 the dynamics of the average of excitations
$\langle N_{k}\rangle=\langle\mathfrak{a}_{k}^{+}\mathfrak{a}_{k}\rangle$. We
set $\omega_{r}=\omega_{1}=\omega_{2}=\omega_{3}=1$, $J_{13}=J_{12}=0.1$ and
$J_{23}=0$ (indicating an open chain). The three populations simultaneously
appear and undergo damped oscillatory behaviour to a steady state of
excitation. This multi-oscillatory regime is due to normal frequencies
$\Omega_{k}(k=1,2,3)$ and $\Gamma$. Besides, the three plots show that the
excitations in particles $a$ and $b$ are the same which implies that the
excitations exchange between them. We see that the excitations satisfy the
inequalities $\langle N_{3}\rangle\leq\langle N_{1,2}\rangle$ all the time.
Furthermore, we notice that as long as $\Gamma$ is large, the stronger the
excitations are. The Schrödinger dynamics is obtained for
$\Gamma\rightarrow\infty$, where the excitations undergo an undamped
oscillatory behavior, and the steady excitations disappear. We also notice
that as $\Gamma$ decreases, the oscillations disappear and the excitations
exponentially revive.

Figure 3: (color online) The effect of decoherence rate $\Gamma$ on the
dynamics of virtual excitations. The parameters
$\omega_{1}=\omega_{2}=\omega_{3}=\omega_{r}=1$, $J_{23}=J_{12}=0.1$, and
$J_{13}=0$ were set. The case where $\Gamma\to+\infty$, corresponds to
Schrödinger dynamics.

###  V.3 Effects of coupling $J_{13}$ and central frequency $\omega_{2}$

To investigate the evolution of excitations under the effects of coupling
$J_{13}$ and the central frequency $\omega_{2}$, we plot the dynamics of
excitations in three particles by varying $J_{13}$ and $\omega_{2}$,
respectively, in Figs. 4 and 5. First, we set
$\omega_{r}=\omega_{1}=\omega_{2}=\omega_{3}=1$, $J_{13}=J_{12}=0.1$ as well
as varying the coupling between particles $a$ and $c$. As expected, because of
the coupling $J_{13}$, the excitations in particle $c$ become more intensive
than those in particles $a$ and $b$. Furthermore, the higher the coupling
$J_{1}3$, the more excitations are significant in particles $a$ and $c$,
implying that excitations are transferred between them. Now, we weakly couple
the particles $(a|b)$ and $(b|c)$, and ultra-strongly couple the particles
$(a|c)$. The frequencies of particles a and c are assumed to be equal to
$\omega_{r}=1$, and we vary the central frequency $\omega_{2}$. As we move
away from resonance, the excitations in particle c become more intense, owing
to $J_{13}$. Additionally, we notice an inverted dynamics for $\langle
N_{1}\rangle$ and $\langle N_{2}\rangle$, which means the transfer of
excitations between them. Decreasing the frequency $\omega_{2}$ to $0.8$, we
observe the decrease of excitations in all particles. Fulfilling the condition
$\omega_{1}+\omega_{3}=4\omega_{2}$, we observe a resonance, which amplifies
the excitations in particles $a$ and $b$, and diminishes them in particle $c$.
It is also worth noting that by lowering the central frequency to
$\omega_{2}=0.3$, we observe an amazing extinction of excitations in particle
$b$ as well as a transfer of excitations between particles $a$ and $c$. These
results show that, when excitations are important in a particle, they will be
less important in at least one of the remaining particles. This indicates that
excitations are not freely generated or annihilated. For this, we elucidate
the polygamy of excitations among the three parties of the whole system.

Figure 4: (color online) The effect of central coupling $J_{13}$ on dynamics
of virtual excitations. The parameters
$\omega_{1}=\omega_{2}=\omega_{3}=\omega_{r}=1$,$J_{23}=J_{12}=0.1$ and
$\Gamma=100$ were set.

Figure 5: (color online) The effect of central frequency $\omega_{2}$ on the
dynamics of virtual excitations. The parameters
$\omega_{1}=\omega_{3}=\omega_{r}=1$,$J_{23}=J_{12}=0.01$, $J_{13}=0.9$ and
$\Gamma=50$ were set.

###  V.4 Virtual excitations and polygamy

To investigate the interplay between excitations in the bi-partitions of three
oscillators. We define the bipartite virtual excitations in modes $k$ and $l$
as follows

$\displaystyle N_{k|l}$ $\displaystyle=$ $\displaystyle\sqrt{\langle
N_{k}\rangle\langle N_{l}\rangle}.$ (86)

This can be seen as a geometric mean of the number of excitations. As a
result, the excitations in the bi-partition vanish if and only if they vanish
at least in one mode. The present choice will be justified later, when the
interconnection between excitations and entanglement is addressed. Similarly,
we define the excitation in the partition $(k|lm)$ by

$\displaystyle N_{k|lm}$ $\displaystyle=$ $\displaystyle\sqrt{\langle
N_{k}\rangle N_{l|m}}.$ (87)

The excitations will be polygamous [43] if they satisfy the following
triangular inequality

$\displaystyle N_{k|lm}\leq N_{k|l}+N_{k|m}.$ (88)

These inequalities constrain the generation and extinction of excitations.
This finding can be used to consider virtual excitations as a kind of quantum
correlation [44, 11, 12]. Furthermore, the virtual excitations can be used to
quantify quantum resources beyond entanglement. Here, to show the interplay
between excitations, we plot in Fig. 6, the dynamics of the trade-off
quantities

$\displaystyle\delta^{i}_{jk}$ $\displaystyle=$ $\displaystyle
N_{i|j}+N_{i|k}-N_{i|jk},\quad k,l,m=a,b,c.$ (89)

We numerically show that excitations are polygamous. This entails that
excitations are not freely generated or annihilated between the several
parties of the whole system. However, they undergo the triangular inequality
given in (88).

Figure 6: (color online) The illustration of polygamy of excitations. The
parameters $\omega_{1}=\omega_{2}=\omega_{3}=\omega_{r}=1$, and $\Gamma=50$
were set.

###  V.5 Entanglement and virtual excitations

We tackle in this paragraph the interconnection between excitations and
entanglement. For this aim, we plot in Fig. 7 the dynamics of entanglement and
excitations. First, we assume an open chain $J_{13}=0$ and ultra-strongly
coupled oscillators $J_{12}=J_{23}=0.1$. Surprisingly, the bi-partition $a|b$
is more entangled than $b|c$. This originates from the number of excitations
in those bi-partitions, i.e., $N_{a|b}\leq N_{b|c}$. Another result that
deserves attention is the entanglement in the bi-partition $(a|c)$. The
oscillators are indeed decoupled at $J_{13}=0$, but they are indirectly
coupled via the oscillator $b$. As a result, the numbers of excitations in the
partitions $a|c$ and $b|c$ are equal, and $E_{ac}=E_{bc}$. Second, we set the
oscillators $a$ and $b$ to $J_{13}=J_{23}=0$ and strongly coupled them. Thus,
the oscillator $c$ does not contain excitations. As a result, the partitions
$j|c$, while $(j\neq c)$ are not entangled and do not contain excitations.
Finally, it is worthwhile to mention that both profiles exhibit similar
dynamics. This points out that excitations can be used as quantifiers of
entanglement [11].

Figure 7: (color online) The interplay between the dynamics of entanglement
and virtual excitations. The parameters
$\omega_{1}=\omega_{2}=\omega_{3}=\omega_{r}=1$, and $\Gamma=50$ were set.

##  VI Concluding remarks

In this work, we analyzed the interconnection between excitations and
entanglement. The system of interest is made up of three non-resonant
oscillators. First of all, by making use of Euler rotation and squeezor
transformations, we have obtained the diagonalized form of the Hamiltonian in
terms of creation and annihilation operators. We have shown that the Milburn
density $\rho$ is Gaussian, and this is based on the unitary Milburn
evolution, the quadratic form of the Hamiltonian and the Gaussian nature of
the initial ground state. Additionally, due to the Gaussian nature of the
quantum density, we have derived the Milburn density in its symplectic
covariance matrix form.

The corresponding covariance matrix is expressed as a poissonian sum of the
product of $26$ symplectic matrices. This finding led us to quantify the
entanglement together with virtual excitations only based on the covariance
matrix. Furthermore, we investigated the resonant case and discovered that
steady values of excitations are dependent on the coupling $J$ and the
resonant frequency $\omega_{r}$. We have also shown that the necessary time to
establish the steady states can be controlled by the normal frequencies
$\Omega_{j}$ and $\Gamma$. In addition, the excitations were analyzed between
all oscillators, and it was found that they exchanged between them. The
redistribution of excitations among partitions was investigated, and it was
discovered that virtual excitations, like quantum discord, violate the
principle of monogamy. And we have shown that excitations are freely generated
or annihilated but constrained by a triangular inequality called the polygamy
inequality.

As a consequence, excitations are polygamous and are not freely distributed
between the parts of the system. This issue can be developed later, in both
continuous and discrete variables, to unify and understand the hierarchy of
quantum correlations via virtual excitations-based measures. Finally, we have
analyzed the interplay between excitations and entanglement. It is found that
excitations are necessary to maintain the pairwise entanglement between
oscillators. As a result, the hierarchy relationship between excitations and
entanglement is demonstrated. Another point worth noting is that as one
approaches the hermitianity point, i.e., $\min(\Omega_{j})\to 0$, the
excitations and entanglement become less susceptible to decoherence.

The present work will not remain at this stage, and therefore, we are willing
to extend it to deal with other issues of the entanglement and related
matters. On the other hand, we hope that our work gives a pulse to a new route
to quantifying and processing quantum information.

## References

  * [1] A. F. Kockum, A. Miranowicz, S. D. Liberato, S. Savasta, and F. Nori, Nat. Rev. Phys. 1, 19 (2019).
  * [2] A. Ridolfo, M. Leib, S. Savasta, and M. J. Hartmann, Phys. Rev. Lett. 109, 193602 (2012).
  * [3] L. Garziano, V. Macrì, R. Stassi, O. Di Stefano, F. Nori, and S. Savasta, Phys. Rev. Lett. 117, 043601 (2016).
  * [4] T. Niemczyk, Nat. Phys. 6, 772 (2010).
  * [5] O. Di Stefano, Phys. Rev. Lett. 122, 030402 (2019).
  * [6] P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, Rev. Mod. Phys. 91, 025005 (2019).
  * [7] S. Terradas-Briansó, C. A. González-Gutiérrez, F. Nori, L. Martín-Moreno, and D. Zueco, arXiv:2205.07915 [quant-ph] (2022).
  * [8] A. Le Boité, M. J. Hwang, H. Nha, and M. B. Plenio, Phys. Rev. A 94, 033827 (2016).
  * [9] S. De Liberato, Phys. Rev. Lett. 112, 016401 (2014).
  * [10] I. Lizuain, J. Casanova, J. J. García-Ripoll, J. G. Muga, and E. Solano, Phys. Rev. A 81, 062131 (2010).
  * [11] R. Hab-arrih, A. Jellal, D. Stefanatos, and A. Merdaci, Quantum Rep. 3, 684 (2021).
  * [12] J-Y. Zhou,Y-H. Zhou, X-L. Yin, J-F. Huang, and J-Q. Liao, Sci. Rep. 10, 12557 (2020).
  * [13] F. J. Garcia-Vidal, C. Ciuti, and T. W. Ebbesen, Science 373, 34244383 (2021).
  * [14] C. M. Wilson, G. Johansson, A. Pourkabirian, M. Simoen, J. R. Johansson, T. Duty, F. Nori, and P. Delsing, Nature 479, 376 (2011).
  * [15] J. R. Johansson, G. Johansson, C. M. Wilson, and F. Nori, Phys. Rev. Lett. 103, 147003 (2009).
  * [16] B. Wang, J-M. Hu, V. Macrí, Z-L. Xiang, and F. Nori, arXiv:2204.08238 [quant-ph] (2022).
  * [17] A. Ridolfo, M. Leib, S. Savasta, and M. J. Hartmann, Phys. Rev. Lett. 109, 193602 (2012).
  * [18] K. R. Brown, C. Ospelkaus, Y. Colombe, A. C. Wilson, D. Leibfried, and D. J. Wineland, Nature 471, 196 (2011).
  * [19] J. Eisert, M. B. Plenio, S. Bose, and J. Hartley, Phys. Rev. Lett. 93, 190402 (2004).
  * [20] A. R. Urzúa, I. Ramos-Prieto, F. Soto-Eguibar, V. Arrizón, and Héctor M. Moya-Cessa, Sci. Rep. 9, 16800 (2019).
  * [21] J. Chen, Z. Li, X-Q Luo, W. Xiong, M. Wang, and H-C. Li, Opt. Express 29, 32639 (2021).
  * [22] G. Giedke, B. Kraus, M. Lewenstein, and J. I. Cirac, Phys. Rev. A 64, 052303 (2001).
  * [23] A. Merdaci and A. Jellal, Phys. Lett. A, 384, 126134 (2020).
  * [24] R. Hab-Arrih, A. Jellal and A. Merdaci, Int. J. Geom. Methods Mod. Phys. 18, 08, 2150120 (2021).
  * [25] G. Manzano, F. Galve, and R. Zambrini, Phys. Rev. A 87, 032114 (2013).
  * [26] X. Deng, C. Tian, X. Su, and C. Xie, Sci. Rep. 7, 44475 (2017).
  * [27] D. E. Bruschi, C. Sabín, and G. S. Paraoanu, Phys. Rev. A 95, 062324 (2017).
  * [28] J. C. Gonzalez-Henao, E. Pugliese, S. Euzzor, S. F. Abdalah, R. Meucci, and J. A. Roversi, Sci. Rep. 5, 13152 (2015).
  * [29] D. Stefanatos, Quantum Sci. Technol. 2, 014003 (2017).
  * [30] R. Hab-arrih, A. Jellal, D. Stefanatos, E. H. El Kinani, Virtual excitations and quantum correlations in ultra-strongly coupled harmonic oscillators under intrinsic decoherence (2022), submitted to JPA.
  * [31] A. Streltsov, G. Adesso, M. Piani, and D. Bruß, Phys. Rev. Lett. 109, 050503 (2012).
  * [32] G. L. Giorgi, Phys. Rev. A 84, 054301 (2011).
  * [33] D. X. Macedo and I. Guedes, J. Math. Phys. 53, 052101 (2012).
  * [34] M. A. Lohe, J. Phys. A: Math. Theor 42, 035307 (2009).
  * [35] G. J. Milburn, Phys. Rev. A 44, 5401 (1991).
  * [36] G. J. Milburn, New J. Phys. 8, 96 (2006).
  * [37] A. R. Urzúa and H. M. Moya-Cessa, Pramana 96, 72 (2022).
  * [38] W. H. Louisell, Quantum statistical properties of radiation (Wiley Classics Library Edition Published, 1990).
  * [39] G. Adesso, S. Ragy, and A. R. Lee, Open Syst. Inf. Dyn. 21, 1440001 (2014).
  * [40] G. Adesso, A. Serafini, and F. Illuminati, Phys. Rev. A 70, 022318 (2004).
  * [41] T. J. Osborne and F. Verstraete, Phys. Rev. Lett. 96, 220503 (2006).
  * [42] V. Coffman, J. Kundu, and W. K. Wootters, Phys. Rev. A 61, 052306 (2000).
  * [43] X-L. Zong, H-H. Yin, W. Song, and Z-L. Cao, Front. Phys. 10, 880560 (2022).
  * [44] R-C. Yang, G. Li, and T-C. Zhang, Quantum Inf. Process. 12, 493 (2013).