\section{\hspace{-3mm} \large{Runtime Energy-Accuracy Optimization}}  \label{sec:problem}

\subsection{Preliminaries}
\vspace{-0.5mm}

We consider human activity monitoring applications implemented on energy-constrained IoT devices.
We denote the period over which the total energy budget is provided as $T_P$, as summarized in Table~\ref{tab:problem_symbols}.
\textit{REAP} computes the energy allocations at runtime with a period of $T_P$, which is set to one hour in our experiments.
If the energy consumption over this period exceeds the amount of harvested energy and remaining battery level,
the device powers down and misses user activity.
Hence, our goal is to maximize the active time and the
expected accuracy over a given period $T_P$.

Suppose that the IoT device can operate at $N$ distinct DPs.
The recognition accuracy achieved by design point $i$ is denoted by $a_i$,
while the corresponding power consumption is given as $P_i$ for $1\leq i \leq N$.
In addition to these design points, we denote the time that the device remains off as $t_{off}$.
Finally, the power consumption during the off period,
which is due to the energy harvesting and the battery charging circuitry, is denoted by $P_{off}$.

\renewcommand{\arraystretch}{0.90
\begin{table}[]
	\caption{\small{Summary of symbols used in the optimization problem.}}
	\label{tab:problem_symbols}
	\vspace{-3mm}
	\begin{tabular}{@{}cc|cc@{}}
		\toprule
		Symbol & Description & Symbol & Description \\ \midrule
		$T_P$ & Activity period & $J(t)$ & Objective function \\ \midrule
		$E_b$ & Energy budget & $t_i$ & Active time of $DP_i$ \\
		\midrule
		$a_i$ & \begin{tabular}[c]{@{}c@{}}Recognition accuracy\\ of
		$DP_i$\end{tabular} & $\alpha$ &
		\begin{tabular}[c]{@{}c@{}}Accuracy-active time\\ trade-off
		parameter\end{tabular} \\ \midrule
		$P_{off}$ & \begin{tabular}[c]{@{}c@{}}Power consumption\\ in the off
		state\end{tabular} & $P_i$ & \begin{tabular}[c]{@{}c@{}}Power
		consumption\\ of $DP_i$\end{tabular} \\ \bottomrule
	\end{tabular}
\vspace{-3mm}
\end{table}
\vspace{-1mm}

\subsection{Optimization Problem Formulation}
In a given activity period, the system may operate at different design points,
resulting in varying levels of active time and accuracy.
Let $t_i$ denote the amount of time DP $i$ is utilized during $T_P$.
The active time of the device is simply given by the sum of the active times of each DP: $\sum_{i=1}^N t_i$.
Likewise, the expected accuracy over the activity period can be expressed as
$E\{a\} = \frac{1}{T_P}\sum_{i=1}^{N} a_it_i$.
The expected accuracy is a useful metric that incorporates both active time and accuracy,
but it does not allow emphasis of one over the other.
Therefore, we define a generalized cost function and solve the following optimization problem:
\vspace{-2mm}
\begin{align}\label{eq:objective}
& \text{$maximize$} & J(t) = \frac{1}{T_P}\sum_{i=1}^{N} a_i^\alpha t_i &   \\
\label{eq:constraint1}
& \text{$subject~to$} &  t_{off} + \sum_{i=1}^{N} t_i &= T_P
\\\label{eq:constraint2}
& \text{$and$} &  P_{off}t_{off} + \sum_{i=1}^{N}P_it_i \hspace{0.5mm} & \leq E_b
\\\label{eq:constraint3}
&& t_i &\geq 0~~~~\hspace{1mm} 0\leq i \leq N
\end{align}
\vspace{-5mm}

\noindent \textbf{Objection function $J(t)$:}
The parameter $\alpha$ in Equation~\ref{eq:objective}
enables a smooth trade-off between the active time and accuracy.
When $\alpha = 1$, the objective function reduces to the expected accuracy.
Similarly, when $\alpha=0$, the objective function reduces to total active time.
In general, the objective function gives a higher weight to the active time when $\alpha < 1$.
In contrast, the design points with higher accuracy are preferred for $\alpha > 1$.

\noindent\textbf{Constraints:}
The constraint given in Equation~\ref{eq:constraint1} states that the
sum of the active times and off period equals to the overall period $T_P$.
Similarly, Equation~\ref{eq:constraint2} specifies the energy budget constraint.
The left-hand side gives the sum of the energy consumed in the off state and active states.
The energy consumption of the $i^\mathrm{th}$ DP is given by the product of its power consumption $P_i$ and the active time allocated to it $t_i$.
Energy budget $E_b$ on the right-hand side is determined by energy allocation
techniques using the expected amount of harvested energy and battery
capacity~\cite{kansal2007power,bhat2017near}.
Finally, Equation~\ref{eq:constraint3} ensures that all active times
are non-negative.

\vspace{-1mm}
\subsection{Runtime Optimization Algorithm}

The solution of the proposed optimization problem
provides the duration of active times for each design points
such that they can collectively maximize the objective function in
Equation~\ref{eq:objective}.
It is important to solve this problem at runtime because the available energy budget is not known at design time.
Furthermore, the importance given to accuracy versus active time (i.e., $\alpha$) may change due to user preferences.

The optimization objective and the constraints in Equations~\ref{eq:objective}, \ref{eq:constraint1} and \ref{eq:constraint2}
are linear in the decision variables $t_i$ and $t_{off}$ for $1\leq i \leq N$.
Therefore, we use a procedure based on the simplex
algorithm~\cite{cormen2009introduction}, as
outlined in Algorithm~\ref{real_algorithm}.
It takes the energy budget $E_b$, Pareto-optimal DPs,
and the maximum number of iteration as input arguments.
The output is a vector with the
values of decision variables $t_i$ and $t_{off}$ for $1\leq i \leq N$ that
maximize the objective value.
We start the optimization process by constructing
a tableau with the initial conditions. The first row of the tableau describes
the objective function, while the other rows describe the constraints.
In each iteration of the procedure, we first find the pivot column by finding the column with the largest value in the last row of
the tableau.
Using the pivot column, we next find the pivot row in the tableau
in Line 8 of the algorithm.
Then, we update the tableau using the pivot column and row. The procedure is
terminated when all the entries in the last row are non-positive. In this case,
the pivot column is set as negative and the
optimal solution is returned.
\textit{Our implementation takes 1.5~ms on our IoT device
prototype with five design points}. We also observe that the algorithm takes
only 8~ms for up to 100 design points.
Since we run the optimization algorithm every hour, it takes
a negligible portion of the activity period and energy budget.

\vspace{-2.5mm}
\begin{algorithm}[h]
	\SetKwData{Left}{left}\SetKwData{This}{this}\SetKwData{Up}{up}
	\SetKwFunction{Union}{Union}\SetKwFunction{FindCompress}{FindCompress}
	\SetKwInOut{Input}{Input}\SetKwInOut{Output}{Output}
	\caption{The \textit{REAP} Procedure} \label{real_algorithm}
	\Input{Design points, energy budget $E_b$, max. iterations}
	\Output{Time allocated to each design point}
	\BlankLine
	\SetAlgoLined
	Initialize the $tableau$ with objective function and constraints  \\
	Add slack variables for inequality constraints \\
	\While {iter $\le$ max. iterations}  {
		$PivotCol \leftarrow$  $findPivotCol(tableau)$\\
		\If{$PivotCol < 0$}{\Return Optimal Solution}
		$PivotRow \leftarrow findPivotRow(tableau, PivotCol)$ \\
		Update the $tableau$ using the $PivotCol$ and $PivotRow$ \\
	}
\vspace{-1mm}
\end{algorithm}

\subsection{Accuracy -- Active Time Trade-off Analysis}

\textcolor{black}{
Next, we analyze how \textit{REAP~} can exploit trade-off between the accuracy and active time using the parameter $\alpha$ in objective function $J(t)$ in Equation~\ref{eq:objective}.
Since Section~\ref{sec:impact_allocated_energy} considered the expected accuracy ($\alpha$= 1),
this section considers $\alpha > 1$, which gives more emphasis for higher accuracy.}

As a representative example, Figure~\ref{fig:alpha2_4} shows the comparison of objective values of the five design points with \textit{REAP~} when $\alpha$ is set to 2.
\textit{REAP~} always achieves higher performance than the lowest energy design DP5, since accuracy is given higher weight.
The difference between \textit{REAP~} and  DP5 increases further as $alpha$ grows.
When the energy budget is less than 6~J, DP4 outperforms all the other DPs,
while \textit{REAP~} successfully matches it.
In contrast, DP1, DP2, and DP3 have a very low performance, since they are mostly in the off state.
When the energy budget exceeds 6~J, there is sufficient energy to provide a higher accuracy, but DP4 cannot exploit it.
Hence, the higher accuracy design points become affordable and start outperforming DP4 one by one.
Notably, \textit{REAP~} consistently outperforms or matches the static DPs, as we have
also observed in Figure~\ref{fig:accuracy_act_time_Ec}.
For example, DP3 is able to provide the same performance as \textit{REAP~} when the energy budget is 6.5~J.
As the energy allocation increases beyond 6.5~J,
\textit{REAP~} starts outperforming DP3 by optimally switching between DP1, DP2 and DP3.
This trend continues until the energy allocation reaches 9.9~J,
beyond which there is sufficient energy to support DP1 alone. Thus, \textit{REAP~} reduces to DP1 in this region.
In summary, \textit{REAP~} exceeds or matches the
performance of any individual DP.

\begin{figure}[t]

	\centering
	\includegraphics[width=0.89\columnwidth]{figures/alpha2_plot.pdf}
	\vspace{-4mm}
	\caption{The objective value $J(t)$ given in
				Equation~\ref{eq:objective} of static design points is normalized over
				$J(t)$ of \textit{REAP} with $\alpha$= 2.}
	\label{fig:alpha2_4}
	\vspace{-4mm}
\end{figure}

\section{Experimental Evaluation} \label{sec:experiments}

\vspace{-0.5mm}
\subsection{Experimental Setup} \label{sec:experimental_setup}
\vspace{-0mm}

\noindent\textbf{IoT device:} We use a custom prototype
based on the TI-Sensortag IoT platform~\cite{ TI_sensortag}
to implement the proposed design points.
The prototype consists of a TI CC2650 MCU, Invensense MPU-9250 motion sensor unit,
a stretch sensor and energy harvesting circuitry.
Sensors are sampled at 100~Hz and the MCU runs at 47~MHz frequency.
Power measurements from the prototype and data from 14 user subject studies
are used to obtain the 24 design points.
This data will be released to the public at \textit{blind-url}.

\noindent\textbf{Energy harvesting data:}
We use the solar radiation data measured by the NREL Solar
Radiation Research Laboratory
to obtain the energy harvesting profile from January 2015 to October
2018~\cite{NREL}.
We use the profile for each hour within this data to generate the energy budget.
These energy budgets are then used to evaluate \textit{REAP} and the static design points in Section~\ref{sec:case_study}.

\input{files/energy_accuracy_cooptimization.tex}

\input{files/activetime_accuracy_tradeoff.tex}

\input{files/optimization_solar_data.tex}

\section{Introduction}
\vspace{-0.5mm}
Wearable low-power internet-of-things (IoT) devices
enable a wide range of health monitoring, activity tracking, and wide area sensing
applications~\cite{liu2017gazelle,lara2013survey,bonomi2009detection}.
These devices must stay on for as long as possible to observe the user activities.
At the same time, they have to provide the maximum quality of service, such activity recognition accuracy.
These two objectives compete with each other since higher accuracy comes at the cost of larger energy consumption.
Since weight and form-factor constraints prohibit large batteries, feasibility
of these devices depends critically on optimizing the energy-accuracy trade-off
optimally at runtime.

\vspace{1mm}
Widely used dynamic power management techniques optimize the power-performance trade-off
by switching between different power states at runtime~\cite{cochran2011pack}.
High-performance states are used to execute computationally heavy
workloads at the expense of larger power consumption.
In contrast, low-performance states are used during light workloads to save
power.
In analogy, energy-accuracy trade-off in self-powered devices can be optimized by utilizing multiple design points.
\textcolor{black}{This is a challenging proposition
since characterizing the accuracy analytically
is much harder than developing power consumption and performance models.
For example, we consider an activity recognition application, where a wearable device infers the user activities,
such as jogging, by processing motion sensor data.
The recognition accuracy is a strong function of the users.
Hence, energy-accuracy optimization of requires
 \textit{user studies} and \textit{optimally chosen design points},
in addition to a \textit{runtime optimization algorithm}
that utilizes multiple design points.}

\vspace{1mm}
This paper presents a Runtime Energy-Accuracy oPtimization framework (\textit{REAP}) for energy-constrained IoT devices.
While our framework is general, we focus on health and activity monitoring applications
where a  wearable device processes sensor inputs to infer user activities.
The recognized activities are sent
to a gateway, such as a smartphone, for further processing.
\textit{REAP~} co-optimizes the accuracy and active time under a tight energy budget.
This optimization is enabled by the following three contributions.

\noindent \textbf{User studies for accuracy evaluation:}
We perform experiments with 14 users to recognize six activities:
\textit{sit}, \textit{stand}, \textit{walk}, \textit{jump}, \textit{drive},
\textit{lie down} and \textit{transitions} among them.
During these experiments, we collect 3-axis accelerometer and stretch sensor data.
We obtain a total of \textcolor{black}{3553 activity windows} from these
experiments.
After labeling, we utilize this data for evaluating the accuracy
of the human activity classifiers used in this work.

\noindent \textbf{Pareto-optimal design points:}
A common baseline in activity monitoring applications is to obtain
a classifier with the highest recognition accuracy.
High accuracy is obtained by using a sophisticated set of sensors,
features, and classification algorithms,
all of which imply a larger energy consumption, hence, lower active time.
Other design points can be obtained by reducing the number of sensors and feature set to save energy.
In turn, the energy savings lead to longer active time under a given harvested energy budget, albeit with lower accuracy.
\textit{To enable this work, we implemented 24 design points (DPs) with varying energy-accuracy trade-offs} on our hardware prototype.
Among them, we choose five Pareto-optimal DPs as our primary designs used at runtime.
We provide detailed execution time and power consumption breakdown for sensing, feature generation and processing steps
for each of these five DPs.

\noindent \textbf{Runtime optimization algorithm:}
Given an energy budget, two fundamental objectives are to maximize the recognition accuracy
and the amount of time the device is on, i.e., the active time.
We first formulate this co-optimization problem assuming that
there are $N$ design points with different energy-accuracy trade-offs.
We define a general objective function that enables us to tune the importance of active time versus recognition accuracy.
Then, we propose an efficient runtime algorithm that determines how much each
DP should be used so as to optimize the accuracy-active time
trade-off.
Our solution reveals the amount of time the device must operate in each of these design points.

\vspace{1mm}
Experimental results using a custom prototype based on TI Sensortag~\cite{TI_sensortag} IoT board
show that \textit{REAP} outperforms all static design points under
a range of energy budget constraints.
\textit{REAP~} achieves both 46\% higher expected accuracy and
66\% longer active time compared to the highest performance DP.
\textit{REAP} also achieves comparable active time to the lowest energy design
points
while providing significantly higher expected accuracy.
This makes \textit{REAP} suitable for use in a wide range of energy harvesting profiles.

\textit{The major contributions of this paper are as follows:}
\begin{itemize}
\vspace{-2mm}
\item A runtime technique to co-optimize the accuracy and active time
of energy-harvesting IoT devices.

\item Pareto-optimal design points with varying energy-accuracy
trade-offs for human activity recognition (HAR).

\item Experiments on a custom prototype with 14 user studies
that show significant improvements both in expected accuracy and active time compared to static design points.

\end{itemize}

\vspace{-1mm}
In the rest, Section~\ref{sec:related_work} reviews the related work.
Section~\ref{sec:problem} and Section~\ref{sec:design_points} present
the accuracy-active time optimization problem and DPs used in this paper, respectively.
Finally, Section~\ref{sec:experiments} presents the experimental results
and Section~\ref{sec:conclusions} concludes the paper.

\section{Conclusions} \label{sec:conclusions}

This paper presented a runtime accuracy-active time optimization technique for
energy-constrained IoT devices.
The proposed approach dynamically chooses design points with different
energy-accuracy trade-offs to co-optimize the accuracy and active time under
energy budget constraints.
To demonstrate the effectiveness in a realistic setting, we implemented a human
activity recognition application on a custom IoT prototype.
We presented five Pareto-optimal design points with different energy-accuracy trade-offs.
We achieve 46\% higher expected accuracy and 66\% longer
active time compared to the highest performance design point,
and 22\% to 29\% higher accuracy than low-power design points without sacrificing the active time.

\vspace{-1mm}

\subsection{Duty Cycle Optimization} \label{sec:optimization}

\section{Related Work} \label{sec:related_work}
\vspace{-0.5mm}

Energy harvesting for IoT devices has recently received significant attention
due to their small form-factor and low capacity batteries~\cite{sudevalayam2011energy,jayakumar2014powering}.
These devices can be broadly categorized into two classes.
The first class of devices rely solely on harvested energy and turn off when no energy is
harvested~\cite{shenck2001energy}
The second class of devices uses a small battery as a backup
to extend the active time~\cite{kansal2007power,vigorito2007adaptive,buchli2014dynamic}.
These approaches manage the power consumption of the device such that the total
energy consumed over a finite horizon is equal to the harvested energy.
This ensures a long device lifetime without battery replacement or manual
charging.
\textcolor{black}{\textit{REAP} is applicable to all devices that operate under
a fixed energy budget}.

Intermittent nature of ambient energy sources necessitates the development of
energy allocation and duty cycling
algorithms~\cite{buchli2014dynamic,shaikh2016energy,vigorito2007adaptive,bhat2017near}.
For example, the work in~\cite{kansal2007power} uses linear programming to determine the duty
cycle of the application as a function of the harvested energy.
Similarly, the algorithm in~\cite{vigorito2007adaptive} uses a linear quadratic
controller to assign the duty cycle of the device while maintaining a set
battery level.
An algorithm for dynamic power and energy management of an energy-harvesting
node for long-term energy-neutral operation is presented
in~\cite{buchli2014dynamic}.
However, these approaches choose between on- and off power states, which leads to sub-optimal operation.
Furthermore, they do not have a notion of accuracy or any concrete application,
unlike our work.

Human activity and health monitoring using wearable devices
have been an active area of research due to their potential benefits to sports,
patients with movement disorders, and
elderly~\cite{jafari2007physical,espay2016technology,lara2013survey,bhat2018online}.

A recent work proposes a wearable system for mobile analysis of running using motion sensors~\cite{liu2017gazelle}.
The authors selectively identify the best sampling points to maintain high accuracy while reducing sensing and analysis energy overheads.
The work in~\cite{jafari2007physical} presents a framework to
detect falls by using a wearable device equipped with accelerometers.
Authors in~\cite{bonomi2009detection} design a classifier that detects physical
activity using a body-worn accelerometer.
While this study offers an accurate classifier for human activity recognition, it cannot sustain
operation under tight energy budget constraints.
Based on this observation, we find Pareto-optimal design points for the HAR application
that offer varying levels of accuracy and energy consumption. Then, we use
these design points to maximize the expected accuracy of HAR.

\fontdimen2\font=0.72e
In summary, we present a unique combination of (1) a
runtime energy-accuracy optimization technique, and
(2) experimental evaluation with five concrete design points for HAR.
\textit{We will release our detailed power--performance characterizations and
user subject data to stimulate research in these areas}.

\section{\hspace{-3mm} \large{Human Activity Recognition Case Study}}\label{sec:design_points}

\definecolor{Gray}{gray}{0.9}
\newcolumntype{g}{>{\columncolor{Gray}}c}
\renewcommand{\arraystretch}{0.8
\begin{table*}[t]
	\footnotesize
	\centering
	\caption{\small{Accuracy, execution time, power and energy consumption of different human activity recognition application design points.}}
	\vspace{-4mm}
	\label{tab:design_points}
	\resizebox{\linewidth}{!}{
		\begin{tabular}{@{}ccg|cccc|cc|gg@{}}
			\toprule
			\multicolumn{3}{c}{Design point description} \vline &
			\multicolumn{4}{c}{MCU exec. time distribution (ms)} \vline

			& \multicolumn{4}{c}{Per activity summary} \\ \midrule
			\begin{tabular}[c]{@{}c@{}}DP \\ no.\end{tabular} & Features &
			\begin{tabular}[c]{@{}c@{}}Accuracy\\ (\%)\end{tabular} &
			\begin{tabular}[c]{@{}c@{}}Accel. \\ features\end{tabular} &
			\begin{tabular}[c]{@{}c@{}}Stretch \\ features\end{tabular} &
			\begin{tabular}[c]{@{}c@{}}NN\\ classifier\end{tabular} & Total &

			\begin{tabular}[c]{@{}c@{}}MCU\\ energy\\ (mJ)\end{tabular} &
			\begin{tabular}[c]{@{}c@{}}Sensor \\ energy\\ (mJ)\end{tabular} &
			\begin{tabular}[c]{@{}c@{}}Energy\\ (mJ)\end{tabular} &
			\begin{tabular}[c]{@{}c@{}}Power\\ (mW)\end{tabular} \\ \midrule
			1 & \begin{tabular}[c]{@{}c@{}} Statistical acceleration,\\ 16-FFT
				stretch\end{tabular} & 94 & 0.83 & 3.83 & 1.05 & 5.71 & 2.38 &
			2.10
			& 4.48 & 2.76 \\ \midrule
			2 & \begin{tabular}[c]{@{}c@{}}Statistical y-axis accel.,\\ 16-FFT
				stretch\end{tabular} & 93 & 0.27 & 3.83 & 1.00 & 5.10  & 2.29 &
			1.43 & 3.72 & 2.30 \\ \midrule
			3 & \begin{tabular}[c]{@{}c@{}}Statistical x- and y-axis \\
				accel. (0.8~s), 16-FFT stretch
			\end{tabular} & 92 & 0.27 & 3.83 & 0.90 & 5.00 & 2.10 &
			0.84 & 2.94 & 1.82 \\ \midrule
			4 & \begin{tabular}[c]{@{}c@{}}Statistical y-axis\\
				accel. (0.6~s), 16-FFT stretch
			\end{tabular} & 90 & 0.14 & 3.83 & 1.00 & 4.97 & 2.09 &
			0.57
			& 2.66 & 1.64 \\ \midrule
			5 & 16-FFT stretch & 76 & 0.00 & 3.83 & 0.88 & 4.71 & 1.85 & 0.08 &
			1.93 & 1.20 \\ \bottomrule
			\vspace{-3mm}
		\end{tabular}
	}
\end{table*}

\textit{REAP} is broadly applicable to energy-harvesting IoT devices that operate under a fixed energy budget.
In order to illustrate the optimization results on a real example,
we employ human activity recognition, i.e., HAR, as a driver application.

\vspace{-2mm}
\subsection{\hspace{-1mm}Background and Baseline Implementation}\label{sec:har_background}

There is a steady increase in the use of wearable and mobile devices
for the  treatment of movement disorders and obesity-related
diseases~\cite{espay2016technology}.
This technology enables data collection while the patients perform their daily activities.
The first step in this effort is to understand what activity the user is
performing at a given time.
For example, the gait quality of the patient cannot be checked unless we know the user is walking.
Therefore, HAR on mobile devices has recently attracted
significant attention~\cite{lara2013survey}.

We implement a HAR application on a custom prototype based on the
TI-Sensortag IoT board~\cite{TI_sensortag} and a passive stretch sensor.
Figure~\ref{fig:har_overview} shows an overview of our HAR application. It
starts with sampling of the accelerometer and stretch sensors. The streaming
sensor data
is fed to the TI-CC2650 MCU to generate the feature vector. The feature vector
is then processed by a parameterized neural network to infer the activity of
the user.
Finally, the inferred activity is transmitted to a host device, such as a
phone, using the Bluetooth Low Energy~(BLE) protocol.

\begin{figure}[t]
	\centering
	\includegraphics[width=1\linewidth]{figures/har_overview.pdf}
	\vspace{-7mm}
	\caption{\small{Overview of the human activity recognition application.}}
	\label{fig:har_overview}
    \vspace{-4mm}
\end{figure}

The energy consumption and recognition accuracy of HAR depends on the types of
sensors used, active time of the sensors, the type of features and the
complexity of the classifier.
Figure~\ref{fig:dp_combinations} shows the trade-off in energy and accuracy of
HAR for different
choices of sensors, features and classifiers. The left side of
Figure~\ref{fig:dp_combinations} shows the different configurations available for the sensors.
For instance, we can use all three axes of the accelerometer or turn off
selected axes to lower the energy consumption.
In the extreme case, we can
turn off the accelerometer to completely eliminate its energy consumption.
Once we choose the configuration of the sensors,
we can choose the sensing period i.e. the time for which sensors are active for each activity duration.
By default, the sensors are turned on during the full activity duration.
Turning off the sensors early, such as after 50\% of the activity
duration, provides energy savings at the cost of missed data points, hence accuracy.
We can also control the complexity of the features to trade-off accuracy and energy consumption. Complex features, such as Fast Fourier Transform~(FFT) and
Discrete Wavelet Transform~(DWT), offer a higher accuracy at the expense of a higher energy consumption. In contrast, statistical features have a lower energy consumption, albeit with a lower accuracy.
Finally, the structure and depth of the NN classifier can be controlled to obtain further energy-accuracy trade-off, as illustrated in Figure~\ref{fig:dp_combinations}.

We exploit this trade-off between energy and accuracy to design 24 different DPs implemented on the TI-Sensortag based prototype,
as described in the following section.

\begin{figure}[h]
	\centering
	\vspace{-2.5mm}
	\includegraphics[width=0.98\linewidth]{figures/dp_tradeoff_table.pdf}
	\vspace{-4mm}
	\caption{The knobs used to obtain design points with different energy-accuracy trade-offs.}
	\vspace{-3.5mm}
	\label{fig:dp_combinations}

\end{figure}

\vspace{-1mm}
\subsection{Pareto-Optimal Design Points}
We design a total of 24 DPs by exploiting the energy-accuracy trade-off illustrated in Figure~\ref{fig:dp_combinations}.
We start by using all the axes of the accelerometer, generating complex
features, and using an NN classifier with 3 hidden layers, which provide the highest recognition accuracy.
Then, we progressively reduce the number of axes of the accelerometer and sensing period  to reduce the energy consumption of the DPs.
We always use the passive stretch sensor in our DPs, since it has a
low energy consumption.
There is a need for detailed accuracy and energy consumption characterization of each DP to obtain the Pareto-optimal design points. To find the accuracy of each
design point, we performed experiments with 14 different users.
We have obtained a total of 3553 activity windows from the experiments and labeled each window with the corresponding activity.
Each DP is designed using 60\% of this data for training, 20\% for validation and the remaining 20\% for testing.

\begin{figure}[t]

	\centering
	\includegraphics[width=0.80\linewidth]{figures/dp_accuracy_vs_energy.pdf}
	\vspace{-4mm}
	\caption{\small

		The energy-accuracy trade-off of
		various design points. The dashed line connects the selected five
		design points.}
	\label{fig:dp_plot}
	\vspace{-6mm}
\end{figure}

All 24 design points are implemented on our prototype
to profile the execution time and measure the power consumption using the test pads on our prototype.
Figure~\ref{fig:dp_plot} shows the recognition accuracy and energy per
activity for each design point. As expected, each DP offers a
unique energy-accuracy trade-off.
For example, DP1 shows the highest accuracy with the highest energy consumption while DP5 shows the lowest
recognition accuracy and energy consumption. However, some
design points do not offer any benefit in the energy-accuracy trade-off. For
example, the design point marked with a red rectangle is dominated by DP2, DP3
and
DP4. Hence, we consider five Pareto-optimal design points shown using black diamonds~(DP1 to DP5) for validating
the proposed REAP algorithm.
\textit{Table~\ref{tab:design_points} summarizes the details of the configuration, accuracy, execution time and energy for five Pareto-optimal DPs}.
Next, we provide a short description of the Pareto-optimal DPs.

\vspace{1mm}
\noindent\textbf{Design Point-1~(DP1):}
DP1 offers the highest accuracy by utilizing all three axes of the
accelerometer for the entire activity window of 1.6~s.
It uses 16-FFT of the stretch sensor data and statistical features of the accelerometer, such as the mean and standard deviations.
DP1 leads to highest accuracy of 94\% at the cost of highest energy consumption
of 4.48~mJ per activity.
The energy break-down in Figure~\ref{fig:dp1_energy} shows that about 47\% of the energy consumption is due to the sensors.
Thus, reducing the sensor activity is an effective mechanism to save energy.

\begin{figure}[t]
	\vspace{-1mm}
	\centering
	\includegraphics[width=0.68\linewidth]{figures/dp1_pie_chart.pdf}
	\vspace{-3mm}
	\caption{\small{\textcolor{black}{Energy consumption distribution of DP1
	over one-hour activity period $T_P$.
	Total energy consumption is 9.9~J.}}}
	\label{fig:dp1_energy}
	\vspace{-3mm}
\end{figure}

\vspace{1mm}
\noindent\textbf{Design Point-2~(DP2):}
DP2 reduces the sensory energy  by utilizing only the y-axis of the
accelerometer along with the stretch sensor.
As depicted in Table~\ref{tab:design_points}, the energy consumption of the sensor is reduced to 1.43~mJ from 2.10~mJ.
It achieves an accuracy of 93\%, which is only 1\% lower than DP1.

\vspace{1mm}
\noindent\textbf{Design Point-3~(DP3):}
As shown in Figure~\ref{fig:dp_combinations}, reducing the sensing period leads to a lower energy consumption.
DP3 exploits this by sampling the x- and y- axes of the accelerometer for 50\% of each activity window, i.e., 0.8~s.
As a result, the energy consumption of the sensor is reduced to 0.84~mJ and the
total energy consumption of DP3 is reduced to 2.94~mJ per activity, while the
recognition accuracy drops to 92\%.

\vspace{1mm}
\noindent\textbf{Design Point-4~(DP4):}
DP4 is similar to DP3, except that sensing period of accelerometer is
further reduced to 40\% (0.6 s).
This reduces the energy consumption of DP4 to 2.66~mJ per activity with
recognition accuracy of 90\%.

\vspace{1mm}
\noindent\textbf{Design Point-5~(DP5):}
DP5 uses only the stretch sensor for data features to minimize energy consumption.
The energy consumption is reduced to 1.93~mJ per activity which is the lowest energy consumption among all our design points. However, it also shows the lowest recognition accuracy of 76\%.

\noindent\textbf{Offloading to a host:}
Finally, we note that the raw sensor data can be directly sent to a host device, such as a smartphone or server, for processing.
To assess the viability of this alternative, we implemented and measured its energy consumption.
Sending the raw sensor data over BLE consumes 5.5 mJ per activity without any
significant increase in the recognition accuracy.
In contrast, transmitting just the recognized  activity consumes only about 0.38 mJ per activity.
Hence, offloading is not an energy-efficient choice.

\subsection{\hspace{-2mm}Expected Accuracy and Active Time Analysis}
\label{sec:impact_allocated_energy}
\vspace{-0.5mm}

We first analyze the results of the proposed optimization approach as a function of the allocated energy over one-hour activity period $T_P$.
In the most energy-constrained scenario,
the minimum energy required to run the energy harvesting and monitoring circuitry is 0.18~J.
In the opposite extreme, 9.9~J energy is sufficient to run DP1,
the most power hungry design point, throughout $T_P$.
Therefore, we sweep the allocated energy starting with 0.18~J,
and find the optimal active time of each DP using the proposed
approach.

Figure~\ref{fig:accuracy_act_time_Ec}(a) shows the expected accuracy ($\alpha = 1$) as function of the energy budget.
The expected accuracy of all the design points approaches to zero
when the energy budget is close to 0.18~J,
since the device is almost always off.
As the energy budget increases (Region 1),
the accuracy of all DPs starts growing
since they can become active.
None of the design points can afford to stay 100\% active under the energy budget in Region 1.
We observe that design points with lowest energy consumption (DP5)
achieve significantly higher accuracy because they can stay in the
active state much longer.
\textit{REAP~} successfully matches or exceeds
the accuracy of DP5 under the most energy constrained scenario.
When the energy budget goes over 4.3~J,
DP5 can remain active throughout the activity period but its recognition
accuracy saturates.
The other DPs benefit from more energy in Region 2, while \textit{REAP~} outperforms all by utilizing them optimally.
At 5~J energy budget, for example, \textit{REAP} utilizes
DP4 42\% of the time and DP5 for 58\% of the time to optimize the expected
accuracy.
Finally, all design points can remain active throughout the activity period
when the energy budget is larger than 9.9~J.
Hence, their accuracy saturates, and \textit{REAP} reduces to DP1 beyond this
point.
In summary, \textit{REAP~} consistently outperforms or matches
the accuracy of all individual DPs
by utilizing \textit{multiple DPs} optimally.

\begin{figure}[t]

	\centering
	\includegraphics[width=0.95\linewidth]{figures/alpha1_v2.pdf}
	\vspace{-4mm}
	\caption{(a) Expected accuracy of \textit{REAP~} and
				design points.
				(b) Active time of each DP normalized to \textit{REAP~}.}
	\vspace{-5mm}
	\label{fig:accuracy_act_time_Ec}
\end{figure}

The active time of each DP \textit{normalized to} \textit{REAP~} is plotted in Figure~\ref{fig:accuracy_act_time_Ec}(b).
DP5 is expected to have the longest active time since it has the least energy consumption.
\textit{REAP} successfully matches its active time in all the regions.
In Region 1, \textit{REAP~} also achieves 2.3$\times$ larger active time
compared to DP1 while providing significantly better accuracy.
\textit{REAP} consistently provides longer active times compared to DP1, DP2 and DP3
until the energy budget becomes large enough
to sustain them throughout the activity period $T_P$.

\subsection{Case Study using Real Solar Energy Data}
\label{sec:case_study}
\vspace{-0.5mm}

In this section, we evaluate \textit{REAP~} under the real solar radiation data measured by NREL Solar Radiation Research Laboratory
at Golden, Colorado.
This data is used to calculate the amount of energy that can be harvested by a flexible solar cell~\cite{sp3-37} on our prototype.
Using the harvested energy budget, we compare the performance of \textit{REAP~}
against the static DPs over an entire month.
Figure~\ref{fig:one_month_plot} shows
the performance of \textit{REAP~} normalized to DP1, DP3, and DP5 as a function of $\alpha$.
Due to space limitation, we plot the DPs with the highest performance (DP1), lowest energy (DP5),
and best trade-off (DP3).
Our gains with respect to DP2 and DP4 are larger than that of DP3.

When active time is emphasized in the objection function ($\alpha=0.5$),
\textit{REAP~} outperforms DP1 by 1.4$\times$--2.2$\times$ with an average improvement of 1.6$\times$ across the month.
DP1 suffers the most in this case as it has the largest energy consumption among all the DPs.
Since accuracy becomes more important with larger $\alpha$,
the improvement of \textit{REAP~} over DP1 reduces.
However, we still obtain 1.1$\times$--1.3$\times$ improvement even for $\alpha=8$.
We observe a similar trend in improvements for \textit{REAP~} for DP3 as well.
The improvement is 1.1$\times$--1.4$\times$ for $\alpha=0.5$, and it
gradually decreases with larger $\alpha$.
The improvements over DP3 are relatively lower, since DP3 offers the best trade-off between
energy consumption and accuracy among our Pareto-optimal design points.

\begin{figure}[t]
	\vspace{-3mm}
	\centering

	\includegraphics[width=1\linewidth]{figures/alpha_up_to_8.pdf}

	\vspace{-5mm}
	\caption{Performance (i.e., J(t) ) achieved by \textit{REAP~} normalized to DP1,
	DP3, and DP5 during the month of September 2015. Error bars represent the
	range of improvement.}
	\label{fig:one_month_plot}

\end{figure}
Finally, we see that the improvements over DP5 follow the opposite trend.
When $\alpha = 0.5$, DP5 is able to match the active time of \textit{REAP~} due to its lower
energy consumption.
However, the performance of DP5 diminishes severely with increasing $\alpha$.
In summary, \textit{REAP~} can
provide a higher performance than any individual design point under any optimization objective.
If the user needs a higher accuracy, \textit{REAP~} can
successfully adapt to new requirements.
This can be utilized by the IoT device to tune
its performance as user needs change.


\section{Introduction}\label{sec1}

For decades,  the classical relativistic  field equations with  stable solitary wave solutions (solions\footnote{According to some well-known references such as  \cite{rajarama}, a solitary wave solution is a soliton when it reappears  without any distortion after collisions. The  stability is essentially    a necessary condition for a solitary wave solution to be a soliton. However, in this paper,  we only adopt  the stability condition for the definition of a soliton solution.}) have drawn the interest of many physicists   \cite{rajarama,Das,lamb,Drazin,TS}. In fact,  soliton solutions  behave like real particles as   they have the non-disperse  localized energy density functions and satisfy the standard relativistic energy-momentum relations. For example, the real nonlinear Klein-Gordon  (RNKG)\cite{mandel,peskin} systems in $1+1$ dimensions with kink (anti-kink) solutions \cite{rajarama,TS,phi45,OV,GH,MM1,MR,JRM1,JRM2,DSG1,DSG2,DSG3,MM2,waz,ana1,
ana2,Kink1,Kink2,Kink3,Kink4,Kink5,Kink6,Kink7,Kink8,Kink9,Kink10,Kink11,Kink12,Kink13,
Kink14,Kink15,Kink16,Kink17}, Skyrme model of baryons \cite{TS,SKrme,SKrme2,SKrme3,SKrme4} and 't Hooft Polyakov model which yields magnetic monopole solutions   \cite{rajarama,TS,toft,pol,MKP,TOF,CTOPO}
 in $3+1$ dimensions are three well-known  systems which yield stable solitary wave solutions or solitons. One should note that, all the three systems mentioned above, yield topological solitons and the topological feature  is the main reason  behind  their stability. With topological  solutions, there are generally  complicated conditions on the boundaries to have a multi particle-like solution. However,  with the non-topological solitary wave solutions,  each arbitrary multi particle-like  solution can be obtained easily  just by adding distinct  far apart  solitary wave solutions together.

There have  been many works on the non-topological solitary  wave solutions  so far  \cite{rajarama,Das,lamb,Drazin,TS,N1,N2,N3,N4,N5,N6,N7,N8,N9,N10, N11,N12,N13,Vak3,Vak4,Vak5,Vak6,Lee3,Scoleman,
R1,R2,R3,R4,R5,R6,R7,R8,VD}.  However, the famous relativistic non-topological solitary  wave solutions are   Q-balls \cite{Vak3,Vak4,Vak5,Vak6,Lee3,Scoleman,R1,R2,R3,R4,R5,R6,R7,R8,VD,JOP,Non,JOP2}.
With  non-topological solitary wave solutions, an  important   criterion for  examining   the stability  is  the classical (or  Vakhitov-Kolokolov) criterion \cite{Vak3,Vak4,Vak5,Vak6,Vak1,Vak2,Vak777}.  The classical stability criterion is based on  examining dynamical equations  when they are  linearized for  small fluctuations above the background of the solitary wave solution. If the linearized equation does not lead to any unstable growing mode, the solitary wave solution  is called  a stable solution classically.     There is another stability criterion  called the energetical stability criterion.   In fact,  a special  solitary wave solution is  energetically stable if any arbitrary variation above its background  leads to an increase in the total energy. In other words, an energetically stable solitary wave solution would be stable against any arbitrary deformation \cite{Derrick,MM4,MMM1,MM5,MMM2}.

A solitary wave solution which is energetically stable would be a single solution among the  other (close) solutions. For example, the   kinks (anti-kinks) as the well-known topological solitary wave solutions  of the real nonlinear Klein-Gordon (RNKG) systems are inevitably energetically stable \cite{rajarama,Derrick}. However, a solitary wave solution, which is classically  stable, is not necessarily an energetically stable solution; or it is not a  single solution with the minimum  energy among the other (close) solutions. For example, some of the Q-ball solutions are classically stable \cite{Vak3,Vak4,Vak5,Vak6,Vak777}, but none of the them are energetically stable \cite{MM5}.  In fact, the energetical  stability criterion  is at a higher level  than the classical stability criterion. In other words,  if a solitary wave solution is  energetically  stable, it would undoubtedly be classically  stable, as well. Moreover, if a solitary wave solution is not  classically stable, it would not be  an energetically stable solution, as well.

In this paper, in line with the previous works \cite{MM4,MMM1,MM5,MMM2},  we   introduce a special  RNKG  model in $1+1$ dimensions with a well-formed non-topological  solitary wave solution which is essentially unstable  \cite{rajarama,Vak3}. But we will show  how adding a proper term  to the original  RNKG Lagrangian density, transforms the special solitary wave solution (SSWS) into  an energetically  stable  object. We call this additional term the “\textit{ stability catalyzer}", because  it behaves as  a massless spook\footnote{We chose  the word “\emph{spook}" so that it will not be confused with words like   “\emph{ghost}" and “\emph{phantom}", which have meaning in the  literature.}  which surrounds the SSWS  and guarantees its energetical stability. In other words,  it prevents any change in the internal structure of  the SSWS, and  leaves  the dominant dynamical equations  and other properties  of the SSWS invariant. It should be noted that, we  consider the model in $1+1$ dimension just for the sake of  simplicity, as it can be extended to $3+1$ dimensions with some modifications.

There is a  parameter $B$  in the stability  catalyzer term, which  leads to more stability for the SSWS at larger values.  In other words,  the larger the values the greater will be the increase in the total energy for any arbitrary small variation above the background of the SSWS.
Hence, if one considers   a system with an extremely  large  value of parameter $B$, then the other  configurations of the fields (which are not very  close to the SSWS and the vacuum state) need extremely  large energies  to be created; meaning that,
the possible solutions  of the  system with the finite energies are only  the  free  far apart SSWSs, as   multi   particle-like solutions.

The present  paper has been organized as follows: In the next section,     we  set up the basic equations for the RNKG systems with a single scalar field and consider a special RNKG model  with a special  non-topological non-vibrational  solitary wave solution which is essentially unstable. In section 3,  we  introduce the   stability catalyzer term, which serves   to introduce an  extended KG system\footnote{The extended KG systems were introduced in Ref.~\cite{MM4}. } with  a single non-topological energetically stable  SSWS. Section 4 presents an in-depth  study of the   stability of the SSWS for any arbitrary small deformations according to the energetical stability criterion.  The last section is devoted to  summary and conclusion.

\section{Single  field RNKG systems in $1+1$ dimensions}\label{sec2}

The simplest form of the real nonlinear Klein Gordon (RNKG) systems in $1+1$ dimensions  can be introduced by the following Lagrangian density:
\begin{equation} \label{lag}
{\cal L}_{o}=\partial^\mu \varphi\partial_\mu \varphi-U(\varphi),
\end{equation}
in which $\varphi$ is a single real scalar field  and $U(\varphi)$ is called  the field potential. Note that,  we set the speed of light to one ($c=1$) throughout the paper for the sake of simplicity. Using the principle of least  action,  the related equation of motion is
\begin{equation} \label{eom}
\Box \varphi=\frac{\partial^{2}\varphi}{\partial t^2}-\frac{\partial^{2}\varphi}{\partial x^2}=-\frac{1}{2}\frac{dU}{d\varphi}.
\end{equation}
Using the Noether's theorem \cite{mandel,peskin}, one can simply obtain the energy-momentum tensor:
\begin{eqnarray} \label{ed}
T^{\mu\nu}=2\partial^{\mu}\varphi\partial^{\nu}\varphi-g^{\mu\nu}{\cal L}_{o},
\end{eqnarray}
 in which $g^{\mu\nu}$  is the Minkowski metric. The $T^{00}$ ($T^{10}=T^{01}$) component of this tensor is the same energy  (momentum) density function, which for the Lagrangian density (\ref{lag}), is simplified   to
\begin{eqnarray} \label{emt}
T^{00}=\varepsilon_{o}(x,t)=\dot{\varphi}^2+\varphi'^2+U(\varphi)\quad (T^{01}=2\dot{\varphi}\varphi'),
\end{eqnarray}
in which the dot and the prime are symbols for time and  space derivatives respectively. The  integration of $T^{00}$ ($T^{10}$) over the whole space yields the same total energy $E$ (momentum $P$) of the system and always remains constant.

In general, there is not a stable  non-topological non-vibrational solitary wave solution  for the RNKG systems in $1+1$ dimensions  \cite{rajarama,Vak3}. For example, if one considers a special  field potential  as follows:
\begin{equation} \label{po}
U(\varphi)= \varphi^4(1-\varphi^2),
\end{equation}
then the equation of motion (\ref{eom}), for a static (non-moving and  non-vibrational) solution $\varphi=\varphi_{o}(x)$,  is simplified   to
\begin{equation} \label{ste}
\frac{d^2 \varphi_{o}}{d x^2}=2\varphi_{o}^3-3\varphi_{o}^5,
\end{equation}
which has a  non-topological solution as follows:
\begin{equation} \label{ss}
\varphi_{o}(x)=\frac{\pm 1}{\sqrt{1+ x^2}}.
\end{equation}
Applying the Lorentz  transformations, the moving   version of this solution (\ref{ss}) can be obtained as well:
\begin{equation} \label{sss}
\varphi_{v}(x,t)=\varphi_{o}(\tilde{x})=\frac{\pm 1}{\sqrt{1+ \tilde{x}^2}},
\end{equation}
where $\tilde{x}=\gamma(x-vt)$, $\gamma=1/\sqrt{1-v^2}$ and $v$ is the velocity. However,  the field  potential (\ref{po}) for $|\varphi|>\frac{\sqrt{6}}{3}$ is  decreasing and for
$|\varphi|>1$ takes negative values. Therefore, the special solitary wave solution (\ref{sss}) is essentially unstable and without  violating  the conservation  energy law, the
effect of any small perturbation,  causes  the profile of the localized solution (\ref{ss}) to blow  up \cite{rajarama}(see Figs.~\ref{pp} and \ref{3}). In \cite{Vak3},  the instability of the non-topological non-vibrational solitary wave solutions  of the RNKG systems in $1+1$ dimensions  are generally  referred to the existence of the growing modes.

\begin{figure}[ht!]
   \centering
   \includegraphics[width=110mm]{pp.jpg}
   \caption{The non-topological static  solitary wave solution (\ref{ss}) of the RNKG system (\ref{po}) is essentially unstable and spontaneously blows  up. This figure shows a continuous representation of all the possible configurations of the spontaneous evolution of the  SSWS (\ref{ss}), which is initially considered to be at rest. It is obtained  from a finite difference method in Matlab for the  SSWS (\ref{ss}) as the initial condition  of the PDE (\ref{eom}).  } \label{pp}
\end{figure}

\begin{figure}[ht!]
   \centering
   \includegraphics[width=160mm]{3.jpg}
   \caption{ Instead of Fig.~\ref{pp}, many configurations of the field at discrete times  in the range $0<t<5$  can be plotted in a 2D figure for better understanding. All  the different configurations have the same energy equaling the rest energy of the SSWS    (\ref{ss}) at $t=0$.} \label{3}
\end{figure}

In general, since the theory is relativistic,  the same    well-known  relativistic relations between the moving and non-moving versions  of any arbitrary solution,  which has a  localized energy density function,   would be obtained,  meaning that:
 \begin{eqnarray} \label{d}
&&E=m = \int_{-\infty}^{+\infty} T^{00} dx=\int_{-\infty}^{+\infty} [\dot{\varphi_{v}}^{2}+\varphi_{v}'^2+U(\varphi_{v})]  dx=\gamma E_{o}=\gamma m_{o}, \\&&
P= \int_{-\infty}^{+\infty} T^{01}dx=2\int_{-\infty}^{+\infty} \dot{\varphi_{v}}\varphi_{v}' dx=\gamma m_{o} \textbf{v}.
\end{eqnarray}
where $E_{o}$ ($m_{o}$) is the same rest energy (mass) of the   solution. For the special solution (\ref{ss}), according to Eq.~(\ref{d}), the related  rest energy is $E_{o}=\pi/4$. Moreover, the width of any arbitrary moving solution is always smaller  than  its non-moving version, exactly according to the Lorentz contraction law.

\section{The stability catalyzer term for a SSWS}\label{sec3}

Here,  we attempt   to find a proper additional term for the original  Lagrangian density (\ref{lag}) that could  guarantee   the energetical  stability of the SSWS (\ref{ss}). However, similar to  a catalyzer,  we   expect that it has no role in the   dominant dynamical equation  and the other properties of the  SSWS (\ref{ss}). In other words, we first expect that the standard Eq.~(\ref{eom}) to remain  the dominant  dynamical equation  only  for the SSWS (\ref{ss}), and second,  the rest energy of the SSWS (\ref{ss}) to be a  minimum   among the energies of the  other (close) solutions. The  other complementary discussions are the same as those sufficiently presented in \cite{MM5}.

 However, we assume   a new extended KG  Lagrangian density  as follows:
\begin{equation} \label{lag2}
{\cal L}={\cal L}_{o}+F=\partial^\mu \varphi\partial_\mu \varphi-U(\varphi)+F,
\end{equation}
where $F$ is the same unknown additional (stability catalyzer) term, which should  be properly identified.
 We expect the new extended Lagrangian density  (\ref{lag2}) to be reduced  to the same original  version (\ref{lag}) just for the SSWS  (\ref{ss}), that is,   the new additional term $F$ should be zero only for the SSWS (\ref{ss}).  Note that,  the new extended system (\ref{lag2}) and the original RNKG system (\ref{lag}) are essentially different relativistic field  systems with   different solutions  except  for the SSWS (\ref{ss}), which is considered  to be a common solution. According to  the standard relativistic  Lagrangian densities  in physics, we expect  the unknown additional  scalar term $F$ to be  a function of the allowed scalars $\varphi$ and $\partial_{\mu}\varphi\partial^{\mu}\varphi$. However,  the new equation  of motion is
\begin{eqnarray} \label{new}
\left[\Box \varphi+\frac{1}{2}\frac{dU}{d\varphi}\right]+\frac{1}{2}\left[ \frac{\partial}{\partial x^{\mu}}\left(\frac{\partial F }{\partial (\partial_{\mu}\varphi)}\right)-\left(\frac{\partial F}{\partial \varphi}\right)  \right]=0.
\end{eqnarray}
For the SSWS (\ref{ss}) to  still remain  a solution of the new equation of motion (\ref{new}) (or  the new  equation of motion (\ref{new}) is reduced  to the same original version (\ref{eom})),  since the first part of this new equation, according to the  same original Eq.~(\ref{eom}), is satisfied automatically, i.e. $\left[\Box \varphi_{o}+\frac{1}{2}\frac{dU}{d\varphi_{o}}\right]=0$, and the functional $F$ is not essentially linear in  ${\cal L}_{o}$, so we conclude that the two distinct terms $\frac{\partial}{\partial x^{\mu}}\left(\frac{\partial F }{\partial (\partial_{\mu}\varphi)}\right)$ and $\frac{\partial F}{\partial \varphi}$ must be  zero independently for the SSWS (\ref{ss}).

To meet all these requirements, one can conclude that  $F$ must be a function of the powers of ${\cal K}$  (${\cal K}^{n}$'s with $n\geq3$), where    ${\cal K}$  is a  special scalar functional
\begin{equation}\label{k}
{\cal K}=\partial_{\mu}\varphi\partial^{\mu}\varphi+U(\varphi)=\dot{\varphi}^2-\varphi'^2+\varphi^{4}(1-\varphi^{2}),
\end{equation}
which is defined  to be zero  when we have a  SSWS (\ref{ss}).  For example, a simple choice   for the functional $F$ is
\begin{equation} \label{FF}
F=B{\cal K}^3,
\end{equation}
where $B$ is  a real positive number. For this special choice  (\ref{FF}), we obtain
\begin{eqnarray} \label{DF}
  && \frac{\partial}{\partial x^{\mu}}\left(\frac{\partial F }{\partial (\partial_{\mu}\varphi)}\right)= B{\cal K}\left[6{\cal K}\partial_{\mu}\partial^{\mu}\varphi+12\partial_{\mu}{\cal K}\partial^{\mu}\varphi \right], \nonumber  \\&&
\frac{\partial F}{\partial \varphi}=\left[3 B{\cal K}^{2}\frac{\partial {\cal K}}{\partial \varphi}\right],\nonumber
\end{eqnarray}
which are both obviously zero for the SSWS (\ref{ss}). In fact, each term on the right hand side of the above relations   contains a power of ${\cal K}$ and hence all are zero for the SSWS (\ref{ss}). Therefore, with this special choice (\ref{FF}), we can be certain  that the previous SSWS (\ref{ss}) is again a solution of the new extended system (\ref{lag2}), and the new dynamical equation (\ref{new}) is reduced to the same original  one (\ref{eom}), as its dominant dynamical equation.

The  energy density functional  of the new extended system (\ref{lag2}) can be obtained  easily
\begin{eqnarray} \label{nte}
T^{00}=\frac{\partial{\cal L}}{\partial \dot{\varphi}}\dot{\varphi}-{\cal L}=\left[\dot{\varphi}^{2}+\varphi'^2+U(\varphi)\right]+\left[B{\cal K}^{2}(6\dot{\varphi}^{2}-{\cal K})\right]=\varepsilon_{o}+\varepsilon_{1}.
\end{eqnarray}
According to Eq.~(\ref{k}), the second part of the energy density  (\ref{nte}) becomes
\begin{equation} \label{e1}
\varepsilon_{1}=B{\cal K}^{2}(5\dot{\varphi}^{2}+\varphi'^2+\varphi^{4}(\varphi^{2}-1)),
\end{equation}
which is zero  for the SSWS (\ref{ss}) and the vacuum state $\varphi=0$. However, it is not   a positive definite function,  because function $\varphi^{4}(\varphi^{2}-1)$ in the range $0<|\varphi|<1$ is negative. Hence, we cannot be certain   about the energetical  stability of the  SSWS (\ref{ss}).

In order to introduce  a new proper  additional term $F$ for which  the energetical stability of the SSWS (\ref{ss}) is properly guaranteed,  we have to use  a new  scalar field $\theta$ which  can be called the  phase field or the catalyzer field. However,  the new  proper additional  term $F$ can be introduced as follows:
\begin{equation} \label{FFF}
F=B\sum_{i=1}^{3}{\cal K}_{i}^3,
\end{equation}
where
 \begin{eqnarray} \label{e5}
  &&{\cal K}_{1}= \varphi^4 \mathbb{S}_{1},\\ \label{e52}&&
{\cal K}_{2}= \varphi^4 \mathbb{S}_{1}+\mathbb{S}_{2}, \\ \label{e53}&&
{\cal K}_{3}= \varphi^4 \mathbb{S}_{1}+\mathbb{S}_{2}+2 \varphi^2 \mathbb{S}_{3},
\end{eqnarray}
and
\begin{eqnarray} \label{sc2}
  &&\mathbb{S}_{1}=\partial_{\mu}\theta\partial^{\mu}\theta-1,\\ \label{sc22}&&
 \mathbb{S}_{2}=\partial_{\mu}\varphi\partial^{\mu}\varphi+\varphi^{4}(1-\varphi^{2}),\\ \label{sc23}&&
 \mathbb{S}_{3}=\partial_{\mu}\varphi\partial^{\mu}\theta.
\end{eqnarray}
In general, since  $\mathbb{S}_{i}$'s   are three independent scalars, it is not possible for them to be zero simultaneously except for  the non-trivial SSWS (\ref{ss}) with $\theta=\omega_{s}t= \pm t$. In other words, $\mathbb{S}_{1}=0$, $\mathbb{S}_{2}=0$ and $\mathbb{S}_{3}=0$ are three independent coupled nonlinear PDE's which do not have any non-trivial  common solutions except for the SSWS (\ref{ss}) with $\theta=\pm t$ (see the Appendix A). In fact, the same result applies to ${\cal K}_{i}$'s,  since ${\cal K}_{i}$'s are three independent linear combinations of the scalars  $\mathbb{S}_{i}$'s,  they are not zero simultaneously  except for  the SSWS (\ref{ss}) with  $\theta=\pm t$. Note that, for a moving version of the SSWS (\ref{sss}), which moves at the velocity of $v$, the  proper phase function $\theta$, for which all $\mathbb{S}_{i}$'s  would be zero simultaneously,  is $\theta=k_{\mu}x^{\mu}$, i.e. the boosted version of  function $\theta= \pm t$, provided
\begin{equation} \label{pro}
k^{\mu}\equiv (k^{0},k^{1})=(\omega,k)=(\omega,\omega v),
\end{equation}
where $\omega=\gamma\omega_{s}$ and $\omega_{s}=\pm 1$.

However, the dynamical equations of motion of the extended KG system  (\ref{lag2}) with the new additional term (\ref{FF}) can be obtained easily as follows:
\begin{eqnarray} \label{geq}
&&\Box \varphi+\dfrac{1}{2}\frac{dU}{d\varphi}+\frac{1}{2}\left[ \frac{\partial}{\partial x^{\mu}}\left(\frac{\partial F }{\partial (\partial_{\mu}\varphi)}\right)-\left(\frac{\partial F}{\partial \varphi}\right)  \right]=0,\\ \label{geq2} &&   \frac{\partial}{\partial x^{\mu}}\left(\frac{\partial F }{\partial (\partial_{\mu}\theta)}\right)=0.
\end{eqnarray}
Again, it is easy to show that all the different first and second derivatives of   $F$ (\ref{FFF}), which were seen in the Eqs.~(25) and (26),   for the SSWS (\ref{ss}) with $\theta=\pm t$, would be zero simultaneously. In other words, for the SSWS (\ref{ss}) with $\theta=\pm t$, Eq.~(26)  is automatically  satisfied  and  Eq.~(25)  is reduced to the same standard original version (\ref{eom}) as the dominant dynamical equation of the free SSWS (\ref{ss}). Note that the SSWS (\ref{ss}) in the new extended system (\ref{lag2}) must be now considered along with a    scalar field $\theta=\pm t$.   However, from here on in this paper,  the non-moving SSWS is as follows:
 \begin{equation} \label{SSWS}
\varphi_{s}(x)=\varphi_{o}(x)=\frac{\pm 1}{\sqrt{1+ x^2}},\quad\quad \theta_{s}(t)=\omega_{s}t=\pm t.
\end{equation}
Hence, the   moving version of the SSWS (\ref{SSWS}) would be
 \begin{equation} \label{SSWS2}
\varphi_{v}(x,t)=\varphi_{s}(\tilde{x})=\frac{\pm 1}{\sqrt{1+ \tilde{x}^2}},\quad\quad \theta_{v}(x,t)=\theta_{s}(\tilde{t})=k_{\mu}x^{\mu}=\omega t-kx.
\end{equation}
where $\tilde{t}=\gamma(t-vx)$ and $\tilde{x}=\gamma(x-vt)$.

The   energy density functional  of  the  extended  KG system  (\ref{lag2}) with the new additional  term (\ref{FFF})  is
\begin{eqnarray} \label{MTE}
&&\varepsilon(x,t)=T^{00}=\frac{\partial{\cal L}}{\partial \dot{\varphi}}\dot{\varphi}+\frac{\partial{\cal L}}{\partial \dot{\theta}}\dot{\theta}-{\cal L}=\varepsilon_{o}+\varepsilon_{1}+\varepsilon_{2}+\varepsilon_{3}=\nonumber\\&&
\quad\quad\quad\left[\dot{\varphi}^{2}+\varphi'^2+U(\varphi)\right]+B\sum_{i=1}^{3}\left[3C_{i}
{\cal K}_{i}^{2}-{\cal K}_{i}^3\right],
\end{eqnarray}
which are divided  into four distinct  parts and
\begin{equation}\label{cof}
C_{i}=\dfrac{\partial{\cal K}_{i}}{\partial \dot{\theta}}\dot{\theta}+\dfrac{\partial{\cal K}_{i}}{\partial \dot{\varphi}}\dot{\varphi}=
\begin{cases}
\quad\quad 2\varphi^4\dot{\theta}^{2} &  i=1
\\
2(\dot{\varphi}^{2}+\varphi^4\dot{\theta}^2) & i=2
\\
2(\dot{\varphi}+\varphi^2\dot{\theta})^{2}
 & i=3.
\end{cases}
 \end{equation}
 After a straightforward calculation, one can obtain:
 \begin{eqnarray} \label{eis}
&&\varepsilon_{1}=B{\cal K}_{1}^2[5\varphi^4\dot{\theta}^2+\varphi^4\theta'^2+\varphi^4]\geq 0,\\ \label{eis2}&&
\varepsilon_{2}=B{\cal K}_{2}^2[5\varphi^4\dot{\theta}^2+5\dot{\varphi}^2+\varphi^4\theta'^2+ \varphi'^2+\varphi^{6}]\geq 0, \\ \label{eis3}&&
\varepsilon_{3}=B{\cal K}_{3}^2[5(\varphi^2\dot{\theta}+\dot{\varphi})^2+(\varphi^2 \theta'+ \varphi')^2+\varphi^{6}]\geq 0.
\end{eqnarray}
All terms in the above relations are now positive definite, therefore all the $\varepsilon_{i}$'s ($i=1,2,3$) are positive definite functions and bounded from below by zero. All $\varepsilon_{i}$'s ($i=1,2,3$) are zero simultaneously just for the trivial vacuum state $\varphi=0$ and the non-trivial SSWS (\ref{SSWS}), just as we expected  from the catalyzer. Now, if parameter  $B$
is considered  to be a  large number, since  at least one of the ${\cal K}_{i}$'s  is a non-zero  function for any other solution,  then at least one of the  $\varepsilon_{i}$'s ($i=1,2,3$), which all contain parameter $B$, would be a  large positive function. It means that for  other solutions (except for the ones which are very close to the vacuum $\varphi=0$), the related  energies are always larger than the rest energy of the SSWS (\ref{SSWS}). Unlike $\varepsilon_{i}$'s ($i=1,2,3$), which are three absolute positive functions and are minimum for the SSWS (\ref{SSWS}), $\varepsilon_{o}$ is not  an absolute positive function and is not a minimum for the SSWS (\ref{SSWS}).
In the next section, we will show  the role of  $\varepsilon_{o}$ in the stability considerations, meaning that if we take an extended system (\ref{lag2}) with a large parameter $B$ (approximately $B>10^2$), it would be  completely ineffective.

\section{the energetically stability of the SSWS}

In general, a solitary wave solution (such as kink and anti-kink solutions of the RNKG  systems)  is  energetically stable if its rest energy  is  at a minimum   among the energies of the other (close) solutions. In other words, for an energetically stable solitary wave solution, any arbitrary deformation (variation)  above the background, leads to an increase in the total energy. In this section, we specifically  examine  the energetical  stability of  the SSWS  (\ref{SSWS}). In fact,   we are going to consider the variation of the total energy versus  any arbitrary small deformation above the background of the SSWS  (\ref{SSWS}), which is at rest. In general, any deformed version of the  SSWS (\ref{SSWS}) can be introduced  as follows:
\begin{equation} \label{so1}
\varphi(x,t)=\varphi_{s}(x)+\delta \varphi(x,t) \quad \textrm{and} \quad \theta(x,t)=\theta_{s}(t)+\delta \theta(x,t),
\end{equation}
where $\delta \varphi$ and $\delta \theta$  are considered as  arbitrary small  functions  of space-time.
Now, if we insert  (\ref{so1}) in   $\varepsilon_{o}(x,t)$ and  keep the terms
up to the second  order of small variation $\delta \varphi$, then it yields
\begin{eqnarray} \label{so3}
&&\varepsilon_{o}(x,t)=\varepsilon_{os}(x)+\delta\varepsilon_{o}(x,t)\approx\varepsilon_{os}(x)+\delta\varepsilon_{o1}(x,t)+
\delta\varepsilon_{o2}(x,t)=\left( \varphi_{s}'^2+U(\varphi_{s})\right)+\nonumber\\&&
\quad\quad 2\left(\varphi_{s}'  (\delta \varphi')+\frac{1}{2}\frac{dU(\varphi_{s})}{d\varphi_{s}}(\delta \varphi)\right)+\left((\delta\dot{\varphi})^2+ (\delta\varphi')^2 +\frac{1}{2}\frac{d^2U(\varphi_{s})}{d\varphi_{s}^2}(\delta \varphi)^2\right)
\end{eqnarray}
where $\varepsilon_{os}$, $\delta\varepsilon_{o1}$ and $\delta\varepsilon_{o2}$   are defined on the right hand side of the above equation, respectively.  $\varepsilon_{os}(x)=\left( \varphi_{s}'^2+U(\varphi_{s})\right)$ is the energy density function of the non-moving SSWS  (\ref{SSWS}).  $\delta\varepsilon_{o1}$ and $\delta\varepsilon_{o2}$ are  functionals  of the  first and second  order of the small  variation $\delta \varphi$, respectively.
Note that, for a non-moving SSWS (\ref{SSWS}), $\dot{\varphi_{s}}=0$, $\theta_{s}'=0$ and $\dot{\theta_{s}}=\omega_{s}=\pm 1$. It is obvious that $\delta\varepsilon_{o1}$, $\delta\varepsilon_{o2}$, and hence  $\delta\varepsilon_{o}$ are not necessarily the positive definite small  functionals.
Now, one can do a similar procedure  for the additional terms $\varepsilon_{i}$'s ($i=1,2,3$).  If  one inserts   a slightly deformed  SSWS (\ref{so1}) in  $\varepsilon_{i}$ ($i=1,2,3$), it yields
\begin{eqnarray} \label{so4}
&&\varepsilon_{i}(x,t)=\varepsilon_{is}+\delta\varepsilon_{i}=\delta\varepsilon_{i}=B[3(C_{is}+\delta C_{i})({\cal K}_{is}+\delta{\cal K}_{i})^{2}-({\cal K}_{is}+\delta{\cal K}_{i})^{3}]=\nonumber\\&&
B[3(C_{is}+\delta C_{i})(\delta{\cal K}_{i})^{2}-(\delta{\cal K}_{i})^{3}]\approx B[3C_{is}(\delta{\cal K}_{i})^{2}-(\delta{\cal K}_{i})^{3}]\approx B[3C_{is}(\delta{\cal K}_{i})^{2}]>0\nonumber\\&&
\end{eqnarray}
in which $\varepsilon_{is}=0$, ${\cal K}_{is}=0$ and $C_{is}=\omega_{s}^2\varphi_{s}^4$ are related   to the SSWS (\ref{SSWS}). According to Eq.~(\ref{so4}), since $C_{i}>0$ (\ref{cof}),
$\delta\varepsilon_{i}$'s  are   positive definite,  as is generally  expected  from   Eqs.~(\ref{eis}), (\ref{eis2}) and (\ref{eis3}).

According to Eqs.~(\ref{e5})-(\ref{sc23}),  keeping up the terms to  the first  order of small variations for the deformed SSWS (\ref{so1}),  we have
\begin{eqnarray} \label{sdd}
&&\delta{\cal K}_{1}\thickapprox 2\omega_{s}\varphi_{s}^4\delta \dot{\theta},\nonumber\\&&
\delta{\cal K}_{2}\thickapprox \delta{\cal K}_{1}-2\varphi_{s}' (\delta \varphi')-2(3\varphi_{s}^5-2\varphi_{s}^3)\delta \varphi,\nonumber\\&&
\delta{\cal K}_{3}\thickapprox \delta{\cal K}_{2}+2\varphi_{s}^2(\omega_{s}\delta\dot{\varphi}-\varphi_{s}'\delta \theta').
\end{eqnarray}
Since $\delta{\cal K}_{i}$'s are linear in the first  order of small variations  $\delta\varphi$, $\delta\theta$ and their derivatives ($\delta \dot{\theta}$, $\delta \varphi'$, $\delta\dot{\varphi}$ and $\delta \theta'$), thus,  according to Eq.~(\ref{so4}),  $\delta\varepsilon_{i}$'s are positive definite linear functions of the second order of small variations and their derivatives,  which   are all multiplied by $B$.

For any arbitrary small variations  $\delta\varphi$ and $\delta\theta$  above the background of a non-moving  SSWS (\ref{SSWS}), the variation of the total   energy  can be calculated  by the integration of $\delta\varepsilon$  over the whole space:
\begin{equation}\label{Eo}
\delta E=\int_{-\infty}^{+\infty} (\delta \varepsilon) ~dx=\int_{-\infty}^{+\infty} (\delta\varepsilon_{o}+\sum_{i=1}^{3}\delta\varepsilon_{i}) ~dx=\sum_{j=0}^{3}\delta \mathbb{E}_{j}.
\end{equation}
To show that the SSWS (\ref{SSWS}) is energetically  stable, we must  prove that $\delta E$ is always positive for any arbitrary small deformation. In other words, if any arbitrary deformation needs external energies to occur, then the SSWS (\ref{SSWS}) is an energetically stable solitary wave solution. Since $\delta\varepsilon_{1}$, $\delta\varepsilon_{2}$ and $\delta\varepsilon_{3}$ are positive definite  functions, then their  integration over the whole space, i.e. $\delta \mathbb{E}_{1}$, $\delta \mathbb{E}_{2}$ and $\delta \mathbb{E}_{3}$,  always leads  to positive  values. Now, let us to focus  on the $\delta \mathbb{E}_{o}$:
\begin{equation}\label{Eoo}
\delta \mathbb{E}_{o}=\delta \mathbb{E}_{o1}+\delta \mathbb{E}_{o2}=\int_{-\infty}^{+\infty}  \delta\varepsilon_{o1} ~dx+\int_{-\infty}^{+\infty} \delta\varepsilon_{o2}~ dx,
\end{equation}
where, $\delta \mathbb{E}_{o1}$   is the contribution of the first  order  of variations in  $\delta \mathbb{E}_{o}$, which  we will show that it would be zero in general.  For the not-deformed  non-moving SSWS (\ref{SSWS}), according to Eq.~(\ref{eom}),  as its  dominant  dynamical equation, we can use  $\varphi''_{s}=\frac{d^2\varphi_{s}}{dx^2}$ instead of $\frac{1}{2}\frac{dU(\varphi_{s})}{d\varphi_{s}}$  to obtain:
\begin{eqnarray} \label{sodf}
&&\delta\varepsilon_{o1}=
  2\left[ \varphi_{s}'  (\delta \varphi')+\frac{1}{2}\frac{dU(\varphi_{s})}{d\varphi_{s}}(\delta \varphi)\right]=  2\left[ \varphi_{s}'  (\delta \varphi')+(\delta  \varphi)\varphi''_{s}\right]=2\frac{d}{dx}\left(\delta\varphi\frac{d\varphi_{s}}{dx}\right).
\end{eqnarray}
Hence, the integration of $\delta\varepsilon_{o1}$ over the  whole space leads to
\begin{equation}\label{sxc}
\int_{-\infty}^{+\infty}\delta\varepsilon_{o1}~dx=\left.2(\delta\varphi\frac{d\varphi_{s}}{dx})\right|_{+\infty} -\left.2(\delta\varphi\frac{d\varphi_{s}}{dx})\right|_{-\infty}=0.
\end{equation}
Note that, $\delta\varphi$ and $\frac{d\varphi_{s}}{dx}$ are zero at $\pm\infty$.
Therefore, the following result is generally valid:
\begin{equation}\label{of}
\delta E=\int_{-\infty}^{+\infty}(\delta\varepsilon_{e}) ~dx=\int_{-\infty}^{+\infty} (\delta\varepsilon_{o2}+\sum_{i=1}^{3}\delta\varepsilon_{i})~dx,
\end{equation}
where $\delta\varepsilon_{e}=\delta\varepsilon-\delta\varepsilon_{o1}=\delta\varepsilon_{o2}+\sum_{i=1}^{3}\delta\varepsilon_{i}$, and can be called the effective variation of the energy density function.
Now, if one can prove that for all arbitrary small variations, $\delta\varepsilon_{e}$ always remains a positive function, then the integration of that would be  always positive as well, and the energetical stability condition is fulfilled.
Note that, $\sum_{i=1}^{3}\delta\varepsilon_{i}$ is essentially positive definite, but  $\delta\varepsilon_{o2}$ is not necessarily  a positive function.

Since $\frac{1}{2}\frac{d^2U(\varphi_{s})}{d\varphi_{s}^2}=-15\varphi_{s}^4+6\varphi_{s}^2>0$ for $|\varphi_{s}|<\frac{\sqrt{10}}{5}$, hence undoubtedly, $\delta\varepsilon_{o2}=(\delta\dot{\varphi})^2+ (\delta\varphi')^2 +\frac{1}{2}\frac{d^2U(\varphi_{s})}{d\varphi_{s}^2}(\delta \varphi)^2$ itself would be positive for the points $x$  that   $|\varphi_{s}(x)|$ is  less than  $\frac{\sqrt{10}}{5}$, and then $\delta\varepsilon_{e}>0$ for such points.
But,  function $\frac{1}{2}\frac{d^2U(\varphi_{s})}{d\varphi_{s}^2}$, for the points $x$ that    $|\varphi_{s}(x)|>\frac{\sqrt{10}}{5}$,  would be negative, and we cannot be certain  that $\delta\varepsilon_{o2}$ (and then  $\delta\varepsilon_{e}$) is always positive. Nevertheless, if one considers a system with a large value of parameter $B$,   we can be certain  that  $\delta\varepsilon_{e}>0$ for all points.  In fact, $|\delta\varepsilon_{o2}|$ is a function    of the second order of $\delta\varphi$, $\delta\varphi'$ and $\delta\dot{\varphi}$  which does  not contain parameter $B$,   but $\delta\varepsilon_{i}$'s ($i=1,2,3$) are also functions   of the second order of variations $\delta \varphi$, $\delta\theta$ and their derivatives  which are multiplied by $B$. Hence, we are certain   that always   $\sum_{i=1}^{3}\delta\varepsilon_{i}\gg |\delta\varepsilon_{o2}|$ or $\delta\varepsilon_{e}\approx \sum_{i=1}^{3}\delta\varepsilon_{i}>0$, provided that  $B$ is  a large number (approximately $B>10^2$). Accordingly, for the arbitrary variations $\delta\varphi$ and $\delta\theta$, $\delta E$ would be always positive and then we ensure that the SSWS (\ref{SSWS}) is  an energetically stable object, meaning that, its energy would be at a minimum among the  other (close) solutions. It should be noted  that, the   theory is relativistic, hence confirming the energetical stability of the SSWS at rest, is  generalized to all moving versions at any arbitrary speed.

To summarize, according to the previous considerations,   for any small deformation above the background of the SSWS (\ref{SSWS}), we finally have
\begin{equation}\label{bbn}
E=E_{s}+\delta E\approx\frac{\pi}{4}+\sum_{i=1}^{3}\int_{-\infty}^{+\infty}
\delta\varepsilon_{i} ~dx=\frac{\pi}{4}+3B\sum_{i=1}^{3}\int_{-\infty}^{+\infty}
[C_{is}(\delta{\cal K}_{i})^{2}]dx,
\end{equation}
where, $E$ is the total energy of the small deformed SSWS (\ref{so1}), and $E_{s}=\pi/4$ is the rest energy of the SSWS (\ref{SSWS}). It is true that for  the hypothetical  solutions (\ref{so1}), which are  close to the SSWS (\ref{SSWS}), the field variations $\delta\varphi$, $\delta\theta$ and hence  $\delta{\cal K}_{i}$'s  are  small, but the term  $[3B C_{is}(\delta{\cal K}_{i})^{2}]$ is not necessarily  small,  because it contains  the large parameter  $B$. Hence, $\delta E$ is not necessarily small as well (see Fig.~\ref{mmmm}).   For any  close solution (\ref{so1}), with two specific small variations  $\delta\varphi$ and $\delta\theta$, there are three  specific $\delta{\cal K}_{i}$'s that are used to obtain the total energy (\ref{bbn}). Since $\delta E$ is proportional to the integration of   $\sum_{i=1}^{3}C_{is}(\delta{\cal K}_{i})^2$, thus  for any arbitrary close solution (\ref{so1}), there are always continuously  closer solutions with smaller $\delta\varphi$ and $\delta\theta$ and hence   smaller $\sum_{i=1}^{3}C_{is}(\delta{\cal K}_{i})^2$,  which leads to smaller $\delta E$. Therefore,  none of the close solutions (\ref{so1}),  are  energetically stable. Note that, the close solutions (\ref{so1}) are those for which  the approximations (\ref{so3})  and (\ref{so4}) are valid.

\begin{figure}[ht!]
   \centering
   \includegraphics[width=160mm]{mmmm.jpg}
   \caption{Variations  of  the total  energy $E$ versus small $\xi$  and different $B$'s at $t=0$ for the SSWS (\ref{SSWS}).  The Figs a-f are related to different variations (\ref{var1})-(\ref{var6}) respectively. Note that for the case $\xi=0$, in all figures,  the total energy is the same rest energy of the SSWS (\ref{SSWS}), i.e. $E(\xi=0)=E_{o}=\frac{\pi}{4}$. } \label{mmmm}
\end{figure}
\begin{figure}[ht!]
   \centering
   \includegraphics[width=110mm]{2.jpg}
   \caption{The variation  of the total energy $E$ versus time or all continuous  profiles which all together form    Fig.~\ref{pp}. Some of these  profiles are shown in   Fig.~\ref{3}. As  expected,  a horizontal line is obtained for  case $B=0$,  that is, the line $E=\pi/4$. } \label{2}
\end{figure}

Numerically, we should  study the stability of the SSWS (\ref{SSWS})  for some    arbitrary small hypothetical deformations.
For example, six  arbitrary slightly deformed  SSWSs  (\ref{so1}) can be introduced as follows:
\begin{equation}\label{var1}
\varphi(x,t)=\varphi_{s}+\delta \varphi=\frac{\pm 1}{\sqrt{1+ x^2}}+\xi~ \exp{(-x^{2})},\quad\quad \theta(x,t)=\omega_{s}t,
\end{equation}
 \begin{equation}\label{var2}
\varphi(x,t)=\varphi_{s}+\delta \varphi=\frac{\pm 1+\xi}{\sqrt{1+ x^2}},\quad\quad \theta(x,t)=\omega_{s}t,
 \end{equation}
\begin{equation}\label{var3}
\varphi(x,t)=\varphi_{s}+\delta \varphi=\frac{\pm 1}{\sqrt{ 1+\xi+ x^2}},\quad\quad \theta(x,t)=\omega_{s}t,
\end{equation}
 \begin{equation}\label{var4}
\varphi(x,t)=\varphi_{s}+\delta \varphi=\frac{\pm 1}{\sqrt{ 1+(1+\xi) x^2}},\quad\quad \theta(x,t)=\omega_{s}t,
\end{equation}
\begin{equation}\label{var5}
\varphi(x,t)=\frac{\pm 1}{\sqrt{1+ x^2}},\quad\quad \theta(x,t)=\theta_{s}+\delta \theta=\omega_{s}t+\xi~ t,
\end{equation}
\begin{equation}\label{var6}
\varphi(x,t)=\frac{\pm 1}{\sqrt{1+ x^2}},\quad\quad \theta(x,t)=\theta_{s}+\delta \theta=\omega_{s}t+\xi~ \exp{(-x^{2})},
\end{equation}
in which $\xi$ is a small parameter, which    can be an  indication  of  the order  of   deformations (variations) for any  kind  of  the small deformations (\ref{var1})-(\ref{var6}). All of the deformed functions (\ref{var1})-(\ref{var6})  turn to the same free non-deformed SSWS (\ref{SSWS}) for  $\xi=0$.  For all arbitrary deformations (\ref{var1})-(\ref{var6}),   Fig.~\ref{mmmm}(a-f) show  how  larger values of parameter  $B$ lead to more   stability. In other words, the larger  is the value of $B$ the greater would be the   increase in the total energy versus $|\xi|$. Figure~\ref{mmmm}(a-c) show  clearly  why the case $B=0$ leads  to an energetically  unstable SSWS. In other words, for the case $B=0$, in Figure.~\ref{mmmm}(a-c), the rest energy of the SSWS, i.e. $E_{s}=E(\xi=0)=\frac{\pi}{4}$, is not a minimum.   Note that, the case $B=0$ is  the same original RNKG system (\ref{lag}) with the same SSWS (\ref{ss}).

Some   hypothetical  deformations for the SSWS (\ref{SSWS})   can be considered  as those that appear in Figs.~\ref{pp} and \ref{3}. According to Fig.~\ref{pp}, the profile    of the initial SSWS (\ref{ss}) at $t=0$, does not  remarkably change at the
time interval $0 < t < 3$. In fact, it changes so slightly  that is not visible in Fig.~\ref{pp}. If one numerically  calculates
the total energy of the deformed SSWS   at different times
$0 < t < 5$,  Fig.~(\ref{2}) is obtained.  It reaffirms that for arbitrary small deformations above the background of the
SSWS (\ref{SSWS}), for example in the range $0 < t < 3$, the larger values of $B$ lead to more stability, that is, the larger the values  the greater will be the increase in the total energy for any arbitrary small variation above the background of the SSWS (\ref{SSWS}).
In general,    the total energy always increases (and  increases more for larger values of $B$) versus the amount of any arbitrary variation above the background of the SSWS (\ref{SSWS}). Although the parameter $B$ can be taken a  large value,  it would not affect the dominant dynamical equation (\ref{eom}) and the observable of the SSWS (\ref{SSWS}).

At the initial times (the times that are close to  $t=0$), a  multi lump (particle-like) solution with different velocities can be easily constructed just by adding distinct  far apart SSWSs (\ref{SSWS2}) together. For example, for $N$ distinct   SSWSs (\ref{SSWS2}) at different velocities $v_{j}$ and different initial positions $a_{j}$,  provided that $a_{j+1}-a_{j}\gg 1$, the multi particle-like solution at initial times  is as follows:
 \begin{equation}\label{SSWS3}
\varphi(x,t)=\sum_{j=1}^{N}\left[\frac{\pm 1}{\sqrt{1+\gamma_{j}^2(x-v_{j}t-a_{j})^2}}\right],
\end{equation}
where $\gamma_{j}=1/\sqrt{1-v_{j}^2}$. Since the phase field $\theta$ for each SSWS (\ref{SSWS2}) depends on its velocity, hence it must change from one to another.   That is to say,   if there are two SSWSs (i.e. $N=2$) with  one of them being at rest ($v_{1}=0$)  and the other  moving ($v_{2}\neq 0$), then the phase field must change from $\theta=\omega_{s}t$ at the position of the first SSWS  to  $\theta=k_{\mu}x^{\mu}=\omega_{2} t-k_{2}x$ at the position of the second SSWS. In the regions between two SSWSs, the scalar field $\varphi$ is almost zero and hence   $\varepsilon$ is almost zero everywhere.  Thus, there is not any rigorous restriction on $\theta$ to be in  the standard forms  $\theta=\omega_{s}t$ and $\theta=k_{\mu}x^{\mu}$ as  the special solutions of the PDE $\mathbb{S}_{1}=0$. In other words,  where the scalar  field   $\varphi$ is almost zero, the phase field $\theta$ is completely free and evolve without any rigorous restriction, i.e. it can change  slowly from $\theta=\omega_{s}t$ to $\theta=k_{\mu}x^{\mu}$  in  the spaces  for which  $\varphi\thickapprox 0$. In fact, for the case $\varphi\thickapprox 0$, it is not necessary    to satisfy   $\mathbb{S}_{i}=0$ or ${\cal K}_{i}=0$ ($i=1,2,3$) simultaneously, because for such situations, all $\varepsilon_{i}$'s ($i=1,2,3$) would automatically be almost zero simultaneously  without any restrictive condition.

The  general dynamical equations (\ref{geq}) and (\ref{geq2}), as two coupled nonlinear PDEs, have infinite solutions. Similar to the  SSWS (\ref{SSWS}), some of these solutions may be stable or there may not have any other stable solution at all. However,
if one considers  a system with a large  parameter $B$, as we have already shown, there is not any energetically stable solution among  the close solutions (\ref{so1}) at all.
But how is it possible  to know if the system has other stable solutions or not? In the present situation, it is by no means our goal to answer this question that can be mathematically important. In general, if one demonstrates that the other possible  stable solutions have very large energies, then physically, it is not an important issue, because they require  very large external energies  to be created and  it is outside the scope of our research. In fact,   except for  the  solutions which are very close\footnote{In order for a better understanding of the matter, we can provide a qualitative definition of  a conventional criterion.  We could designate  “\textit{very close solutions}" for those “\textit{close solutions}" (\ref{so1}) whose energy, for example for the case $B=10^{20}$, is less than $10^{3}$, or the ones for which  the magnitude of the variations $\delta\varphi$ and $\delta\theta$ is approximately less than $10^{-9}$.}  to the free far apart SSWSs (\ref{SSWS3}) and trivial vacuum state $\varphi=0$, the other (possible stable) solutions  can not be physically simply created. To put it differently, for any (possible stable) solution, it is not possible for all ${\cal K}_{i}$'s to be zero simultaneously, thus  at least one of  the $\varepsilon_{i}$'s, which contain the  large  parameter  $B$,  is a non-zero large function and then  the   energy is much  larger than the rest energy of a SSWS (\ref{SSWS}). Accordingly, the other (possible stable) solutions  need large  external  energies to be created, which is very unlikely to occur  physically. For example, if $B=10^{20}$, for a hypothetical small deformation (\ref{var2}) with  $\xi=\pm 10^{-5}$, the total energy is in the order of $10^{11}$, that is, the required  external energy must be of the order  $10^{11}$ to create  a small deformed SSWS (\ref{var2}) with  $\xi=\pm 10^{-5}$!

Furthermore, for other (possible stable) solutions, the additional  term $F$ is no longer zero  and is a large functional.  Thus, in the coupled PDEs (\ref{geq}) and (\ref{geq2}), the terms $\Box \varphi$ and $\frac{1}{2}\frac{dU}{d\varphi}$ are very small  compared with  the terms which contain $F$. In other words, the general   dynamical equations (\ref{geq}) and (\ref{geq2}), for the other solutions, which are not very close to the free far apart SSWSs (\ref{SSWS3}) and the vacuum state $\varphi=0$, are reduced to
\begin{eqnarray} \label{dfg}
&& \frac{\partial}{\partial x^{\mu}}\left(\frac{\partial F }{\partial (\partial_{\mu}\varphi)}\right)-\left(\frac{\partial F}{\partial \varphi}\right)  =0,\\ \label{dfg2} &&   \frac{\partial}{\partial x^{\mu}}\left(\frac{\partial F }{\partial (\partial_{\mu}\theta)}\right)=0,
\end{eqnarray}
as their dominant dynamical equations. In terms of functional ${\cal K}_{i}$'s, since $F=B\sum_{i=1}^{3}{\cal K}_{i}^{3}$, equivalently Eqs.~(\ref{dfg}) and (\ref{dfg2}),  turn into
\begin{eqnarray} \label{dlo}
&&\sum_{i=1}^{3} \left[2{\cal K}_{i}(\partial_{\mu}{\cal K}_{i})   \frac{\partial{\cal K}_{i}}{\partial(\partial_{\mu}\varphi)}    +   {\cal K}_{i}^2\partial_{\mu}\left(\frac{\partial{\cal K}_{i}}{\partial(\partial_{\mu}\varphi)}\right)    -    {\cal K}_{i}^2\frac{\partial{\cal K}_{i}}{\partial \varphi}   \right]=0,~~~~~~~\\&&\label{dlo1}
\sum_{i=1}^{3} \left[2{\cal K}_{i}(\partial_{\mu}{\cal K}_{i})   \frac{\partial{\cal K}_{i}}{\partial(\partial_{\mu}\theta)}    +   {\cal K}_{i}^2\partial_{\mu}\left(\frac{\partial{\cal K}_{i}}{\partial(\partial_{\mu}\theta)}\right)       \right]=0.
\end{eqnarray}
The important  point is that the parameter $B$ is again  ineffective in these equations. In other words, any solution of the coupled PDEs (\ref{dlo}) and (\ref{dlo1}) would be approximately  a solution of the general dynamical equations (\ref{geq}) and (\ref{geq2}) as well, provided that $B$ is considered to be a large number.

For the  far apart free SSWSs (\ref{SSWS3}), the dominant dynamical equation is the same simple equation (\ref{eom}),  and the role of the catalyzer term (\ref{FFF}) is approximately zero. However, when they get close to each other and
their profiles change slightly, the role of the stability catalyzer term  $F$ becomes important and strongly opposes a closer approach and causes  more change in their  profiles.
In this situation, the dominant dynamical equations are the same general forms  (\ref{geq}) and (\ref{geq2}).
Accordingly, as  expected  in the collision between the SSWSs (\ref{SSWS3}), they reappear  after collision processes without any significant  distortion.
For example,  for two SSWSs (\ref{SSWS3}) which  are initialized to collide with each other at the same speed $v=1-10^{-11}$, that is very close to the speed of light, the  kinetic energy of each one  is $K=E-E_{s}\thickapprox 2.2\times 10^{5}$. Now, if $B=10^{20}$, for a hypothetical deformation like (\ref{var2}), such  kinetic energy can only cause a  variation in the order of  $\xi\approx10^{-8}$ for each SSWS, that is, they reappear after the collision without any significant  distortion.  In general, for two free SSWSs (\ref{SSWS3}), while they are far apart, the total energy is finite, but when they come close together to interact, their profiles can   not have  notable  deformations, because their  initial kinetic energy must be a  huge  value  to  generate a  remarkable deformation for each SSWS, which is not  physically simple to be provided.

Based on all that has been said so far,  it is obvious that  if one considers  a system with an extremely   large value of $B$ (for example $B=10^{20}$ or even more),  the other  configurations of the fields $\varphi$ and $\theta$,   which are not very  close to the free far apart SSWSs (\ref{SSWS3}) and the vacuum state $\varphi=0$, require  extreme  energy to be created.  From a  physical  point of view, this issue can be interesting, as it classically explains how a system  leads to many identical particles with  specific characteristics. In other words, if one considers this system as a real physical system,   it is not possible to provide  an extremely large   external energy at a special place  for creating the other (possible stable) configurations of the fields. Thus, the only non-trivial stable configurations  of the fields  with the finite energies would be  any number of the  far apart SSWSs (\ref{SSWS3}) as a multi particle-like solution.
Similar to the  quantum field theory, the free far apart SSWSs (\ref{SSWS3})   can be  classically called the  quanta of the  system. Again, it should be noted  that, it does not matter to us whether the system has other possible  stable solutions.
What is important to us here is that the only non-trivial stable solutions with  finite energies are any number of the free far apart SSWSs (\ref{SSWS3}), as many identical particles with the specific characteristics,   which can be interesting for physicists.

\section{Summary and conclusion}\label{sec6}

In this paper, we introduced an extended  Klein-Gordon  system as an example (\ref{lag2}), which analytically yields an energetically stable solitary wave solution (\ref{SSWS}). In other words, it leads to a soliton solution. The new Lagrangian density (\ref{lag2}) is composed of two distinct parts, first, the original  part which is a known standard RNKG system (\ref{lag}), and second, an additional part, which can be called  the stability catalyzer term (\ref{FFF}). The original standard RNKG  Lagrangian density (\ref{lag}) is introduced for a single scalar field $\varphi$. But,  to introduce a proper   stability catalyzer term, it is necessary to use a different scalar field  $\theta$ (phase field) along with the original scalar field  $\varphi$, meaning that the stability catalyzer term is a functional of  $\varphi$ and $\theta$ simultaneously. The role of the  stability catalyzer term  seems as a massless spook   which surrounds the SSWS (\ref{SSWS}) and   guarantees the  stability of the SSWS (\ref{SSWS}). The stability catalyzer term  has no role in    the dominant   dynamical equation  and the other properties  of the SSWS (\ref{SSWS}), meaning that the general  dynamical equations (\ref{geq}) are reduced to the same   known standard RNKG version (\ref{eom}) just for the SSWS (\ref{SSWS}).  However, it guarantees   the energetical stability of the SSWS (\ref{SSWS}). Therefore,
any arbitrary small deformation above the background of  the SSWS (\ref{SSWS}) leads to an increase in the total energy. In other words, the rest energy of the SSWS is at a  minimum among the other solutions of the new extended KG system (\ref{lag2}) except for the ones which are very close to the vacuum state $\varphi=0$.

There is a parameter $B$ in the stability catalyzer term, the larger values of which, leads to more stability; meaning that, the larger the values the greater will be the increase in the total energy for any arbitrary small variation above the background of the SSWS (\ref{SSWS}). Hence,  considering a system with an extremely large value of $B$, leads to a classical  system  only with multi particle-like solutions. In fact,  the other   solutions of the system, which are not very close to the free far apart  SSWSs (\ref{SSWS3})  and the vacuum state $\varphi=0$, require infinite   amounts of energy to be created, and  are not physically possible to occur.  Thus, physically the possible stable solutions of the system with the finite energies, are either  any number of the far apart SSWSs (as any number of identical free particles) or the trivial  vacuum state. In other words, the SSWS (\ref{SSWS}) can be considered as the quantum of this classical system.

\section*{Acknowledgement}

The author  wishes to express his appreciation to the Persian Gulf University Research Council for their constant support.