\section{Introduction}
The knowledge of a person's spatial whereabouts is an important component of a ubiquitous system. In the past years, several researchers have proposed various methods for automatic collection of a user's spatial trajectories using smartphone applications~\cite{Chon14,Kjaergaard11,Stylianou15} running in the background. Given appropriate algorithms for processing this spatial data in order to extract useful information such as stay-points and user patterns, important applications emerge. Examples of such applications include predicting a user's movement~\cite{Herder14}, tourism~\cite{Horanont14}, extracting individual life patterns~\cite{Ye09} and  recommending locations and activities~\cite{Zheng10a}.

In this work, we are interested in the user's stay-points. A \textit{stay-point} is a location where the user has stopped for some amount of time. This could be the user's home or work where the duration is large or could be the user's children's school where a daily transportation is done. Even though the stay-point identification problem has been previously addressed by several researchers, the problem was mainly tackled as a point clustering problem where usually existing clustering techniques were used. An important aspect of previous work, was the identification of the maximum distance between two points such that they are considered to be in the same location and the minimum time between them such that they are labeled as stay-points~\cite{Ashbrook03, Kang05, Palma08, Zheng10, Horanont14, Ranvier15}. In contrast, recent work~\cite{Zhou14, Yadav14, Thomason15} attempts to reduce the dependency on these thresholds but introduces other thresholds.  Another problem is that confusion is generated to the reader as there is no agreement on the values of distance and time thresholds because their values vary significantly among different works. In addition, despite the fact that the exhibited performance of previous work was quite good, demonstrating stay-point identification rates of over $90\%$, this performance is mainly achieved after experimental fine tuning of the different threshold values such that they fit very well their datasets.

\begin{figure}[t]
  \centering\begin{tabular}{c}
     \begin{minipage}{6cm}
     \centering
     \includegraphics[width=6cm]{19-2-16.png}
    \end{minipage}
\\
     \begin{minipage}{6cm}
     \centering
     \includegraphics[width=6cm]{19-2-16-spatial-curve-stay-points.png}
    \end{minipage}

\end{tabular}
\caption{Method overview. We transform the user's trajectory for a single day to a 2D spatial curve. Using extrema extraction, we identify the user's stay-points; these are depicted using rectangles.}
\label{fig:methodOverview}
\end{figure}
It is evident that an appropriate method that reduces significantly the dependence on experimentally deduced threshold values or a method that requires little or no fine tuning does not exist. In contrast, we propose a different perspective towards stay-point identification as depicted in figure~\ref{fig:methodOverview}. Instead of directly processing the recorded raw spatial trajectories, we first transform each daily spatial trajectory to a two-dimensional ($2D$) discrete spatial curve. Now the problem of stay-point identification is transformed from a thresholding problem to an extrema extraction problem. Since, the spatial curve models displacement over time and it's first derivative models speed over time, one can intuitively see it as unifying distance and time thresholds into a single threshold. However this new threshold required to identify an extremum has a well-posed numerical value.
Among others, we show that the proposed method produces good results on both sparse and dense spatial trajectory data, does not depend in any way on the uniformity of the data points and it is stable over large positioning errors that raw trajectory data exhibit.

This article has two main contributions:
\begin{itemize}
\item We propose a novel method for stay-point identification. This method tackles the problem from a completely different perspective. In this same time it eliminates dependency on threshold values that are experimentally deduced and are good only on specific datasets.
\item The proposed method is very stable over different trajectory densities, works well on uniformly and non-uniformly spaced raw data points.
\end{itemize}

\section{Related Work}
\label{sec:relatedwork}
Past work on stay-point extraction can be divided into two groups. The work that depends on time and distance thresholds~\cite{Ashbrook03,Hariharan04, Kang05, Palma08, Zheng09, Zheng10,  Nishida14, Horanont14, Ranvier15, Pavan15} and the work that attempts to eliminate this dependency~\cite{Zhou14, Yadav14, Thomason15}.

One of the first methods was developed by Ashbrook and Starner~\cite{Ashbrook03} that proposed that stay-points are points that have stay time of at least $10$ minutes and used the k-means clustering algorithm in order to discover the optimal distance that reveals stay-points. The next year, Hariharan et al.~\cite{Hariharan04} used the same thresholds but they proposed that their values should be tuned according to the application. Nevertheless, they show results for distance and time threshold of $50$ meters and $10$ minutes, respectively. Kang et al.~\cite{Kang05} suggest that a good value for the distance threshold is $30-50m$ and for time is $6$ minutes for long stays and proposed a time threshold of $100-150$ seconds for quick stays. Even though Palma et al.~\cite{Palma08} still use a distance threshold they try to compute it by computing the quantile function - the inverse of the cumulative distribution function - of the gaussian curve for the arithmetic mean and standard deviation of the list of distances between consecutive points. The quantile function and a time threshold value, proposed to be set to two minutes, are used to improve the DBSCAN clustering algorithm for the purpose of stay-point identification. Later, Zheng et al.~\cite{Zheng10,Zheng09} proposed to identify stay-points using a time threshold of $20$ minutes and a distance threshold of $200$ meters. Also, Nishida et al.~\cite{Nishida14} propose a $20$ meter distance threshold and $5$ minute time threshold. Horahont et al.~\cite{Horanont14} consider the distance threshold of $196$ meters and the time threshold of $14$ minutes to be good as identified out of experiments with $15$ people. Finally, Pavan et al.~\cite{Pavan15} decide that for stay-point labeling the distance threshold is $50$ meters and the time threshold is $50$ seconds. They also introduce a speed threshold defined as the median of all speeds reached within the trajectory data.

In contrast, other researchers~\cite{Zhou14, Yadav14, Thomason15} proposed different approaches towards eliminating this dependency. Zhou et al.~\cite{Zhou14} in order to detect a stay-point introduce another threshold denoted as $\mathbf{K}$. This is used in a time stamp clustering algorithm to pre-process the trajectory data such that a classical clustering algorithm such as k-means and DBSCAN produce better results. In addition, the authors experiment with a time threshold.

Yadav et al.~\cite{Yadav14} record Cell Ids and use a graph based clustering algorithm in order to identify stay points. The graph is built such that the cell ids are the vertices and two cell ids are incident if their time difference is less than $\alpha$; this ensures successive cell ids are connected. Then they cluster vertices, within the same cell id, using their edge weight and vertex degree which are controlled by two arbitrary parameters $\eta$ and $\eta'$, respectively. Each cluster within a Cell represents a stay-point.

Thomason et al.~\cite{Thomason15} propose a gradient-based visit extractor algorithm. This algorithm works as follows: insert all points into a buffer until the user has moved more than some threshold (computed via the gradient) or the time difference between the first point and the last point in the buffer is greater than another time threshold; then if there is some duration between the first and last points in the buffer, then a stay-point is identified. The threshold for the gradient, that controls the distance the user moved, depends on two parameters, which are experimentally tuned per dataset.

It is obvious that there is no consensus on the number of thresholds and their the best values. In the previous work, the values for time threshold range from $50$ seconds to $20$ minutes for long stays and the values for distance threshold range from $20$ meters to $200$ meters. Even though these numbers came out of experimental results with groups of people, these results are confusing when studying, evaluating or even attempting to implement a previous method. In addition, we believe that the stay-point identification problem can be more effectively tackled as a curve extremum problem than a clustering problem for reasons including reduction of the two main thresholds of time and distance to a single more intuitive threshold and the simplification of implementation.
Furthermore, we show in our experimental results (see table \ref{table:inflectionPointDurations}) that these hard threshold values lead to false identification of stay-points.

\section{Proposed Method}
\label{sec:proposedMethod}
We propose to transform the problem from a spatial to a geometry problem such that the stay-points become curve extrema. In this section, we provide the problem definition, describe the problem transformation and explain how we identify the stay-points.

\subsection{The Problem}
\label{subsec:theProblem}

\textbf{Objective}. Given a spatial trajectory (of a person) for a single day, the objective is to identify all the stay-points of this person.

\textbf{Input}. The spatial trajectory is represented as the set $S=\{loc_i=(t_i,\mathbf{l}_i)|\,i=1...n\}$, where  $\mathbf{l}_i=(latitude,longitude)$ are geographical coordinates and $t_i$ (timestamp) is the date/time that $\mathbf{l}_i$ is recorded.

\textbf{Output}.  The stay-points, denoted as $\mathbf{s}_i$, which are the locations the user has stayed for some amount of time.

\subsection{Problem Transformation}
\label{subsec:problemTransformation}
As already mentioned, we tackle the problem of stay-point identification under a totally different perspective with respect to the previous work. The first step is the transformation of the spatial trajectory to a spatial curve as it is depicted in figure~\ref{fig:methodOverview}.

The \textit{spatial curve} is a $2D$ increasing trajectory curve $C$ composed of points  $\mathbf{p}_i=(x_i,y_i)$, where $x_i$ is the elapsed time since midnight (i.e. the beginning of the day) in minutes and $y_i$ is the total displacement, in kilometers, until minute $x_i$.

The transformation is achieved via a bijective map $M$ (eq.~\ref{eq:map}), which is composed of the functions $f$ (eq.~\ref{eq:f},~\ref{eq:f_explained}) and $g$ (eq.~\ref{eq:g},~\ref{eq:g_explained}). The function $f$ (eq.~\ref{eq:f_explained}) transforms the timestamp $t_i$ to the $x_i$-coordinate of a point and the function $g$ (eq.~\ref{eq:g_explained}) transforms the position $\mathbf{l}_i$ to the $y_i$-coordinate of a point.

\begin{gather}
\label{eq:map}
M: (t_i,\mathbf{l}_i)\rightarrow (x_i,y_i) \\
\label{eq:f}
f:t_i\rightarrow x_i\\
\label{eq:g}
g:\mathbf{l}_i\rightarrow y_i \\
\label{eq:f_explained}
\left. f(t_i)=t_i-midnight\right\vert_{\mbox{in minutes}}\\
\label{eq:g_explained}
g(l_i)=\sum_{i=2}^{n} distance(l_i,l_{i-1})\mbox{, }g(l_1)=0
\end{gather}

The function $g$ transforms the geographical coordinate $\mathbf{l}_i$ to the cumulative distance from the first geographical coordinate $\mathbf{l}_1$.

\begin{figure}[h]
  \centering\begin{tabular}{c}
     \begin{minipage}{7cm}
     \centering
     \includegraphics[width=7cm]{19-2-16-stay-points-inflection-points.png}
    \end{minipage}
\\
(a)\\
     \begin{minipage}{7cm}
     \centering
     \includegraphics[width=7cm]{19-2-16-der1-stay-points-inflection-points.png}
   \end{minipage}
  \\
 (b) \\
     \begin{minipage}{7cm}
     \centering
     \includegraphics[width=7cm]{19-2-16-der2-stay-points-inflection-points.png}
    \end{minipage}
\\
(c)
\end{tabular}
\caption{Spatial curve overlayed with inflection points (diamonds) and stay-points (rectangles). (a) spatial curve, (b) $1^{st}$ derivative spatial curve, (c) $2^{nd}$ derivative spatial curve.}
\label{fig:spatialCurveWithDerivatives}
\end{figure}
\subsection{Stay-point Identification}
\label{subsec:stay-pointIdentification}
Given the resulting $2D$ spatial curve, we are in place to identify the daily stay-points. We achieve this task by performing extrema extraction on this curve.

\textbf{Extrema extraction}. Let's see what a stay-point looks like on a spatial curve. We remind that for every point $(x_i,y_i)$ of the spatial curve, $x_i$ is the minute, elapsed from midnight, that the cumulative distance $y_i$ is recorded. It is also necessary to remind that the spatial curve is an increasing curve. We claim that \textit{\textbf{the stay-points are the local minima of the $1^{st}$ derivative of the spatial curve}}; the rationale follows.

The first derivative of the spatial curve $C'(x,y)=\Delta y/\Delta x$ provides the speed of movement. So, when the user stops the speed is zero or  ideally $C'(x,y)=0$. An observation we make is that the $1^{st}$ derivative values are always greater than or equal to zero (fig.~\ref{fig:spatialCurveWithDerivatives}b) because the spatial curve is an increasing one (fig.~\ref{fig:spatialCurveWithDerivatives}a). Consequently it is easy to see that all stay-points are local minima of the $1^{st}$ derivative curve. These local minima can be identified using the $2^{nd}$ derivative (fig.~\ref{fig:spatialCurveWithDerivatives}c) check, which is $C''(x,y)=\Delta^2 y/\Delta x^2=0$. However, in practice the second derivative would never be zero. Instead, a numerically stable widely used approach is to identify these local minima by checking for zero-crossings of the $2^{nd}$ derivative~\cite{Haralick84}. The  $2^{nd}$ derivative will cross zero, from positive to negative, indicating deceleration until the stop. When the user starts moving again, the  $2^{nd}$ derivative will cross zero again, indicating acceleration due to the movement.

A drawback of the $2^{nd}$ derivative zero-crossing check is that even though this check is adequate for identifying extrema of the spatial curve, it also identifies inflection points of the curve. \textit{Inflection points} are stationary points where the sign of curve's curvature changes (or the curve's concavity changes). Given what the spatial curve represents, inflection points may indicate short stops, e.g. due to traffic jams.

Generally it is easy to decide between an extremum or an inflection point by just checking the concavity of the curve or by thresholding on the $1^{st}$ derivative value because the $1^{st}$ derivative is in the vicinity of zero ($C'(x,y)\approx 0$) only for extrema points and conversely it is far larger than zero for inflection points. However, in this application where positioning errors are abundant and the spatial curve is not smooth, these approaches cannot be applied without some error filtering. As smoothing the curve is costly and using hard thresholding iterates the problem of having to identify a good threshold, we decide to use soft thresholding by using a confidence value, as explained later.

\textbf{Stay-point region.} It is evident that the zero-crossings of the $2^{nd}$ derivative return a pair of points ($\mathbf{p}_i,\mathbf{p}_j)$, where $i<=j$. Even though we may get a pair of points we could also get a single point, if $i=j$. In addition, we could get more points (when $i-1<j$), i.e. a stay-point region $[i, j]$. This is completely normal as it depends on the number of locations returned by the trajectory tracking technique when the user is not moving. Considering the scenario that the user has slowed down significantly, let's say due to traffic, closeby the destination, before completely stopping, it is possible that this region contains outlier points. Similarly, it is possible that the whole region corresponds to inflection points.

\textbf{Extremum or Inflection?} Given a pair of points ($\mathbf{p}_i,\mathbf{p}_j)$ that corresponds to consecutive zero-crossings, the aim is to compute a confidence value that provides insight into whether the corresponding region $[i,j]$ is an extremum region.

First we compute the confidence value $C(\mathbf{p}_w)$ for each point $\mathbf{p}_w$, $i \le w \le j$ according to equation \ref{eq:confidence}, where $\mathbf{p}^{'}_{w}$ is the $1^{st}$ derivative on point $\mathbf{p}_w$, $\mathbf{p}^{'}_{min}$ is the minimum $1^{st}$ derivative over all points of the curve. Equation~\ref{eq:confidence} sets the confidence to $100\%$ when the first derivative is less than $e$. When the first derivative's value is increasing the confidence decreases until it becomes zero. Our experiments showed that a good value for $e$ is $0.05$ as this doesn't produce false minima. However, $e$ could be even smaller. A smaller $e$ can be used to reduce error due to noise, if such exists, by making it stricter to identify a stay-point.

The second step is to compute the confidence $C(\mathbf{p}_i,\mathbf{p}_j)$ for the whole region (see eq.~\ref{eq:totalConfidence}). The confidence value for the region $[i, j]$ is the confidence of the subregion $[k, w]$ (see eq.~\ref{eq:partialConfidence}) with maximum confidence value.

\begin{equation}
\begin{split}
C(\mathbf{p}_w) &=100\mbox{, if }\mathbf{p}_w^{'}<=e\\
                           &=100-\frac{\mathbf{p}^{'}_w-e}{\mathbf{p}^{'}_{min}}\mbox{, if }\mathbf{p}^{'}_w>e\\
                           &=0\mbox{, if }\frac{\mathbf{p}^{'}_w-e}{\mathbf{p}^{'}_{min}}>100
\end{split}
\label{eq:confidence}
\end{equation}

\begin{equation}
C(\mathbf{p}_k,\mathbf{p}_w)=\frac{\sum_{i=k}^w C(\mathbf{p}_i)}{w-k+1}
\label{eq:partialConfidence}
\end{equation}

\begin{equation}
C(\mathbf{p}_i,\mathbf{p}_j)=\max_{i\le k\le w\le j} C(\mathbf{p}_k,\mathbf{p}_w)
\label{eq:totalConfidence}
\end{equation}

After the computation of the confidence $C(\mathbf{p}_i,\mathbf{p}_j)$ we classify the region $[k, w]$ as a stay-point region when $C(\mathbf{p}_i,\mathbf{p}_j)\ge 80$. When the confidence falls in the range $[60,80)$ the region is considered a candidate stay-point region and we discard regions with lower confidence values. A candidate stay-point region may become a stay-point if the user identifies it as such or if an analysis that aims to identify user stay-point patterns indicates that it is a stay-point. Such an analysis is out of the scope of the current work.

\subsection{Computational Complexity}
The proposed method requires the computation of the $1^{st}$ and $2^{nd}$ derivatives and the computation of a maximum confidence value for each candidate region. We compute the derivatives using central differences with $2^{nd}$ order accuracy. Hence the derivatives can be computed concurrently. The confidence value requires $p^{'}_{min}$, which is the minimum $1^{st}$ derivative over all points of the curve. Similarly this can be obtained during the computation of the $1^{st}$ derivative. In addition, the maximum confidence value for a region can be computed in linear time as computationally this would mean to compute the confidence value for the region and then trim from left and right all points such that confidence increases. Therefore the computational complexity of the proposed method is $O(N)$, where $N$ is the number of points of the spatial curve.

In addition, replacing $\mathbf{p}^{'}_{min}$ with the constant $0.001$, which numerically is a good one, relaxes the need to have all $1^{st}$ derivatives prior to computing $\mathbf{p}^{'}_{min}$. The advantage is that the proposed method becomes real-time as it is able to process trajectory points as they arrive and in this way identify stay-points as they appear.

\section{Experimental Evaluation}
\label{sec:evaluation}
The objective of the evaluation is to examine if the proposed work performs in line to the theoretical analysis and identify factors that may affect performance.

To this end, we have collected daily spatial trajectories for a period of $28$ days. This data was collected using different spatial tracking techniques. We executed the proposed method on the data produced by these techniques and collected all automatically generated stay-points with their confidence values. In order to make a proper evaluation, we compared them with the ground truth, i.e. the actual stay-points.

It is important to mention that we don't perform an evaluation on publicly available datasets such as GeoLife~\cite{Zheng10} because the ground truth is not provided. As a result it is not possible to properly evaluate the performance of our method.

\subsection{Data Acquisition}
\label{subsec:dataAcquisition}
The iOS operating system via the coreLocation framework provides two services for spatial tracking: (a) the standard location service and (b) the significant location service. Service (a) relies on the GPS sensor, generates uniformly spaced data with various densities. Service (b) relies on the cell towers, generates non-uniformly spaced, sparse data, as it provides a new location at least every $500m$ or $5$ minutes.

We implemented three applications that use the standard location service in order to acquire the trajectory path using three different densities: $100$, $250$ and $500$ meters. We denote these as $SLS-100$, $SLS-250$ and $SLS-500$. In addition, we implemented one application that uses a combined approach~\cite{Stylianou15}, which is denoted as \textit{hybrid}. This approach uses primarily the significant location service for spatial tracking. Whenever a new location is acquired it schedules the standard location service to power up in one minute in order to capture extra locations; this ensures that if a stay-point occurred the correct position of this point is captured. Examples of trajectories generated by two different techniques for the same day are provided in figure~\ref{fig:example}. We quickly observe the differences in data density and uniformity.

The trajectory points were acquired using these four different iOS applications running on an iphone 6S and an iphone 5. Both devices run the iOS $10.x$ operating system. The iPhone 6S was used to run the $SLS-500$ and the hybrid and the iPhone 5 was used to run the $SLS-100$ and $SLS-250$. In total, we acquired trajectories for $60$ days and $28$ unique days. We acquired trajectories for $8$ days using all four techniques, $8$ days using two techniques ($SLS-500$, hybrid) and $12$ days using only the hybrid approach.

It is important to mention that the acquired data do not go through any noise removal process for the purposes of increasing accuracy rates. This means that the trajectory data exhibit noise from $10m$, which is normal, to about $1000$ meters, which is extreme. The underlying reason is that we want to demonstrate, via the experiments, robustness to noise.

\begin{figure}[h]
  \centering\begin{tabular}{c}
     \begin{minipage}{8cm}
     \centering
     \includegraphics[width=7cm]{gps.png}
    \end{minipage}
\\
     \begin{minipage}{8cm}
     \centering
     \includegraphics[width=7cm]{SLUandGPS.png}
   \end{minipage}
\end{tabular}
\caption{Trajectory tracking examples for a single day. Top: SLS-500. Bottom: Hybrid.}
\label{fig:example}
\end{figure}

\subsection{Ground Truth}
\label{subsection:groundtruth}
During the evaluation the subject manually logged in the end of each day its stay-points including approximate arrival and departure times.  Figure~\ref{fig:typical_day} shows a state diagram for the subject's stay-points for a single day. The subject was moving in urban and suburban areas by car and was concurrently using the two iPhones to automatically track its daily whereabouts.

The ground truth stay-points are utilized in two ways: (a) we manually check whether these were captured in the raw trajectory data and (b) we compare them to the automatically identified stay-points in order to evaluate the accuracy of the proposed method.

\begin{figure}[h]
  \centering\begin{tabular}{c}
     \begin{minipage}{7cm}
     \centering
     \includegraphics[width=7cm]{typical_wednesday.png}
   \end{minipage}
\end{tabular}
\caption{One typical Wednesday.}
\label{fig:typical_day}
\end{figure}

\subsection{Results}
\label{subsection:comparison}
In order to evaluate the accuracy of the proposed method, we compared the stay-points extracted from the raw trajectory data, acquired using four different techniques, to the ground truth stay-points. The comparison aims to reveal the accuracy for automatic stay-point identification, the accuracy for inflection point identification, the percentage of false positive stay-points. In addition, we compared the automatically approximated stay-point duration to the actual duration and also examined whether a time threshold for identifying stay-points could exhibit similar accuracy to the proposed method.

\begin{table}[htb]
\centering
\caption{Summary of results. It includes the number of trajectory days, average number of trajectory points per day, the success rate and the percentage of false positives.}
\begin{tabular}{|l | p{1cm}| p{1.5cm}| p{1.5cm}| p{1.5cm}|}
\hline
         & \# of stay-points & Points per day & Success rate (\%) & False positive (\%) \\ \hline
Hybrid   & 180         & 48             & 86                & 2.8                 \\ \hline
SLS-100 & 56        & 97             & 98                & 13                  \\ \hline
SLS-250 & 56        & 78             & 84                & 9                   \\ \hline
SLS-500 & 108         & 54             & 60                & 5                   \\ \hline
\end{tabular}
\label{table:summaryResults}
\end{table}

Table~\ref{table:summaryResults} shows the overall success rate of stay-point identification for all data acquisition techniques. The best accuracy is achieved by the $SLS-100$ technique with about $98\%$ success rate, and the second best accuracy by the hybrid approach with $86\%$ success rate.  In addition, we observe that our proposed method exhibits a very few number of false positives -these are falsely identified stay-points-, which can be as low as $2.8\%$.
Since the hybrid method produces the least number of data points, we believe that the reduced data density reduces the number of false positives and that the stay-point identification accuracy is not directly related to the density of the data points.

We also observe that given adequate data points the results could reach up to nearly $100\%$ accuracy rates. However, we notice that the increased number of data points increases the possibility for false positives as these increase up to $13\%$. The results show that the false positives and the missed stay-points are not repeating over the same days (e.g. all Tuesdays). Therefore, a good method that infers the user's repeating patterns will not include most of false positives and will increase the success rate for stay-points due to inference of the missing ones. After studying the results, we conclude that false positives occurred mainly due to urban traffic and very rarely due to other reasons. The increasing number of data points produced by the techniques $SLS-500$, $SLS-250$ and $SLS-100$, captured urban traffic with increasing accuracy. In one such case, the author's (driver) car was stationary for more $15$ minutes. Unfortunately, such cases are quite difficult to identify automatically without incorporating a learning method into the process.

It is as important to investigate the reasons for failing to identify some stay-points. According to table~\ref{table:missingStaypoints}, several stay-points were not identified by the hybrid method because (a) there were no trajectory points captured in their vicinity ($3.3\%$), (b) there were very few trajectory points in the vicinity of stay-points ($5.6\%$), (c) the stay was very short ($2.2\%$) of about a minute or so or (d) stay-points were identified but with lower than $80\%$ confidence ($1.7\%$). If we assume non-existence of reasons due to trajectory tracking techniques then we can deduce that the success rate of our proposed method is at least $90\%$, for the hybrid method. This hypothesis could theoretically increase the success rate for $SLS-250$ and $SLS-500$ to about $85\%$. Given the fact that the proposed method executed on raw trajectory data, we believe that the identification rate is quite satisfactory and can increase due to pre-processing. Pre-processing can help to reduce the effect of erroneous data as the trajectory data, in some cases, exhibited positioning errors from $5m$, which is normal, to more than $1000m$, which is abnormal, both in urban and suburban areas.

In addition, the results summarized in both tables~\ref{table:summaryResults} and~\ref{table:missingStaypoints} show that the proposed formula for computing the confidence value has very good results. First, only very few stay-points are not identified due to low confidence values; most of these are candidate stay-points. Second, the number of false positives due to high confidence can be as low as $2.8\%$. We note that the main reason for increased false positives e.g. in $SLS-100$ was not the confidence value, as explained before. Third, the proposed method assigned correct low confidence value to a total of $204$ inflection points, acquired by all four techniques, which are $50\%$ of the total number of stay-points.

\begin{table}[htb]
\centering
\caption{The reasons some stay-points were not identified.}
\label{table:missingStaypoints}
\begin{tabular}{|l | p{1cm}| p{1.5cm}| p{1.5cm}| p{1.5cm}|}
\hline
                               & Short stay (\%)           & Few points (\%)           & No points (\%)            & Low confidence (\%)       \\ \hline
Hybrid  & 2.2 & 5.6&3.3 & 1.7 \\ \hline
SLS-250 & 1.8 & 0   & 7.3&1 \\ \hline
SLS-500                       & 0.9                      & 3.7                      & 21.3                     & 0                        \\ \hline
\end{tabular}
\end{table}

Further, we have evaluated the accuracy of the automatically computed duration by the two best performing trajectory tracking techniques for all the common stay-points and days by comparing them to the actual duration. Given the stay-point region $[k,w]$, the duration of stay was estimated by considering the time to distance from $\mathbf{p}_{k-1}$ to $\mathbf{p}_k$ and from $\mathbf{p}_{w}$ to $\mathbf{p}_{w+1}$ in addition to the time between $\mathbf{p}_k$, $\mathbf{p}_w$. The time to distance was estimated assuming traveling speed of $50Km/h$. The actual duration ranges from $5$ minutes to several hours. Table~\ref{table:durationComparison} summarizes the comparison. Both trajectory tracking techniques over-estimated the duration by about $19\%$ and $13\%$ on average. These deviations are on average about $12.5$ minutes for both techniques. However, the standard deviation shows that the deviations from the average for the $SLS-100$ technique are smaller than the hybrid technique. This makes sense because $SLS-100$ produces twice the number of data points to the hybrid technique.

\begin{table}[htb]
\centering
\caption{Comparisons of the estimated duration for 43 stay-points using trajectory data from two different techniques. The actual total duration is 4194 minutes.}
\label{table:durationComparison}
\begin{tabular}{|l|p{1.5cm}|p{1cm}|p{1cm}|p{1cm}|p{1cm}|}
\hline
 &       & \multicolumn{4}{c|}{Duration (deviation from actual)}                       \\ \cline{2-6}
                  & Total Minutes & Average (\%) & STD (\%) & Average (min) & STD (min) \\ \hline
Hybrid            & 4517  & 19           & 30        & 12.6          & 20.2      \\ \hline
SLS-100           & 4294  & 13           & 17        & 12.5          & 16.7      \\ \hline
\end{tabular}
\end{table}

Finally given the estimated duration, we examined whether a simple threshold could be adequate as most of the previous work suggests. For instance, one could claim that since, according to our experiments, the minimum actual duration is $5$ minutes then it would be adequate to just use a time threshold of $5$ minutes in order to identify a stay-point. To this end, we have again examined the duration of all the inflection points for the best performing techniques. We remind that an inflection point is an extremum point (or region) that is assigned a confidence value less than $60\%$.

Our findings, which are summarized in table~\ref{table:inflectionPointDurations}, reveal that the minimum, maximum and average duration for the hybrid technique are $4$,  $49$ and $9$ minutes, respectively. In addition, the minimum, maximum and average duration for the $SLS-100$ technique are $2$, $24$ and $4$ minutes, respectively. The aforementioned results demonstrate that there is a significant gap between the minimum and maximum duration. Probably the most interesting outcome is that a time threshold of $5$ minutes~\cite{Nishida14} would label $30$ and $16$ inflection points as stay-points for the hybrid and $SLS-100$ techniques, respectively. Consequently, the accuracy rates would have been significantly smaller as false positives will increase $6$ times and $2.5$ times for the hybrid and $SLS-100$ techniques, respectively.

\begin{table}[htb]
\centering
\caption{The duration of inflection points for the Hybrid and SLS-100.}
\label{table:inflectionPointDurations}
\begin{tabular}{|l|l|l|l|l|l|}
\hline
 & \#of Points & \multicolumn{4}{c|}{Duration} \\ \cline{3-6}
                  &                              & Min     & Max    & Average   & $\ge 5'$ \\ \hline
Hybrid            & 32                           & 4       & 49     & 9   &   30  \\ \hline
SLS-100           & 54                           & 2       & 24     & 4   &  16   \\ \hline
\end{tabular}
\end{table}

\section{Conclusion and Future Work}
We have presented a novel method that identifies stay-points using a geometry-based approach. In a nutshell, our method transforms the initial trajectory path to a time series curve such that the stay-point identification problem becomes an extrema extraction problem. The proposed method avoids using experimentally deduced distance, time or other threshold values. In contrast, it uses a confidence formula in order to robustly distinguish stay-points from stationary points. It is also runs in linear time and can run in real time.

We evaluated our method on raw trajectory data, for $28$ unique days, acquired using four different techniques with varying densities.
The results demonstrate that the proposed method has $85-98\%$ accuracy and it is insensitive to the trajectory point uniformity, density and noise. The confidence value, computed for each stay-point region, correctly identifies all stationary points and generates very few false positive stay-points.

In the future, we aim to further improve the proposed method's accuracy by investigating various ways for either improving the quality of the raw trajectory data or by improving current trajectory tracking techniques. A way to do either is to incorporate a machine learning method. In addition, we plan to extend our method to $3D$ such that we can identify user patterns using $3D$ feature extraction.

\bibliographystyle{plain}


\section{Introduction}

Quantum systems are generally never completely isolated but interact with an environment with which they may exchange energy and other physical observable quantities. Their properties are naturally affected by their coupling to the external world. The understanding and control of the influence of an environment on a given physical system is of crucial importance in different fields of physics and in technological applications and led to a large amount of investigations ~\cite{reb,aba,kos,gol,cor}.
Among the quantities of interest energy exchange between a system and its environment is of prime importance. The total system composed of the considered system and its environment is closed. It conserves the physical observables, in particular the total energy it contains. However, the existence of an interaction between its two parts leads to possible exchanges between them, energy and other observables can be transfered between them in both directions. In the present work we study the conditions under which this tranfer occurs.

In order to investigate the process a number of different approaches have been developed
\cite{bel,bag,fl1,fl2,esp}. They make use of the cumulant method as well as other developments ~\cite{car,sch,def}. In the following we shall use the cumulant method in order to work out the energy exchange and the speed at which this exchange takes place~\cite{gua}.

In previous studies ~\cite{kr1,kr2,kr3} we defined criteria which allow to classify open quantum systems with respect to their behaviour in the presence of an environment. In order to do so we used a general formulation which relies on the examination of the properties of the density operator of the system and its environment. The dynamical behaviour of this operator is governed by the structure of the Hamiltonians of the system, its environment and the coupling Hamiltonian which acts between them. The method introduces a general form of the total density operator and avoids the determination of explicit solutions of the equation of motion which governs the time evolution of the system. We follow this formalism in order to work out different cases of physical interest concerning the energy exchange process.\\

The analysis is presented in the following way. In section 2 we define the energy exchange and its rate starting from the characteristic function which generates these quantities in terms of its first and second moment. Section 3 introduces a general formal expression of the density operator at the initial time and the structure of the total Hamiltonian of the system, the environment and their interaction. The expressions of the energy exchange and the exchange rate are explicitly written out. In section 4 we analyze and work out two different cases. They exemplify the role of the interaction  between the system and the environment which may or may not induce the time divisibility property of the open system. Conclusions are drawn in section 5.

\section{Energy transfer between S and E: the cumulant approach}

Consider a system $S$ coupled to an environment $E$. The time dependent density operator of the total system is $\hat \rho_{SE}(t)$ and $S$ and $E$ are coupled by an interaction $\hat H_{SE}$. The interaction may generate an energy exchange between the two parts. This quantity can be worked out by means of the cumulant method ~\cite{gua} which was developed in a series of works, first for closed driven systems ~\cite{esp}, later extended to open systems, see f.i.~\cite{fl1,fl2} and further approaches quoted in ~\cite{gua}.

Define the modified density operator

\ba
\rho^{\sim}_{SE}(t,0)=Tr_{E}[{\hat U_{\eta/2}(t,0)\hat \rho_{SE}(0)\hat U^{+}_{-\eta/2}(t,0)}]
\label{eq1}
\end{eqnarray}
where $Tr_{E}$ is the trace over the states in $E$ space and
\ba
\hat U_{\eta}(t,0)=\exp(+i\eta \hat H_{SE})\hat U(t,0)\exp(-i\eta \hat H_{SE})
\label{eq2}
\end{eqnarray}
Here $\hat U(t,0)$ is the time evolution operator of the interacting total system $S+E$.

The characteristic function obtained from the generating function reads
\ba
\chi^{(\eta)}(t)=Tr_{S}(\rho^{\sim}_{SE}(t,0))
\label{eq3}
\end{eqnarray}

From $\chi^{(\eta)}(t)$ one derives the energy exchange between the environment and the system

\ba
\Delta E(t)=\frac{d \chi^{(\eta)}(t)}{d (i\eta)}_{| \eta=0}
\label{eq4}
\end{eqnarray}
The speed at which the energy flows between the system and its environment is given by

\ba
V_{E}(t)=\frac{\partial \dot \chi^{(\eta)}(t)}{\partial (i\eta)}_{|\eta=0}
\label{eq5}
\end{eqnarray}
where $\dot \chi^{(\eta)}(t)$ is the time derivative of $\chi^{\eta}(t)$. If the energy flows from the sytem to the environment $V_{E}(t)$ is positive and negative if the flow is reversed.

\section{Energy flow and speed of energy exchange}

\subsection{The density operator}

At time $t=0$ we choose the system to stay in a mixed state

\ba
|\psi(0)\rangle=\sum_{{i_k}}c_{i_k}|i_k\rangle
\label{eq6}
\end{eqnarray}
where the normalized states $\{|i_{k}\rangle\}$ are eigenstates of the Hamiltonian $\hat H_{S}$. The environment is described in terms of its density matrix chosen as
$\{ |\alpha\rangle d_{\alpha,\alpha}\langle \alpha|\}$ where $\{d_{\alpha,\alpha}\}$ are the statistical weights of the density matrix in a diagonal basis of states $\{\alpha\}$.

Given these bases of states in $S$ space and in $E$ space the density operator of the total system at
time $t=0$ is written as~\cite{buz}

\ba
\hat \rho_{SE}(0)=\hat \rho_{S}(0) \otimes \hat \rho_{E}(0)
\notag\\
\hat \rho_{S}(0)=\sum_{k,l}|i_{k}\rangle c_{i_{k}}c^{*}_{i_{l}}\langle i_{l}|
\notag\\
\hat \rho_{E}(0)=\sum_{\alpha}|\alpha \rangle d_{\alpha,\alpha} \langle \alpha|
\label{eq7}
\end{eqnarray}
The density operator $\hat \rho_{SE}$ describes a system in the absence of an interaction
$\hat H_{SE}$ at $t=0$, hence in the absence of an initial entanglement between $S$ and $E$.

The total Hamiltonian $\hat H$ reads

\ba
\hat H=\hat H_{S}+\hat H_{E}+\hat H_{SE}
\label{eq8}
\end{eqnarray}

and
\ba
\hat H_{S}|i_{k}\rangle=\epsilon_{i_{k}}|i_{k}\rangle
\notag\\
\hat H_{E}|\gamma\rangle=E_{\gamma}|\gamma\rangle
\label{eq9}
\end{eqnarray}
where $\{|i_{k}\rangle\}$ and $\{|\gamma\rangle\}$ are the eigenvector bases of $\hat H _{S}$ and $\hat H_{E}$. At time $t$ the evolution of $S+E$ is given by

\ba
\hat\rho_{S+E}(t)= \hat U(t)\hat \rho_{S+E}(0) \hat U^{+}(t)
\label{eq10}
\end{eqnarray}
where $\hat U(t,0)=e^{-i\hat H t}$ is the evolution operator of $S+E$.

\subsection{Explicit expressions of the energy exchange and the speed of the flow}

Using the definitions given in Eqs. (4) and (5) the energy transfer and speed of the energy flow read

\ba
\Delta E(t)=\frac{1}{2}Tr_{S}Tr_{E}\{[\hat H_{E},\hat U(t,0)]\hat\rho_{S+E}(0)\hat U^{+}(t,0)\}+h.c.
\label{eq11}
\end{eqnarray}
and

\ba
V_{E}(t)=\frac{1}{2}Tr_{S}Tr_{E}\{[\hat H_{E},\frac{d}{dt}\hat U(t,0)]\hat\rho_{S+E}(0)\hat U^{+}(t,0)\}+h.c.
\label{eq12}
\end{eqnarray}

Developing these expressions in the bases of states given above they read

\ba
\Delta E(t)=\frac{1}{2}\sum_{j,\gamma}\sum_{i_{1},i_{2},\gamma_{1}}c_{i_{1}} c{*}_{i_{2}}
\langle j \gamma|[\hat H_{E},e^{-i\hat H t}]|i_{1} \gamma_{1}\rangle
\notag\\
d_{\gamma_{1},\gamma_{1}}\langle i_{2} \gamma_{1}|e^{+i\hat H t}| j \gamma \rangle +h.c.
\label{eq13}
\end{eqnarray}
and

\ba
V_{E}(t)=\frac{(-i)}{2} \sum_{j,\gamma}  \sum_{\alpha_{1},i_{1},i_{2}}
c_{i_{1}} c^{*}_{i_{2}} d_{\alpha_{1},\alpha_{1}}
\notag\\
\sum_{j_{1}, \gamma_{1}}\{\langle j \gamma|[\hat H_{E},\hat H_{SE}]|j_{1}, \gamma_{1}\rangle
\langle j_{1}, \gamma_{1}|e^{-i\hat H t}|i_{1}\alpha_{1}\rangle
\langle i_{2}\alpha_{1}||e^{+i\hat H t}|j \gamma \rangle
\notag\\
+ \langle j \gamma|\hat H|j_{1} \gamma_{1}\rangle
\langle j_{1} \gamma_{1}|[\hat H_{E},e^{-i\hat H t}]|i_{1} \alpha_{1}\rangle
\langle i_{2} \alpha_{1}|e^{+i\hat H t}|j \gamma \rangle
\notag\\
- \langle j \gamma|\hat H e^{-i\hat H t}|i_{1} \alpha_{1}\rangle
\langle i_{2}\alpha_{1}||[\hat H_{E},  e^{+i\hat H t}]|j \gamma \rangle \}
+h.c.
\label{eq14}
\end{eqnarray}

\section{Properties of the interaction Hamiltonian}

We use now the general expressions of the energy transfer and its exchange rate in order to test the role of the interaction Hamiltonian in this process. Since $S$ and $E$ are distinct different physical systems their Hamiltonians verify the commutation relation $[\hat H_{S},\hat H_{E}]=0$. Former work ~\cite{kr1} has shown that one may consider two cases which are of special interest:\\

(a) $[\hat H_{E},\hat H_{SE}]=0$\\

(b) $[\hat H_{S},\hat H_{SE}]=0$

\subsection{Case (a)}

It has been shown elsewhere ~\cite{kr1,kr2} that if $H_{E}$ and $H_{SE}$ commute the evolution of the system $S$ is characterized by the divisibility property which is a specific property of Markovian systems.
Since the Hamiltonian of the system $S$  commutes in practice with $\hat H_{E}$ it follows that $\hat H_{E}$ commutes with the whole Hamiltonian $\hat H$, hence also $\hat U(t,0)$ and its derivatives. Going back to the expressions of $\Delta E(t)$ and $V_{E}(t)$ it comes out that $\Delta E(t)=V_{E}(t)=0$.

There is no energy exchange between $S$ and $E$ in this case. The physical explanation is the following: the divisibility property imposes that the environment stays in a fixed state at a fixed energy which blocks any possible transfer of energy between the two parts of the total system $S+E$. This is due to the fact that $\hat H_{SE}$ is diagonal in $E$ space, hence the considered state $|\gamma \rangle$ stays the same over any time interval.

\subsection{Evolution of the energy: an impurity immersed in a bosonic condensate}

We introduce the Hamiltonian of a fermionic impurity interacting with a Bose-Einstein condensate ~\cite{chr}. It is given by $\hat H= \hat H_{S}+\hat H_{E}+\hat H_{SE}$ where

\begin{center}
\ba
\hat H_{S}=\sum_{\vec k}\epsilon_{\vec k}c^{+}_{\vec k}c_{\vec k}
\label{eq15}
\end{eqnarray}
\end{center}
\begin{center}
\ba
\hat H_{E}=\sum_{\vec k}e_{\vec k}a^{+}_{\vec k}a_{\vec k}+\frac{1}{2V}\sum_{\vec k_{1}\vec k_{2}
\vec q}V_{B}(\vec q)(\vec q)a^{+}_{\vec k_{1}}a^{+}_{\vec k_{1}+\vec q}a^{+}_{\vec k_{2}-\vec q}
a_{\vec k_{2}}a_{\vec k_{1}}
\label{eq16}
\end{eqnarray}
\end{center}
\begin{center}
\ba
\hat H_{SE}=\frac{1}{V}\sum_{\vec k_{3}\vec k_{4}\vec q}c^{+}_{\vec k_{3}+\vec q}c_{\vec k_{4}}
a^{+}_{\vec k_{4}-\vec q}a_{\vec k_{3}}
\label{eq17}
\end{eqnarray}
\end{center}
where $[c,c^{+}]$ and $[a,a^{+}]$ are fermion and boson annihilation and creation operators.

We consider the case where the momentum transfer $\vec q=0$. Then one expects that there is no energy exchange between $S$ and $E$. This is indeed so since a simple calculation shows that
$[\hat H_{E},\hat H_{SE}]=0$ in this case.

\subsection{Case (b)}

Start from the general expression of $\Delta E(t)$ given by Eq.(13).

Since $[\hat H_{S},\hat H_{SE}]=0$ all the matrix elements are diagonal in $S$ space and the expression takes the form

\ba
\Delta E(t)=\sum_{j}|c_{j}|^{2}\sum_{\gamma,\gamma_{1}}(E_{\gamma}-E_{\gamma_{1}})
\langle j \gamma|e^{-it \hat H}|j \gamma_{1}\rangle d_{\gamma_{1}\gamma_{1}}
\langle j \gamma_{1}|e^{+it \hat H}|j \gamma \rangle
\label{eq18}
\end{eqnarray}

In order to follow the evolution of $\Delta E(t)$ in time we determine the expression of $V_{E}(t)$.
A somewhat lengthy but straightforward calculation leads to the following expression:

\ba
V_{E}(t)=V^{(1)}_{E}(t)+ c.c. +V^{(2)}_{E}(t)+ c.c.
\label{eq19}
\end{eqnarray}
where
\ba
V^{(1)}_{E}(t)=(-i)\sum_{j}|c_{j}|^{2}\sum_{\gamma,\gamma_{1}}(E_{\gamma}-E_{\gamma_{1}})
\langle j \gamma|\hat H_{SE}|j \gamma_{1}\rangle
\notag\\
\sum_{\gamma_{2}}\langle j \gamma_{1}|e^{-it \hat H}|j \gamma_{2}\rangle
d_{\gamma_{2},\gamma_{2}} \langle j \gamma_{2}|e^{+it \hat H}|j \gamma \rangle
\label{eq20}
\end{eqnarray}
and $V^{(1)*}_{E}(t)$ its complex conjugate. It comes out that the c.c. $V^{(1)*}_{E}(t)=V^{(1)}_{E}(t)$ which means that $V^{(1)}_{E}(t)$ is real.
The second term reads

\ba
V_{E}^{(2)}(t)= (-i)\sum_{j}|c_{j}|^{2}\sum_{\gamma \gamma_{1}}(E_{\gamma}+\epsilon_{j})
\sum_{\gamma_{2}}(E_{\gamma}-E_{\gamma_{2}})
\notag\\
\langle j \gamma_{1}|e^{-it \hat H}|j \gamma_{2}\rangle
d_{\gamma_{2},\gamma_{2}}\langle j \gamma_{2}|e^{+it \hat H}|j \gamma \rangle
\label{eq21}
\end{eqnarray}

It is easy to see that $V_{E}^{(2)}(t) +c.c.=0$, hence $V_{E}(t)=2V_{E}^{(1)}(t)$.

For $t=0$

\ba
\Delta E(0)= \sum_{j}|c_{j}|^{2}\sum_{\gamma,\gamma_{1}}(E_{\gamma}-E_{\gamma_{1}})
\langle j \gamma|j \gamma_{1}\rangle d_{\gamma_{1}\gamma_{1}}
\langle j \gamma_{1}|j \gamma \rangle=
\notag\\
\sum_{j}|c_{j}|^{2}\sum_{\gamma,\gamma_{1}}(E_{\gamma}-E_{\gamma_{1}})\delta_{\gamma \gamma_{1}}
\label{eq22}
\end{eqnarray}
Hence $\Delta E(0)=0$ which could have been anticipated from the symmetry property of the expression of $\Delta E(t)$.\\

But contrary to case $(a)$ the energy tranfer is now different from zero. It varies with time hence $\Delta E(t)$ increases or de decreases as a function of the sign of $V_{E}(t)$. This is due to the fact that now many channels can open in $E$ space because $\hat H_{SE}$is no longer diagonal in this space.

\subsection{Evolution of the energy: two examples}

In order to illustrate this case we develop two models on which we exemplify the time dependence of the energy transfer between a system and its environment when $[\hat H_{S},\hat H_{SE}]=0$.

\begin{itemize}

\item {\bf First example}: we consider a model consisting of a 2-level state system $E$,
$[|\gamma\rangle=|1\rangle,|2\rangle]$. The Hamiltonian $\hat H$ decomposes into two parts
$\hat H_{0}=\hat H_{S}+\hat H_{SE}$ and $\hat H_{E}$. We consider the case where
\ba
[\hat H_{0},[\hat H_{0},\hat H_{E}]]=[\hat H_{E},[\hat H_{0}, \hat H_{E}]]=0
\label{eq23}
\end{eqnarray}
and $[\hat H_{0},\hat H_{E}]=c1$ where $c$ is a number. Then

\ba
e^{\hat H_{0}+\hat H_{E}}=e^{\hat H_{0}}e^{\hat H_{E}}e^{c/2}
\label{eq24}
\end{eqnarray}
and

\ba
e^{i(\hat H_{0}+\hat H_{E)}t}=e^{i\hat H_{0}t}e^{i\hat H_{E}t}e^{ct^{2}/2}
\label{eq25}
\end{eqnarray}

The quantities which enter the expressions which follow are defined in Appendix A.

The expressions of $\Delta E(t)$ and $V_{E}(t)$ read

\ba
\Delta E(t)=e^{ct^{2}}\Delta_{12}\sum_{j}|c_{j}|^{2}(d_{22}-d_{11})[a_{j}^{(12)2}(t)+b_{j}^{(12)2}(t)]
\label{eq26}
\end{eqnarray}
where $\Delta_{12}=E_{1}-E_{2}$, $d_{11},d_{22}$ the weights of the states in $E$ space and

\ba
V_{E}(t)=2e^{ct^{2}} \Delta_{12}\sum_{j} |c_{j}|^{2}[I^{(j)}_{12}Re(\langle 2|\hat \Omega_{j}t)|1\rangle + R^{(j)}_{12}
Im (\langle 2|\hat \Omega_{j}(t)|1\rangle]
\label{eq27}
\end{eqnarray}
with

\ba
Re\langle 2|\hat \Omega_{j}(t)|1\rangle=a_{j}^{11}(t)[a_{j}^{21}(t)\cos(\Delta_{12}t)+b_{j}^{21}(t)
\sin(\Delta_{12}t)]d_{11}
\notag\\
+a_{j}^{22}(t)a^{21}(t)d_{22}
\notag\\
Im\langle 2|\hat \Omega_{j}(t)|1\rangle=a_{j}{11}(t)
[-b_{j}^{21}(t)\cos(\Delta_{12}t)+a_{j}^{21}(t)\sin(\Delta_{12}t)]d_{11}
\notag\\
+b_{j}^{21}(t)a_{j}^{22}(t)d_{22}
\label{eq28}
\end{eqnarray}

Both $\Delta E(t)$ and $V_{E}(t$ are oscillating functions of time. $\Delta E(t)$ keeps a fixed sign depending on the sign of $\Delta_{12}$, $V_{E}(t)$ may change sign with time. The energy and the speed of the energy transfer decays to zero for real $c$ real and negative.\\

\item
{\bf Second example}: we consider the Hamiltonian $\hat H=\hat H_{S}+\hat H_{E}+\hat H_{SE}$ which governs the coupling of a phonon field to the electron in the BCS theory of superconductivity.

The total Hamiltonian of the electron-phonon system reads

\begin{center}
\ba
\hat H_{S}=\sum_{\vec k_{1}}\epsilon_{\vec k_{1}}c^{+}_{\vec k_{1}}c_{\vec k_{1}}
\label{eq29}
\end{eqnarray}
\end{center}
\begin{center}
\ba
\hat H_{E}=\sum_{\vec q}\hbar \omega_{\vec q}a^{+}_{\vec q}a_{\vec q}
\label{eq30}
\end{eqnarray}
\end{center}
\begin{center}
\ba
\hat H_{SE}=V_{ph-e}(\vec q)\sum_{\vec k_{2}\vec q}(a^{+}_{-\vec q}+a_{\vec q})
c^{+}_{\vec k_{2}+\vec q}c_{\vec k_{2}}
\label{eq31}
\end{eqnarray}
\end{center}
where $V_{ph-e}(\vec q)$ is the phonon-electron interaction, $[a,a^{+}]$ and $[c,c^{+}]$ are phonon and electron annihilation and creation operators. We consider the case where the phonons evolve in the zero mode $\vec q=0$.

Then $\hat H_{E}=\hbar \omega_{0}a^{+}_{0}a_{0}$ and
$\hat H_{SE}=V_{ph-e}(0)\sum_{\vec k_{2}}(a^{+}_{0}+a_{0})c^{+}_{\vec k_{2}}c_{\vec k_{2}}$.

The density operator of the total system $S+E$ at $t=0$ is given by the expression
\ba
\hat \rho^{(n_{0},n'_{0})}_{(\vec k,\vec k')}=\frac{1}{2\pi(n_{0}!n'_{0}!)^{1/2}}|\vec k n_{0}\rangle    \langle \vec k' n'_{0}|
\label{eq32}
\end{eqnarray}
Working out the commutation relation between $\hat H_{S}$ and $\hat H_{SE}$ leads to
$[\hat H_{S},\hat H_{SE}]=0$, hence the evolution operator can be written as
\ba
e^{-i\hat Ht}=e^{-i\hat H_{S}t}e^{-i(\hat H_{E}+\hat H_{SE})t}
\label{eq33}
\end{eqnarray}
and in explicit form

\ba
e^{-i\hat Ht}=e^{-it\sum_{k=1}^{\nu}\epsilon_{\vec k}c^{+}_{\vec k}c_{\vec k}}
e^{-i[\omega_{0} a^{+}_{0}a_{0}+\nu V_{ph-e}(0)(a^{+}_{0}+a_{0})]t}
\label{eq34}
\end{eqnarray}
where $\nu$ is the number of electron states. The quantity of interest concerns the time dependence of $\Delta E(t)$ and $V_{E}(t)$ given by Eqs.(18) and (19), hence the time evolution of the matrix elements of $e^{-i\hat Ht}$

\ba
E^{(n_{0},n'_{0})}_{\vec k \vec k}(t)=\langle \vec k n_{0}|e^{-it\sum_{k=1}^{\nu}\epsilon_{\vec k}c^{+}_{\vec k}c_{\vec k}}e^{-i[\omega_{0} a^{+}_{0}a_{0}+\nu V_{ph-e}(0)(a^{+}_{0}+a_{0})]t}|\vec k n'_{0}\rangle
\label{eq35}
\end{eqnarray}
These matrix elements and their hermitic conjugates which enter the expressions of $\Delta E(t)$ and $V_{E}(t)$ can be worked out explicitly using the Zassenhaus development~\cite{za}, see Appendix B and~\cite{ca}.

They read
\ba
E^{(n_{0},n'_{0})}_{\vec k \vec k}(t)=N^{-1}e^{-i\epsilon_{k}t}e^{-i\omega_{0}n_{0}t}  F^{(n_{0},n'_{0})}(t)
\label{36}
\end{eqnarray}
with
\ba
F^{(n_{0},n'_{0})}(t)= \sum_{n_{2}\leq n_{0},n_{2}\leq n_{3}}
\sum_{n_{4}\leq n_{3},n_{4}\leq n'_{0}}(-i)^{n_{0}+n_{3}}(-1)^{n'_{0}+n_{2}-n_{4}}
\notag\\
\frac{n_{0}!n'_{0}!(n_{2}!)^{2}(n_{3}!)^{2}[\alpha(t)^{n_{0}+n_{3}-2n_{2}}]
[\zeta(t)^{n'_{0}+n_{3}-2n_{4}}]}{(n_{2})^{2} (n_{4})^{2}(n_{0}-n_{2})!(n_{3}-n_{4})!
(n_{3}-n_{2})!(n'_{0}-n_{4})!}e^{\Psi(t)}
\label{37}
\end{eqnarray}
and the normalization factor $N= 2\pi(n_{0}!n'_{0}!)^{1/2}$. The functions $\alpha(t)$, $\zeta(t)$
and $\Psi(t)$ read

\ba
\alpha(t)=\frac{\nu V_{ph-e}(0)\sin\omega_{0}t}{\omega_{0}}
\label{eq38}
\end{eqnarray}

\ba
\zeta(t)=\frac{\omega_{0}[1-\cos\nu V_{ph-e}(0)t]}{\nu V_{ph-e}(0)}
\label{eq39}
\end{eqnarray}

\ba
\Psi(t)=-\frac{1}{2}[\frac{\nu^{2}V^{2}_{ph-e}(0)\sin^{2}(\omega_{0} t)}{\omega_{0}^{2}}+\frac{\omega_{0}^{2}(1-\cos\nu V_{ph-e}(0)t)^{2}}{\nu^{2}V_{ph-e}^{2}}]
\label{eq40}
\end{eqnarray}

As one can see from these expressions $E^{(n_{0},n'_{0})}_{\vec k \vec k}(t)$ are oscillating functions of time which leads to the conclusion that the energy transfer and the transfer velocity oscillate continuously and stay finite over any interval of time.

\end{itemize}

\section{Conclusions, remarks}

In the present work we used a cumulant approach~\cite{gua} in order to study the energy transfer between a system and its environment and the rate at which this transfer evolves in time.

If the Hamiltonian of the environment commutes with the interaction between the system and the environment
there results an absence of energy transfer between the two subsystems. This can be explained in the following way. As already seen in former work the commutation property is a sufficient condition for divisibility in the time behaviour of an open system, one of the properties which characterize Markov processes ~\cite{sti,riv}.
In this case it has been shown~\cite{kr1,kr2,kr3} that the environment keeps in the energy state in which it was at the origin of time and stays there over any interval of time. Hence the energy of the environment is blocked in a fixed state so that it cannot feed the system and does not receive energy from it. Time delays are correlated with the possibility of the environment to jump between different states in closed systems as well as in open ones (see f.i.~\cite{dod,def1,def2} and references quoted in there). This is the case  in non-Markovian systems ~\cite{zhe,hai,ren,gua1}.

The experimental realization of the absence of energy transfer may be obtained under different conditions:

\begin{itemize}

\item the strength of the interaction can be chosen such that it keeps very weak and hence does not allow any possible jump to another level in the case of a discrete environment spectrum.

\item the temperature of the environment is kept close to zero so that the ground state is the only accessible state.

\item the commutation relation between the environment and the interaction is rigorously verified which is the case discussed in the present work.

\end{itemize}

We considered also the case where the system and the interaction Hamiltonians commute with each other.
In this case energy can flow from the environment to the system and back, there is no blocking effect
coming from the environment since several states are at hand and the energy exchange can take place.
In order to illustrate this situation we worked out two model systems corresponding to different physical situations. In the first case the energy exchange oscillates but stops exponentially, in the second case which models the electron-phonon system the energy exchange and the exchange speed oscillate in time and never go to zero. This is so because the system behaves coherently as it has been shown in ref.~\cite{kr3,lid}.

\section{Appendix A}

We define the following matrix elements which enter the expressions of $\Delta E(t)$ and $V_{E}(t)$ given by Eqs.(16-18)

\ba
\langle j 1|e^{-i\hat Ht}|j 1\rangle=e^{ct^{2}/2}e^{-i\epsilon_{j}t}e^{-iE_{1}t}a_{j}^{11}(t)
\notag\\
\langle j 2|e^{+i\hat Ht}|j 2\rangle=e^{ct^{2}/2}e^{+i\epsilon_{j}t}e^{+iE_{2}t}a_{j}^{22}(t)
\notag\\
\langle j 1|e^{+i\hat Ht}|j 2\rangle=e^{ct^{2}/2}e^{+i\epsilon_{j}t}e^{+iE_{1}t}
(a_{j}^{12}(t)+ib_{j}^{12}(t))
\notag\\
\langle j 2|e^{-i\hat Ht}|j 1\rangle=e^{ct^{2}/2}e^{-i\epsilon_{j}t}e^{-iE_{2}t}
(a_{j}^{21}(t)-ib_{j}^{21}(t))
\label{eq41}
\end{eqnarray}
where $\epsilon_{j}$ is the eigenvalue of state $|j\rangle$ and $|E_{k}\rangle,(k=1,2)$ the eigenvalues of states the states $|\gamma\rangle$
\ba
a_{j}^{11}(t)=Re(\langle j 1|e^{-i \hat H_{SE}t}|j 1 \rangle)
\notag\\
a_{j}^{22}(t)=Re(\langle j 2|e^{+i \hat H_{SE}t}|j 2 \rangle)
\notag\\
a_{j}^{12}(t)=Re(\langle j 1|e^{+i \hat H_{SE}t}|j 2 \rangle)
\notag\\
a_{j}^{21}(t)=Re(\langle j 2|e^{-i \hat H_{SE}t}|j 1 \rangle)
\notag\\
b_{j}^{12}(t)=Im(\langle j 1|e^{+i \hat H_{SE}t}|j 2 \rangle)
\notag\\
b_{j}^{21}(t)=Im(\langle j 2|e^{-i \hat H_{SE}t}|j 1 \rangle)
\label{eq42}
\end{eqnarray}
Introduce also

\ba
R^{(j)}_{12}=Re(\langle j 1|\hat H_{SE}|j 2 \rangle)
\notag\\
I^{(j)}_{12}=Im(\langle j 1|\hat H_{SE}|j 2 \rangle)
\label{eq43}
\end{eqnarray}
and

\ba
\hat \Omega_{j}(t)=\sum_{\gamma}e^{-i\hat H t}|j \gamma\rangle d_{\gamma \gamma}\langle j \gamma|
e^{+i\hat H t}
\label{eq44}
\end{eqnarray}

Using these defintions the calculation leads to expressions (26) and (27) of $\Delta E(t)$ and $V_{E}(t)$ in the text.

\section{Appendix B: the Zassenhaus development}

If
$X=-i(t-t_{0})(\hat H_{S}+\hat H_{E})$ and $Y=-i(t-t_{0})\hat H_{SE}$

\ba
e^{X+Y}=e^{X}\otimes e^{Y}\otimes e^{-c_{2}(X,Y)/2!}\otimes e^{-c_{3}(X,Y)/3!}\otimes e^{-c_{4}(X,Y)/4!}...
\label{eq45}
\end{eqnarray}
where

\begin{center}
$c_{2}(X,Y)=[X,Y]$\\
$c_{3}(X,Y)=2[[X,Y],Y]+[[X,Y],X]$\\
$c_{4}(X,Y)=c_{3}(X,Y)+3[[[X,Y],Y],Y]+[[[X,Y],X],Y]+[[X,Y],[X,Y]$\\
\end{center}

The series has an infinite number of term which can be generated iteratively in a straightforward way ~\cite{ca}. If $[X,Y]=0$ the truncation at the  third term  leads to the factorisation of the $X$ and the $Y$ contribution. If $[X,Y]=c$ where $c$ is a c-number the expression corresponds to the well-known Baker-Campbell-Hausdorff formula.