\section{Introduction}
Process mining is a fast growing discipline that combines knowledge and techniques from data mining, process modeling, and process model analysis~\cite{Aalst2016}. Process mining techniques concern the analysis of events that are logged during process execution, where event records contain information on what was done, by whom, for whom, where, when, etc. Events are grouped into cases (process instances), e.g. per patient for a hospital log, or per insurance claim for an insurance company. \emph{Process discovery} plays an important role in process mining, focusing on extracting interpretable models of processes from event logs. One of the attributes of the events is usually used as its label and its values become transition/activity labels in the process models generated by process discovery algorithms.\looseness=-1

The scope of process mining have broadened in recent years from business process management to other application domains, one of them being analysis of events of human behavior with data originating from sensors in smart home environments \cite{Sztyler2015,Tax2016a,Tax2016b}. Table \ref{tab:example_log} shows an example of such an event log. Events in the event log are generated by e.g. motion sensors placed in the home, power sensors placed on appliances, open/close sensors placed on closets and cabinets, etc. Particularly challenging in applying process mining in this application domain is the extraction of meaningful event labels that allow for discovery of insightful process models. Simply using the sensor that generates an event (the \emph{sensor} column in Table \ref{tab:example_log}) as event label is shown to produce non-informative process models that overgeneralize the event log and allow for too much behavior \cite{Tax2016a}. Abstracting sensor-level events into events at the level of human activity (e.g. \emph{eating}, \emph{sleeping, etc.}) using techniques closely related to techniques used in the activity recognition field helps to discover more behaviorally more constrained and insightful process models \cite{Tax2016b}, but applicability of this approach relies on the availability of a reliable diary of human behavior at the activity level, which is often just impossible to obtain.\looseness=-1

\begin{table}[t]
\centering\scriptsize
\caption{An example of an event log from a smart home environment.}
\scalebox{0.7}{
	\begin{tabular}{ccccc}
		\toprule
		Id & Timestamp & Address & Sensor & Sensor value\\
		\midrule
		1 & \textcolor{gray}{03/11/2015} 04:59:54 & \textcolor{gray}{Mountain Rd. 7} & Motion sensor - Bedroom  &  1\\
		2 & \textcolor{gray}{03/11/2015} 06:04:36 & \textcolor{gray}{Mountain Rd. 7} & Motion sensor - Bedroom & 1\\
		3 & \textcolor{gray}{03/11/2015} 08:45:12 & \textcolor{gray}{Mountain Rd. 7} & Motion sensor - Living room &  1\\
		4 & \textcolor{gray}{03/11/2015} 09:10:10 & \textcolor{gray}{Mountain Rd. 7} & Motion sensor - Kitchen & 1\\
		5 & \textcolor{gray}{03/11/2015} 09:12:01 & \textcolor{gray}{Mountain Rd. 7} & Power sensor - Water cooker & 1200\\
		6 & \textcolor{gray}{03/11/2015} 09:15:45 & \textcolor{gray}{Mountain Rd. 7} & Power sensor - Water cooker & 0\\
		\dots & \textcolor{gray}{03/11/2015} \dots & \textcolor{gray}{Mountain Rd. 7} & \dots & \dots \\
		\midrule
		7 & \textcolor{gray}{03/12/2015} 01:01:23 & \textcolor{gray}{Mountain Rd. 7} & Motion sensor - Bedroom &  1\\
		8 & \textcolor{gray}{03/12/2015} 03:13:14 & \textcolor{gray}{Mountain Rd. 7} & Motion sensor - Bedroom &  1\\
		9 & \textcolor{gray}{03/12/2015} 07:24:57 & \textcolor{gray}{Mountain Rd. 7} & Motion sensor - Bedroom &  1\\
		10 & \textcolor{gray}{03/12/2015} 08:34:02 & \textcolor{gray}{Mountain Rd. 7} & Motion sensor - Bedroom &  1\\
		11 & \textcolor{gray}{03/12/2015} 09:12:00 & \textcolor{gray}{Mountain Rd. 7} & Motion sensor - Living room & 1\\
		\dots & \textcolor{gray}{03/12/2015} \dots & \textcolor{gray}{Mountain Rd. 7} & \dots & \dots \\
		\midrule
		12 & \textcolor{gray}{03/14/2015} 03:41:46 & \textcolor{gray}{Mountain Rd. 7} & Motion sensor - Bedroom & 1\\
		13 & \textcolor{gray}{03/14/2015} 05:00:17 & \textcolor{gray}{Mountain Rd. 7} & Motion sensor - Bedroom & 1\\
		14 & \textcolor{gray}{03/14/2015} 08:52:32 & \textcolor{gray}{Mountain Rd. 7} & Motion sensor - Bedroom & 1\\
		15 & \textcolor{gray}{03/14/2015} 09:30:54 & \textcolor{gray}{Mountain Rd. 7} & Motion sensor - Living room & 1\\
		16 & \textcolor{gray}{03/14/2015} 09:35:25 & \textcolor{gray}{Mountain Rd. 7} & Power sensor - TV & 160 \\
		17 & \textcolor{gray}{03/14/2015} 10:27:37 & \textcolor{gray}{Mountain Rd. 7} & Power sensor - TV & 0 \\
		\dots & \textcolor{gray}{03/14/2015} \dots & \textcolor{gray}{Mountain Rd. 7} & \dots & \dots \\
		\midrule
		\dots&\dots&\dots&\dots&\dots\\
		\bottomrule
	\end{tabular}}
	\label{tab:example_log}
	\vspace{-0.4cm}
\end{table}

In our earlier work \cite{Tax2016a} we showed that better process models can be discovered by taking the name of the sensor that generated the event as a starting point for the event label and then refining these labels using information on the time within the day at which the event occurred. The refinements used in \cite{Tax2016a} were based on domain knowledge, and not identified automatically from the data. In this paper, we aim at automatic generation of semantically interpretable label refinements that can be explained to the user, by basing label refinements on data attributes of events. We explore methods to bring parts of the timestamp information to the event label in an intelligent and fully automated way, with the end goal of discovering behaviorally more precise and therefore more insightful process models.\looseness=-1

We start by introducing basic concepts and notations used in this paper in Section \ref{sec:preliminaries}. In Section \ref{sec:approach}, we introduce a framework for the generation of event labels refinements based on the time attribute. In Section \ref{sec:case_study}, we apply this framework on a real life smart home data set and show the effect of the refined event labels on process discovery.  We continue by describing related work in Section \ref{sec:related_work} and conclude in Section \ref{sec:conclusion}.\looseness=-1

\section{Preliminaries}
\label{sec:preliminaries}
In this section we introduce the notions related to event logs and relabeling functions for traces and then define the notions of refinements and abstractions. We also introduce the Petri net process model notation.

We use the usual sequence definition, and denote a sequence by listing its elements, e.g. we write $\langle a_1,a_2,\dots,a_{n} \rangle$ for a (finite) sequence $s:\{1,\dots,n\}\to A$ of elements from some alphabet $A$, where $s(i)=a_i$ for any $i \in \{1,\dots,n\}$; $|s|$ denotes the length of sequence $s$; $s_1 s_2$ denotes the concatenation of sequences $s_1$ and $s_2$. A \emph{language} $\Lan$ over an alphabet $A$ is a set of sequences over $A$. $\Lan^p$ is the prefix closure of a language $\Lan$ (with $\Lan\subseteq\Lan^p$).

An event is the most elementary element of an event log. Let $\mathcal{I}$ be a set of event identifiers, and $\mathcal{A}_1 \times \dots \times \mathcal{A}_{n}$ be an attribute domain consisting of $n$ attributes (e.g. timestamp, resource, activity name, cost, etc.). An event is a tuple $e=(i,a_1,\dots,a_{n})$, with $i\in\mathcal{I}$ and $(a_1, \dots, a_{n})\in \mathcal{A}_1 \times \dots \times \mathcal{A}_{n}$. The \emph{event label} of an event is the attribute set $(a_1\dots,a_n)$; $e_i$, and $e_a$ respectively denote the identifier and label of event $e$. The timestamp attribute of an event is denoted by $a_t$. $\mathcal{E}=\mathcal{I}\times \mathcal{A}_1 \times \dots \times \mathcal{A}_{n}$ is a universe of events over $\mathcal{A}_1, \dots, \mathcal{A}_{n}$.
The rows of Table \ref{tab:example_log} are events from an event universe over the event attributes \emph{timestamp}, \emph{sensor}, \emph{address}, and \emph{sensor value}.

Events are often considered in the context of other events. We call $E\subseteq\mathcal{E}$ an \emph{event set} if $E$ does not contain any events with the same event identifier. The events in Table \ref{tab:example_log} together form an event set. A \emph{trace} $\sigma$ is a finite sequence formed by the events from an event set $E\subseteq{\mathcal{E}}$ that respects the time ordering of events, i.e. for all $k,m\in\N$, $1\leq k < m \leq |E|$, we have: $\sigma(k)_t\leq \sigma(m)_t$. We define the \emph{universe of traces} over event universe $\mathcal{E}$, denoted $\Sigma(\mathcal{E})$, as the set of all possible traces over $\mathcal{E}$. We omit $\mathcal{E}$ in $\Sigma(\mathcal{E})$ and use the shorter notation $\Sigma$ when the event universe is clear from the context.\looseness=-1

Often it is useful to partition an event set into smaller sets in which events belong together according to some criterion. We might for example be interested in discovering the typical behavior within a household over the course of a day. In order to do so, we can e.g.~group together events with the same \emph{address} and the same day-part of the \emph{timestamp}, as indicated by the horizontal lines in Table \ref{tab:example_log}. For each of these event sets, we can construct a trace; time stamps define the ordering of events within the trace. For events of a trace having the same time stamps, an arbitrary ordering can be chosen within a trace.

An \emph{event partitioning function} is a function $ep: \mathcal{E} \to T_{id}$ that defines the partitioning of an arbitrary set of events $E\subseteq\mathcal{E}$ from a given event universe $\mathcal{E}$ into event sets $E_1,\ldots,E_j,\ldots$ where each $E_j$ is the maximal subset of $E$ such that for any $e_1, e_2\in E_j$, $ep(e_1)= ep(e_2)$; the value of $ep$ shared by all the elements of $E_j$ defines the value of the \emph{trace attribute} $T_{id}$. Note that multidimensional trace attributes are also possible, i.e. a combination of the name of the person performing the event activity and the date of the event, so that every trace contains activities of one person during one day. The event sets obtained by applying an event partitioning can be transformed into traces (respecting the time ordering of events).\looseness=-1

An event log $L$ is a finite set of traces $L \subseteq \Sigma(\mathcal{E})$. $A_L\subseteq\mathcal{A}_1 \times \dots \times \mathcal{A}_{n}$ denotes the \emph{alphabet of event labels} that occur in log $L$. The traces of a log are often transformed before doing further analysis: very detailed but not necessarily informative event descriptions are transformed into some \emph{informative} and \emph{repeatable} labels. For the labels of the log in Table~\ref{tab:example_log}, the sensor values could be abstracted to \textit{on}, and \textit{off} or labels can be redefined to a subset of the event attributes, e.g. leaving the sensor values out completely. Next to that, if the event partitioning function maps each event from Table~\ref{tab:example_log} to its address and the day-part of the timestamp, these attributes (indicated in gray) become the trace attribute and can safely be removed from individual events.

After this relabeling step, some traces of the log can become identically labeled (the event id's would still be different). The information about the number of occurrences of a sequence of labels in an event log is highly relevant for process mining, since it allows differentiating between the mainstream behavior of a process (frequently occurring behavioral patterns) and exceptional behavior.

Let $\E, \E'$ be an event universe. A function $l: \E \to \E'$ is an \emph{event relabeling function}. A relabeling function can be used to obtain more useful event labels than the full set of event attribute values. We lift $l$ to event logs. Let $\E,\E_1,\E_2$ be event universes with $\E,\E_1,\E_2$ being pairwise different. Let $l_1: \E \to \E_1$ and $l_2: \E \to \E_2$ be event relabeling functions. Relabeling function $l_1$ is a \emph{refinement} of relabeling function $l_2$, denoted by $l_1\preceq l_2$, iff $\forall_{e_1,e_2\in \E}:l_1(e_1)=l_1(e_2)\implies l_2(e_1)=l_2(e_2)$; $l_2$ is then called an \emph{abstraction} of $l_1$.

The goal of process discovery is to discover a process model that represents the behavior seen in an event log. A frequently used process modeling notation in the process mining field is the Petri net~\cite{Reisig1998}. Petri nets are directed bipartite graphs consisting of transitions and places, connected by arcs. Transitions represent activities, while places represent the enabling conditions of transitions. Labels are assigned to transitions to indicate the type of activity that they model. A special label $\tau$ is used to represent invisible transitions, which are only used for routing purposes and not recorded in the execution log.

A \emph{labeled Petri net} $N=\langle P,T,F,A_M,\ell\rangle$ is a tuple where $P$ is a finite set of places, $T$ is a finite set of transitions such that $P \cap T = \emptyset$,  $F \subseteq (P \times T) \cup (T \times P)$ is a set of directed arcs, called the flow relation, $A_M$ is an alphabet of labels representing activities, with $\tau \notin A_M$ being a label representing  invisible events, and $\ell:T\rightarrow A_M\cup \{\tau\}$ is a labeling function that assigns a label to each transition. For a node $n \in (P \cup T)$ we use $\bullet n$ and $n \bullet$ to denote the set of input and output nodes of $n$, defined as $\bullet n =\{n|(n',n)\in F\}$ and $n \bullet =\{n|(n,n')\in F\}$. An example of a Petri net can be seen in Figure \ref{fig:example_petri_net}, where circles represent places and squares represent transitions. Gray transitions with smaller width represent $\tau$ transitions.\looseness=-1

\begin{figure}[t]
	\centering
	\scalebox{0.67}{
		\begin{tikzpicture}
		[node distance=1.11cm,
		on grid,>=stealth',
		bend angle=20,
		auto,
		every place/.style= {minimum size=4mm},
		every transition/.style = {minimum size = 8.5mm}

		]
		\node [place, tokens = 1] (p1) [label=below:$p_1$]{};
		\node [transition] (b) [label=below:$t_1$, right = of p1] {$A$}
		edge [pre] node[auto] {} (p1);
		\node [place] (p2) [label=below:$p_2$, right = of b]{}
		edge [pre] node[auto] {} (b);
		\node [transition] (c) [label=below:$t_3$, below right = of p2] {$C$}
		edge [pre] node[auto] {} (p2);
		\node [transition] (c2) [label=below:$t_2$, above right = of p2] {$B$}
		edge [pre] node[auto] {} (p2);
		\node [place] (p3) [label=below:$p_3$, above right = of c]{}
		edge [pre] node[auto] {} (c)
		edge [pre] node[auto] {} (c2);
		\node [transition] (d) [label=below:$t_4$, right = of p3] {$D$}
		edge [pre] node[auto] {} (p3);
		\node [place] (p4) [label=above:$p_4$, below right = of d]{}
		edge [pre] node[auto] {} (d);
		\node [transition] (e) [label=above:$t_5$, right = of p4] {$F$}
		edge [pre] node[auto] {} (p4);
		\node [place] (p6) [label=above:$p_6$, right = of e]{}
		edge [pre] node[auto] {} (e);
		\node [place] (p5) [label=above:$p_5$, above right = of d]{}
		edge [pre] node[auto] {} (d);
		\node [transition] (f) [label=above:$t_6$, right = of p5] {$E$}
		edge [pre] node[auto] {} (p5);
		\node [place] (p7) [label=above:$p_7$, right = of f]{}
		edge [pre] node[auto] {} (f);
		\node [transition] (g) [label=below:$t_7$, minimum width=3mm, fill=lightgray, below right = of p7] {}
		edge [pre] node[auto] {} (p7)
		edge [pre] node[auto] {} (p6);
		\node [place] (p8) [label=below:$p_8$, right = of g]{}
		edge [pre] node[auto] {} (g);

		\end{tikzpicture}}
	\caption{An example Petri net.}
	\label{fig:example_petri_net}
	\vspace{-0.25cm}
\end{figure}

A state of a Petri net is defined by its \emph{marking} $M \in \mathbb{N}^{P}$ being a multiset of places. A marking is graphically denoted by putting $M(p)$ tokens on each place $p\in P$. A pair $(N,M)$ is called a marked Petri net. State changes occur through transition firings. A transition $t$ is enabled (can fire) in a given marking $M$ if each input place $p\in \pre{t}$ contains at least one token. Once a transition fires, one token is removed from each input place of $t$  and one token is added to each output place of $t$, leading to a new marking. An \emph{accepting Petri net} is a 3-tuple $(N,M_i,M_f)$ with $N$ a labeled Petri net, $M_i$ an initial marking, and $M_f$ a set of final markings.
Many process modeling notations, including accepting Petri nets, have formal executional semantics and a model defines a \emph{language of accepting traces} $\Lan$. For the Petri net in Figure \ref{fig:example_petri_net}, the language of accepting traces is $\{\langle A,B,D,E,F\rangle,\langle A,B,D,F,E\rangle,\langle A,C,D,E,F\rangle,\langle A,C,D,F,E\rangle\}$.

\section{A Framework for Time-based Label Refinements}
\label{sec:approach}
To generate potential label refinements for every label based on time we take a clustering based approach by identifying dense areas in time space for each label. The time part of the timestamps consists of values between $\textit{00:00:00}$ and $\textit{23:59:59}$, equivalent to the timestamp attribute from Table \ref{tab:example_log} with the day-part of the timestamp removed. This timestamp can be transformed into a real number hourfloat representation in interval $[0,24)$. We chose to apply soft clustering (also referred to as fuzzy clustering), which has the benefit of assigning to each data point a likelihood of belonging to each cluster. A well-known approach to soft clustering is based on the combination of the Expectation-Maximization (EM) algorithm with mixture models, which are probability distributions consisting of multiple components of the same probability distribution. Each component in the mixture represents one cluster and the probability of a data point belonging to that cluster is the probability that this cluster generated that data point. The EM algorithm is used to obtain a maximum likelihood estimate of the mixture model parameters, i.e. the parameters of the probability distributions in the mixture.\looseness=-1

\begin{figure}[t]
	\centering
	\includegraphics[width=0.675\linewidth]{gmm_plates_cupboard}
	\caption{The histogram representation and a Gaussian Mixture Model fitted to timestamps values of the plates cupboard sensor in the van Kasteren data set \cite{Kasteren2008}.\looseness=-1}
	\label{fig:problem_data_set}
	\vspace{-0.2cm}
\end{figure}

 A well-known example of a mixture model is the Gaussian Mixture Model (GMM), where the components in the mixture distributions are normal distributions. The data space of time is, however, non-euclidean: it has a circular nature, e.g. $23.99$ is closer to $0$ than to $23$. This circular nature of the data space introduces problems for GMMs, as shown in Figure \ref{fig:problem_data_set}. The GMM fitted to the timestamps of the sensor events consists of two components, one with the mean at $9.05$  and one with a mean at $20$. The histogram representation of the same data shows that some events happened just after midnight, which is actually closer on the clock to $20$ than to $9.05$. The GMM however is unaware of the circularity of the clock, which results  in the mixture model that seems inappropriate when visually comparing with the histogram. The field of circular statistics (also referred to as directional statistics), concerns analysis of such circular data spaces (cf. \cite{Mardia2009}).\looseness=-1

Here, we introduce a framework for generating refinements of event labels based on time attributes using techniques from the field of circular statistics. This framework consists of three stages:
\begin{description}
	\item[Data-model pre-fitting stage]{A known problem with many clustering techniques is that they return clusters even when the data should not be clustered. In this stage we assess how many clusters the events of a sensor type contain.}
	\item[Data-model fitting stage]{In this stage we cluster the events of a sensor type by timestamp using a mixture consisting of components that take into account the circularity of the data.}
	\item[Data-model post-fitting stage]{In this stage the quality of the label refinements is assessed from both a cluster quality perspective and a process model (event ordering statistics) perspective.}
\end{description}

\subsection{Data-model pre-fitting stage}
We now describe a test for uniformity, a test for unimodality, and a method to select the number of clusters in the data.
\subsubsection{Uniformity Check - Rao's Spacing Test}
Rao's spacing test~\cite{Rao1976} tests the uniformity of the timestamps of the events from a sensor around the circular clock. This test is based on the idea that uniform circular data is distributed evenly around the circle, and $n$ observations are separated from each other $\frac{360}{n}$ degrees. The null hypothesis is that the data is uniform around the circle.

Given $n$ successive observations $f_1,\ldots,f_n$, either clockwise or counterclockwise, the test statistics $U$ for Rao's Spacing Test is defined as $U = \frac{1}{2}\sum_{i = 1}^{n}\mid T_i - \lambda \mid$, where $\lambda = \frac{360^\circ}{n}$, $T_i = f_{i + 1} - f_{i}$ for $1 \le i \le n - 1$ and $T_n = (360^\circ - f_n) + f_1$.

\subsubsection{Unimodality Check - Hartigan's Dip Test}
Hartigan's dip tests~\cite{Hartigan1985} the null hypothesis that the data follows a unimodal distribution on a circle. When the null hypothesis can be rejected, we know that the distribution of the data is at least bimodal. Hartigan's dip test measures the maximum difference between the the empirical distribution function and the unimodal distribution function that minimizes that maximum difference.

\subsubsection{Number of Component Selection - Bayesian Information Criterion}
The Bayesian Information Criterion (BIC)~\cite{Schwarz1978} introduces a penalty for the number of model parameters to the evaluation of a mixture model. Adding a component to a mixture model increases the number of parameters of the mixture with the number of parameters of the distribution of the added component. The likelihood of the data given the model can only increase by adding extra components, adding the BIC penalty results in a trade-off between number of components and the likelihood of the data given the mixture model. BIC is formally defined as $\textit{BIC} = -2 * ln\hat{L} + k * ln(n)$, where $\hat{L}$ is a maximized value for the data likelihood, $n$ is the sample size, and $k$ is the number of parameters to be estimated. A lower BIC value indicates a better model. We start with $1$ component, and iteratively increase from $k$ to $k+1$ components as long as the decrease in BIC is larger than 10, which is the threshold for decisive evidence of high BIC \cite{Kass1995}.

\subsection{Data-model fitting stage}
We cluster events generated by one sensor using a mixture model consisting of components of the von Mises distribution, which is a circular version of the normal distribution. This technique is based on the approach of Banerjee et al. \cite{Banerjee2005}, who introduce a clustering method based on a mixture of von Mises-Fisher distribution components, which is a generalization of the $2$-dimensional von Mises distribution to $n$-dimensional spheres. A probability density function for a von Mises distribution with mean direction $\mu$ and concentration parameter $\kappa$ is defined as $\textit{pdf}(\theta \mid \mu, \kappa) = \frac{1}{2\pi I_0(\kappa)}e^{\kappa\cos(\theta - \mu)}$, where mean $\mu$ and data point $\theta$ are expressed in radians on the circle, such that $0 \le \theta \le 2\pi, ~0 \le \mu \le 2\pi, ~\kappa \ge 0$.
\textit{$I_0$} represents the modified Bessel function of order 0, defined as $I_0(k) = \frac{1}{2\pi}\int_0^{2\pi} e^{\kappa\cos(\theta)}d\theta$. As $\kappa$ approaches 0, the distribution becomes uniform around the circle. As $\kappa$ increases, the distribution becomes relatively concentrated around the mean $\mu$ and the von Mises distribution starts to approximate a normal distribution. We fit a mixture model of von Mises components using the package movMF \cite{Hornik2014} provided in R.

\subsection{Data-model post-fitting stage}
After fitting a mixture of von Mises distributions to the sensor events, we perform a goodness-of-fit test to check whether the data could have been generated from this distribution. We describe the Watson $U^2$ statistic~\cite{Watson1962}, a goodness-of-fit assessment based on hypothesis testing.
The Watson $U^2$ statistic measures the discrepancy between the cumulative distribution function $F(\theta)$ and the empirical distribution function $F_n(\theta)$ of some sample $\theta$ drawn from some population and is defined as $U^2 = n\int_0^{2\pi} \Big[ F_n(\theta) - F(\theta) - \int_0^{2\pi} \big\{ F_n(\phi) - F(\phi) \big\} dF(\phi) \Big]^2 dF(\theta)$.

Furthermore we assess the quality of refining the event label into a new label for each cluster from a process perspective using the label refinement evaluation method described in \cite{Tax2016a}. This method tests whether the log statistics that are used in many process discovery algorithms become significantly more deterministic by applying the label refinement.

\section{Case Study}
\label{sec:case_study}

\begin{figure}[t]
	\centering
	\begin{minipage}[t]{.5\textwidth}
	\centering
	\includegraphics[width=\linewidth]{bic_bedroom_door}
	\caption{BIC values for different numbers of components in the mixture model.}
	\label{fig:num_component_selection_kasteren}
	\end{minipage}\hfill
	\begin{minipage}[t]{.4\textwidth}
		\centering
		\vspace{-3.9cm}
		\captionof{table}{Estimated parameters for a mixture of von Mises components for bedroom door sensor events.\looseness=-1}
		\label{tab:bedroom_components}
		\vspace{0.5cm}
		\begin{tabular}{c|ccc}
			\toprule
			Cluster & $\alpha$ & $\mu$ (radii) & $\kappa$ \\
			\midrule
			Cluster 1 & 0.76 & 2.05 & 3.85\\
			Cluster 2 & 0.24 & 5.94 & 1.56\\
			\bottomrule
		\end{tabular}
	\end{minipage}
	\vspace{-0.2cm}
\end{figure}

\begin{figure}[t]
	\centering
	\subfloat[]{
		\includegraphics[width=0.98\linewidth]{van_kasteren_original_3}
		\label{fig:kasteren_original}
	}\\
	\subfloat[]{
		\hspace{-0.025\linewidth}
		\includegraphics[width=1.0125\linewidth]{van_kasteren_refined_3}
		\label{fig:kasteren_refined}
	}
	\caption{Process models discovered on the van Kasteren data with sensor-level labels (a) and refined labels (b) with the Inductive Miner infrequent (20\% filtering) \cite{Leemans2013}.}
	\vspace{-0.2cm}
\end{figure}

\begin{figure}
	\centering
	\hspace{-0.025\linewidth}
	\includegraphics[width=1.0125\linewidth]{van_kasteren_refined_sec}
	\caption{Inductive Miner infrequent (20\% filtering) \cite{Leemans2013} result after a second label refinement.}
	\label{fig:kasteren_refined2}
	\vspace{-0.2cm}
\end{figure}

We show the results of our time-based label refinements approach on the real life smart home data set described in van Kasteren et al. \cite{Kasteren2008}. The van Kasteren data set consists of 1285 events divided over fourteen different sensors. We segment in days from midnight to midnight to define cases. Figure \ref{fig:kasteren_original} shows the process model discovered on this event log with the Inductive Miner infrequent \cite{Leemans2013} with 20\% filtering, which discovers a process model that describes the most frequent 80\% of behavior in the log. Note that this process model overgeneralises allowing too much behaviour. At the beginning a (possibly repeated) choice is made between five transitions. At the end of the process, the model allows any sequence over the alphabet of  five activities, where each activity occurs at least once.

We illustrate our proof of concept by applying the framework to the \emph{bedroom door} sensor. Rao's spacing test results in a test statistic of $241.0$ with $152.5$ being the critical value for significance level $0.01$, indicating that we can reject the null hypothesis of a uniformly distributed set of \emph{bedroom door} timestamps. Hartigan's dip test results in a p-value of $3.95\times10^{-4}$, indicating that we can reject the null hypothesis that there is only one cluster in the \emph{bedroom door} data. Figure \ref{fig:num_component_selection_kasteren} shows the BIC values for different numbers of components in the model. The figure indicates that there are two clusters in the data, as this corresponds to the lowest BIC value. Table \ref{tab:bedroom_components} shows the mean and $\kappa$ parameters of the two clusters found by optimizing the von Mises mixture model with the EM algorithm. A value of $0=2\pi$ radii equals midnight. After applying the von Mises mixture model to the \emph{bedroom door} events and assigning each event to the maximum likelihood cluster we obtain a time range of [3.08-10.44] for cluster 1 and a time range of [17.06-0.88] for cluster 2. The Watson $U^2$ test results in a test statistic of $0.368$ and $0.392$ for cluster 1 and 2 respectively with a critical value of $0.141$ for a $0.01$ significance level, indicating that the data is likely to be generated by the two von Mises distributions found. The label refinement evaluation method \cite{Tax2016a} finds statistically significant differences between the events from the two \emph{bedroom door} clusters with regard to their control-flow relations with other activities in the log for 10 other activities using the significance level of $0.01$, indicating that the two clusters are different from a control-flow perspective. Figure \ref{fig:kasteren_refined} shows the process model discovered with the Inductive Miner infrequent with 20\% filtering after applying this label refinement to the van Kasteren event log. The process model still overgeneralizes in general, but the label refinement does help restricting the behavior, as it shows that the evening \emph{bedroom door} events are succeeded by one or more events of type \emph{groceries cupboard}, \emph{freezer}, \emph{cups cupboard}, \emph{fridge}, \emph{plates cupboard}, or \emph{pans cupboard}, while the morning \emph{bedroom door} events are followed by one or more \emph{frontdoor} events. It seems that this person generally goes to the bedroom in-between coming home from work and starting to cook. The loop of the \emph{frontdoor} events could be caused by the person leaving the house in the morning for work, resulting in no logged events until the person comes home again by opening the \emph{frontdoor}. Note that in Figure \ref{fig:kasteren_original} \emph{bedroom door} and \emph{frontdoor} events can occur an arbitrary number of times in any order. Figure \ref{fig:kasteren_original} furthermore does not allow for the \emph{bedroom door} to occur before the whole block of kitchen-located events at the beginning of the net. \looseness=-1

Label refinements can be applied iteratively. Figure \ref{fig:kasteren_refined2} shows the effect of a second label refinement step, where \emph{Plates cupboard} using the same methodology is refined into two labels, representing time ranges [7.98-14.02] and [16.05-0.92] respectively. This refinement shows the additional insight that the evening version of the \emph{Plates cupboard} occurs in directly before or after the microwave.

\section{Related Work}
\label{sec:related_work}
Refining event labels in the event log is closely related to the task of mining process models with duplicate activities, in which the resulting process model can contain multiple transitions/nodes with the same label. From the point of view of the behavior allowed by a process model, it makes no difference whether a process model is discovered on an event log with refined labels, or whether a process model is discovered with duplicate activities such that each transition/node of the duplicate activity precisely covers one versions of the refined label. The first process discovery algorithm capable of discovering duplicate tasks was proposed by Herbst and Karagiannis in 2004 \cite{Herbst2004}, after which many others have been proposed, including the Genetic Miner \cite{DeMedeiros2007}, the Evolutionary Tree Miner \cite{Buijs2012}, the $\alpha^*$-algorithm \cite{Li2007}, the $\alpha^\#$-algorithm \cite{Gu2008}, the EnhancedWFMiner \cite{Folino2009}, and a simulated annealing based algorithm \cite{Song2008}. An alternative approach has been proposed by V\'{a}zques-Barreiros \cite{Vazquez-Barreiros2015} et al., who describe a local search based approach to repair a process model to include duplicate activities, starting from an event log and a process model without duplicate activities. Existing work on mining models with duplicate activities all base their duplicate activities on how well the event log fits the process model, and do not try to find any semantic difference between the multiple versions of the activities in the form of data attribute differences.\looseness=-1

The work that is closest to our work is the work by Lu et al. \cite{Lu2016}, who describe an approach to pre-process an event log by refining event labels with the goal of discovering a process model with duplicate activities. The method proposed by Lu et al., however, does not base the relabelings on data attributes of those events but instead bases them solely on the control flow context, leaving uncertainty whether two events relabeled differently are actually semantically different.

Another area of related work is data-aware process mining, where the aim is to discover rules with regard to data attributes of events that decide decision points in the process. De Leoni and van der Aalst \cite{DeLeoni2013} proposed a method that discovers data guards for decision points in the process based on alignments and decision tree learning. This approach relies on the discovery of a behaviorally well-fitting process model from the original event log. When only overgeneralizing process models (i.e. allowing for too much behavior) can be discovered from an event log, the correct decision points might not be present in the discovered process model at all, resulting in this approach not being able to discover the data dependencies that are in the event log. Our label refinements use data attributes prior to process discovery to enable discover more behaviorally constrained process models by bringing parts of the event attribute space to the event label. \looseness=-1

\section{Conclusion \& Future Work}
\label{sec:conclusion}
We have proposed a framework based on techniques from the field of circular statistics to refine event labels automatically based on their timestamp attribute. We have shown through a proof of concept on a real life event log that this framework can be used to discover label refinements that allow for discovery of more insightful and behaviorally more specific process models.
An interesting area of future work is to explore the use of other types of event data attributes to refine labels, e.g. power values of sensors. A next research step would be to explore label refinements based on multiple data attributes combined. This would bring challenge of clustering on partially circular and partially euclidean data spaces.\looseness=-1

\bibliographystyle{acm}


\section{Introduction}
The thermodynamic state of gas in cosmic web filaments has important
implications for observations and theoretical predictions. It will
affect halos which are fed by those filaments, as well as small halos
that accrete onto those filaments as part of the cosmic hierarchical
growth.

Spherical virialization of gas in halos has been a prediction of
galaxy formation models for decades. In particular, it has been shown
\citep{ro77,silk77,binney77,white91} that a comparison between cooling
times and ages of galactic halos can predict the transition from
galaxies to groups and clusters. \citet[hereafter BD03]{bd03} derived
a stability criterion against gravitational collapse of the gas in the
presence of significant cooling. They find that for halos below
$M_{crit}\simeq 10^{12} \ensuremath{M_{\odot}} $ a hot gaseous halo is not expected to
form, and gas will free-fall until it reaches the disk, at which point
it will stop, radiating its kinetic energy abruptly at that
point. This has been confirmed in multiple hydrodynamical simulations
\citep[e.g.\@ ][]{keres05,ocvirk08,cafg11} and successfully reproduces
star forming galaxies at high-z \citep{db06,dekel09} and the
color-magnitude bi-modality
\citep{db08,croton06,cattaneo06}. Observational
indications of this scenario are
gradually accumulating \citep[e.g.\@][]{dijkstra09,kimm10,martin15}.

In this letter we derive a criterion for the stability of virial shocks around filaments and sheets that form the cosmic web,
analogous to the BD03 criterion for halos. Following \citet[hereafter
  FG84]{fg84}, we construct self-similar density profiles of filaments.  We apply our stability analysis
to these profiles to identify filaments around which a stable virial shock is
expected to form. This criterion is translated to a more useful form
by identifying which halos are expected to be fed by these
filaments. We find that filament instability influences a
large portion of halos in the universe throughout cosmic age.

The stability of filaments has been addressed before, numerically
\citep{harford11} and analytically \citep{freundlich14,breysse14}, but
without taking into account cooling, and by analyzing the stability of
initially static filaments, ignoring the effects of the shock at the
filaments' edge.

In \S~2 we derive the stability criterion for the existence of
virialized gas in 1,2 and 3 dimensional collapse. In \S~3 we
relate our local criterion to cosmic filaments according to
the self-similar solutions of FG84. In \S~4 we relate
these filaments to typical halo masses that will likely be fed by
them. In \S~5 we summarize and conclude.

\section{virial shock stability in spherical, cylindrical and planar geometry}
\label {sec:analytics}
The analysis in BD03 was performed for spherical accretion in the presence of cooling. In this section we generalize that derivation for
infall onto spherical,
cylindrical (onto filaments), and planar  (onto sheets or
disks) objects.

The ideal gas equation of state (EoS) is
\begin{equation}
P=(\gamma-1)\rho e,
\label{eq:ideal}
\end{equation}
with
$\gamma=\left(\frac{\partial~ {\rm ln}P}{\partial~
    \rm{ln}\rho}\right)_s$ the adiabatic index, and $\rho,e,P$ the density, internal energy
and pressure respectively.
The adiabatic index measures the ``stiffness'' of the EoS, or the adiabatic pressure response to compression. By analogy, we define the effective EoS index, $\gamma_{\rm eff},$ of a parcel
of gas undergoing compression as its pressure response  along a Lagrangian trajectory:
\begin{equation}
\gamma_{\rm eff} \equiv
\frac{d~\rm{ln}P}{d~\rm{ln}\rho}=\frac{\rho}{P}\frac{\dot{P}}{\dot{\rho}},
\label{eq:geff}
\end{equation}
with upper dot implying a Lagrangian time derivative.
\refereebf{Following BD03, this form is transformed into}
\begin{equation}
\gamma_{\rm eff}=\gamma-\frac{\rho}{\dot{\rho}}\bigg{\slash} \frac{e}{q}.
\label{eq:geff2}
\end{equation}
\refereebf{with $q$ related to the \citet{sutherland93} cooling rate according to}
\begin{equation}
q=\rho\Lambda_{\rm cool}(T,Z)=\rho(N_a^2
\chi^2/\mu^2)\Lambda_{\rm mic}(T),
\label{eq:lambda}
\end{equation}
\refereebf{
where ${N_a}/{\mu}$ is the
  average number of molecules per unit mass and $\chi$ is the number
  of electrons per particle. This parametrization allows for altered ionization state of the gas \citep[see ][]{cantalupo10}.
 This equation illustrates how the effective stiffness of the gas is set by two competing timescales: the compression time, $\frac{\rho}{\dot{\rho}}$, and the cooling time, $\frac{e}{q}$. If radiative losses are significant during a compression time, we get effective softening, $\gamma_{\rm eff}<\gamma$. In the absence of cooling, we recover $\gamma_{\rm eff}=\gamma$.}

The compression time for post-shock infalling gas depends on the velocity and geometry of the infall through the Lagrangian continuity equation,
\begin{equation}
\dot{\rho}=-\rho \nabla \cdot u=-\rho\left((n-1)\frac{u}{r}+\frac{\partial
  u}{\partial r}\right),
\label{eq:nabla}
\end{equation}
with $u$ being the velocity of the gas, and $n$ the dimensionality of
the infall: $n=1,2,3$ corresponds to planar, filamentary and spherical collapse respectively. Assuming that the post-shock flow is
homologous,
\begin{equation}
u=\frac{u_s}{r_s}r
\label{eq:homo}
\end{equation}
with $r_s$ the shock radius and $u_s$ the velocity directly below the shock, we can apply \equ{nabla} to show  that the post-shock gas must contract uniformly, independent of $r$,
\begin{equation}
\frac{\dot{\rho}}{\rho}=-n\frac{u_s}{r_s}={\rm Const}.
\label{eq:rhodot}
\end{equation}
This homologous behavior is seen in 1D simulations for the spherical
case \citep{bd03}, and is a standard assumption in hydrostatic
stability calculations. It also roughly matches self-similar
solutions of gaseous spherical infall \citep{bertschinger85}.

We consider a quasi-static configuration where the shock radius, $r_s$, and the post-shock velocity profile \equ{homo} are approximately constant. We test the stability of this configuration by assuming that it is initially hydrostatic,
\begin{equation}
\ddot{r}=-\frac{1}{\rho}\nabla P+a_g=0
\label{eq:rddot}
\end{equation}
($a_g$ being the gravitational acceleration)
but has an inward post-shock velocity, and checking the sign of
the auxiliary
acceleration or ``jerk'', $\dddot r$, that forms due to this motion.
A positive jerk
will instigate outwards motion implying a stable configuration,
while a negative $\dddot{r}$ will lead to collapse, implying instability.\footnote{The hydrostatic condition, \equ{rddot}, assumes the acceleration implied by the homologous velocity profile can be considered negligible compared to the deceleration across the shock. This is a good approximation for strong shocks (see \S~3). We avoid the inclusion of a homologous acceleration term because it complicates the derivation considerably but yields only an insignificant quantitative correction.}

The gravitational acceleration depends on the dimensionality of
the potential well. Since the perturbation is Lagrangian, the
mass enclosed below the gas parcel is constant, so
\begin{equation}
a_g=-Ar^{1-n},
\label{eq:ag}
\end{equation}
with $A$ a positive constant\footnote{for spherical, cylindrical and
  planar configuration the gravitational acceleration is
  $-GM/r^2,-2Gl/r,-2\pi G \Sigma$ respectively, with $G,M,l,\Sigma$ constants.}.
We convert the $\nabla=\frac{d}{dr}$ operator in \equ{rddot} to a Lagrangian
mass derivative, by noting that a mass element is related to a spatial
differential according to:
\begin{equation}
dm=B\rho r^{n-1}dr,
\label{eq:dm}
\end{equation}
with $B$ a positive constant\footnote{Note that the units of $dm$ here
  depend on the dimensionality of the infall, $n$: it is mass for
  spherical infall, mass per unit length for cylindrical infall
  and mass per unit area for planar infall.}.
Plugging \cref{eq:ag,eq:dm} into \equ{rddot} we find
\begin{equation}
\ddot{r}=-Br^{n-1}\frac{dP}{dm}-Ar^{1-n}=0.
\label{eq:rddot2}
\end{equation}
From here on we shall denote all mass derivatives with $'$.

The rate of change in the acceleration is the time derivative
of \equ{rddot2},
\begin{equation}
\dddot{r}=-B(n-1)r^{n-2}uP'-Br^{n-1}\dot{P'}-A(1-n)r^{-n}u,
\end{equation}
noting that $u=\dot{r}.$
Eliminating the last term by use of \equ{rddot2}, collecting terms,
and exchanging the spatial and time derivative of $P$
one gets:
\begin{equation}
\dddot{r}=-Br^{n-1}[2(n-1)\frac{u}{r}P'+(\dot{P})'].
\label{eq:dddr}
\end{equation}

The calculation of $(\dot{P})'$ is somewhat lengthy. We first derive
an expression for $\dot{P}$:
\begin{align}
  \dot{P}&=\gamma_{\rm eff}\frac{\dot{\rho}}{\rho}P=\gamma\frac{\dot{\rho}}{\rho}P-\frac{q}{e}P
  \nonumber\\
  &=\gamma\frac{\dot{\rho}}{\rho}P-\frac{\rho P
  \Lambda_{\rm cool}}{e}=\gamma\frac{\dot{\rho}}{\rho}P-(\gamma-1)\rho^2\Lambda_{\rm cool},
\label{eq:pdot2}
\end{align}
with the first equality due to \equ{geff}, the second equality to
\equ{geff2}, third to \equ{lambda} and fourth to \equ{ideal}.
By \equ{rhodot} we note that the term $\frac{\dot{\rho}}{\rho}$ is
independent of the spatial derivative. We further assume that the
cooling function, $\Lambda_{\rm cool},$ does not change significantly due to the
change of temperature of the Lagrangian mass element \footnote{This assumption is reasonable for $\Lambda_{\rm cool}(T,Z),$ except near $10^4K,$ and is necessary for an analytic solution to be possible. We neglect it here, at the risk of a slight error near the lower boundary of the unstable regime.}, so $\Lambda_{\rm cool}$ can
also be taken out of the derivative. We assume that nearby mass elements directly below the shock start
with the same thermodynamic conditions (i.e.\@ they lie on the same adiabat). Combined with the definition in \equ{geff}, this indicates that
\begin{equation}
\frac{P'}{\rho'}=\gamma_{\rm eff}\frac{P}{\rho}.
\label{eq:ptag}
\end{equation}
Differentiating \equ{pdot2} then yields:
\begin{align}
&(\dot{P})'= \label{eq:pdottag}\\
&\gamma \frac{\dot{\rho}}{\rho}P'-2(\gamma-1)\rho
\Lambda_{\rm cool}\rho'=\gamma
\frac{\dot{\rho}}{\rho}P'-2(\gamma-1)q\frac{1}{\gamma_{\rm eff}}\frac{\rho}{P}P'=\nonumber\\
&\gamma
  \frac{\dot{\rho}}{\rho}P'-\frac{2}{\gamma_{\rm eff}}\:\frac{q}{e}P'=\gamma
  \frac{\dot{\rho}}{\rho}P'-\frac{2}{\gamma_{\rm eff}}(\gamma-\gamma_{\rm eff})\frac{\dot{\rho}}{\rho}P'.\nonumber
\end{align}
The second equality is due to \equ{ptag}, the third by reverse use of
\cref{eq:ideal,eq:lambda}, and the fourth by separating
$\frac{q}{e}$ from \equ{geff2}.

Inserting \equ{pdottag} into \equ{dddr}, and converting
$\frac{\dot{\rho}}{\rho}$ and $\frac{u}{r}$ to $\frac{u_s}{r_s}$
according to \cref{eq:homo,eq:rhodot}, we finally get:
\begin{equation}
\dddot{r}=-Br^{n-1}\frac{u_s}{r_s}P'[2(n-1)-n\gamma+\frac{2n}{\gamma_{\rm eff}}(\gamma-\gamma_{\rm eff})].
\end{equation}
$u_s<0$ because below the standing shock there is inwards velocity, and
$P'<0$ so that the pressure force, $-\nabla P,$ is positive to balance gravitation in the quasi-hydrostatic halo, so the factor before the square brackets is always negative. A
positive jerk, or stability, thus occurs when
\begin{equation}
2(n-1)-n\gamma+\frac{2n}{\gamma_{\rm eff}}(\gamma-\gamma_{\rm eff})<0.
\label{eq:stability}
\end{equation}
We first note that in the absence of cooling $\gamma_{\rm eff}=\gamma$ and
the stability condition reduces to:
\begin{equation}
\gamma>2-\frac{2}{n}.
\end{equation}
For the spherical case ($n=3$) the stability criterion is
$\gamma>\frac{4}{3},$ recovering a well known result. For filamentary
accretion ($n=2$) of adiabatically collapsing gas stability is gained when
$\gamma>1$, and for planar collapse $(n=1)$ when $\gamma>0.$

When cooling is present the stability criterion is:
\begin{equation}
\gamma_{\rm eff}>\frac{2n\gamma}{n\gamma+2},
\end{equation}
which, for monoatomic gas ($\gamma=\frac{5}{3}$) is
$\frac{10}{7},~\frac{5}{4},~\frac{10}{11}$ for spherical, cylindrical
and planar collapse respectively.

We note that for $\gamma=\frac{5}{3}$ the critical value of
$\gamma_{\rm eff}$ in the presence of cooling is always somewhat larger
than the adiabatic critical value. Hence, monoatomic gas that is
cooling with some local $\gamma_{\rm eff},$ (with entropy and energy decreasing due to cooling) is always less stable
than an adiabatic gas with a softened EoS
$\gamma=\gamma_{\rm eff}$.

\section{Stability of cosmological filaments}
\label {sec:selfsimilar}

\begin{figure}
\begin{center}
\includegraphics[width=3.3in]{fig1.eps}
\caption{\label{fig:massnorm} Mass profiles of self-similar
  cylindrical accretion normalized to the virial radius and virial
  mass, as defined by the first shell crossing, for various initial
  perturbation power-laws coefficients $\epsilon$.}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=3.3in]{fig2.eps}
\caption{\label{fig:ml} Self-similar DM
  trajectory and mass profile of an $\epsilon=0.8$ perturbation. At
  time $t_*$ a specific shell with mass $M_*$ and radius $x_*$ is at
  turnaround. The unitless time is $\tau=t/t_*$ and unitless radius is
  $\lambda=x/x_*$\refereebf{ (for consistency, all variables are named according to FG84)}. {\it Red (left axis):} the self-similar shell
  trajectory with $\lambda$ normalized to the current turnaround
  radius at each time $\Lambda(\tau).$ {\it Blue (right axis):} the
  self-similar mass profile in units of $M_*.$ The black vertical line
  corresponds to the radius at which shell crossing first occurs, and
  is evident in the mass profile as well as in the trajectory.
}
\end{center}
\end{figure}

Filaments grow (in girth) by accreting gas from their
surroundings. The infall geometry is primarily cylindrical, although
most of the gas is channeled along sheets, and becomes spherical in
the vicinity of halos. The flow parallel to the filament below and
above the shock is continuous, and can be factored out locally with a
proper shift of the frame of reference. Dark matter (DM) and gas
accrete onto filaments together. Separation between gas and DM occurs
when gas becomes thermalized and is decelerated due to its pressure.
We wish to determine where and how this gas is thermalized,
particularly in the presence of cooling, which softens the effective
EoS of the gas. To do so, we must connect global properties of
filaments: their mass per unit length and their density and velocity
profiles, to local cooling and contraction rates that determine
the stability of the gas.

An ideal framework for connecting the large-scale properties of
filaments to local conditions is the self-similar solutions of
FG84.  We do not present here the derivation and results, and refer
the reader to the original paper. Following FG84, we numerically solve
for the self-consistent density profile and trajectories of infalling cylindrical
DM shells, starting from an initial mass
perturbation within an Einstein de-Sitter universe. In this framework,
a filament is characterized by its mass per unit length, and by its
initial perturbation:
\begin{align}
M_{fil}(r)&=M_{fil}^0(r)\left(1+\frac{\delta M_{fil}(r)}{M_{fil}^0(r)}\right)\nonumber\\
M_{fil}^0(r)&\equiv\pi r^2\rho_u,
\end{align}
with $\rho_u$ the universal density at the initial time. The perturbation is defined as a function of the unperturbed mass:
\begin{equation}
\frac{\delta M_{fil}}{M_{fil}^0}=\left(\frac{M_{fil}^0}{M_*}\right)^{-\epsilon},
\end{equation}
with $M_*$ a reference mass.
 $\epsilon$ varies between $0$ and $1$, where $1$ corresponds to the
most localized perturbation that still grows with $M_{fil}^0$, and $0$
to a long range perturbation for which the density of the perturbation
is still decreasing with $M_{fil}^0.$ \fig{massnorm} shows the
resulting mass profiles (normalized to the virial radius and the
virial mass) for various values of $\epsilon.$ It is evident that the
profile within the virial radius of the filament depend weakly on
$\epsilon.$ The visible discontinuities in gradient correspond to
caustics in DM shells as they turn around consecutively. The outermost
caustic, or ``first shell crossing'' is defined as the virial
radius. The mass profiles in \fig{massnorm} are normalized to that
radius.

\fig{ml} shows the self-similar trajectory of a DM shell and the
filament's mass profile. The trajectory's radius is normalized to the
turnaround radius at each time, so before turnaround, at $\tau<1,$ the
spatial coordinate $\lambda/\Lambda>1.$ The mass profile as a function
of that same spatial variable is also present. The vertical black line
corresponds to the event of first crossing (the virial
radius), and is roughly where the virial shock will occur for a
gaseous shell.

From the self-similar solution we extract the infall velocity and
density at every radius. The infalling velocity and density are for
DM trajectories and correspond to infalling gaseous shells only before
they pass through the virial shock, at which point their pressure
becomes significant and their trajectories diverge from those of the
DM. The stability of the post-shock gas depends on the local compression
rate, density and temperature for the {\it post-}shocked gas. These
values are approximated from the {\it pre}-shocked ones by use of the
strong shock approximation, that is valid as long as the pre-shocked
velocity is much larger than the pre-shocked speed of sound
($c_{s}^0$). Assuming that the $c_s^0\lesssim 10 \mathrm{km\,s^{-1}},$ we show later
  (\cref{fig:gmx_0.2,fig:gmx_0.99}), that this
  approximation goes from being marginally satisfied for
  the smallest filaments to being fully justified for the large
filaments. Using the full shock conditions will not change the results
significantly, and requires knowledge of the thermodynamic state of
the cold gas. As the effective EoS becomes softer, the virial shock
ceases to expand, and starts to collapse. At its
critical state, we expect the shock to be at rest. Using the strong
shock approximation and assuming the shock is at rest, the post-shock
values are:
\begin{align}
&\rho_1=\frac{\gamma+1}{\gamma-1}\rho_0\\
&u_1=\frac{\gamma-1}{\gamma+1}u_0\nonumber\\
&e_1=\frac{2u_0^2}{(\gamma+1)^2}\nonumber\\
&T_1=\frac{\mu(\gamma-1)}{N_ak_B}e_1,\nonumber
\end{align}
with subscript $0$ and $1$ denoting pre-shocked and post-shocked
variables respectively, $\rho_0,u_0$ given from the numerical solution
of the self-similar collapse, and $e_1,T_1$ are the internal energy
and temperature of the post-shock variables. $\mu=0.61$ is the mean
molecular weight for primordial, fully ionized gas, and $N_a,k_B$ are
Avogadro's number and Boltzmann's constant. The
conversion from DM density to gas density is achieved by multiplying
the density by a universal baryonic fraction ($f_b=0.17$ throughout
this work). This value is reasonable as long as the DM and gas flow
together, i.e.\@ for pre-shocked gas.

\begin{figure}
\begin{center}
\includegraphics[width=3.3in]{fig3.eps}
\caption{\label{fig:gmx_0.2} Stability of $\epsilon=0.2$ perturbation
  as a function of filament mass and radius with respect to the
  filamentary virial radius. {\it colormap:} $\gamma_{\rm eff}$ - white color corresponds to $\gamma_{eff}=\gamma=5/3.$ All
  values below $1.25$ are unstable. {\it contours:} the infall
  velocity.}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=3.3in]{fig4.eps}
\caption{\label{fig:gmx_0.99} Same as \cref{fig:gmx_0.2}, but for $\epsilon=0.99$}
\end{center}
\end{figure}

For a filament characterized by $M$ and $\epsilon$ we calculate
the post-shock values for every radius below $r_{\rm vir},$ and use
\equ{rhodot} and \equ{lambda} inserted into \equ{geff2} to calculate
$\gamma_{\rm eff}$.  \cref{fig:gmx_0.2,fig:gmx_0.99}
show $\gamma_{\rm eff}$ as a function of $M_{fil}$ and $r/R_{\rm
  vir},$ with $R_{\rm vir}$ defined as the radius of first shell crossing of the self-similar solutions (see \cref{fig:ml}). In regions where $\gamma_{\rm eff}$ drops below the threshold
for stability, $\gamma_{\rm eff}<\gamma_{crit}=1.25$, the filament
cannot sustain a virial shock. Moreover, if a shock were to form in regions
where $\gamma_{\rm eff}<0,$ the post-shock pressure would decrease
even as it contracts. The mass range for which the filament
is unstable at least at some radius for both values of $\epsilon$ is
between $10^{11}-10^{13} \ensuremath{M_{\odot}} \, \ensuremath{{\rm Mpc}} ^{-1}.$ In these plots, radii of density caustics leave horizontal features, and peaks in the cooling curve create features parallel to infall velocities. For each radius, the sharp drop in $\gamma_{eff}$ as mass exceeds the lower threshold for instability is due to the post-shock temperature exceeding $10^4K,$ where the cooling rate grows by many orders of magnitude.

\section{Implication for accretion onto halos}
\label {sec:halos}
\begin{figure}
\begin{center}
\includegraphics[width=3.3in]{fig5.eps}
\caption{\label{fig:stab_0.2} Stability of filaments falling into
  halos as a function of halo mass and redshift, for $\epsilon=0.2.$ {\it colormap:} the stable filament fraction (see text).}
\end{center}
\end{figure}
We now wish to relate the filament masses (per unit length) to the typical masses of halos fed by such filaments. A
full analysis of the filament distribution that accrete onto certain
halos requires cosmological N-body simulations. For simplicity,
we choose an alternative avenue that approximates the relation between
halos and their filaments. \citet{danovich12} study the filamentary
nature of mass accretion onto high redshift galaxies. They find that
typically, halos that originate from high-$\sigma$ peaks in the
initial perturbation accrete most ($f_{acc}\simeq70\%$) of their mass
in filaments, out of which, $95\%$ originates from the combined flow
in the $3$ largest filaments. Although the work analyzes high redshift
galaxies, we expect the same to be true for low-redshift clusters, who are also high-$\sigma$ peaks. Using these values, we estimate the
typical accretion rate through each filament as
\begin{equation}
\dot{M}_{acc}=\frac{f_{acc}}{3}\dot{M}_{halo},
\end{equation}
with $\dot{M}_{acc},\dot{M}_{halo}$ the mass flow rate of gas within a
filament and halo accretion rate
respectively. For the halo accretion rate we use the fit from
\citet{neistein06},
\begin{equation}
\frac{\dot{M}_{halo}}{M_{halo}}= \left(\frac{0.03}{Gyr} \right) \left(1+z\right)^{2.5},
\end{equation}
with $z$ the redshift. Finally, the mass of the filament per unit
length is related to the flow rate by
$M_{fil}=\dot{M}_{acc}/v_{\rm vir}(M_{halo},z)$, using the halo virial
velocity, $v_{\rm vir},$ as an estimate for filaments' velocity as they
accrete onto halos. A typical filament feeding a high-$\sigma$ peak halo of mass
$M_{halo}$ will thus have an estimated mass of:
\begin{equation}
M_{fil}\simeq \frac{f_{acc}}{3} \left(\frac{0.03}{Gyr} \right) \left(1+z\right)^{2.5}\frac{M_{halo}}{v_{\rm vir}}.
\end{equation}

In \fig{stab_0.2} we use the inverse of this transformation to show the stable fraction of a filament, $f_{stable},$ as a function of the halo mass and redshift:
\begin{equation}
f_{stable}\equiv \frac{1}{R_{\rm vir}}\int_0^{R_{\rm vir}}\Theta\left[\gamma_{\rm eff}(r)-1.25\right]dr,
\end{equation}
with $R_{\rm vir}$ the virial radius of the filament, and $\Theta$ the Heaviside function.
From \cref{fig:stab_0.2} it is evident that for $10^{10}M_{\odot}\lesssim M_{halo}\lesssim 10^{13}M_{\odot}$ at $z=3$, as well as for $10^{12}M_{\odot}\lesssim M_{halo}\lesssim 10^{15}M_{\odot}$ at $z=0$ halos are expected to be fed by filaments that are not in hydrostatic stability.

\section{Summary and Discussion}
\label {sec:summary}
We have shown that in the presence of significant cooling, the
accretion process of gas onto cosmic-web structures will not always
proceed according to the standard virialization scenario of the
infall-heating-cooling sequence.
The analysis shown
here can be applied for accretion onto spherical halos, cylindrical filaments and planar sheets.
For filament of $10^{11}-10^{13} \ensuremath{M_{\odot}} \, \ensuremath{{\rm Mpc}} ^{-1}$, we show that
gas is expected to fall without ever passing a shock, resulting in
dense, thin filaments with low entropy. This is in complete analogy
to spherical cold accretion onto halos that have been shown in BD03
and demonstrated in observations and simulations.

Using a simplified toy model for the relation between halo mass and
redshift to typical filaments that feed it, we show that throughout
cosmic history galaxies and clusters are affected by that
instability. In particular, high-z star forming galaxies
($M_{halo}>10^{10} \ensuremath{M_{\odot}} $ at $z=3$), and low redshift groups and
clusters ($10^{12} \ensuremath{M_{\odot}} \lesssim M_{halo}\lesssim 10^{15} \ensuremath{M_{\odot}} $ at
$z=0$) will be fed by filaments for which the gas is unstable.

The process that eventually stops the infall is still
unclear, and we postulate that it is either angular momentum support
from an original helicity of the filament, or by reduction of the
cooling rate due to self-shielding of the gas.  A prediction of this
work is thus that filament gas, in the non-stable regime, will be highly
rotating and angular momentum supported. Both processes are hard to identify
in simulations, and have not been examined so far. In their absence,
gas in simulations will flow towards the center of the filament until
it approaches the numerical gravitational smoothing length, at
which point the force will diminish. This
indicates that the density and entropy of gas in unstable filaments are a numerical artifact and will not converge to the right values. This
problem will be examined in future work.

\refereebf{The lack of virialized gas in filaments is expected to significantly affect the outcome of galaxies falling onto the filament, and of halos fed by the filament. Halos falling onto filaments are expected to looe gas through ram pressure stripping, and to enrich the filament with metals. Both these processes will be suppressed when galaxies fall into filaments with no stable atmosphere. Penetration of cosmic-web gas directly to galaxies affects the ISM state, and the gas available for star formation and AGNs, as well as their feedback efficiencies. \citet{mandelker16} analyze the Kelvin-Helmholtz stability of supersonic filaments. They find that filaments lose stability via bulk modes, that correspond to standing waves reflecting through edges of the filament. These results do not account for the effects of gravitational attraction towards the center of the filament, and to angular momentum support, both expected to stabilize the filament further. These effects will be addressed in future work.

Observationally, the temperature of the filament could affect its detectability through Lyman$-\alpha$ absorption \citep{narayanan10,wakker15} and emission \citep{martin15}. The temperature of the filaments will also affect the soft X-ray background and the total amount of gas in the ``warm phase'' \citep{cen99,dave01}. All these effects are left for analysis in future work.}

\acknowledgments We thank Oliver Hahn for useful discussions. YB and DP have been supported by ISF grant 1059/14.


\section*{Abstract}{\small
Recent $\gamma$-ray observations show that middle aged supernova remnants (SNRs) interacting with molecular clouds (MCs) can be sources of both GeV and TeV emission. Based on the MC association, two scenarios have been proposed to explain the observed $\gamma$-ray
emission. In one, energetic cosmic ray (CR) particles escape from the SNR
and then illuminate nearby MCs, producing $\gamma$-ray emission, while the other involves direct
interaction between the SNR and MC. In the direct interaction scenario, re-acceleration of pre-existing CRs in the ambient medium is investigated while particles injected from the thermal pool are neglected in view of the slow shock speeds in middle aged SNRs. However, standard diffusive shock acceleration (DSA) theory produces a steady state particle spectrum that is too flat compared to observations, which suggests that the high energy part of the observed spectrum has not yet reached a steady state. We derive a time dependent DSA solution in the test particle limit for re-acceleration of pre-existing CRs case and show that it is capable of reproducing the observed $\gamma$-ray emission in SNRs like  IC 443 and W44, in the context of a MC interaction model. We also provide a simple physical picture to understand the time dependent DSA spectrum. A spatially averaged diffusion coefficient around the SNR can be estimated through fitting the $\gamma$-ray spectrum. The spatially averaged diffusion coefficient in middle aged SNRs like IC 443 and W44 is estimated to be $\sim 10^{25}\rm cm^2s^{-1}$ at $\sim \rm 1GeV$, which is between the Bohm limit and interstellar value.

\vspace{10mm}
\normalsize}
\end{minipage}

\section{Introduction}
Recent observations from both space-based GeV observatories and ground-based TeV observatories reveal $\gamma$-ray emission from several middle aged supernova remnants (SNRs) which are interacting with molecular clouds (MCs), e.g. W44 \citep{Abdo10a,Giuliani11,Uchiyama12}, IC443 \citep{Albert07,Acciari09,Abdo10b}, W28 \citep{Aharonian08,Abdo10c,Hanabata14} and W51C \citep{Abdo09,Aleksic12}. The observed $\gamma$-ray emission from these objects shows a smooth transition from the GeV band to the TeV band (if detected) and they all peak in the GeV band. The characteristic $\pi^0$-decay signature identified in IC 443 and W44 \citep{Giuliani11,Ackermann13} provides possible direct evidence for cosmic ray (CR) proton acceleration in SNRs, making middle aged SNRs detected in $\gamma$-rays an important class of objects for understanding  CR acceleration in SNRs.

So far two scenarios have been proposed to explain the observed spatial correlation between the $\gamma$-ray emission region and the MC region with hadronic origin: one is the escaping scenario \citep{Gabici09,Fujita09,LC10,Ohira11}, which focuses on the CR particles that escaped from the remnant, and the other one is the direct interaction scenario \citep{bykov00,Uchiyama10,Inoue10,TC14,TC15,Cardillo16}, which studies the CR particles accelerated at the remnant. In the escaping scenario, CR particles having escaped from the remnant together with the pre-existing CR in the ambient medium illuminate the nearby MCs, producing the $\gamma$-ray emission. A collision between the SNR and MC is not necessarily required and the resulting  $\gamma$-ray emission can be external to the remnant. A good example for this scenario might be W28. In the direct interaction scenario, a collision between the remnant and MC occurs, as indicated by both observation and theory. Shocked clouds are identified in middle aged SNRs and are spatially coincident with the $\gamma$-ray emission region  \citep[e.g.,][]{Abdo10a,Abdo10b,uchiyama11,Nicholas12}. Besides, the MC interaction creates a region with high density which then becomes an ideal site for $\pi^0$-decay emission.

Diffusive shock acceleration (DSA) is believed to be the particle acceleration mechanism in most astrophysical environments involving shock waves \citep[e.g.,][]{Bell78,B&E87}. The theory naturally produces a power law spectrum of energetic particles in a steady state which is close to the observed CR spectrum after taking into account  propagation effects. DSA is considered to be the most plausible particle acceleration mechanism at the SNR shock front for both scenarios.

Here, we limit our discussion to the direct interaction scenario. In this scenario, re-acceleration of pre-existing CRs is taken into account, while particle injection through the thermal pool is neglected in view of the slow radiative shock \citep{B&C82,Uchiyama10,TC14}. If this picture is correct, then it implies a transition of seed particles for CR acceleration in SNRs from  injected seed particles to pre-existing CRs.  When does this transition occur and how would it affect the accelerated CR spectrum in SNRs could be interesting problems for future study.

The standard DSA theory produces too flat a steady state particle spectrum compared to that indicated by observations, so it has been suggested that the observed particle spectrum in the energy range of interest has not reached a steady state yet. An exponential cutoff has been explored in \cite{Uchiyama10} and \cite{TC14} as the high energy cutoff for  MC interaction models. It is found that for remnants like IC 443, W28 and W51C with TeV detection, the exponential cutoff tends to underproduce the TeV emission. In \cite{TC15}, we developed a time dependent DSA solution in the test particle limit for the case of re-acceleration of pre-existing CRs and show that the time dependent solution in combination with the MC interaction model in \cite{TC14} is capable of explaining the observed $\gamma $-ray emission in remnants like W44 and IC 443 from GeV to TeV bands.

In the following section, we first discuss  the properties of the time dependent DSA solution  and then provide a simple physical picture to understand the spectral shape of the time dependent DSA solution.

\section{Time dependent DSA solution}
Assume that the pre-existing CRs around the SNR have the same spectrum as that observed locally,
which shows a broken power law. Based on calculations in \cite{TC15} for the time dependent DSA spectrum, there is a critical momentum below which the spectral shape is consistent with the steady state spectrum, and above which the spectral shape simply follows the input CR spectrum with possible hardening due to energy dependent diffusion (see figure 1-4 in \cite{TC15} for details).

In the Appendix of \cite{TC15}, we provided a simple phenomenological way to understand the  time dependent DSA spectrum for energy independent diffusion. Here, we discuss the physical picture in more detail and then discuss  possible implications of such a picture. In DSA, energetic particles bounce back and forth around the shock front due to self generated magnetic turbulence \citep[e.g.][]{Bell78}. Every time a particle moves across the shock discontinuity it gains some energy and every time the particle enters the downstream region it has chance to escape from the DSA site. According to \cite{Drury83}, the momentum gain rate for a particle undergoing DSA is
\begin{equation}
\frac{dp}{dt}=\frac{U_1-U_2}{3}\, \left( \frac{\kappa_1}{U_1} +\frac{\kappa_2}{U_2}\right)^{-1}  p,
\label{energy_gain}
\end{equation}
where $U$ is the flow velocity, $\kappa$ is the diffusion coefficient and subscripts 1 and 2 are for physical quantities in the upstream and downstream regions, respectively.
The probability for a particle to leave the acceleration site at downstream region is given by $P=4U_2/v$ \citep{Drury83}.

For a particle with initial momentum $p_i$, after $n$ cycles of acceleration the particle momentum becomes (See e.g. \cite{Drury83} for detail)
\begin{equation}
p_n\sim \prod^n_{k=1}\left[1+\frac{4(U_1-U_2)}{3v_k}\right]p_i,
\end{equation}
leading to
\begin{equation}
 {\rm ln}(p_n/p_i)\sim\frac{4(U_1-U_2)}{3}\sum^n_{k=1}\frac{1}{v_k}.
\end{equation}
The probability for a particle to stay at the acceleration site after $n$ cycles of back and forth motion is
\begin{equation}
P_n\sim\prod^n_{k=1} \left( 1-\frac{4U_2}{v_k}\right),
\end{equation}
leading to
\begin{equation}
\ln P_n  \sim  -4U_2\sum^n_{k=1}\frac{1}{v_k}
	\approx -\frac{3U_2}{U_1-U_2} {\rm ln}(p_n/p_i).
	\label{escape_probability}
\end{equation}
For a strong shock $P_n\approx p_i/p_n$.

Next consider the acceleration of a group of particles with momentum spectrum $N(p)$. Based on the conservation of particle number, the new particle spectrum $\overline{N}(p)$ is related to the input particle spectrum $N(p)$ by
 \begin{equation}
\overline{N}(p_n)dp_n=N(p_i)P_ndp_i,
\label{conservation}
\end{equation}
where $p_i$ and $p_n$ are the particle momentum before and after the DSA respectively. $P_n$ is the probability for a particle to stay at the DSA site after $n$ cycles of acceleration and equals  $p_i/p_n$ according to eq.\  \ref{escape_probability}. For constant $U$ and $\kappa$, the energy gain rate $dp/dt$ given by eq.\  \ref{energy_gain} is proportional to the particle momentum $p$, which implies that $dp_n/dp_i=p_n/p_i$. Now eq.\  \ref{conservation} becomes $\overline{N}(p_n)p_n^2\approx N(p_i)p_i^2$, which means in the particular $\log[N(p)p^2]-\log (p)$ plane the whole DSA process simply works like a horizontal shift in the $\log(p)$ axis. The amount of shift is determined by
\begin{equation}
{\rm ln}\left( \frac{p_n}{p_i}\right)=\frac{U_1-U_2}{3}\, \left( \frac{\kappa_1}{U_1} +\frac{\kappa_2}{U_2}\right)^{-1}  t,
\label{shift}
\end{equation}
which is linearly proportional to $t$ in the $\log(p)$ axis.

In the case of re-acceleration of pre-existing CRs in a SNR, energetic particles with the same input CR spectrum are being injected into the system all the time. Particles injected at different times will be shifted by different amounts in the $\log(p)$ axis due to different acceleration times, as shown in Fig \ref{schematic}. As a result, the accumulated particle spectrum at the shock front is determined by the sum of a group of input CR spectra with different amounts of shift. Since the amount of shift in $\log(p)$ axis is linearly proportional to $t$ according to eq.\  \ref{shift}, all the shifted spectra have the same weight for the sum in $\log[N(p)p^2]-\log (p)$ plane.
Now we can understand the time dependent DSA spectrum presented in \cite{TC15}. Because our input CR spectrum follows roughly a broken power law, the accumulated CR spectrum at the shock front  contains three parts according to the above discussion. In the low energy and high energy parts the accumulated particle spectrum maintains the shallow and steep power law in the input CR spectrum, respectively, as the shifted spectra all share the same power law index. At intermediate energies the accumulated particle spectrum shows a plateau which is mainly due to the break in the input CR spectrum.

\section{Discussion}
If the observed $\gamma$-ray emission from those middle aged SNRs interacting with MCs is the result of time dependent DSA,  according to eq.\  \ref{shift} the extension of the plateau region in the $\gamma$-ray spectrum reflects the logarithmic ratio $\ln(p_f/p_i)$ between the initial and final momentum. Thus it can be used to estimate the spatially averaged diffusion coefficient around the SNRs. In SNRs, the spatially averaged diffusion coefficient in the upstream region $\kappa_1$ is much larger than that in the downstream region $\kappa_2$, as $\kappa_1$ gradually approaches the interstellar value when particles move far away from the shock front. Hence, according to eq.\ \ref{shift} we have
\begin{equation}
\kappa_1 \approx\frac{(U_1-U_2)U_1t}{3\ln (p_f/p_i)}\approx 10^{25}{\rm cm^2s^{-1}}\frac{r-1}{r\ln(p_f/p_i)}\left(\frac{U_1}{100\rm ~km~s^{-1}}\right)^2\left(\frac{t}{10^4~ \rm yrs}\right),
\end{equation}
where $r$ is the shock compression ratio and $r=4$ for strong shocks. In middle aged SNRs $\ln(p_f/p_i)$ is of order  unity based on the observed $\gamma$-ray spectrum. Hence the spatially averaged diffusion coefficient in these middle aged SNRs at the energy range around the plateau region ($\sim \rm GeV$) is approximately $10^{25}\rm cm^2s^{-1}$. The number is close to the value derived in \cite{TC15} through detailed fitting and is found to be between the Bohm limit and the interstellar value.
\begin{figure}[htb]
\center
\includegraphics[width=0.8\textwidth]{schematic.pdf}
\caption{Schematic figure for the time dependent DSA solution.}
\label{schematic}
\end{figure}

\small
\section*{Acknowledgments}
We are very grateful to P. Boumis and the meeting organizers for arranging a very stimulating conference.