\section{Introduction}
\label{sec:intro}
In the past decade, a novel class of visual sensors, known as Dynamic Vision Sensors (DVS), has attracted a lot of attention from the academia as well as the industry. This is fueled by its exceptional performance compared to conventional image sensors in terms of temporal resolution, latency, dynamic range and power consumption. Such breakthrough is achieved via a fundamentally different principle of operation. Conventional image sensors output a synchronous sequence of intensity frames, whereas the DVS outputs a sequence of events asynchronously, whereby each event indicates a change in intensity at a particular pixel location and time instant.
\par

\begin{figure}[t]
    \centering
        \begin{subfigure}{0.48\linewidth}
    		\centering
			\includegraphics[width=1\linewidth]{{figures/aps_slider_depth_49}}
			\caption{}
			\label{bei:fig:slider_depth_aps}
		\end{subfigure}
		\begin{subfigure}{0.48\linewidth}
            \centering
            \includegraphics[width=1\linewidth]{{figures/canny_0.15_0.5_slider_depth_49}}
            \caption{}
            \label{bei:fig:slider_depth_canny}
		\end{subfigure}
		\begin{subfigure}{0.48\linewidth}
            \centering
            \includegraphics[width=1\linewidth]{{figures/sofea_bei_slider_depth_49}}
            \caption{}
            \label{bei:fig:slider_depth_sofea_bei}
		\end{subfigure}
	    \begin{subfigure}{0.48\linewidth}
            \centering
            \includegraphics[width=1\linewidth]{{figures/1895_sPixel.png}}
            \caption{}
            \label{bei:fig:slider_depth_mueggler2015}
	    \end{subfigure}
    \caption{(a). An event camera is moved on a slider in front of objects at different depth \cite{Mueggler2017}; (b) Edge map obtained via Canny edge detection is visualized for reference; (c) Scene-adaptive optical flow estimation is subsequently employed for lifetime estimation to produce crisp event stream images; (d) Superevents extracted using a fully convolutional network on the event stream image (best viewed on monitor).}
    \label{fig:intro_fig}
\end{figure}

Due to the above-mentioned unique advantages, event cameras have been trialed for a wide range of high-level vision tasks, namely object segmentation \cite{evseg2019}, motion segmentation \cite{Mitrokhin_2020_CVPR}, pattern recognition \cite{Ramesh2019, Ramesh2018}, gesture recognition \cite{eddd2020, chen2020}, space situational awareness \cite{bagchi2020event, chin2019star}, to name a few. Despite these advances in addressing these problems using conventional or deep learning methods that have made use of novel neural architecture design and training schemes suited for event-based vision, previous works have ignored mid-level representation that are critical in bridging gap between the sparse, disjoint information in the event stream and downstream visual tasks. As these novel cameras respond to structure and texture of the scene while not capturing intensity information, it is challenging to obtain groupings of events, i.e. mid-level representations, as done in superpixel extraction for image data \cite{stutz2018superpixels}.
\par
In computer vision, superpixels are perceptually similar pixels grouped together to effectively reduce the number of image primitives for subsequent processing. This mid-level representation has been widely adopted for scene understanding\cite{wang2019revisiting} via detecting and segmenting salient regions \cite{zhao2019pyramid, zhao2015saliency, mostajabi2015feedforward}, and also for high-level tasks such as object detection \cite{li2018contour, shu2013improving} and tracking \cite{feng2019dynamic, yeo2017superpixel}. However, these ideas are yet to be explored in the event-based domain for mid-level feature extraction. One main reason is possibly that, the standard convolution operation in convolutional neural networks (CNNs) is defined on a regular image grid with intensity information. While a few attempts have been made to modify deep architectures \cite{asynconvet2018, asyncconvnet2020}, performing convolution over aysnchronous data remains challenging.
\par
In this paper, we introduce a fully convolutional network (FCN) based superevents extraction algorithm that is in the spirit of \cite{yang2020sup}. The input to the FCN is a Binary Event Image (BEI) obtained via lifetime augmentation \cite{Mueggler2015}. In this regard, we make use of the open-source event-based flow estimation method SOFEA \cite{Low_2020_CVPR_Workshops}. It is worth noting that the recent works targeting super-resolution of the event stream\cite{mohd2020eventsr,wang2020eventsr} is very different to the task considered in this paper. Super-resolution aims to reconstruct high resolution representations directly from the event stream and is unrelated to superevents (or superpixels) mid-level feature extraction.
\par
The primary contributions of this paper are:
\begin{enumerate}
\item A scene-adaptive noise rejection method for event-based optical flow estimation (\autoref{bei:subsec:application}) and a robust lifetime augmentation method for accurate BEI estimation (\autoref{bei:sec:lifetime}).
\item A scene-adaptive BEI rendering time interval which approximately normalizes for the expected apparent displacement magnitude of features with apparent motion (\autoref{bei:sec:render}).
\item A superevents extraction pipeline for event cameras leveraging on the above contributions using a fully convolutional network.
\end{enumerate}

\section{Related Work}
\label{sec:related}
As previously mentioned, the output of a DVS is a stream of asynchronous events, as opposed to a set of images captured at a constant rate by a standard camera. Due to this unconventional output, the general design strategy of event-based algorithms has been to either implicitly or explicitly buffer a certain number of past events to perform a given task. The most straight-forward method is to accumulate events over a fixed time interval. For instance, \cite{Authors, Authorsbmvc19, Kogler2009} used contrived time intervals of $5-66 \mathit{ms}$ to estimate binary event images for object tracking and stereo vision. Nevertheless, the motion-dependent sensing aspect of the DVS is a major hindrance for choosing the optimal time interval to prevent motion blur or a low density event image.
\par
Another way to obtain an event stream representation is to implicitly choose a time interval by accumulating a fixed number of past events \cite{Ramesh2019, gallego2016event, Rebecq2016}. However, the critical issue is choosing the number of events, which unfortunately depends on the scene complexity and spatial resolution of the DVS itself. \cite{Liu2018} proposed a method to accumulate past events until the density of events in a local spatial neighborhood reaches a predefined threshold. It adaptively selects the number of events to buffer by considering both scene complexity and spatial resolution. With that said, all the methods above accumulate all events in a particular time interval without correcting for apparent motion.
\par
A more principled approach is known as the ``motion correction''. In \cite{Gallego2018, Zhu2017}, uniform optical flow is first assumed for all events within a small spatio-temporal neighborhood. Using this assumption, an event image is then estimated by warping all events in the spatio-temporal neighborhood to a given reference time, according to the set of optical flow parameters that best describe the apparent motion of events. The optimal set of parameters is searched using iterative methods, therefore it is more computationally expensive compared to all the above methods. The uniform optical flow assumption is replaced with a spatially uniform Bézier curve motion assumption in \cite{Seok2020}. In general, a relatively large time interval (or a large number of events) is required for accurate estimation of the relative intensity change magnitude, due to inherent noise and latency of the DVS.
\par
Due to the limitations of the above event stream representation strategies, this paper adopts the concept of ``lifetime'' introduced by \cite{Mueggler2015} to represent how long will it take for an event to trigger a new associated event in its neighboring pixel location, due to persistent apparent motion. Hence, by activating pixels of the binary event image according to the estimated lifetime and position of events, a binary event image can be defined for any time instant in the past. Lifetime estimation itself can be done using the so-called ``Local Plane Fit''' \cite{Benosman2014}, which is one of the most widely adopted approaches for optical flow estimation. Normal flow is estimated by fitting a plane in a local spatio-temporal neighborhood of the Surface of Active Events (SAE), which describes the timestamp of the most recent event for each pixel location, about the incoming event. The single-shot optical flow estimation algorithm (SOFEA) \cite{Low_2020_CVPR_Workshops} considers selected events in the local spatio-temporal neighborhood with the same polarity, and thus is heavily robust to noise.
\par
Nonetheless, the limitations of SOFEA make it rather non-ideal for direct BEI estimation. One of the main issues is that SOFEA suffers from low event density in select scenes, mainly due to refractory filtering, which implies that a significant amount of events are lost in the process. In other words, the refractory period prevents updating the SAE if an event has arrived within the refractory period. This is undesirable as the estimated BEI (similar to an edge map) is then not representative of the scene. The use of a refractory filter also limits the maximum speed of scene dynamics. Apart from that, the validity of the local planar assumption is dependent on the scene dynamics and the local scene structure. Thus, events associated to high acceleration motion or moving edges with large curvatures (e.g. a corner) will most likely be regarded as noise by the goodness-of-fit noise rejection method, as it is evaluated with respect to a planar model. Thus, we propose a pipeline for BEI estimation and lifetime augmentation, termed as SOFEA-BEI, to deal with the shortcomings of vanilla SOFEA.

\section{SOFEA-BEI}

\label{sec:bei}
An event $\boldsymbol{e}$ can be characterized by a 4-tuple $(x, y, t, p)$, where $x, y$ represents the spatial location or position of the event in pixel coordinates, $t$ represents the time of occurrence of the event and $p \in \{  -1, +1 \}$ represents the polarity of the event. On top of that, the event stream can be defined as a sequence of events $(\boldsymbol{e}_n)_{n \in \mathbb{N}}$ such that $i<j \implies t_i \leq t_j$, where $\boldsymbol{e}_i = (x_i, y_i, t_i, p_i), \boldsymbol{e}_j = (x_j, y_j, t_j, p_j)$.

\par
Disregarding event polarity, each event carries information in the three-dimensional $xyt$ spacetime domain. Consequently, the sequence of most recent events over the entire spatial domain, which is a subset of the event stream $(\boldsymbol{e}_n)_{n \in \mathbb{N}}$, can be simply described by a spatio-temporal surface $\Sigma_e$. The spatio-temporal surface $\Sigma_e$ is better known as the Surface of Active Events (SAE), as coined in \cite{Benosman2014}. The mathematical definition of the SAE is given by:
\begin{equation}
	\begin{split}
		\Sigma_e: \mathbb{R}^2 &\to\mathbb{R} \\
		\boldsymbol{p} &\mapsto \Sigma_e( \boldsymbol{p} ) = t
	\end{split}
\end{equation}
where $\boldsymbol{p} = \begin{bmatrix}x & y\end{bmatrix}^\top$.

\par
A simple expression for $\Sigma_e$ can be obtained by approximating the SAE as a plane in a local spatial neighborhood. With the approximation above, the local planar SAE can be sufficiently described by a set of plane parameters $\boldsymbol{\Pi} = \begin{bmatrix} a & b & c & d \end{bmatrix}^\top$. An interesting property of the SAE is that its spatial gradient $\nabla \Sigma_e(\boldsymbol{p}) = \begin{bmatrix} \frac{\partial \Sigma_e(\boldsymbol{p})}{\partial x} & \frac{\partial \Sigma_e(\boldsymbol{p})}{\partial y} \end{bmatrix}^\top$ encodes information of normal flow associated to the active events. For a local planar SAE, $\nabla \Sigma_e(\boldsymbol{p})$ can be further defined as $\begin{bmatrix} -\frac{a}{c} & -\frac{b}{c} \end{bmatrix}^\top$. SOFEA improved upon existing local plane fitting methods by using an efficient greedy algorithm to optimally select the spatially neighboring events associated to the incoming event for accurate non-iterative local SAE plane fitting. SOFEA also incorporated a mathematically simpler description of the local planar SAE given by the additional constraint that the incoming event lies on it.
\par
The refractory filter originally incorporated in SOFEA is removed from SOFEA-BEI so as to retain a high event density, and also to eliminate its introduced limit on the speed of scene dynamics. We then post-process the lower accuracy optical flow estimates to obtain highly satisfactory lifetimes for the purpose of BEI rendering. In particular, post-processing is performed via extended lifetime augmentation and lifetime reset. Apart from that, the goodness-of-fit noise rejection method in SOFEA is replaced with a novel scene-adaptive noise rejection method. In order to reduce visual redundancy between successive BEI renders, a motion-normalized BEI rendering time interval is also introduced. \autoref{bei:fig:sofea_bei_pipeline} illustrates the overall pipeline of SOFEA-BEI.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.55\linewidth]{figures/sofea_bei_pipeline}
    \caption{Pipeline of SOFEA for Binary Event Image Estimation (SOFEA-BEI).}
    \label{bei:fig:sofea_bei_pipeline}
\end{figure}

\subsection{Scene-adaptive Noise Rejection}
\label{bei:subsec:application}
The proposed noise rejection method for flow estimation does not depend on the goodness-of-fit criteria\cite{Aung2018} for the plane fitting, i.e. check whether the fitted plane predicts timestamps of associated events within a tolerance level. Instead, we utilize an estimated upper bound of the true underlying distribution of optical flow speed, to implicitly assess the approximate upper bound likelihood of a given optical flow speed estimate. More precisely, our method rejects events with optical flow speed estimates that fall in the tails of the approximate probability upper bound, as it is highly improbable to observe such a value. Such an approach is possible due to the fact that spurious noise events and erroneous  flow estimates generally have an extreme flow speed associated to it. Hence, this simple noise rejection method is adaptive to the scene dynamics and independent of the local scene structure with practical bounds of the true optical flow speed.

\subsection{Extended Lifetime Augmentation and Reset}
\label{bei:sec:lifetime}

\par
As a variant of the Local Plane Fit method, SOFEA estimates the normal flow of an event by performing plane fitting on its corresponding local SAE. Thus, the global SAE can be approximated as a piecewise-continuous set of planes, with each plane defined for a distinct $1 \times 1$ pixel spatial neighbourhood. Theoretically, a BEI with infinite spatial resolution can then be estimated by evaluating the SAE at a particular time instant. Nevertheless, this BEI estimation approach requires highly accurate normal flow estimates, as well as high computational cost for the evaluation of each local time slice. Hence, we adopt the lifetime augmentation approach \cite{Mueggler2015}, with some refinements, for relatively accurate and efficient BEI estimation, even with low accuracy normal flow estimates.

\par
Ideally, the lifetime of an event signifies the amount of time required for an associated (new) event to be generated in its neighboring pixel location, due to apparent motion of their common corresponding feature. Suppose lifetimes of events are known, BEI can be rendered by simply activating specific pixels, given by the position of events, for a duration of their augmented lifetimes. However, event lifetimes depend on the continuous-time trajectory and velocity of the apparent motion, which are not trivial to estimate. Thus, the lifetime of an event is approximated as the maximum time required for the event to ``travel'' $1$ pixel, supposing the apparent motion velocity of the associated feature is constant, or equivalently the normal flow of the event is constant. Note that the lifetime is defined as the maximum time because in Local Plane Fit methods, only normal flow is estimated. As the normal flow  speed always underestimates the optical flow speed, only the maximum time can be estimated. Mathematically, the lifetime, $\tau$ of an event is given by:
\begin{equation}
	\tau = \frac{1}{\| \boldsymbol{v_\perp}  \|} = \| \nabla \Sigma_e(\boldsymbol{p}) \|
    \label{bei:eq:lifetime}
\end{equation}
In other words, $\tau$ is merely given by the spatial gradient magnitude of the associated local SAE.
\par
Nevertheless, the scene representation accuracy of the estimated BEI depends greatly on the accuracy of the estimated lifetimes, which in turn relies on the normal flow estimation accuracy. Apart from that, the constant velocity assumption may be violated in scenes with high acceleration or deceleration, resulting in an overestimation or underestimation of lifetimes. In order to enhance the quality of BEI estimates in general scenes and also when low accuracy normal flow estimates are used, we propose to augment each event with an extended lifetime, $\tau_\mathit{ext, bei}$ instead, given by:
\begin{equation}
	\tau_\mathit{ext, BEI} = \kappa_\mathit{ext, BEI} \ \tau, \ \kappa_\mathit{ext, BEI} \geq 1
    \label{bei:eq:ext_lifetime_bei}
\end{equation}
By doing so, it can be observed that edges in the BEI are approximately $\kappa_\mathit{ext, bei}$ pixels thick, due to the overestimation of lifetime by an approximate factor of $\kappa_\mathit{ext, BEI}$. To alleviate that, we also perform edge thinning by resetting (i.e. set to zero) the extended lifetime of the neighboring event, which is in the opposite flow direction of a given event, because the neighboring event generally indicates the ``source'' of the given event. As the optimal value of $\kappa_\mathit{ext, bei}$ might be different for accurate scene representation of the rendered BEI, and for accurate estimation of the scene complexity, given by $\widetilde{C} (t_j)$, an event can also be optionally augmented with a second extended lifetime for scene complexity estimation, as follows:
\begin{equation}
	\tau_\mathit{ext, \widetilde{C}} = \kappa_\mathit{ext, \widetilde{C}} \ \tau, \ \kappa_\mathit{ext, \widetilde{C}} \geq 1
    \label{bei:eq:ext_lifetime_C_tilde}
\end{equation}
With both lifetime extension and lifetime reset, the eventual effective lifetime, $\tau_\mathit{eff, BEI}$ or $\tau_\mathit{eff, \widetilde{C}}$ gives a relatively better estimate of the true lifetime. This is due to the larger error tolerance of normal flow estimates, given by the closed-loop nature of the lifetime estimation process. This allows for a crisp and accurate scene representation of the estimated BEI, which is visually similar to an edge map. Isolated active pixels in the estimated BEI can also be optionally filtered out, as noise in a BEI manifests as impulse noise.

\subsection{Scene-adaptive BEI Rendering}
\label{bei:sec:render}
Although \autoref{bei:sec:lifetime} allows for continuous-time BEI rendering, it is computationally infeasible and generally unnecessary, as a new BEI rendering is usually only required whenever it differs significantly from the most recent BEI rendered. Thus, a BEI rendering time interval is defined to govern the rendering time instants. Similar to how standard cameras capture intensity images at a fixed rate, it is rather reasonable to set it as a constant. Nevertheless, the amount of redundancy increases as the scene dynamics decreases. Thus, we propose a scene-adaptive BEI rendering time interval that approximately normalizes for the expected apparent displacement magnitude of features with apparent motion. In particular, it approximately normalizes for the expected number of events generated in an event substream.
\subsubsection{Event Generation Model}
\label{bei:subsec:model}
In this subsection, we will first devise a simplified model to describe how the event stream is formed by numerous features with apparent motion to assist scene-adaptive BEI rendering:
\begin{enumerate}
\item $c(t)$ is a deterministic function that defines the number of features with apparent motion (also known as the scene complexity), for all time $t$. A feature is defined as a unique surface point in the scene, within the field-of-view of the DVS, that causes events to be generated if it has apparent motion. \label{bei:def:c}
\item $\boldsymbol{v}_i(t) = \begin{bmatrix} \|\boldsymbol{v}_i(t)\| & \angle \boldsymbol{v}_i(t) \end{bmatrix}^\top$ is a deterministic vector-valued function that defines the optical flow of events associated to feature $i \in \{ \ 1, 2, \ldots, c(t)\ \}$ with apparent motion, for all time $t$, assuming infinite spatial resolution of the image sensor (hence events are generated in continuous-time). In particular, $\|\boldsymbol{v}_i(t)\|$ gives the speed and $\angle \boldsymbol{v}_i(t)$ gives the orientation of the optical flow. $\boldsymbol{v}_i(t)$ is also equivalent to the apparent motion velocity of feature $i \in \{ \ 1, 2, \ldots, c(t)\ \}$ at time $t$. \label{bei:def:v}
\item A short time interval, $\mathcal{I}(t) = ( \ t-\Delta T(t), t \ ]$ is defined for all time $t$, where $\Delta T(t)$ is a deterministic function that gives the length of the time interval. \label{bei:def:i}
\item $f_i(t)$ is a deterministic function that defines the event generation frequency associated to feature $i \in \{ \ 1, 2, \ldots, c(t)\ \}$ with apparent motion, for all time $t$. \label{bei:def:f}
\end{enumerate}
\par
With the definitions above, the assumptions are:
\begin{enumerate}
\item $c(t')$ and $\boldsymbol{v}_i(t'), \forall i \in \{ \ 1, 2, \ldots, c(t')\ \}$ have the same respective values for all time $t'$ in $\mathcal{I}(t)$. In other words, $c(t)$ and $\boldsymbol{v}_i(t), \forall i \in \{ \ 1, 2, \ldots, c(t)\ \}$ are constant throughout $\mathcal{I}(t)$. \label{bei:assmp:constant_c_v}
\item Given that $\boldsymbol{v}_i(t')$ has the same value for all time $t'$ in $\mathcal{I}(t)$, $f_i(t')$ also has the same value. In other words, if $\boldsymbol{v}_i(t)$ is constant throughout $\mathcal{I}(t)$, it is also the case for $f_i(t)$. Implicitly, this assumes that all pixels have the same contrast threshold, which could be negligibly violated in real DVS data due to random hardware mismatch. This assumption also disregards the subtle variation of $f_i(t')$ if the apparent motion direction of feature $i$, or equivalently $\angle \boldsymbol{v}_i(t')$, is non-Ordinal or non-Cardinal with respect to the DVS pixel grid. \label{bei:assmp:constant_f}
\end{enumerate}
Assumptions \ref{bei:assmp:constant_c_v} and \ref{bei:assmp:constant_f} can be combined to yield the following: $c(t'), \boldsymbol{v}_i(t')$ and $f_i(t')$, $\forall i \in \{ \ 1, 2, \ldots, c(t')\ \}$ have the same respective values for all time $t'$ in $\mathcal{I}(t)$. As a result, the event generation frequency of the event stream, given by:
\begin{equation}
	f(t') = \sum\limits_{i=1}^{c(t')} f_i(t')
    \label{bei:eq:f}
\end{equation}
also has the same value for all time $t'$ in $\mathcal{I}(t)$. Also, it can be observed that a direct relationship exists between $\|\boldsymbol{v}_i(t)\|$ and $f_i(t)$. In fact, it is given by:
\begin{equation}
	\|\boldsymbol{v}_i(t)\| = \alpha(\angle \boldsymbol{v}_i(t)) \ f_i(t)
    \label{bei:eq:v_f_relation}
\end{equation}
where $\alpha(\theta)$ is a deterministic function that specifies the conversion from frequency to pixel units depending on the orientation of optical flow.
\subsection{Scene-adaptive Rendering}
Scene-adaptive motion normalization can be achieved between successive BEI renders by using the expected number of events generated in an event substream within a small temporal window. To formally derive the scene-adaptive BEI rendering time interval, we first present the following definitions:
\begin{enumerate}
\item A short time interval, $\mathcal{I}_\mathit{BEI} (t) = [ \ t, \ t + \Delta T_\mathit{BEI} (t) \ )$ is defined for all time $t$, where $\Delta T_\mathit{BEI} (t)$ is a deterministic function that gives the length of the time interval (the scene-adaptive BEI rendering time interval).
\item $\mathbf{s}_i(t)$ is a deterministic vector-valued function that defines the apparent displacement of feature $i \in \{ \ 1, 2, \ldots, c(t)\ \}$ with apparent motion from the start to the end of the time interval $\mathcal{I}_\mathit{BEI} (t)$ (i.e. between time $t$ and $t + \Delta T_\mathit{BEI} (t)$). Mathematically, $\mathbf{s}_i(t)$ is given by:
\begin{equation}
	\mathbf{s}_i(t) = \int_t^{t + \Delta T_\mathit{BEI} (t)}
    	\begin{bmatrix}
			\boldsymbol{v}_{i, x}(t') \\
            		  \boldsymbol{v}_{i, y}(t')
		\end{bmatrix}
        	  dt'
              \label{bei:eq:s}
\end{equation}
\item $\{ \ S(t),\ t \in [ \ 0, \infty \ ) \ \}$ is a continuous stochastic process, whereby it is equally likely for $S(t)$ to give any of $\|\boldsymbol{s}_i(t)\|, \ i \in \{ \ 1, 2, \ldots, c(t)\ \}$. In other words, $S(t)$ describes the distribution of $\|\boldsymbol{s}_i(t)\|$. With that, the  expected value of $S(t)$ is given by: \label{bei:def:S}
\begin{equation}
	\mathbb{E} [ S(t) ] = \frac{1}{c(t)} \sum\limits_{i=1}^{c(t)} \|\boldsymbol{s}_i(t)\|
	\label{bei:eq:exp_S}
\end{equation}

\item $\{ \ F(t),\ t \in [ \ 0, \infty \ ) \ \}$ is a continuous stochastic process, whereby it is equally likely for $F(t)$ to give any of $f_i(t), \ i \in \{ \ 1, 2, \ldots, c(t)\ \}$. In other words, $F(t)$ describes the distribution of $f_i(t)$. Thus, the expected value of $F(t)$ is given by (\autoref{bei:eq:f}):
\begin{equation}
	\mathbb{E} [ F(t) ] = \frac{1}{c(t)} \sum\limits_{i=1}^{c(t)} f_i(t) = \frac{f(t)}{c(t)}
    		  \label{bei:eq:exp_F}
\end{equation}
\label{bei:def:F}
\item $\{ \ V(t),\ t \in [ \ 0, \infty \ ) \ \}$ is a continuous stochastic process, whereby it is equally likely for $V(t)$ to give any of $\|\boldsymbol{v}_i(t)\|, \ i \in \{ \ 1, 2, \ldots, c(t)\ \}$. In other words, $V(t)$ describes the distribution of $\|\boldsymbol{v}_i(t)\|$. With that, the  expected value of $V(t)$ is given by (\autoref{bei:eq:v_f_relation}):
\begin{equation}
	\mathbb{E} [ V(t) ] = \frac{1}{c(t)} \sum\limits_{i=1}^{c(t)} \|\boldsymbol{v}_i(t)\| = \frac{1}{c(t)} \sum\limits_{i=1}^{c(t)} \alpha(\angle \boldsymbol{v}_i(t)) \ f_i(t)
	\label{bei:eq:exp_V}
\end{equation}
Because $\alpha(\angle \boldsymbol{v}_i(t)) \in [ \ 1, \sqrt{2} \ ]$, the following inequality regarding $\mathbb{E} [ V(t) ]$ and $\mathbb{E} [ F(t) ]$ holds:
\begin{equation}
	\mathbb{E} [ F(t) ] \leq \mathbb{E} [ V(t) ] \leq \sqrt{2} \ \mathbb{E} [ F(t) ]
	\label{bei:eq:exp_v_exp_f_ineq}
\end{equation}
\end{enumerate}
Along with the above definitions, we impose similar assumptions, as given by Assumptions \ref{bei:assmp:constant_c_v} and \ref{bei:assmp:constant_f} in \autoref{bei:subsec:model}, but on the time interval $\mathcal{I}_\mathit{BEI} (t)$ instead. Thus, \autoref{bei:eq:s} can be simplified to:
\begin{equation}
	\mathbf{s}_i(t) =
    	\begin{bmatrix}
			\boldsymbol{v}_{i, x}(t') \\
            		  \boldsymbol{v}_{i, y}(t')
		\end{bmatrix}
        	  \Delta T_\mathit{BEI} (t)
              \label{bei:eq:s_simple}
\end{equation}
Besides, \autoref{bei:eq:exp_S} can also be simplified to give:
\begin{equation}
	\mathbb{E} [ S(t) ] = \frac{1}{c(t)} \sum\limits_{i=1}^{c(t)} \|\boldsymbol{v}_i(t)\| \ \Delta T_\mathit{BEI} (t)
	\label{bei:eq:exp_S_simple} = \mathbb{E} [ V(t) ] \ \Delta T_\mathit{BEI} (t)
\end{equation}
Substituting \autoref{bei:eq:exp_v_exp_f_ineq} into the above equation, the following bounds can be imposed on $\mathbb{E} [ S(t) ]$:
\begin{equation}
	\mathbb{E} [ F(t) ] \ \Delta T_\mathit{BEI} (t) \leq \mathbb{E} [ S(t) ] \leq \sqrt{2} \ \mathbb{E} [ F(t) ] \ \Delta T_\mathit{BEI} (t)
	\label{bei:eq:exp_s_exp_f_ineq}
\end{equation}
Note that $\mathbb{E} [ F(t) ] \ \Delta T_\mathit{BEI} (t)$ actually gives the expected number of events generated in an event substream within $\mathcal{I}_\mathit{BEI} (t)$ (proof omitted). Similar to the case of $\Delta T (t)$, the choice of $\Delta T(t)_\mathit{BEI}$ is arbitrary, despite being deterministic. Therefore, we suppose that $\Delta T(t)_\mathit{BEI}$ is defined such that it coincides with:
\begin{equation}
	\Delta T(t)_\mathit{BEI} = \frac{2}{1 + \sqrt{2}} \cdot \frac{ \overline{s} }{ \hat{F} (t) }
    \label{bei:eq:delta_t_bei}
\end{equation}
in the first realization of $\hat{F} (t)$, where $\overline{s}$ is a constant representing the desired value of $\mathbb{E} [ S(t) ]$. Substituting \autoref{bei:eq:delta_t_bei} into \autoref{bei:eq:exp_s_exp_f_ineq} gives the following:
\begin{equation}
	\frac{2}{1 + \sqrt{2}} \cdot \frac{ \mathbb{E} [ F(t) ] }{ \hat{F} (t) } \cdot \overline{s} \leq \mathbb{E} [ S(t) ] \leq \frac{2  \sqrt{2}}{1 + \sqrt{2}} \cdot \frac{ \mathbb{E} [ F(t) ] }{ \hat{F} (t) } \cdot \overline{s}
	\label{bei:eq:exp_s_exp_f_ineq_expanded}
\end{equation}
Thus, $\Delta T(t)_\mathit{BEI}$ is defined in such a way that $\overline{s}$ is the midpoint of the above bounds on $\mathbb{E} [ S(t) ]$. Equivalently, it is also defined such that $\frac{2}{1 + \sqrt{2}} \overline{s}$ is the desired expected number of events generated in an event substream within $\mathcal{I}_\mathit{BEI} (t)$ (i.e. suppose the mean squared error of $\hat{F}(t)$ is kept small relative to $\mathbb{E} [ F(t) ]$, $\mathbb{E} [ F(t) ] \ \Delta T_\mathit{BEI} (t) \approx \frac{2}{1 + \sqrt{2}} \overline{s}$). Therefore, motion normalization can be achieved between successive BEI renders by using $\Delta T(t)_\mathit{BEI}$ as the rendering interval.

\begin{figure}[t]
    \centering
		\includegraphics[width=1\linewidth]{{figures/network_diagram}}
    \caption{Schematic representation of the FCN used for extracting superevents from SOFEA-BEI}
    \label{fig:network_fig}
\end{figure}
\section{SOFEA-BEI-FCN Superevents}
\label{sec:superevents}
In this section, we introduce our CNN-based superevents segmentation method using the output of SOFEA-BEI. We make use of the architecture proposed in \cite{yang2020sup}, which has been shown to perform well on RGB compared to several state-of-the-art methods. Superpixels as defined for standard images are supposed to delineate the edges, and thus binary images containing edge-like structures can also be well-segmented by the fully convolutional neural network architecture proposed in \cite{yang2020sup}.
\par
The model in \cite{yang2020sup} divides a given image into a number of grids with preset dimensions. These grids form the initial superpixels. The model then makes use of a series of convolutions followed by one of deconvolutions in order to achieve the desired segmentation, as illustrated in Figure \ref{fig:network_fig}. Further, features from earlier layers are concatenated to those from deeper layers via skip connections. The final output is an assignment map

\begin{equation}
	\begin{split}
		g: \mathbb{Z}^{H \times W} &\to N_r \\
		\boldsymbol{r} &\mapsto g( \boldsymbol{r} ) = n
	\end{split}
\end{equation}
where $\boldsymbol{r} = (y, x)$ and $n$ denotes the superpixel that each pixel has been assigned to. Also, $N_r$ denotes the nine superpixels in the neighbourhood of the given pixel. A lot of the simplicity in \cite{yang2020sup} is the direct result of this part of their formulation: any given pixel can only be assigned to one of the nine superpixels in its immediate neighbourhood. This results in compact superpixels and oversegmenation to an extent that is greater than what is observed in traditional superpixel segmentation methods. \cite{yang2020sup} justify this by arguing that since the superpixels thus obtained will eventually be used in downstream tasks, the latter may group perceptually similar superpixels if such a need arises. Hence, this formulation does not compromise on the quality of superpixels for use in downstream tasks, while also helping to decrease the number of training and inference time computations.
\par
For the sake of better visualization in Figure \ref{fig:network_fig}, superevents with little to no events have been merged. In practice, this was achieved by computing the ratio of the number of active events to that of the total number of pixels inside a given superevent. They were merged together whenever this ratio turned out to be less than a pre-determined threshold (set to 0.05 in this work).

\begin{figure}[t]
	\centering

	\begin{subfigure}{0.38\linewidth}
    	\centering
		\includegraphics[width=1\linewidth]{{figures/aps_dynamic_6dof_405}}
		\caption{APS intensity image}
		\label{bei:fig:dynamic_6dof_aps}
	\end{subfigure} \hspace{1em}
	\begin{subfigure}{0.38\linewidth}
    	\centering
        \includegraphics[width=1\linewidth]{{figures/canny_0.2_0.6_dynamic_6dof_405}}
		\caption{Edge map via \cite{Canny1986}}
		\label{bei:fig:dynamic_6dof_canny}
	\end{subfigure}

    \begin{subfigure}{0.38\linewidth}
    	\centering
		\includegraphics[width=1\linewidth]{{figures/sofea_bei_dynamic_6dof_405_ext1_noreset}}
		\caption{$\kappa_\mathit{ext, BEI} = 1$, no reset}
		\label{bei:fig:dynamic_6dof_sofea_bei_ext1_noreset}
	\end{subfigure} \hspace{1em}
	\begin{subfigure}{0.38\linewidth}
    	\centering
        \includegraphics[width=1\linewidth]{{figures/sofea_bei_dynamic_6dof_405_ext1_reset}}
		\caption{$\kappa_\mathit{ext, BEI} = 1$, with reset}
		\label{bei:fig:dynamic_6dof_sofea_bei_ext1_reset}
	\end{subfigure}

	\begin{subfigure}{0.38\linewidth}
    	\centering
		\includegraphics[width=1\linewidth]{{figures/sofea_bei_dynamic_6dof_405_ext3_noreset}}
		\caption{$\kappa_\mathit{ext, BEI} = 3$, no reset}
		\label{bei:fig:dynamic_6dof_sofea_bei_ext3_noreset}
	\end{subfigure} \hspace{1em}
	\begin{subfigure}{0.38\linewidth}
    	\centering
        \includegraphics[width=1\linewidth]{{figures/sofea_bei_dynamic_6dof_405_ext3_reset}}
		\caption{$\kappa_\mathit{ext, BEI} = 3$, with reset}
		\label{bei:fig:dynamic_6dof_sofea_bei_ext3_reset}
	\end{subfigure}

	\begin{subfigure}[t]{0.38\linewidth}
    	\centering
        \includegraphics[width=1\linewidth]{{figures/sofea_bei_flow_wheel_dynamic_6dof_405_ext3_reset}}
		\caption{$\kappa_\mathit{ext, BEI} = 3$, with reset, augmented with normal flow orientation}
		\label{bei:fig:dynamic_6dof_sofea_bei_flow_ext3_reset}
	\end{subfigure} \hspace{1em}
	\begin{subfigure}[t]{0.38\linewidth}
    	\centering
		\includegraphics[width=1\linewidth]{{figures/sofea_bei_dynamic_6dof_405_ext5_reset}}
		\caption{$\kappa_\mathit{ext, BEI} = 5$, with reset}
		\label{bei:fig:dynamic_6dof_sofea_bei_ext5_reset}
	\end{subfigure}

	\caption{Comparison of BEI estimated using SOFEA-BEI, under various extended lifetime augmentation and lifetime reset configurations, on the {\fontfamily{pcr}\selectfont dynamic\_6dof} sequence.}
	\label{bei:fig:lifetime_analysis}
\end{figure}

\begin{figure}[t]
	\centering

    \begin{minipage}[t]{0.29\linewidth}
		\begin{subfigure}{1\linewidth}
    		\centering
			\includegraphics[width=1\linewidth]{{figures/aps_slider_far_70}}
			\caption{}
			\label{bei:fig:slider_far_aps}
		\end{subfigure}

		\begin{subfigure}{1\linewidth}
            \centering
            \includegraphics[width=1\linewidth]{{figures/canny_0.2_0.5_slider_far_70}}
            \caption{}
            \label{bei:fig:slider_far_canny}
		\end{subfigure}

		\begin{subfigure}{1\linewidth}
            \centering
            \includegraphics[width=1\linewidth]{{figures/sofea_bei_slider_far_70}}
            \caption{}
            \label{bei:fig:slider_far_sofea_bei}
		\end{subfigure}

		\begin{subfigure}{1\linewidth}
            \centering
            \includegraphics[width=1\linewidth]{{figures/mueggler2015_slider_far_70}}
            \caption{}
            \label{bei:fig:slider_far_mueggler2015}
		\end{subfigure}
    \end{minipage} \hspace{1em}
    \begin{minipage}[t]{0.29\linewidth}
		\begin{subfigure}{1\linewidth}
    		\centering
			\includegraphics[width=1\linewidth]{{figures/aps_urban_100}}
			\caption{}
			\label{bei:fig:urban_aps}
		\end{subfigure}

		\begin{subfigure}{1\linewidth}
            \centering
            \includegraphics[width=1\linewidth]{{figures/canny_0.2_0.5_urban_100}}
            \caption{}
            \label{bei:fig:urban_canny}
		\end{subfigure}

		\begin{subfigure}{1\linewidth}
            \centering
            \includegraphics[width=1\linewidth]{{figures/sofea_bei_urban_100_bei}}
            \caption{}
            \label{bei:fig:urban_sofea_bei}
		\end{subfigure}

		\begin{subfigure}{1\linewidth}
            \centering
            \includegraphics[width=1\linewidth]{{figures/mueggler2015_urban_100}}
            \caption{}
            \label{bei:fig:urban_mueggler2015}
		\end{subfigure}
    \end{minipage} \hspace{1em}
    \begin{minipage}[t]{0.29\linewidth}
		\begin{subfigure}{1\linewidth}
    		\centering
			\includegraphics[width=1\linewidth]{{figures/aps_dynamic_rotation_1000}}
			\caption{}
			\label{bei:fig:dynamic_rotation_aps}
		\end{subfigure}

		\begin{subfigure}{1\linewidth}
            \centering
            \includegraphics[width=1\linewidth]{{figures/canny_0.07_dynamic_rotation_1000}}
            \caption{}
            \label{bei:fig:dynamic_rotation_canny}
		\end{subfigure}

		\begin{subfigure}{1\linewidth}
            \centering
            \includegraphics[width=1\linewidth]{{figures/sofea_bei_dynamic_rotation_1000}}
            \caption{}
            \label{bei:fig:dynamic_rotation_sofea_bei}
		\end{subfigure}

		\begin{subfigure}{1\linewidth}
            \centering
            \includegraphics[width=1\linewidth]{{figures/mueggler2015_dynamic_rotation_1000}}
            \caption{}
            \label{bei:fig:dynamic_rotation_mueggler2015}
		\end{subfigure}
    \end{minipage}

	\caption[BEI scene representation accuracy comparison between SOFEA-BEI and \cite{Mueggler2015} on the {\fontfamily{pcr}\selectfont slider\_far}, {\fontfamily{pcr}\selectfont urban} and {\fontfamily{pcr}\selectfont dynamic\_rotation} sequences.]{BEI scene representation accuracy comparison between SOFEA-BEI (3\textsuperscript{rd} row) and \cite{Mueggler2015} (4\textsuperscript{th} row) on the {\fontfamily{pcr}\selectfont slider\_far} (1\textsuperscript{st} column), {\fontfamily{pcr}\selectfont urban} (2\textsuperscript{nd} column) and {\fontfamily{pcr}\selectfont dynamic\_rotation} (3\textsuperscript{rd} column) sequences. APS intensity images (1\textsuperscript{st} row) and edge maps obtained via Canny edge detection (2\textsuperscript{nd} row) are also visualized.}
	\label{bei:fig:bei_comparison_2}
\end{figure}
\section{Experiments}

\label{sec:experiments}
We use sequences from the event camera dataset \cite{Mueggler2017} for testing SOFEA-BEI and subsequently the superevents extraction using the FCN. The superevents concept is introduced for the first time using event cameras, to the best of our knowledge. Thus, evaluation is focused on comparing SOFEA-BEI-FCN to SOFEA-BEI-SLIC \cite{slicachanta2012} and APS-FCN with simultaneously recorded intensity frames from the publicly available dataset (upper baseline). The ground truth segmentation results for SOFEA-BEI are obtained from the corresponding APS-Canny edge maps for the calculation of boundary precision and recall \cite{yang2020sup} with the natural assumption that APS frames have more information for segmentation compared to BEI. Next we qualitatively evaluate SOFEA-BEI results comapared to vanilla SOFEA.

\begin{figure*}[t]
	\centering

	\begin{subfigure}{0.3\linewidth}
    	\centering
		\includegraphics[width=1\linewidth]{{figures/slider_depth_superevents}}
		\caption{{\fontfamily{pcr}\selectfont slider\_depth}}
		\label{bei:fig:slider_depth}
	\end{subfigure}
		\begin{subfigure}{0.31\linewidth}
    	\centering
		\includegraphics[width=1\linewidth]{{figures/shapestrans_superevents}}
		\caption{{\fontfamily{pcr}\selectfont shapes\_translation}}
		\label{bei:fig:slider_trans}
	\end{subfigure}
   \begin{subfigure}{0.33\linewidth}
    	\centering
		\includegraphics[width=1\linewidth]{{figures/dynrot_superevents}}
		\caption{{\fontfamily{pcr}\selectfont dynamic\_rotation}}
		\label{bei:fig:dyn_rot}
	\end{subfigure}
		\caption{Comparison of SOFEA-BEI superevents, using FCN and SLIC methods, to APS-FCN superpixels.}
	\label{bei:fig:superevents_analysis}
\end{figure*}
\subsection{Extended Lifetime Augmentation \& Reset}
\label{bei:subsec:exp_ext_lifetime_augm_reset}
We carry out a qualitative analysis on the effectiveness of extended lifetime augmentation and reset in SOFEA-BEI for enhanced BEI estimation accuracy via the {\fontfamily{pcr}\selectfont dynamic\_6dof} sequence \cite{Mueggler2017}. \autoref{bei:fig:lifetime_analysis} illustrates the BEI estimated using SOFEA-BEI, under various extended lifetime augmentation and lifetime reset configurations, on the {\fontfamily{pcr}\selectfont dynamic\_6dof} sequence. An APS intensity image of the scene is shown in \autoref{bei:fig:dynamic_6dof_aps}, and its corresponding edge map generated using Canny edge detection is given in \autoref{bei:fig:dynamic_6dof_canny}. Since an edge map of the APS image gives a good representation of a BEI, we consider \autoref{bei:fig:dynamic_6dof_canny} as the reference for the BEI estimates.  It is worth noting that there is no fixed ``sampling interval'' for the whole pipeline, as SOFEA-BEI is an event-based method adapted to dynamically handle object/camera motion.
\par
\autoref{bei:fig:dynamic_6dof_sofea_bei_ext1_noreset} illustrates the BEI estimated using the original lifetime augmentation approach \cite{Mueggler2015} (equivalently, without lifetime extension or $\kappa_\mathit{ext, BEI} = 1$) without performing lifetime reset, while \autoref{bei:fig:dynamic_6dof_sofea_bei_ext1_reset} shows the BEI estimated with lifetime reset done  for edge thinning. It can be observed that edges in both BEI estimates are generally one pixel thick, as opposed to BEI estimated using the direct accumulation approach (\autoref{fig:intro_fig}). However, several pixel thick edges still exist in \autoref{bei:fig:dynamic_6dof_sofea_bei_ext1_noreset} due to the overestimation of lifetime. This issue is almost completely rectified after lifetime reset is incorporated, which produces a more crisp representation of the scene.
\par
On the other hand, extended lifetime augmentation is performed in \autoref{bei:fig:dynamic_6dof_sofea_bei_ext3_reset} to \autoref{bei:fig:dynamic_6dof_sofea_bei_ext5_reset}. Observe that edges in \autoref{bei:fig:dynamic_6dof_sofea_bei_ext3_noreset} are approximately three pixels thick with $\kappa_\mathit{ext, BEI} = 3$ and lifetime reset not done. The incorporation of lifetime reset in \autoref{bei:fig:dynamic_6dof_sofea_bei_ext3_reset} introduces a substantial improvement in terms of scene representation accuracy, when compared to \autoref{bei:fig:dynamic_6dof_sofea_bei_ext1_reset}. Note that the amount of detail captured is now similar to that of the edge map. Nevertheless, the amount of noise increases as well, due to events with extended lifetimes not being reset. Fortunately, there is still an apparent increase in the Signal-to-Noise Ratio (SNR). \autoref{bei:fig:dynamic_6dof_sofea_bei_flow_ext3_reset} augments \autoref{bei:fig:dynamic_6dof_sofea_bei_ext3_reset} with normal flow orientation estimates, given by the color of the active pixels according to the color wheel displayed. Observe that active pixels associated to the same edge have a common color, which corresponds to the direction perpendicular to the edge. This is desirable for an effective lifetime reset. \autoref{bei:fig:dynamic_6dof_sofea_bei_ext5_reset} illustrates the BEI estimated with $\kappa_\mathit{ext, BEI} = 5$ and lifetime reset. Note that the amount of detail captured is almost the same as that of $\kappa_\mathit{ext, BEI} = 3$, despite a significant increase of noise until the point where feature tracks can be clearly seen. This suggests that $\kappa_\mathit{ext, BEI} = 3$ and lifetime reset is sufficient for SOFEA-BEI to achieve accurate BEI estimation.
\par
To validate SOFEA-BEI, \cite{Mueggler2015} is employed as the benchmark to evaluate its performance. In addition, edge maps generated using Canny edge detection on the APS intensity image also serves as a ground truth reference for BEI estimates. The qualitative comparison is performed with \cite{Mueggler2015} on three sequences from \cite{Mueggler2017}. This collection of sequences covers a wide range of scenes with varying complexity, dynamic range and motion profiles. The details captured in a BEI need not be captured in the edge map, due to limitations of the edge detection, motion blur in the APS image, and high dynamic range of the scene.
\par
\autoref{bei:fig:bei_comparison_2} illustrates the comparison on the three mentioned sequences. In general, it can be observed that SOFEA-BEI estimates achieve remarkable scene representation, similar to the corresponding edge maps. BEI estimates of \cite{Mueggler2015} however are generally sparse in nature, as it fails to capture the abundance of detail in the scene. This is heavily reflected in \autoref{bei:fig:slider_far_mueggler2015} and \autoref{bei:fig:urban_mueggler2015}. As noted earlier, there exist details not captured in the edge map, but reflected in the SOFEA-BEI estimates. For instance, the concentric rings of the dartboard in the {\fontfamily{pcr}\selectfont slider\_depth} sequence, the top left edges of the building in the {\fontfamily{pcr}\selectfont urban} sequence and the details of objects in the {\fontfamily{pcr}\selectfont dynamic\_rotation} sequence. This is made possible due to the high dynamic range of the event stream and accurate SOFEA-BEI estimates.

\subsection{Superevents}
\label{subsec:supeventsresults}
\autoref{bei:fig:superevents_analysis} illustrates the boundary precision and recall for different sequences, which measures the boundary adherence of the superevents/pixels according to the ground truth with a 3$\times$3 tolerance window as done in earlier works \cite{yang2020sup}. The performance of the APS frames using FCN is much better compared to its BEI counterpart using the event stream, owing to the richness of the intensity information. This does not imply that the event stream is lacking information, but rather the way it is used in this work as a BEI. To investigate whether the lower performance is due to the use of the pre-trained network from \cite{yang2020sup}, the popular hand-crafted SLIC \cite{slicachanta2012} was used on BEI to extract superevents. This doubt is ruled out as SOFEA-BEI-SLIC has generally lower performance, especially in the shapes dataset where the edge maps form a strong ground truth. A short clip of the superevents can be viewed here\footnote{\url{https://tinyurl.com/hv3bkz5h}}. This video shows good boundary separation between the superevents at different depths, especially when a new object closer to the camera enters the field-of-view.
\par
To the best of our knowledge, we have introduced the first mid-level perceptual grouping of events using event-based flow estimation that does not limit the rate at which superevents are extracted. In theory, superevents can be extracted for every flow estimate using the proposed scheme although the scene-adaptive BEI is better suited for generating superevents regardless of scene and camera dynamics. Overall, superevents have several potential applications, including native semantic segmentation, generating templates for object tracking, depth estimation, feature extraction, etc., and showcasing them is beyond the scope of this paper.

\section{Conclusion}
\label{sec:conc}
This paper proposed the first superevents extraction method for event cameras, leveraging on a robust binary event image generation scheme and fully convolutional networks. This opens up a wide range of possibilities in the neuromorphic vision domain to aggregate sparse, disjoint visual information into more meaningful intermediate representations. We evaluated the superevents on the event camera dataset and highlighted the gap in boundary scores between the intensity and events modalities. Future work to address this gap using completely asynchronous processing shall be explored.

\bibliographystyle{ACM-Reference-Format}


\section{Introduction}
To search for the CP and phase boundary in the QCD phase diagram, RHIC has undertaken its first phase BES(Beam Energy Scan) Program\cite{b2,b2a,b2b,b2c,b2d}. Since the moments of the distributions of conserved quantities, for example net-baryon number, in the relativistic heavy ion collisions are related to the correlation length $\xi$ of the system\cite{b3,b3a}, they are believed to be good signatures of the QCD phase transition and its CP.

It is known that the moments of the baryon number are proportional to the various order baryon-number susceptibilities, and their relations are shown in \cite{b4}. It can be seen that when relating the susceptibilities to the moments a volume term appears. In order to cancel the volume term, the products of the moments, $S\sigma$ and $\kappa\sigma^2$, are constructed as the experimental observables. The results in RHIC of these observables show a centrality and energy dependence\cite{b5}. When studying the quark numbers at finite temperature and chemical potential, it is found that the quark-number density is determined by the corresponding dressed quark propagator only\cite{b6,b6a}. Then by calculating the derivatives of the quark-number density with respect to $\mu$, we can calculate the quark-number susceptibilities at finite temperature and chemical potential.

In this paper, we obtain the dressed quark propagator in the Nambu-Jona-Lasinio(NJL) model. The NJL model\cite{b7,b7a,b7b}, as a widely adopted phenomenological model of QCD, is used to study the dynamical chiral symmetry breaking(DCSB) and the interaction between hadrons. The NJL model keeps the basic symmetries of QCD and simplifies the interactions in QCD to four-body interactions.
\section{Moments by NJL Model}
\label{two}

The commonly used Largrangian of the two-flavor NJL model is
\begin{equation}
\mathcal{L}=\bar{\psi}(i\gamma^{\mu}\partial_{\mu}-m)\psi+G[(\bar{\psi}\psi)^2+(\bar{\psi}i\gamma_5\vec{\tau}\psi)^2],
\label{1}
\end{equation}
where $m=m_u=m_d$, $\psi=(\psi_u,\psi_d)^T$ and $G$ is the effective coupling strength of the four-point quark interaction. With the mean field approximation of Eq.{\ref{1}}, the effective quark mass $M$ is
\begin{equation}
M=m-2G<\bar{\psi}\psi>,
\label{2}
\end{equation}
where the quark condensate at finite temperature($T$) and chemical potential($\mu$) is defined as
\begin{eqnarray}
\begin{split}
<\bar{\psi}\psi>=&-TN_cN_f\sum^{+\infty}_{-\infty}\int^{\Lambda}\frac{d^3p}{(2\pi)^3}\frac{4M}{p^2+M^2+\widetilde{\omega}^2_n}\\
=&-\frac{N_cN_fM}{2\pi^2}\int^{\Lambda}dp\frac{p^2}{\sqrt{p^2+M^2}}(th(\frac{\sqrt{p^2+M^2}+\mu}{2T})+th(\frac{\sqrt{p^2+M^2}-\mu}{2T})),
\end{split}
\end{eqnarray}
where $N_c=3$ is the number of color and $N_f=2$ is the number of flavor, $\widetilde{\omega}_n=2\pi nT+i\mu$ and the NJL model three-momentum non-covariant cut-off is $\Lambda$. In this paper, we adopt a varying coupling strength $G$\cite{b8,b8a,b8b,b8c,b9}, which is used to reproduce the Lattice results at finite temperature\cite{b9}:
\begin{equation}
G=G_1+G_2<\bar{\psi}\psi>,
\label{g2}
\end{equation}
where $G$  in the NJL model represents the effective gluon propagator, $G_2<\bar{\psi}\psi>$ reflects the contribution of the two-quark condensate to the gluon propagator and $G_1$ reflects all the other condensate contributions. Since the quark propagator and gluon propagator are coupled with each other by QCD and the quark propagator in Nambu and Wigner phase are very different, then the corresponding gluon propagators in these two phases should be different too. At the same time Lattice results have shown that the gluon propagator, rather than a constant, evolves with temperature. That is to say a constant $G$ does not meet these requirements. In \cite{b8c} they investigate how to extract the feedback of quark from gluon propagator and get Eq.{\ref{g2}}.The parameter set used in this paper is $m=5.6 MeV$, $\Lambda=587.9 MeV$, $G_1=5.564\times10^{-6} MeV^{-2}$ and $G_2=-3.16\times10^{-14} MeV^{-5}$, which is proved to be successful in fitting the quark condensate to the Lattice results at finite temperature\cite{b9}. Then by Eq.{\ref{2}} we could get the effective quark mass $M$ as a function of $T$ and $\mu$. Since the quark number density of $N_c$ color and $N_f$ flavor is\cite{b6,b6a}
\begin{eqnarray}
\begin{split}
n_q=&<\psi^+\psi>=-TN_cN_f\sum^{+\infty}_{-\infty}\int\frac{d^3p}{(2\pi)^3}T_r[G(p,\mu)\gamma_4]\\
=&TN_cN_f\sum^{+\infty}_{-\infty}\int\frac{d^3p}{(2\pi)^3}\frac{4i\widetilde{\omega}_n}{p^2+M^2+\widetilde{\omega}^2_n}\\
=&\frac{TN_cN_f}{\pi^2}\int dp\cdot p^2\sum^{+\infty}_{-\infty}(\frac{1}{\epsilon-(i\omega_n-\mu)}-\frac{1}{\epsilon+(i\omega_n-\mu)})\\
=&\frac{N_cN_f}{\pi^2}\int dp\cdot p^2(\frac{1}{e^{\beta(\epsilon-\mu)}+1}-\frac{1}{e^{\beta(\epsilon+\mu)}+1}),
\end{split}
\end{eqnarray}
where $\epsilon=\sqrt{p^2+M^2}$. Then the quark number density of a single color and flavor is
\begin{equation}
n=\frac{1}{\pi^2}\int dp\cdot p^2(\frac{1}{e^{\beta(\epsilon-\mu)}+1}-\frac{1}{e^{\beta(\epsilon+\mu)}+1}).
\label{4}
\end{equation}
The variance of baryon number density is
\begin{eqnarray}
\begin{split}
\sigma^2 &=<(n_{B}-<n_{B}>)^2>\\
&=\frac{1}{3^2}<(n_{q}-<n_{q}>)^2>\\
&=\frac{1}{3^2}<[(n_u+n_d)-<n_u+n_d>]^2>\\
&=\frac{1}{3^2}<(n_u-<n_u>)^2>+<(n_d-n_d)^2>\\
&=\frac{N_f}{3^2}<(n_{f}-<n_{f}>)^2>.
\label{3}
\end{split}
\end{eqnarray}
where $n_f$ represents the quark number density of a single flavor. Since the up quarks and down quarks are independent, as in \cite{b9a}, we have $<(n_{u}-<n_{u}>)(n_{d}-<n_{d}>)>=0$ here. And $n_f=n_1+n_2+n_3$, where $n_1$, $n_2$ and $n_3$ represent the quark number density of three different color respectively. Since the confinement nature of QCD, $<n_in_j>=<n^2>$ ($i,j=1,2,3$). Then the variance of $n_f$ is
\begin{eqnarray}
\begin{split}
<(n_{f}-<n_{f}>)^2>&=<[(n_1+n_2+n_3)-<n_1+n_2+n_3>]^2>\\
&=N_c^2<(n-<n>)^2>,
\end{split}
\end{eqnarray}
and $<(n-<n>)^2>$ could be calculate in theory by the relation $<(n-<n>)^2>=T\frac{\partial n}{\partial \mu}=T\cdot\chi^{(2)}$. Then with Eq.{\ref{3}}, we could get
\begin{eqnarray}
\begin{split}
\sigma^2 &=\frac{N_f}{3^2}<(n_{f}-<n_{f}>)^2>\\
&=\frac{N_fN_c^2}{3^2}<(n-<n>)^2>\\
&=N_f<(n-<n>)^2>=N_f\cdot\chi^{(2)}.
\end{split}
\end{eqnarray}
Similarly, we could get
\begin{eqnarray}
\begin{split}
S\sigma &=\frac{<(n_{B}-<n_{B}>)^3>}{<(n_{B}-<n_{B}>)^2>}=\frac{N_f<(n-<n>)^3>}{N_f<(n-<n>)^2>}\\
&=\frac{T^2\chi^{(3)}}{T\chi^{(2)}}=T\frac{\chi^{(3)}}{\chi^{(2)}},
\end{split}
\end{eqnarray}
and
\begin{eqnarray}
\begin{split}
\kappa\sigma^2 &=\frac{<(n_{B}-<n_{B}>)^4>-3<(n_{B}-<n_{B}>)^3>}{<(n_{B}-<n_{B}>)^2>}\\
&=\frac{<(n-<n>)^4>-3<(n-<n>)^3>}{<(n-<n>)^2>}\\
&=T^2\frac{\chi^{(4)}}{\chi^{(2)}}.
\end{split}
\end{eqnarray}
where $\frac{\partial^2 n}{\partial \mu^2}=\chi^{(3)}$ and $\frac{\partial^3 n}{\partial \mu^3}=\chi^{(4)}$. That is to say by calculating the high order derivatives of Eq.{\ref{4}} with respect to $\mu$, we could get the experimental observables $S\sigma$ and $\kappa\sigma^2$ in the relativistic heavy ion collisions.

\section{Results}
\label{three}
Before comparing the $S\sigma$ and $\kappa\sigma^2$ results of the NJL model with the experimental data in the relativistic heavy ion collisions, we should find the correspondence of the freeze-out temperature $T$ and the quark chemical potential $\mu$ in the NJL results to the collision energies $\sqrt{S_{NN}}$ in the relativistic heavy ion collisions.

\begin{figure}
\centering
\includegraphics[width=1\linewidth]{njl.eps}
\caption{Solid lines represent the upper and the lower limit of the NJL model results. Data points\cite{b11,b11a} are the experimental results of $S\sigma$ and $\kappa\sigma^2$($0-5\%$, $30-40\%$ and $60-80\%$ centralities).}
\label{njl}
\end{figure}

Firstly, we adopt a parameterized scheme by J. Cleymans and H. Osechler, etc in Ref.\cite{b10}. In the paper they proposed that $T=a-b\mu_B^2-c\mu_B^4$ and $\mu_B=d/(1+e\sqrt{S_{NN}})$, where $a=0.166\pm0.002 GeV$,$b=0.139\pm0.016 GeV^{-1}$,$c=0.053\pm0.021 GeV^{-3}$,$d=1.308\pm0.028 GeV$, $e=0.273\pm0.008 GeV^{-1}$ and $\mu_B=3\mu$ is the baryon chemical potential. Since each parameter(a,b,c,d,e) has the value range, we demonstrate the upper and the lower limit of the NJL model results in Fig.\ref{njl} as solid lines. At the same time, the $S\sigma$ and $\kappa\sigma^2$ for net-proton distributions\cite{b11,b11a} of three different centralities($0-5\%$, $30-40\%$ and $70-80\%$) at RHIC are shown as filled squares, triangles and circles respectively.

It is shown that the upper and the lower limit of the NJL model in Fig.\ref{njl} demonstrate an obvious fluctuation when the collision energies $\sqrt{S_{NN}}$ is less than $10 GeV$. And it is reduced when $\sqrt{S_{NN}}$ is more than $20 GeV$. That is to say, the NJL results is sensitive to the parameter of $T$ and $\mu$ when $\sqrt{S_{NN}}$ is less than $10 GeV$. Meanwhile, the NJL results are comprehensively less than the experimental data of both $S\sigma$ and $\kappa\sigma^2$, except the data of $\kappa\sigma^2$ when $\sqrt{S_{NN}}$ is more than $30 GeV$.
\begin{figure}
\centering
\includegraphics[width=1.0\linewidth]{m1.eps}
\caption{$S\sigma$ of the NJL model and the experiments. The NJL results demonstrate as empty squares, triangles and circles($0-5\%$, $30-40\%$ and $60-80\%$ centralities respectively). And experimental data in Ref.\cite{b11,b11a} are shown as the filled patterns($0-5\%$, $30-40\%$ and $70-80\%$ centralities).}
\label{m1}
\end{figure}

Secondly, the correspondence of the freeze-out $T$ and $\mu$ in the NJL results to the collision energies $\sqrt{S_{NN}}$ comes from the experimental results\cite{b12}. Since, in Ref.\cite{b12}, the correspondence is shown in three different centralities($0-5\%$, $30-40\%$ and $60-80\%$), then the NJL results of $S\sigma$ and $\kappa\sigma^2$ are also calculated in the three centralities and demonstrate as empty squares, triangles and circles in Fig.\ref{m1} and Fig.\ref{m2} respectively. And as in Fig.\ref{njl}, the experimental results in Ref.\cite{b11,b11a} are shown as filled patterns. It could be seen that the NJL result fit the experimental data better at large $\sqrt{S_{NN}}$ in Fig.\ref{m1} and Fig.\ref{m2}. And just as in Fig.\ref{njl}, the NJL results fluctuate obviously when the collision energies $\sqrt{S_{NN}}$ is less than $20 GeV$.
\begin{figure}
\centering
\includegraphics[width=1.0\linewidth]{m2.eps}
\caption{$\kappa\sigma^2$ of the NJL model and the experiments. The NJL results demonstrate as empty squares, triangles and circles($0-5\%$, $30-40\%$ and $60-80\%$ centralities respectively). And experimental data in Ref.\cite{b11} are shown as the filled patterns($0-5\%$, $30-40\%$ and $70-80\%$ centralities).}
\label{m2}
\end{figure}

Now it might be concluded that the $S\sigma$ and $\kappa\sigma^2$ results of two-flavor NJL model fit the experiments better at large collision energies and show an obvious fluctuation at small collision energies. And this fluctuation reflect that the NJL model is sensitive to the parameters of $T$ and $\mu$ at small collision energies.

\section{Summary}
\label{four}

The moments of the distributions of conserved quantities, for example net-baryon number, in the relativistic heavy ion collisions are related to the correlation length $\xi$ of the system\cite{b3}, they are believed to be the good signatures of QCD phase transition and its CP. In order to cancel the volume term, the products of the moments, $S\sigma$ and $\kappa\sigma^2$, are constructed as the experimental observables. Since the moments of the baryon number are proportional to the various order baryon-number susceptibilities and  the quark-number density is determined by the corresponding dressed quark propagator only\cite{b6}, then by a reasonable dressed quark propagator we can calculate the moments of quark number at finite temperature and chemical potential.

We obtain the dressed quark propagator from the NJL model\cite{b7}, which is widely adopted to study the dynamical chiral symmetry breaking(DCSB) and the interaction between hadrons. When corresponding the freeze-out temperature $T$ and the quark chemical potential $\mu$ in the NJL results to the collision energies $\sqrt{S_{NN}}$ in the experiments, we choose the parameterized scheme by J. Cleymans and H. Osechler, etc and the RHIC results\cite{b12} at the same time. It is found that the $S\sigma$ and $\kappa\sigma^2$ results of two-flavor NJL model fit the experiments better at large collision energies and show an obvious fluctuation at small collision energies. And this fluctuation reflect that the NJL model is sensitive to the parameters of $T$ and $\mu$ at small collision energies.

\section*{Acknowledgements}

This work is supported by the University Natural Science Foundation of JiangSu Province China (under Grants No.17KJB140003 ).

\begin{appendix}
\section{}
The quark-number density susceptibilities($\chi^{(2)}$, $\chi^{(3)}$ and $\chi^{(4)}$) of a single color and flavor are:
\begin{equation}
\chi^{(2)}=\frac{\partial n}{\partial \mu}=-\frac{1}{\pi^2}\int dp\cdot p^2(\frac{e^s\frac{\partial s}{\partial \mu}}{(e^s+1)^2}-\frac{e^{s'}\frac{\partial s'}{\partial \mu}}{(e^{s'}+1)^2}),
\end{equation}

\begin{eqnarray}
\begin{split}
\chi^{(3)}=\frac{\partial^2 n}{\partial \mu^2}=&-\frac{1}{\pi^2}\int dp\cdot p^2(\frac{e^s(\frac{\partial^2 s}{\partial \mu^2}+(\frac{\partial s}{\partial \mu})^2)+e^{2s}(\frac{\partial^2 s}{\partial \mu^2}-(\frac{\partial s}{\partial \mu})^2)}{(e^s+1)^3}\\
&-\frac{e^{s'}(\frac{\partial^2 s'}{\partial \mu^2}+(\frac{\partial s'}{\partial \mu})^2)+e^{2s'}(\frac{\partial^2 s'}{\partial \mu^2}-(\frac{\partial s'}{\partial \mu})^2)}{(e^{s'}+1)^3}),
\end{split}
\end{eqnarray}

\begin{eqnarray}
\begin{split}
\chi^{(4)}=\frac{\partial^3 n}{\partial \mu^3}=&-\frac{1}{\pi^2}\int dp\cdot p^2(\frac{1}{(e^s+1)^4}(e^s(\frac{\partial^3 s}{\partial \mu^3}+3\frac{\partial s}{\partial \mu}\frac{\partial^2 s}{\partial \mu^2}+(\frac{\partial s}{\partial \mu})^3)\\
&+e^{2s}(2\frac{\partial^3 s}{\partial \mu^3}-4(\frac{\partial s}{\partial \mu})^3)+e^{3s}(\frac{\partial^3 s}{\partial \mu^3}-3\frac{\partial s}{\partial \mu}\frac{\partial^2 s}{\partial \mu^2}+(\frac{\partial s}{\partial \mu})^3))\\
&-\frac{1}{(e^{s'}+1)^4}(e^{s'}(\frac{\partial^3 s'}{\partial \mu^3}+3\frac{\partial s'}{\partial \mu}\frac{\partial^2 s'}{\partial \mu^2}+(\frac{\partial s'}{\partial \mu})^3)\\
&+e^{2s'}(2\frac{\partial^3 s'}{\partial \mu^3}-4(\frac{\partial s'}{\partial \mu})^3)+e^{3s'}(\frac{\partial^3 s'}{\partial \mu^3}-3\frac{\partial s'}{\partial \mu}\frac{\partial^2 s'}{\partial \mu^2}+(\frac{\partial s'}{\partial \mu})^3))),
\end{split}
\end{eqnarray}

where $s=(\epsilon-\mu)/T$ and $s'=(\epsilon+\mu)/T$, that is to say
\begin{eqnarray}
\begin{split}
\frac{\partial s}{\partial \mu}=\frac{1}{T}(\frac{\partial \epsilon}{\partial \mu}-1)=\frac{1}{T}(\frac{M}{\epsilon}\frac{\partial M}{\partial \mu}-1)\\
\frac{\partial s'}{\partial \mu}=\frac{1}{T}(\frac{\partial \epsilon}{\partial \mu}+1)=\frac{1}{T}(\frac{M}{\epsilon}\frac{\partial M}{\partial \mu}+1),
\end{split}
\end{eqnarray}

\begin{equation}
\frac{\partial^2 s}{\partial \mu^2}=\frac{\partial^2 s'}{\partial \mu^2}=\frac{1}{T}\frac{\partial^2 \epsilon}{\partial \mu^2}=\frac{1}{T}(\frac{1}{\epsilon}((\frac{\partial M}{\partial \mu})^2+M\frac{\partial^2 M}{\partial \mu^2})-\frac{1}{\epsilon^3}(M\frac{\partial M}{\partial \mu})^2),
\end{equation}

\begin{eqnarray}
\begin{split}
\frac{\partial^3 s}{\partial \mu^3}&=\frac{\partial^3 s'}{\partial \mu^3}=\frac{1}{T}\frac{\partial^3 \epsilon}{\partial \mu^3}\\
&=\frac{1}{T}(\frac{1}{\epsilon}(3\frac{\partial^2 M}{\partial \mu^2}\frac{\partial M}{\partial \mu}+M\frac{\partial^3 M}{\partial \mu^3})-\frac{3M}{\epsilon^3}\frac{\partial M}{\partial \mu}((\frac{\partial M}{\partial \mu})^2+M\frac{\partial^2 M}{\partial \mu^2})+\frac{3}{\epsilon^5}(M\frac{\partial M}{\partial \mu})^3),
\end{split}
\end{eqnarray}
and the value of $\frac{\partial M}{\partial \mu}$, $\frac{\partial^2 M}{\partial \mu^2}$ and $\frac{\partial^3 M}{\partial \mu^3}$ are obtained from the iterative equations
\begin{equation}
\frac{\partial M}{\partial \mu}=-2G_1\frac{\partial <\bar{\psi}\psi>}{\partial \mu}-4G_2<\bar{\psi}\psi>\frac{\partial <\bar{\psi}\psi>}{\partial \mu},
\end{equation}
\begin{equation}
\frac{\partial^2 M}{\partial \mu^2}=-2G_1\frac{\partial^2 <\bar{\psi}\psi>}{\partial \mu^2}-4G_2<\bar{\psi}\psi>\frac{\partial^2 <\bar{\psi}\psi>}{\partial \mu^2}-4G_2(\frac{\partial <\bar{\psi}\psi>}{\partial \mu})^2,
\end{equation}
\begin{equation}
\frac{\partial^3 M}{\partial \mu^3}=-2G_1\frac{\partial^3 <\bar{\psi}\psi>}{\partial \mu^3}-4G_2<\bar{\psi}\psi>\frac{\partial^3 <\bar{\psi}\psi>}{\partial \mu^3}-12G_2\frac{\partial <\bar{\psi}\psi>}{\partial \mu}\frac{\partial^2 <\bar{\psi}\psi>}{\partial \mu^2}.
\end{equation}

\end{appendix}


\section{Introduction}
\label{sect:intro}
Almost all stars ever formed were likely born in clustered environments \citep[see review by][and references therein]{Zwart10}. Stars inside a cluster are subjected to two-body relaxation which over time leads to more-massive stars sinking to the central regions of the cluster while less-massive stars occupy a larger volume, many of which gradually evaporate to supply the field population of a galaxy. This effect is called mass segregation and has been observed in both galactic globular \citep{hoerner57} and open \citep{Mathieu84} clusters, investigated both through the slope of the mass-function of the cluster members in the core and outskirts \citep[e.g.][]{Bhattacharya17b}, as well as through the radial distribution of massive objects, such as blue stragglers \citep[e.g.][]{Bhattacharya19,Vaidya20,Rao21}. Alongside this internal effect, Galactic potential perturbs star clusters leading to the formation of tidal structures. Such tidal structures are observed in both globular clusters (e.g. Palomar~5 - \citealt{Odenkirchen01}; NGC~5466 - \citealt{Belokurov06}) and open clusters (e.g. Berkeley~17 - \citealt{Chen04,Bhattacharya17b}; NGC~6791 - \citealt{Dalessandro15}).

Galactic open clusters, in particular, are more abundant in the disc where they are additionally subjected to external gravitational effects such as interactions with giant molecular clouds and passage through spiral arms \citep{Spitzer58}. From an observational standpoint, the identification of open cluster morphology, despite the use of novel techniques \citep[e.g.][]{Bhattacharya17a}, has long since relied on stellar density determination against strong field contamination in photometric observations. Low-density features such as tidal tails have thus been particularly difficult to identify with little success \citep{Chen04,Dalessandro15}. Reliable astrometry (proper motion and parallax) from Gaia Data Release 2 \citep[Gaia DR2;][]{Gaia18} and recently Gaia Early Data Release 3 \citep[Gaia EDR3;][]{Gaia20} has allowed for improved field-star decontamination of many nearby clusters through novel membership determination techniques \citep{Agarwal21}. This has resulted in the recent discovery of a few open clusters {with extended stellar coronae \citep[e.g.][]{Carrera19,Meingast21} as well as those with} tidal features -- Hyades: \citet{Roser19}; Coma~Berenices: \citet{Tang19}; {Ruprecht 147: \citet{Yeh19};} Praesepe: \citet{Roser19b}; NGC~2506: \citet{Gao20}; M~67: \citet{Gao20b}; UBC~274: \citet{Castro20}; {Alpha Persei: \citet{Nikiforova20}}.

NGC~752 (cluster center at RA=29.156 deg and DEC=37.809 deg; \citealt{Agarwal21}) is an intermediate age ($\sim1.5$ Gyr old; \citealt{Agueros18}), nearby ($\sim450$ pc) open cluster which was among the clusters whose membership determination was carried out by \citet{Agarwal21} using Gaia DR2 data with their novel membership determination algorithm, ML-MOC. Within their search radius of 2.5$^{\circ}$ around its center, \citet{Agarwal21} found hints of tidal tails in NGC~752. \citet{Hu21} also identify the elongated morphology of the peripheral regions of NGC~752 from the Gaia DR2 data using the members identified by \citet{cantat18}.

In this paper, we report the discovery of the tidal tails of NGC~752 within a large spatial region around the cluster. The data and membership determination are presented in Section~\ref{sect:data}. The analysis of the mass function and tidal structure of NGC~752 is presented in Section~\ref{sect:analysis}. Finally, we discuss {the orbit of the cluster and its dissolution} in Section~\ref{sect:disc}.

\begin{figure*}
	\includegraphics[width=\textwidth]{mlmoc.png}
    \caption{{The spatial, proper-motion and parallax distributions of the \textit{Sample sources} and the member stars identified by ML-MOC.}}
    \label{fig:mlmoc}
\end{figure*}

\begin{table*}
\caption{{Gaia EDR3 members identified in NGC~752 using ML-MOC. The full table will be available in the online version of the journal.}}
\centering
\adjustbox{max width=\textwidth}{
\begin{tabular}{ccccccccccc}
  \hline
    GaiaID & RA & DEC & $\omega$ & $\mu_{\alpha *}$ & $\mu_{\delta}$ & G & BP & RP & RV & $\rm p_{memb}$\\
     & deg & deg & mas & mas/yr & mas/yr & mag & mag & mag & km/s & \\
  \hline
342914075859711744 & 29.23240573 & 37.79978973 & 2.288 & 9.803 & -11.878 & 11.1 & 11.3 & 10.75 & 6.91 & 0.98 \\
342912284857245952 & 29.15366539 & 37.75349579 & 2.272 & 9.892 & -12.043 & 13.64 & 14.04 & 13.07 & -- & 0.97 \\
342914556896034176 & 29.22108601 & 37.8692465 & 2.261 & 9.462 & -11.507 & 10.15 & 10.38 & 9.77 & 18.0 & 0.92 \\
342915411593428736 & 29.16345994 & 37.86139137 & 2.267 & 9.721 & -11.886 & 10.29 & 10.48 & 9.96 & 8.23 & 0.98 \\
342912254793586688 & 29.13259369 & 37.74945795 & 2.141 & 9.794 & -11.62 & 16.2 & 17.08 & 15.27 & -- & 0.68 \\
342913178210782080 & 29.11199587 & 37.79182891 & 2.231 & 9.854 & -12.06 & 17.83 & 19.17 & 16.72 & -- & 0.95 \\
342909845316110336 & 29.31580848 & 37.73702508 & 2.363 & 9.887 & -11.726 & 18.3 & 19.57 & 17.16 & -- & 0.91 \\
342920501130956800 & 29.26040912 & 37.88539331 & 2.235 & 9.815 & -11.984 & 11.21 & 11.43 & 10.84 & 5.1 & 0.97 \\
342920501130956928 & 29.26070082 & 37.88616266 & 2.219 & 10.006 & -11.428 & 13.9 & 14.2 & 13.24 & -- & 0.92 \\
342921257045010432 & 29.22252061 & 37.8972341 & 2.169 & 9.499 & -11.925 & 14.92 & 15.52 & 14.18 & -- & 0.9 \\
  \hline
\end{tabular}
\label{table:data}
}
\end{table*}

\begin{figure*}
	\includegraphics[width=\textwidth]{spat.png}
    \caption{Spatial distribution of the identified members in NGC~752. (a) Coloured by their $\rm p_{memb}$ values. (b) Coloured by their RV values where available and open grey circles where unavailable.  (c) Matched previously-identified members. (d) Coloured by their $\rm A_{V}$ values.}
    \label{fig:spat}
\end{figure*}

\section{Data and membership selection}
\label{sect:data}
Gaia EDR3 provides the positions, trigonometric parallaxes ($\rm\omega$), and proper motions (PM) as well as photometry in three broad-band filters (G, BP, and RP) for almost 1.5 billion stars. It improves upon the Gaia DR2 with the proper motion uncertainty reduced on average by nearly a factor of two while parallax measurements are more accurate by $\sim$20\% \citep{Gaia20}. This allows for reliable photometry down to G$\sim$21 mag. Gaia EDR3 photometry is essentially complete down to G=17 mag although the completeness may be lower in regions of high stellar density such as those of globular clusters \citep{Fabricius20}. Radial velocity measurements in Gaia EDR3 are nearly the same as in Gaia DR2. Further details on Gaia EDR3 are available at \url{https://www.cosmos.esa.int/web/gaia/earlydr3}.

The membership selection in NGC~752 from the Gaia EDR3 data-set is carried out  in a region of radius 5$^{\circ}$ around its center. {This sample is termed \textit{All sources}.} We use ML-MOC \citep{Agarwal21} to identify cluster members using the proper motion and parallax information. While details of the membership selection are described in \citet{Agarwal21}, we briefly state it here. ML-MOC is based on the k-nearest neighbour algorithm \citep[kNN;][]{kNN} and the Gaussian mixture model \citep[GMM;][]{GMM}. It identifies cluster members in the PM--$\rm\omega$ parameter space, independent of the spatial density of the cluster, thereby allowing for the identification of faint extended spatial structures such as tidal tails. As a first step to the membership selection, kNN is utilised to remove the obvious field stars such that the remaining member candidates have more cluster members than field stars. {This is done by applying kNN to those stars within the central region of radius 1$^{\circ}$ to calculate a broad range of PM and $\omega$ values for the cluster members. For NGC~752, we find the range of PM in RA, $\mu_{\alpha *}=$ 8.372 -- 11.372 mas/yr, PM in DEC, $\mu_{\delta *}=$ -13.948 -- -10.348 mas/yr, and $\omega=$ 1.552 -- 5.199. The stars in \textit{All sources} which fall in the aforementioned parameter ranges are termed \textit{Sample sources}. Their spatial, PM and $\omega$ distributions are shown in Figure~\ref{fig:mlmoc}. } {For the second step, a three dimensional GMM is used in the PM--$\omega$ parameter space of the \textit{Sample sources}} to distinguish between the cluster and field members, also assigning a membership probability ($\rm p_{memb}$). Those sources having $\rm p_{memb}\geq0.6$ are considered as high probability members, also including all the radial velocity members. A small number of stars having $\rm 0.2 \leq p_{memb}\leq0.6$ are also considered as members if their $\rm\omega$ values lie within the range of $\rm\omega$ values specified by the members that have $\rm p_{memb}\geq0.8$. {This results in the identification of 282 members in NGC~752 whose spatial, PM and $\omega$ distributions are also shown in Figure~\ref{fig:mlmoc}.} Note the extended tails are clearly identified in the spatial distribution. {The identified members with their Gaia EDR3 IDs, astrometry, photometry, radial velocity (RV) and $\rm p_{memb}$ have been tabulated in Table~\ref{table:data}.}

In Figure~\ref{fig:spat}a, the marked $\rm p_{memb}$ values show that the {extended} tails are present for the high probability members showing their identification as being robust. The mean distance to the cluster is found to be 443.8 pc \citep[see][]{Agarwal21} consistent with literature values \citep{Agueros18}.  The cluster members with RV data from Gaia EDR3 are shown in Figure~\ref{fig:spat}b, having <RV>=7.88 km/s and $\rm\sigma_{RV}=5.03$ km/s. Few cluster members with RV data are also identified among the {extended} features. \citet{Agueros18} had identified 258 members in NGC~752 from photometry and PM from the Tycho-Gaia Astrometric Solution (TGAS) catalog \citep{Gaia16}. Furthermore, they had consolidated membership determinations from previous studies \citep[primarily][]{Daniel94,Mermilliod98} finding 32 literature members as non-members. 144 of the stars identified by \citep{Agueros18} as members have counterparts within 1$\arcsec$ of our identified members while 14 of the literature members classified as non-members also have counterparts within 1$\arcsec$ of our identified members. The culling of non-members in our membership sample as compared to those in the literature is testament to the improved accuracy of the Gaia EDR3 astrometry and our robust membership determination from ML-MOC. Figure~\ref{fig:spat}c shows that the literature members are concentrated in the central parts of NGC~752 as the {extended} tails were not surveyed in previous spatially-restrictive studies whose membership determinations were guided by cluster radial density.

\section{Analysis}
\label{sect:analysis}

\subsection{Reddening}
\label{sect:red}
We utilise the interstellar reddening, $E(B-V)$, map provided by \citet{SF11} available through the NASA/IPAC Infrared Science Archive\footnote{\url{https://irsa.ipac.caltech.edu/applications/DUST}} to obtain the line-of-sight extinction ($\rm A_{V} = 3.1 \times E(B-V)$; \citealt{cardelli89}) towards each identified cluster member in NGC~752. Figure~\ref{fig:spat}d shows the spatial distribution of the cluster members coloured by their $\rm A_{V}$. The mean extinction $\rm <A_{V}>=0.15$mag. The extinction for each cluster member is converted to the Gaia photometric system following the extinction law by \citet{Wang19}. There is indeed some evidence of differential reddening but only a small range of extinction values is covered.

\begin{figure*}
	\includegraphics[width=\textwidth]{cmd.png}
    \caption{[Left] CMD of all the identified members in NGC~752, coloured by their stellar masses estimated from the PARSEC isochrone (black). The dashed grey line shows the $\sim$100\% completeness limit of the Gaia EDR3 data. [Right] The Mass function derived from the PARSEC iscochrones for the entire cluster. The dashed black line shows the fitted mass function where the Gaia EDR3 data is almost complete. The Poissonian uncertainties have been marked.}
    \label{fig:cmd}
\end{figure*}

\subsection{Color Magnitude Diagram and Mass Function}
\label{sect:cmd}
Correcting for extinction and reddening for each individual cluster member, we obtain the de-reddened color magnitude diagram (CMD) of the identified members shown in the left panel of Figure~\ref{fig:cmd}. The PARSEC stellar evolution isochrone\footnote{\url{http://stev.oapd.inaf.it/cgi-bin/cmd}} \citep{Bressan12,Chen14} for the Gaia filters has been plotted for the obtained distance corresponding to an age, t~=~1.5 Gyr and $\rm [M/H] = -0.04$, in line with fitted literature values of age (t~$=1.45\pm0.05$ Gyr) from multi-wavelength photometry by \citet{Twarog15} and metallicity ($\rm <[M/H]>=-0.07\pm0.04$) from spectroscopy of the red giant stars by \citet{Topcu15}.

\begin{figure}
	\includegraphics[width=\columnwidth]{spat_tidal.png}
    \caption{Spatial distribution of the identified members in NGC~752 with the intra-tidal (blue), extra-tidal (red), Tail~1 (encircled in yellow) and Tail~2 (encircled in green) members marked. }
    \label{fig:spat_tidal}
\end{figure}

Stellar mass of each star was estimated from the PARSEC isochrone and is shown in the left panel of Figure~\ref{fig:cmd}. The total observed mass of the cluster members is 284 M$_{\sun}$. The mass function of all the identified members in NGC~752 is shown in the right panel of Figure~\ref{fig:cmd}. The mass function is fitted above the stellar mass value ($\approx0.5$ M$_{\sun}$) corresponding to the Gaia EDR3 completeness limit of G=17 mag. Using the relation log(dN/dM)~=~-(1~+~$\chi$)~$\times$~log(M)~+~constant, where dN represents the number of stars in a mass bin dM with central mass M, we obtain the slope of the mass function for NGC~752, $\chi=-1.26\pm0.07$. Such a negative value of $\chi$, as compared to its value of $\chi=1.37$ as derived by \citep{salpeter55} for the solar neighbourhood, is indicative of a dissolving cluster that has undergone significant mass segregation. Extrapolating the fitted mass function to the stellar hydrogen burning limit of 0.08 M$_{\sun}$, we obtain the total mass of the cluster, M$\rm_{C}=297\pm10$ M$_{\sun}$.

\begin{figure*}
	\includegraphics[width=0.9\textwidth]{mf_tidal.png}
    \caption{[Left] De-reddened CMD of all the identified members in NGC~752 in the four different regions. The black line shows the PARSEC isochrone corresponding to the age and metallicity of NGC~752. [Right] The Mass function derived from the PARSEC iscochrones for each region. The dashed black line shows the fitted mass function where the Gaia EDR3 data is almost complete. The Poissonian uncertainties have been marked.}
    \label{fig:mf_tidal}
\end{figure*}

\subsection{Tidal Structure}
\label{sect:tidal}
While extended tails are clearly visible in the spatial distribution of NGC~752 (Figure~\ref{fig:spat}), we need to estimate the tidal radius (r$\rm_{t}$) of the cluster in order to ascertain the {extra-tidal} nature of these tails. As estimating the tidal radius from fitting the King function \citep{King66} to the radial density of the cluster \citep[as carried out by][]{Agarwal21} is not appropriate for a dissolving cluster with extended structures, we calculate r$\rm_{t}$ following \citet{Pinfield98} as:
\begin{equation}
    $$\rm r_{t}^{3} = \rm \frac{G M_{C}}{2(A-B)^{2}}$$
\end{equation}
where G is the gravitational constant, and A and B are the
Oort constants (A$=15.3\pm0.4$ km/s/kpc, B$=-11.9\pm0.4$ km/s/kpc; \citealt{Bovy17}). We calculate r$\rm_{t}=1.2123^{\circ}=9.52$~pc at the distance of NGC 752. The identified members within and beyond r$\rm_{t}$ are classified as intra-tidal and extra-tidal members respectively  (Figure~\ref{fig:spat_tidal}). Amongst the extra-tidal members, we further note those stars belonging to either of the two tails (1 \& 2; see Figure~\ref{fig:spat_tidal}). At the distance of NGC~752, tail 1 and 2 extend out to $\sim$34 and $\sim$36.5 pc respectively from the center of the cluster.

We further obtain the de-reddened CMDs for the four sub-regions of NGC~752 (intra-tidal, extra-tidal, Tail 1 and Tail 2) as shown in the left panel of Figure~\ref{fig:mf_tidal}. The mass functions for each of the sub-regions are obtained and fitted as shown in the right panel of Figure~\ref{fig:mf_tidal}. The fitted slopes of the mass functions are $\chi=-1.44\pm0.12$, $-0.69\pm0.49$, $-0.87\pm0.57$ and $-0.63\pm0.64$ for the intra-tidal, extra-tidal, Tail 1 and Tail 2 respectively. The mass function slope is most negative in the intra-tidal region, as expected in case of mass segregation, while it is relatively less-steep outside r$\rm_{t}$ despite the large uncertainty.

\section{Discussion}
\label{sect:disc}

\begin{figure*}
	\includegraphics[width=\textwidth]{orbit.png}
    \caption{{Tri-dimensional projections in the galactocentric coordinate system of the identified members in NGC~752 with the intra-tidal (blue), extra-tidal (red), Tail~1 (encircled in yellow) and Tail~2 (encircled in green) members marked. The red arrow follows the cluster orbit. }}
    \label{fig:orbit}
\end{figure*}

\subsection{Cluster orbit}
\label{sect:orbit}

{To further validate the tidal origin of the extended tails, we explore their position with respect to the cluster orbit in a galactocentric coordinate system, similar to that carried out by \citet{Tang19} for Coma Berenices. The Galactocentric coordinates have a positive x direction pointing from the position of the Sun projected to the Galactic midplane to the Galactic center (GC); the y-axis points toward \textit{l}=90$^{\circ}$ and the z-axis roughly points toward \textit{b}=90$^{\circ}$.  Distance to the GC is 8.3 kpc \citep{Reid04}. Figure~\ref{fig:orbit} shows the X-Y (left), Y-Z (middle) and X-Z (right) projections of the cluster members.}

{We compute the average position of all member stars as the cluster mean position giving X=-8595.62 pc, Y=276.26 pc, Z=-152.36 pc, and its mean space velocity (for only those sources having RV measurements) as U=-5.13 km/s, V=226.44 km/s, W=-12.03 km/s. The aforementioned are utilised to calculate the cluster orbit using the python package galpy \citep{Bovy15} which accounts for solar motion and includes the Milky Way potential model \textit{MWPotential2014} (A fitted Galactic potential model incorporating the bulge, disc and halo). Orbits are integrated 100 Myr backward and 100 Myr forward in time.}

{Figure~\ref{fig:orbit} shows the orbit motion of NGC~752 in the three galactocentric projections. The distribution of the extended tails along the cluster orbit is most prominent in the Y-Z plane, similar to that found by \citet{Roser19}, \citet{Tang19} and \citet{Yeh19} for their studied clusters. This is a testament to the tidal nature of the extended tails of NGC~752 with Tails 1 \& 2 being the leading and trailing tails respectively.} Through reliable membership determination up to large radii in conjunction with reliable astrometry from Gaia EDR3, NGC~752 is now added to a growing list of dissolving open clusters exhibiting tidal features.

\subsection{Initial cluster mass}
\label{sect:inimass}
The mass-segregation of NGC~752, evident from the negative value of its mass function slope, and extended tidal features mark NGC~752 as a telltale sign of its dissolution. Environmental effects like tidal shocks due to close interactions with giant molecular clouds, spiral arms, the Galactic disc and, in
general, interactions with the Galactic tidal field or internal dynamical effects like two-body relaxation lead to mass-loss in an open cluster \citep{Dalessandro15}.

\citet{Lamers05} provide an analytical prescription to estimate the mass lost by an open cluster due to the Galactic tidal field. The initial mass ($\rm M_{i}$) is estimated with the following relation:
\begin{equation}
    $$\rm M_{i} = \rm \bigg[\bigg(\frac{M_{C}}{M_{\sun}}\bigg) ^{\gamma} + \frac{\gamma t}{t_{0}} \bigg]^{\frac{1}{\gamma}} [1-q_{ev}(t)]^{-1}$$
    \label{eq:mini}
\end{equation}
where t$\rm_{0}$ is the dissolution time-scale parameter. By comparing the distribution of mass and age of OCs in the solar neighbourhood with theoretical predictions, \citet{Lamers05} obtained t$\rm_{0}=3.3_{-1.0}^{1.4}$ Myr. $\rm \gamma$ is a dimensionless index which depends on the cluster initial
density distribution. We adopt $\rm \gamma=0.62$ which is the typical value for open clusters \citep{Lamers05, Dalessandro15}. The function q$\rm_{ev}$(t) describes the mass-loss due to stellar evolution and can be approximated by the following analytical formula:
\begin{equation}
    $$\rm q_{ev}(t) = \rm (log_{10}(t)-a)^{b}+c$$
    \label{eq:qev}
\end{equation}
where a = $7.03$, b = $0.26$ and c = $-1.80$ for the metallicity of NGC 752 \citep[see Table 1 in][]{Lamers05}. Plugging the values in to Equations~\ref{eq:mini}~\&~~\ref{eq:qev}, we find $\rm M_{i}=11087_{-4666}^{8922}$, i.e,  $\rm M_{i}~=~0.64~-2~\times~10^{4}~M_{\sun}$. Thus, NGC~752 is a descendent of a Young Massive Cluster (which have stellar masses $\gtrsim 10^{4} \rm M_{\sun}$) such as those observed in the Milky Way and other galaxies \citep[see review by][]{Zwart10} that has lost 95.2--98.5\% of its mass to the Galactic field.

\section*{Acknowledgements}
{This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). We thank the anonymous referee for their comments.} This research made use of Astropy-- a community-developed core Python package for Astronomy \citep{Astropy13}, SciPy \citep{scipy}, NumPy \citep{numpy} and Matplotlib \citep{matplotlib}. This research also made use of NASA’s Astrophysics Data System (ADS\footnote{\url{https://ui.adsabs.harvard.edu}}).

\section*{Data Availability}
The data underlying this article are publicly available at \url{https://archives.esac.esa.int/gaia}. {The 282 identified members of NGC~752 along with their Gaia EDR3 identifiers and photometric and astrometric data have been presented in Table~\ref{table:data}}.

\bibliographystyle{mnras}