\section{Introduction} Global Value Numbering is an important static analysis to detect equivalent expressions in a program. Equivalences are detected by assigning \emph{value numbers} to expressions. Two expressions are assigned the same value number if they could be detected as equivalent. The seminal work on GVN by Kildall \cite{Kildall1973} detects all \emph{Herbrand equivalences} \cite{Ruething1999} in non-SSA form of programs using the powerful concept of \emph{structuring} but takes exponential time. Efforts were made to improve on efficiency in detecting equivalences. However the algorithms are either as precise as Kildall's \cite{Saleena2014} or efficient \cite{Ruething1999, Alpern1988, Gulwani2007} but not both.

The strive for combining precision with efficiency has motivated our work in this area. We propose an iterative data-flow analysis GVN algorithm to detect redundancies in SSA form of programs that is precise as Kildall's and efficient (i.e.\,take only polynomial time). As in a data-flow analysis problem, the central challenge is to define a \emph{join} operation to detect all equivalences at a join point in polynomial time such that any later occurrences of redundant expressions could be detected. We introduce the novel concept of \emph{value $\phi$-function} for this purpose.

\section{Terminology}
\paragraph{Program Representation} Input to our algorithm is the Control Flow Graph (CFG) representation of a program in SSA. The graph has empty \emph{entry} and \emph{exit} blocks. Other blocks contain assignment statements of the form $x = e$, where $e$ is an expression which is either a constant, a variable, or of the form $x \oplus y$ such that $x$ and $y$ are constants or variables and $\oplus$ is a generic binary operator. An expression can also be of the form $\phi_{k}(x, y)$, called \emph{$\phi$-functions}, where $x$ and $y$ are variables and $k$ is the block in which it appears. We assume a block can have at most two predecessors and a block with exactly two predecessors is called \emph{join} block. The input and output points of a block are called \emph{in} and \emph{out} points, respectively, of the block. The \emph{in} point of a join block is called \emph{join point}. We may omit the subscript $k$ in $\phi_{k}$ when the join block is clear from the context. In the CFGs we draw, $\phi$-functions appear in join blocks. But for clarity in explaining some of our concepts we assume $\phi$-functions are transformed to \emph{copy} statements and appended to appropriate predecessors of the join block.

\paragraph{Equivalence} Two expressions $e_{1}$ and $e_{2}$ are \emph{equivalent}, denoted $e_{1} \equiv e_{2}$, if they will have the same value whenever they are executed. Two expressions in a path are said to be \emph{equivalent in the path} if they are equivalent in that path. We detect only Herbrand equivalences \cite{Ruething1999} which is equivalence among expressions with same operators and corresponding operands being equivalent.

\section{Basic Concept}
Our main goal is to detect equivalences with a view to detecting redundancies in a program in polynomial time. We introduce the concept of \emph{value $\phi$-function} for the purpose which is explained in this section followed by our method to detect redundancies.

\subsection{Value $\phi$-function} Consider the simple code segment in Fig.\,1(a). Here irrespective of the path taken $x_{1}+y_{1}$ is equivalent to $a_{1}+b_{1}$. In terms of the variables being assigned to, we can say $z_{1}$ is equivalent to same variable $c_{1}$.
\begin{figure}[ht]
        \centering
        \begin{subfigure}[b]{0.3\textwidth}
                \includegraphics[width=17mm, height=15mm]{dc-figure0.eps}
                \caption{Linear program}
                \label{fig:1a}
        \end{subfigure}
        \begin{subfigure}[b]{0.4\textwidth}
                \includegraphics[width=55mm, height=20mm]{dc-figure1.eps}
                \caption{Program with branches}
                \label{fig:1b}
        \end{subfigure}
        \caption{Concept of value $\phi$-function}\label{fig:vpf}
\end{figure}
Now consider the code segment in Fig.\,1(b). Depending on the path taken expression $x_{3}+y_{3}$ is equivalent to either $x_{1}+y_{1}$ or $x_{2}+y_{2}$. In terms of the variables being assigned to, we can say $w_{3}$ is equivalent to merge of different variables -- $p_{1}$ and $q_{2}$. Inspired by the notion of $\phi$-function, we can say $w_{3}$ is equivalent to $\phi(p_{1}, q_{2})$. This notion of $\phi$-function is an extended notion of $\phi$-function as seen in the literature. In the literature, a $\phi$-function has different subscripted versions of the same non-SSA variable, say $\phi(x_{1}, x_{2})$. To express such equivalences, we introduce the concept of \emph{value $\phi$-function} similar to the concept of \emph{value expression} \cite{Saleena2014}.

\paragraph{Value $\phi$-function} A \emph{value $\phi$-function} is an abstraction of a set of equivalent $\phi$-functions (including the extended notion of $\phi$-function). Let $v_{i}$, $v_{j}$ be value numbers and \emph{vpf} be a value $\phi$-function. Then $\phi_{k}(v_{i}, v_{j})$, $\phi_{k}(\emph{vpf}, v_{j})$, $\phi_{k}(v_{i}, \emph{vpf})$, and $\phi_{k}(\emph{vpf}, \emph{vpf})$ are \emph{value $\phi$-functions}.

\paragraph{Partition} A partition at a point represents equivalences that hold in the paths to the point. An equivalence class in the partition has a value number and elements like variables, constant, and value expression. It is also annotated with a value $\phi$-function when necessary. The notation for a partition is similar to that in \cite{Saleena2014} except that a class can be annotated with value $\phi$-function.

\section{Proposed Method}
Using the concept of value $\phi$-function we propose an iterative data-flow analysis algorithm to compute equivalences at each point in the program. The two main tasks in this algorithm are \emph{join} operation and \emph{transfer function}:

\subsection{\emph{Join} operation.} A \emph{join} operation detects equivalences that are common in all paths to a join point. The join is conceptually a class-wise intersection of input partitions. Let $C_{1}$ and $C_{2}$ be two classes, one from each input partition. If the classes have same value number then the resulting class $C$  is intersection of $C_{1}$ and $C_{2}$. If the classes have different value numbers, say $v_{1}$ and $v_{2}$ respectively, then common equivalences are found by intersection of $C_{1}$ and $C_{2}$. The common equivalences obtained are actually a merge of different variables, which is indicated by the difference in value numbers and hence class $C$ is annotated with $\phi(v_{1}, v_{2})$. Now if the classes have different value expressions, say $v_{m}+v_{n}$ and $v_{p}+v_{q}$ respectively, the value expressions may be merged to form a resultant value expression say $v_{i}+v_{j}$. Value expressions $v_{m}+v_{n}$ and $v_{p}+v_{q}$ are merged to get $v_{i}+v_{j}$ by recursively merging classes of $v_{m}$ and $v_{p}$ to get class of $v_{i}$ and classes of $v_{n}$ and $v_{q}$ to get class of $v_{j}$ \cite{Saleena2014}. But merging the value expressions can lead to exponential growth of resulting partition \cite{Gulwani2007}. We do not merge different value expressions now instead merge them at a point where an expression represented by $v_{i}+v_{j}$ actually occurs in the program. This merge is achieved simply by detecting equivalence of $v_{i}+v_{j}$ with $\phi(v_{1}, v_{2})$ and is done during application of transfer function.

\paragraph{Example} Let us now consolidate the concept of join using an example. Consider the case of applying join on partitions $P_{1} = \{v_{1}, x_{1}, x_{3} | v_{2}, y_{1}, y_{3}, v_{1}+1 | v_{3}, z_{1}, z_{3}\}$ and $P_{2} = \{v_{4}, x_{2}, x_{3} | v_{5}, y_{2}, y_{3} | v_{6}, z_{2}, z_{3}, v_{4}+1 \}$. In the classes with value numbers $v_{1}$ in $P_{1}$ and $v_{4}$ in $P_{2}$ there is only one common variable $x_{3}$ and this will appear in a class in the resulting partition $P_{3}$. Since the two classes in $P_{1}$ and $P_{2}$ have different value numbers $v_{1}$ and $v_{4}$, respectively, the resulting class is annotated with value $\phi$-function $\phi(v_{1}, v_{4})$. The class is assigned a new value number, say $v_{7}$. The resulting class is $| v_{7}, x_{3} : \phi(v_{1}, v_{4}) |$. Now consider the classes with value numbers $v_{2}$ in $P_{1}$ and $v_{6}$ in $P_{2}$.  There are no obvious common equivalences in the classes and we don't merge the different value expressions now. Hence no new class is created. Similar strategies are adopted in detecting common equivalences in other pairs of classes one each from $P_{1}$ and $P_{2}$. The resulting partition $P_{3}$ is $\{v_{7}, x_{3} : \phi(v_{1}, v_{4}) | v_{8}, y_{3} : \phi(v_{2}, v_{5}) | v_{9}, z_{3} : \phi(v_{3}, v_{6})\}$.

 \begin{codebox}
 \Procname{$\proc{Join}(P_{1}, P_{2})$}
  \li $P \gets \{ \}$
  \li \For each pair of classes $C_{i} \in P_{1}$ and $C_{j} \in P_{2}$
  \li    \Do
             $C_{k} \gets C_{i} \cap C_{j}$ \>\>\>\>\>\>\> \small \Comment set intersection
  \li        \If $C_{k} \neq \{\}$ and $C_{k}$ does not have value number
  \li        \Then
                 $C_{k} \gets C_{k} \cup \{v_{k}, \phi_{b}(v_{i}, v_{j})\}$ \>\>\>\>\> \small \Comment $v_{k}$ is new value number
  \zi   \>\>\>\>\small \Comment $v_{i} \in C_{i}$, $v_{j} \in C_{j}, $ $b$ is join block
             \EndIf
  \li    $P \gets P \cup C_{k}$ \>\>\>\>\>\>\> \small \Comment Ignore when $C_{k}$ is empty
  \li \Return $P$
\end{codebox}
Note: We define special partition $\top$ such that $\proc{Join}(\top, P) \gets P \gets \proc{Join}(\top, P)$. We assume $\phi$-functions in a join block are transformed to copies and appended to appropriate predecessors of join block.

\subsection{Transfer Function.}\label{trfn} Given a partition $PIN_{s}$, that represents equivalences at \emph{in} point of a statement $s: x \gets e$ the transfer function computes equivalences at its \emph{out} point, denoted $POUT_{s}$. Let \emph{ve} be the value expression of $e$ computed using $PIN_{s}$. If \emph{ve} is present in a class in $PIN_{s}$, then $x$ is just inserted into corresponding class in $POUT_{s}$. Otherwise the transfer function checks whether $e$ could be expressed as a merge of variables represented by a value $\phi$-function \emph{vpf} (as illustrated below). If it is present in a class in $PIN_{s}$ then $x$, \emph{ve} are inserted into corresponding class in $POUT_{s}$. Else a new class is created in $POUT_{s}$ with new value number and $x$, \emph{ve}, \emph{vpf} are inserted into it.

For an example, consider processing the statement $w_{3} = x_{3}+ y_{3}$ as shown in code segment in Fig.\,2.
\begin{figure}[ht]
 \centering
 \includegraphics[width=115mm, height=30mm]{dc-figure2.eps}
\caption{Concept of Transfer Function}\label{fig:tf}
\end{figure}
Since value expression $v_{7}+v_{8}$ of $x_{3}+y_{3}$ is not in $PIN_{3}$, the transfer function proceeds to check whether $x_{3}+y_{3}$ is actually a merge of variables as follows:\\
$x_{3}+y_{3} \equiv v_{7}+v_{8} \equiv \phi(v_{1}, v_{4})+\phi(v_{2}, v_{5}) \equiv \phi(v_{1}+v_{2}, v_{4}+v_{5}) \equiv \phi(v_{3}, v_{6}) $.\\
This implies $x_{3}+y_{3}$ is actually a merge of variables, here $p_{1}$ and $q_{2}$. Since neither $v_{7}+v_{8}$ nor $\phi(v_{3}, v_{6})$ are present in $PIN_{3}$, a new class is created in $POUT_{3}$ with new value number say $v_{9}$ and $w_{3}$, $v_{7}+v_{8}$, and $\phi(v_{3}, v_{6})$ are inserted into it. The classes in $PIN_{3}$ are inserted as such into $POUT_{3}$. The resulting partition $POUT_{3}$ is $\{v_{7}, x_{3} : \phi(v_{1}, v_{4}) | v_{8}, y_{3} : \phi(v_{2}, v_{5}) | v_{9}, w_{3}, v_{7}+v_{8} : \phi(v_{3}, v_{6}) \}$.
\begin{codebox}
 \Procname{$\proc{transferFunction}(x \gets e, PIN_{s})$}
 \li $POUT_{s} \gets PIN_{s}$
 \li $C_{i} \gets C_{i} - \{x\}$	\>\>\>\>\>\>\>\> \small \Comment $x \in C_{i}$, a class in $POUT_{s}$
 \li \emph{ve} $\gets \proc{valueExpr}(e)$
 \li \emph{vpf} $\gets \proc{valuePhiFunc}(ve, PIN_{s})$ \>\>\>\>\>\>\>\> \small \Comment can be NULL
 \li \If \emph{ve} or \emph{vpf} is in a class $C_{i}$ in $POUT_{s}$ \>\>\>\>\>\>\>\> \small \Comment ignore \emph{vpf} when NULL
 \li \Then
       $C_{i} \gets C_{i} \cup \{x, ve\}$  \>\>\>\>\>\> \small \Comment set union
 \li \Else
       $POUT_{s} \gets POUT_{s} \cup \{v_{n}, x, ve :$ \emph{vpf}$\}$ \>\>\>\>\>\> \small \Comment $v_{n}$ is new value number
     \EndIf
 \li \Return $POUT_{s}$
\end{codebox}
The \proc{valuePhiFunc} is a recursive algorithm to compute value $\phi$-function corresponding to input value expression when possible else it returns NULL.

\subsection{Detect Redundancies.} Given partition $POUT$ at \emph{out} of statement $x = e$, expression $e$ is detected to be redundant if there exists a variable in the class of $x$ in $POUT$, other than $x$, or the class of $x$ in $POUT$ is annotated with value $\phi$-function. In the example code in Fig.\,2, consider the case of checking whether $x_{3}+y_{3}$ in the last statement $w_{3} = x_{3}+y_{3}$ is redundant. In the class of $w_{3}$ in $POUT_{3}$ (computed in previous subsection) there are no variables other than $w_{3}$. However the class is annotated with a value $\phi$-function. Hence the expression $x_{3}+y_{3}$ is detected to be redundant.

\begin{theorem}
 Two program expressions are equivalent at a point iff the iterative data-flow analysis algorithm detects their equivalence.
\end{theorem}
\begin{proof}
 This can be proved by induction on the length of a path in a program. \qed
\end{proof}

\section{Complexity Analysis}
Let there be $n$ expressions in a program. The two main operations in this iterative algorithm are join and transfer function. By definitions of \proc{Join} and \proc{transferFunction} a partition can have $O(n)$ classes. If there are $j$ join points, the total time taken by all the join operations in an iteration is $O(n.j)$. The transfer function involves constructing and then looking up for value expression or value $\phi$-function in the input partition. The transfer function of a statement takes $O(n.j)$ time. In an iteration total time taken by transfer functions is $O(n^{2}.j)$. Thus the time taken by all the joins and transfer functions in an iteration is $O(n^{2}.j)$. In the worst case the iterative analysis takes $n$ iterations and hence the total time taken by the analysis is $O(n^{3}.j)$.

\section{Conclusion}
We presented GVN algorithm using the novel concept of value $\phi$-function which made the algorithm precise and efficient.

\bibliographystyle{splncs03}


\section{Introduction}
Simple and productive methods to investigate dynamical features
of solids and molecules are offered by Time-dependent Density-functional Theory (TDDFT) \cite{GM12}.
TDDFT embodies many concepts and formal exact results, but its core
is the 1-1 correspondence \cite{RG84} between time-dependent (TD) external potentials and TD electronic densities, provided the
initial state of the system is given. Every observable of the system is expressed as a TD density-functional.
The calculation of the density is carried out by solving the TD Kohn-Sham (KS) equations \cite{P78}, which are
single-particle Schr\"odinger equations that require an approximation
to the exchange-correlation (XC) potential, a density-functional.

The adiabatic local density approximation (ALDA) \cite{ZS80} to the TD XC potential
is the simplest, useful approximation to study the dynamics of atoms and solids.
However, when applied to molecules,
ALDA often yields
unphysical results, for example, atoms with fractional charges at dissociation, underestimated charge-transfer excitation
energies, missing double excitations, among others.  Two decades of research have shown that it is
challenging to enhance the performance
of ALDA while preserving computational simplicity and elegance.

The TD KS equations describe all of the electrons as part
of a single entity, imposing a limit on the number of atoms that can be simulated
in a reasonable amount of time. This limit can be increased by dividing a molecule into fragments
and performing inexpensive calculations on each individual fragment. Several approximated methods to
investigate the electron dynamics of molecules as composed of fragments are available \cite{CW04,CFK07,N10,P13}.
These methods typically assign a set of TD single-particle Schr\"odinger
equations (not necessarily TD KS equations) to every fragment in the molecule. The fragment electrons are subject to the usual interactions as if
they were in the presence of the fragment nuclei only; and, an extra potential, usually fragment-dependent, accounting
for the interaction between the fragments, is added.
Successful applications to the calculation of solvachromatic shifts \cite{N07,PJ12} and excitation
energy of monomers \cite{CFK07} have been reported.

A rigorous extension of TDDFT for molecules made of chemical fragments was presented in Ref. \cite{MJW13}.
In this extension, a molecule is divided into fragments, each a set of atoms. Every fragment is assigned an initial state, and a
Hamiltonian including a global, auxiliary potential, termed partition potential, which enforces
that the total electronic density is the true TD electronic density of the molecule. We proved
that the partition potential is uniquely determined by the TD electronic density of the system, and that
it can thus be expressed (and approximated) as a density-functional. The linear response and an extension to consider electromagnetic
fields are presented in Ref. \cite{MW14}.
The Hamiltonians used in Refs. \cite{MJW13} and \cite{MW14}, and the aforementioned
approximated methods, are particle-conserving, i.e.,
the average number of electrons in a fragment is time-independent; this restriction is eliminated in this work.

In this paper we first introduce a density-inversion method that improves that of Ref. \cite{MJW13}
and allows us to estimate the partition potential. Using a model for the errors, we
discuss, in close connection with standard TDDFT, the
origin of uncertainties of the partition potential in regions far from the molecule.
We extend our fragment-based TDDFT \cite{MJW13} to allow for variable numbers of
electrons in each fragment, while preserving the uniqueness of observables as density-functionals.
The formalism introduced in this paper serves as a theoretical foundation for the development of
methods accounting for electronic excitations and
electron-transfer processes, without sacrificing explicit use of the electronic density and computational efficiency.

\section{Fragment-based TDDFT}
\subsection{Formulation}
An electron in a fragment, labeled $\alpha$, is subject to a 1-body external potential,
denoted as $v_{\alpha}$. For example, $v_{\alpha}({\bf r})=\sum_{i \in I_{\alpha}} -Z_i/|{\bf r}-{\bf R}_i|$;
$I_{\alpha}$ is a set of labels for the nuclei in fragment $\alpha$.
We assign to each fragment in the molecule a Hamiltonian that includes an auxiliary potential,
the partition potential $v_{\mr{p}}$:
\begin{equation}
\hat{H}_{\alpha}[v_{\mr{p}}](t)=\hat{H}_{\alpha}^0+\int \mr{d}^3\mb{r}~ \hat{n}(\mb{r})v_{\mr{p}}(\mb{r}t)~,
\end{equation}
where $\hat{H}_{\alpha}^0=\hat{T}+\hat{W}+\int \mr{d}^3\mb{r}~\hat{n}(\mb{r})v_{\alpha}(\mb{r})$, $\hat{T}$ and $\hat{W}$
are the kinetic, and coulombic repulsion energy operators, respectively, and $\hat{n}({\bf r})$ is the density operator.
This Hamiltonian does not include external driving forces other than those due to
the nuclei of fragment $\alpha$.
TD displacement of the positions of the nuclei can be described by introducing a time-dependent Hamiltonian
where $v_{\alpha}$ is replaced by the corresponding TD fragment-potential, $\sum_{i\in I_{\alpha}} -Z_i/|{\bf r}-{\bf R}_i(t)|$.

The state of a fragment is described by the evolution of the ket $|\psi_{\alpha}[v_{\mr{p}}](t)\rangle$ in Fock space, which
satisfies the TD Schr\"odinger equation (atomic units are used throughout):
\begin{equation}\label{firstevol}
i\partial_t|\psi_{\alpha}[v_{\mr{p}}](t)\rangle=\hat{H}_{\alpha}[v_{\mr{p}}](t)|\psi_{\alpha}[v_{\mr{p}}](t)\rangle~,
\end{equation}
where
\begin{equation}
|\psi_{\alpha}(t)\rangle=\sum_M \nu_{\alpha,M}|\psi_{\alpha,M}(t)\rangle~.
\end{equation}
$\{\psi_{\alpha,M}\}$ are kets corresponding to states with integer number of particles
and $\{\nu_{\alpha,M}\}$ are the weight amplitudes of those states.
Kets with different number of electrons are orthogonal, $\langle \psi_{\alpha,M}|\psi_{\alpha,M'}\rangle=0~~,M\neq M'$.
The total density is given as
\begin{equation}
n(\mb{r}t)=\sum_{\alpha} n_{\alpha}(\mb{r}t)~,
\label{density_constraint}
\end{equation}
which defines $v_{\mr{p}}$ when
$
n_{\alpha}(\mb{r}t)=
\langle\psi_{\alpha}(t)|\hat{n}(\mb{r})|\psi_{\alpha}(t)\rangle~,
$
as proven in Ref. \cite{MJW13}: Given a set $\{\psi_{\alpha,0},\,v_{\alpha}\}$,
two potentials $v_{\mr{p}}$ and $v_{\mr{p}}'$ that differ by more than a TD constant
{\slshape cannot give rise to the same density}. A corollary of this theorem is that
there is a TD density-functional that, when evaluated at a given TD electronic density, gives
the corresponding TD partition potential \footnote{This theorem and corollary were recently used
in Ref. \cite{HLPC14} to propose an inversion method in the context of embedding potential-functional theory.}.

The partition potential represents the TD electronic density
of the supermolecule, and is decomposed as follows \cite{MW14}:
\begin{equation}
v_{\mr{p}}({\bf r} t)=v_{\mr{G}}({\bf r} t)+v_{\mr{d}}({\bf r} t)~.
\end{equation}
$v_{\mr{G}}$ is a ``gluing potential'', accounting for the correlation between the fragments, and
$v_{\mr{d}}$ is the driving potential the molecule is subject to (e.g. laser field). The gluing potential yields
the shape of the potential such that the TD electronic density is recovered.
The gluing potential satisfies \cite{MW14}:
\begin{equation}
\begin{split}\label{gmotion}
\frac{1}{\mr{i}}\nabla\cdot n({\bf r} t)&\nabla v_{\mr{G}}({\bf r} t)=\langle\psi(t)|[\hat{H}^0,\nabla\cdot\hat{\mb{j}}({\bf r})]|
\psi(t)\rangle\\
        &-\sum_{\alpha}\langle\psi_{\alpha}(t)|[\hat{H}^0_{\alpha},\nabla\cdot\hat{\mb{j}}({\bf r})]|\psi_{\alpha}(t)\rangle~.
\end{split}
\end{equation}
The terms on the right-hand-side of Eq.(\ref{gmotion}) are TD density-functionals.
Approximating these terms and solving the
resulting differential equation yields an estimate of the gluing potential. Another
route to approximating $v_{\mr{G}}$ is
by assuming that the system evolves adiabatically from its ground states, driven
by a very slowly-varying field. In such case, the potential $v_{\mr{G}}$ is obtained from
the adiabatic approximation in ground-state Partition DFT \cite{EBCW10,MW12}:
\begin{equation}
v_{\mr{G}}^{\mr{Ad}}[n(t)]=v_{\mr{p}}^{\mr{Ad}}[n(t)]-v^{\mr{HK}}[n(t)]~,
\end{equation}
where $v^{\mr{HK}}[n(t)]$ is the external perturbation the interacting electrons are subject to in their
ground-state in order to yield the density $n(\mb{r}t)$ (the uniqueness of $v^{\mr{HK}}$ follows from the Hohenberg-Kohn theorem).
The partition potential, $v_{\mr{p}}^{\mr{Ad}}[n(t)]$, is the Lagrange multiplier required
to solve the minimization:
\begin{equation}
\min_{\{\psi_{\alpha}\}\rightarrow n(t)} \sum_{\alpha}\lket{\psi_{\alpha}}\hat{H}_{\alpha}^0\rket{\psi_{\alpha}}~.
\end{equation}
under the constraint of Eq.(\ref{density_constraint}).
The Lagrange multiplier for this problem is unique, up to an arbitrary constant \cite{CW06}.

The TD partition KS equations are:
\begin{equation}
\begin{split}
\mr{i}\partial_t &\phi_{i,\alpha}({\bf r},t)=\Big(-\frac{1}{2}\nabla^2+v_{\mr{Hxc}}[n_{\alpha}]({\bf r},t)
\\ &
+v_{\alpha}({\bf r})+v_{\mr{G}}[n]({\bf r},t)+
        v_{\mr{d}}({\bf r},t)\Big)\phi_{i,\alpha}({\bf r},t)~.
\end{split}
\end{equation}
The density can be obtained by means of:
$
n({\bf r},t)=\sum_{i\alpha} f_{i\alpha}|\phi_{i\alpha}({\bf r},t)|^2~,
$
where
$\{f_{i\alpha}\}$ are the (time-independent) occupation numbers, chosen from a proper ensemble \cite{MJW13}.

\section{Classical Interpretation of the Partition Potential}
Consider a system composed of a single massive particle, and a bath
made of particles much smaller than the massive one. All the particles
in the particle+bath system interact via a potential, $U_{\mr{int}}$, of the form
$\sum_{i>j} u_{ij}(\mb{r}_i-\mb{r}_j)$.
The evolution of the subsystem particle, labeled $\mr{S}$,
is dictated by Eq. (\ref{firstevol}). There is no partitioning
of the external potential because the particle and the bath are subject, implicitly, to a
macroscopic, external, confining potential. The Hamiltonian of the subsystem particle is
thus $\hat{H}_{\mr{S}}=\hat{T}+\int \mr{d}^3\mb{r}~\hat{n}(\mb{r})v_{\mr{p}}(\mb{r}t)$, and the
Hamiltonian of the bath is $\hat{H}_{\mr{B}}=\hat{T}+\hat{W}+\int \mr{d}^3\mb{r}~\hat{n}(\mb{r})v_{\mr{p}}(\mb{r}t)$.
The average position of the particle
is $\bar{{\bf r}}_{\mr{S}}(t)=\int \mr{d}^3{\bf r}~{\bf r}~ n_{\mr{S}}(\mb{r},t)$, where
$n_{\mr{S}}(\mb{r},t)=|\psi(\mb{r},t)|^2$.

By the Ehrenfest theorem and
correspondence principle we have
\begin{equation}
m_{\mr{S}}\frac{\mr{d}^2 \bar{{\bf r}}_{\mr{S}}}{\mr{d} t^2}=\mb{F}_{\mr{p},\mr{S}}(t)~,
\end{equation}
where $\mb{F}_{\mr{p},\mr{S}}(t)=-\int \mr{d}^3\mb{r}~n_{\mr{S}}({\bf r},t)\nabla v_{\mr{p}}(\mb{r},t)$,
and $m_{\mr{S}}$ is the mass of the particle.
Comparison with the equation of motion of the real system indicates that, in the classical limit,
$(\nabla v_{\mr{p}})(\bar{r}_{\mr{S}}(t))=(\nabla_{\mb{r}_{\mr{S}}} U_{\mr{int}})
(\bar{r}_{\mr{S}}(t),\bar{\mb{r}}_{\mr{B}}(t))$, where $\bar{\mb{r}}_{\mr{B}}$ is the average coordinate of the bath.
The partition force at the position of the particle is simply the force exerted on the particle by the bath.

As the mass of the subsystem particle is increased, the density tends to a Dirac
distribution. It follows from the above result and Eq. (\ref{gmotion}) that the shape of
the partition potential for any point but that of the particles
is indefinite. However, for given initial momenta and coordinates of the particles and bath
(this condition replaces the requirement that the initial wavefunction be specified), the trajectory
of the momenta of the total system is in a one-to-one correspondence with the trajectory of partition
forces exerted on each particle. Furthermore, if the assumptions of Langevin dynamics are applicable, the partition force on
the massive particle can be interpreted as $\mb{F}_{\mr{p},\mr{S}}(t)=-\gamma \mb{v}_{\mr{S}}(t)+\mb{F}_{\mr{ran}}(t)$,
where $\mb{v}_{\mr{S}}(t)=\mr{d} \bar{{\bf r}}_{\mr{S}}/\mr{d} t$, $\gamma$ is the friction coefficient, and $\mb{F}_{\mr{ran}}$
is the random force.

The gradient of $v_{\mr{p}}$ can {\slshape only} be known at the position of the particles,
nowhere else. In the original proof by Runge and Gross \cite{RG84}, and its extension for
partition potentials \cite{MJW13}, it was found that
potentials whose gradients grow rapidly in regions distant from the molecule
are not covered. This result, as illustrated in Section 4, is related
to the uncertainty in the estimation of the partition potential near the boundaries of the
simulation box.

\section{Numerial TD Potentials}
TDDFT, in which our formulation is built upon, concerns itself with the simplification of the problem:
\begin{equation}
(\mr{i}\partial_t -\hat{H}^{\lambda}[v](t))|\psi(t)\rangle=0,\quad \rket{\psi(0)}=\rket{\psi_0}~,
\end{equation}
where
\begin{equation}
\hat{H}^{\lambda}[v](t)=\hat{T}+\lambda \hat{W}+\int \mr{d}^3 \mb{r}~ \hat{n}(\mb{r})v(\mb{r}t)~.
\end{equation}
Runge and Gross \cite{RG84} showed that if $v$ is Taylor-expandable around $t=0$ and does not have physical
anomalies in the boundaries, then $v$ determines $n$ uniquely up to a TD constant in the potential;
this theorem (recently used to formulate quantum control in TDDFT \cite{CWG12}) can be extended to include non-analytic potentials \cite{RvL11}.
Let us denote the RG map as $\Lambda_{\psi_0}^{\lambda}$; thus, $n(t)=\Lambda_{\psi_0}^{\lambda}[v](t)$.
The operator $\hat{W}$ can stand for different types of electron-electron interactions, such
as screened coulombic repulsion. If $\lambda=0$, then the electrons are non-interacting.

Suppose a well behaved density, $n^{\mr{ref}}$, and an initial state, $\psi_0$, are known. If
$v_1$ and $v_0$ exist, where $v_{\lambda}(t)=(\Lambda_{\psi_0}^{\lambda})^{-1}[n^{\mr{ref}}](t)$,
then the Hartree-exchange-correlation potential for the system, by definition, is $v_{\mr{HXC}}=v_0-v_1$.
Using the exact map $\Lambda_{\psi_0}$ would require solving the TD Schr\"odinger equation, which is what one wants to
avoid in practical calculations.

For the development of functionals, exploration of the map
$\Lambda_{\psi_0}^{\lambda}$ is fruitful; this map could be investigated by solving the problem
$n^{\mr{ref}}(t)-\Lambda_{\psi_0}^{\lambda}[v](t)=0$, which is a root-finding problem.
The first order response of the density
for some perturbation $\delta v$ is $\delta n(\mb{r}t)=\int \mr{d}^3\mb{r}\mr{d} t'~\chi^{-1}(\mb{r}t,\mb{r}'t')
\delta v(\mb{r}'t')$. The response function $\chi^{-1}$ should decay in the asymptotic regions because
large perturbations of $v$ in those regions have a small response in $n$. This relation between
densities and potentials introduces instabilities in root-finding algorithms aimed at reproducing the density
corresponding to a given external potential.
In the ground-state case, the instabilities could be eliminated by enforcing satisfaction of eigenvalue
constraints. For three dimensional applications, in general, capturing the asymptotic regions is difficult when
Gaussian basis sets are used because they do not have the correct asymptotic behavior.

Instead of looking for the exact root, one can solve a minimization problem:
\begin{equation}
\min_{v\in \mc{V}}\int_0^T \lVert n^{\mr{ref}}(s)-\Lambda^{\lambda}_{\psi_0}[v](s)\rVert^2~ \mr{d} s~,
\end{equation}
where $\lVert f \rVert$ is a suitable spatial norm, and $\mc{V}$ is a space of physical potentials.
The quantities $n^{\mr{ref}}(t)$ and
$\langle \psi[v](t)|\hat{n}(\mb{r})|\psi[v](t)\rangle$ need to be approximated. Let us write
$n^{\mr{ref}}(t)-\Lambda_{\psi_0}^{\lambda}[v](t)=\tilde{n}^{\mr{ref}}(t)-\tilde{\Lambda}_{\psi_0}^{\lambda}[v](t)+
\epsilon[n^{\mr{ref}},v]$. $\tilde{n}^{\mr{ref}}(t)$ is the approximation to $n^{\mr{ref}}(t)$ and
$\tilde{\Lambda}_{\psi_0}^{\lambda}[v]$ is the approximation to $\Lambda_{\psi_0}^{\lambda}$.
If $v^*$ is the exact potential representing $n^{\mr{ref}}$, then the problem becomes $\tilde{n}^{\mr{ref}}=
\tilde{\Lambda}_{\psi_0}^{\lambda}[v^*]+\epsilon$. Because we cannot use exact methods to determine
$n^{\mr{ref}}$ and $\Lambda_{\psi_0}^{\lambda}$, we assume that $\epsilon$ is a random function of the space-time
coordinates. Moreover, one would
expect that $\tilde{n}^{\mr{ref}}$ and $\Lambda_{\psi_0}^{\lambda}$ have smooth timespace gradients, and that
$\epsilon$ displays autocorrelation because the spacing between points is arbitrarily small.

\subsection{Estimation of the Partition Potential}
Let $\mc{V}_{\mr{p}}$ be a space of TD partition potentials, and $\mc{D}$ a space of TD densities and define
the map:
\begin{equation}
\Lambda_{S_0}:\mc{V}_{\mr{p}}\rightarrow\mc{D}~,
\end{equation}
where $S_0=\{\psi_{\alpha,0},v_{\alpha}\}$. For a given TD partition potential, the density
is obtained by evaluation of the above map at the given partition potential; in other words,
$n(t)=\Lambda_{S_0}[v_{\mr{p}}](t)$. This map depends on the history of the partition potential,
i.e., it has memory dependence \cite{MJW13}.

Let $v_{\mr{p}}^*$ be the true partition potential. We assume that, due to numerical
errors, the estimation to the reference density $\tilde{n}^{\mr{ref}}$ is of the
form $\tilde{n}^{\mr{ref}}=\tilde{\Lambda}_{S_0}+\epsilon$, where $\epsilon$
is a {\slshape random function}.
To estimate the partition potential corresponding to $\tilde{n}^{\mr{ref}}$ we minimize:
\begin{equation}
\lVert e[v_{\mr{p}}]\rVert^2_{\mu}=\lVert \tilde{n}^{\mr{ref}}-\tilde{\Lambda}_{S_0}[v_{\mr{p}}]\rVert^2_{\mu}~,
\end{equation}
where $\mr{d}\mu(\mb{r},t)$ is the measure.

Since $\epsilon$ is a function, its probability density function (PDF) is a {\slshape functional}.
The PDF depends on parameters, denoted collectively as $\Theta$, and the PDF itself as $D([\epsilon]|\Theta)$.
The probability that $\epsilon$ is observed in a set $\mc{U}$ is given by the path integral:
\begin{equation}
P(\epsilon\in \mc{U}|\Theta)=\int_{\mc{U}}\mr{d} m_{\mr{L}}[\epsilon]~ D([\epsilon]|\Theta)~,
\end{equation}
where the measure over the space of errors is $m_{\mr{L}}$. The traditional methods
of non-linear regression can be applied to estimate the best parameters of the distribution, $\Theta^*$,
for a given set of observations. Then a Taylor expansion in terms of the parameters can be used
to generate their PDF, which can then be used to estimate the error in the parameters.
In this case, the parameters are the variance and the partition potential.

In the next section, we will expand the partition potential in a spline basis set.
The parameters are the values of the partition potential at the knots, and they
are correlated: A perturbation of the partition potential at
one knot affects the response of the density in other knots. Hence, we must employ
a model of correlated errors. Finding the correct model is a demanding task
beyond the scope of this work. For this reason, we choose a biased model based
on the following observations: i) A measure of the error of the form
$
\int \mr{d}^3{\bf r}\mr{d} t~(\tilde{n}^{\mr{ref}}({\bf r},t)-\tilde{\Lambda}_{S_0}[v_{\mr{p}}])^2
$
suffers of autocorrelation. ii) Far from the molecule, the partition potential
has little influence on the density. iii) Estimating the density is not sufficient;
its spatio-temporal gradient is an important quantity. An error measure accounting
for these observations is:
\begin{equation}\label{err_fun}
\lVert e[v_{\mr{p}}]\rVert^2_{\mu}=
\int \mr{d}\mu({\bf r},t)\{ |\nabla e({\bf r},t)|^2+(\partial_t e({\bf r},t))^2\}~.
\end{equation}
Based on ii), we choose a measure of the form $\mr{d}\mu({\bf r},t)=\mr{d}^3{\bf r}\mr{d} t~\sum_i \tilde{n}^{\mr{ref}}({\bf r}_i,t)\delta({\bf r}-{\bf r}_i)$.
Where $\{{\bf r}_i\}$ are points selected in such a way that $|\nabla e|^2+(\partial_t e)^2$ resembles
a $\chi^2$-distribution.
To apply this measure of error in the next section, we transform
the above measure into a vector norm.
Then, the resultant distribution is expanded
in terms of the partition potential evaluated at the knots, and asymptotic analysis
is applied \cite{SW89}, leading to the random variables required to reproduce the density
within a small error tolerance.

\subsection{1d Electron in a Double-well Potential}\label{texample}
Let us revisit the example considered in Ref. \cite{MJW13}:
A one dimensional ``electron" in a double well potential.
The TD partition equations are:
\begin{equation}
\mr{i}\partial_t \phi_{\alpha}(x,t)=\Big(-\frac{1}{2}\partial^2_x+v_{\alpha}(x)+v_{\mr{p}}(x,t)\Big)\phi_{\alpha}(x,t)~,
\end{equation}
where $\alpha=\mr{L},~\mr{R}$, standing for left and right wells.
The potentials are $v_{\alpha}(x)=v_0/\sqrt{(x-x_{\alpha})^2+a}$; the
parameters are: $v_0=-1$, $x_{\mr{R}}-x_{\mr{L}}=4$, and $a=1$.
The density is obtained by averaging over the orbital densities of both wells:
\begin{equation}
n(x,t)=\frac{1}{2}|\phi_{\mr{L}}(x,t)|^2+\frac{1}{2}|\phi_{\mr{R}}(x,t)|^2~.
\end{equation}
Suppose that the supermolecule evolves from the ground state driven by a monochromatic laser. The
evolution of the system is thus dictated by the solution of:
\begin{equation}
\mr{i}\partial_t\psi(x,t)=\Big(-\frac{1}{2}\partial^2_x+v(x)+v_{\mr{d}}(x,t)\Big)\psi(x,t)~,
\end{equation}
where $v_{\mr{d}}(xt)=Ex\sin \omega t$, and the external potential is $v=v_{\mr{L}}+v_{\mr{R}}$.
The density obtained from the above evolution equation is $n^{\mr{ref}}(xt)=|\psi(xt)|^2$,
which is the target density we wish to represent.
\begin{figure}[tbh!]
\centering
\subfigure[$v$ and $v_{\mr{p}}^0$]{
\includegraphics[scale=0.20]{ini_vp.eps}}
\subfigure[Initial Densities]{
\includegraphics[scale=0.20]{fragden_t_0.eps}}
\\
\subfigure[$v_{\mr{p}}(t=1.0)$]{
\includegraphics[scale=0.20]{vp_t_10.eps}}
\subfigure[$n(t=1.0)$]{
\includegraphics[scale=0.20]{fragden_t_10.eps}}\\
\subfigure[$v_{\mr{p}}(t=8.0)$]{
\includegraphics[scale=0.20]{vp_t_80.eps}}
\subfigure[$n(t=8.0)$]{
\includegraphics[scale=0.20]{fragden_t_80.eps}}
\caption{Snapshots of the partition potential. In a), solid line: Initial partition potential,
dashed line: Left fragment external potential, dashed-dotted line: Right fragment external potential.
In b), d), and f), solid lines: Left electronic density, dashed lines: Right electronic
density. In d) and f) the dashed-dotted line is the total density.
\end{figure}

The laser parameters are $\omega=0.3$, $E_0=0.05$.
We propagate the states of the system using the Crank-Nicholson method;
time step is 0.1, box length is 20, spatial step is 0.2, and total
propagation time is 10 units.
The partition potential is represented in a smoothed, cubic, spline basis set with 22 knots equally spaced in the box. The
initial partition potential is estimated by minimizing the error using sequential quadratic programming (other
useful methods employ the density error directly \cite{RABB15}).
To obtain the initial densities, the error functional of Eq. (\ref{err_fun}) (with the time-dependency dropped)
is minimized.
First the problem $(-1/2\partial_x^2 +v_{\alpha}+v_{\mr{p}}^0)\phi_{n,\alpha}=\epsilon_{n,\alpha}\phi_{n,\alpha}$
is solved (with the finite-difference method) for both wells with some estimation of $v_{\mr{p}}^0$. Then,
the density is compared with that of the system of reference in order
to propose the next estimation in the iterative procedure of sequential quadratic programming;
the constraints are conservation of charge and chemical potential (HOMO) equalization \cite{MW14}.
Figure 1.a. shows the initial partition potential and external potentials of each well.
The initial fragment densities that add up to the ground-state density of the supermolecule
are displayed in Figure 1.b.

The estimation of the evolution of $v_{\mr{p}}$ is determined by using the
step-by-step scheme proposed in Ref. \cite{MJW13}, and the norm of Eq. (\ref{err_fun}):
For example, for $\Delta t=0.1$, the error norm can be approximated as
\begin{equation}
\Delta t \sum_i \{\tilde{n}^{\mr{ref}}(x_i,\Delta t)[(\partial_x e(x_i,\Delta t))^2+(\partial_t e(x_i,\Delta t))^2]\}~.
\end{equation}
The numerical value of the above quantity depends on the value of $v_p$ at the spline knots
and at $t=\Delta t/2$. Hence, this quantity is varied until the above function is minimized and
the total density of the fragments match that of the true system. The procedure
is repeated similarly for the rest of the propagation.
Figure 1.c. shows the partition potential at $t=1.0$;
it is localized in the intermediate region between the fragments.
The fragment electronic densities (Figure 1.d) are also well localized at $t=1.0$.
Because in absence of the partition potential the fragment-densities would just be localized
around their wells, the partition potential must be such that it induces the transfer
of charge from the right fragment into the left fragment (Figure 1.e). However, as we note in Figure 1.f,
the charge transfer in this case is represented by the spreading of the right electronic density
into the left one, not by a change in the fragment populations. Two observations: i) If one were to assign a grid
that is fine enough around the center of the wells and then coarser as one moves away from the wells,
then to describe the density spreading, the grid would need to be time-dependent.
ii) The partition potential must induce the charge transfer and act like a ``spoon''.

\begin{figure}[tbh!]
\centering
\includegraphics[scale=0.35]{errors.eps}
\caption{Error-bar plot of the partition potential at $t=6.2$. Solid line:
Interpolated optimal potential, circles: Interpolation knots.}
\end{figure}

The result of the error estimation in the partition potential at $t=6.2$ is shown in
Figure 2. As expected, in the boundary regions, the error is quite significant, and the error bars extend well beyond the plotting range.
This implies that the shape of the potential
in these regions is not reliable. All space-time points obeying causality are coupled.
For example, variation of a single knot in the boundary affects its neighbors,
introducing large gradient fluctuations.
Thus, the error can spread to regions where the density is non-negligible.
This can cause instabilities in the minimization procedure if a norm such as
$\int \mr{d}^3\mb{r}\mr{d} t~ (\tilde{n}^{\mr{ref}}(\mb{r},t)-n(\mb{r},t))^2$ is employed. For this reason
we recommend the use of derivatives of the density to measure the error.

\section{Variable Occupation Numbers}
Let us assign variable electron-occupation numbers to the fragments.
First, divide the total propagation time into blocks:$[0,\tau)\cup [\tau,2\tau)\cup\ldots\cup [(m-1)\tau,m\tau)$,
where $m\tau$ is the final time of the propagation, and let
\begin{equation}
X_{\alpha}=\{|\xi_{\alpha}^0\rangle,|\xi_{\alpha}^1\rangle,\ldots,|\xi_{\alpha}^{m-1}\rangle\}~,
\end{equation}
be a set of instantaneous kets, in Fock space, for fragment $\alpha$.
At a single time $t=k\tau$, the following minimization is performed to obtain
the set of kets describing the density of the fragmented molecule:
\begin{equation}\label{pdftmin}
\begin{split}
\{|\xi_{\alpha}^k\rangle\}_{\alpha=1}^{N_{\mr{frag}}}=
\arg &\min\Big\{\sum_{\alpha}\langle\psi_{\alpha}^k|
\hat{H}_{\alpha}^0|\psi_{\alpha}^k\rangle ~~\mr{s.t.}\\& \sum_{\alpha}\langle \psi_{\alpha}^k|\hat{n}(\mb{r})|
\psi_{\alpha}^k\rangle= n(\mb{r},k\tau),~\forall \mb{r}\Big\}~,
\end{split}
\end{equation}
The occupation numbers of fragment $\alpha$ are formally expressed as
$ |\nu_{\alpha,M}(k\tau)|^2=|\langle \xi^k_{\alpha,M}|\xi^k_{\alpha}\rangle|^2~.$
Here, $|\xi^k_{\alpha,M}\rangle$ is an optimal ket (obtained from solving the right hand side of Eq. (\ref{pdftmin}))
for fragment $\alpha$ with $M$ electrons.
These numbers are density-functionals.

The evolution operator of fragment $\alpha$ is:
$
\hat{U}_{\alpha}[v_{\mr{p}}](t_1,t_2)=\mc{T}\exp(-\mr{i}\int_{t_0}^{t_1}\mr{d} s~\hat{H}_{\alpha}[v_{\mr{p}}](s))~.
$
Introduce the displaced set of kets:
\begin{equation}
\begin{split}
\tilde{X}_{\alpha}&=\{\hat{U}_{\alpha}(\tau,0)|\xi_{\alpha}^0\rangle,\\&\hat{U}_{\alpha}(2\tau,\tau)|\xi_{\alpha}^1\rangle,\ldots
    \hat{U}_{\alpha}(m\tau,(m-1)\tau)|\xi_{\alpha}^{m-1}\rangle\}~.
\end{split}
\end{equation}
Now let us define the following dyadic product: $(X_{\alpha}\tilde{X}^{\dagger}_{\alpha})(k)=|\xi_{\alpha}^k\rangle\langle \tilde{\xi}_{\alpha}^{k-1}|$.
The symbol $X_{\alpha}\tilde{X}_{\alpha}^{\dagger}$ is the set of dyadic products where the $k$-th component
is the dyadic product between the ket at the beginning of the $k$-th block and the
displaced ket from the $k-1$-th block.
Now, let $\hat{B}_{\alpha}$ be the TD operator:
\begin{equation}
\begin{split}
\hat{B}_{\alpha}(t)&=(W_{\tau}*\ln X_{\alpha}\tilde{X}^{\dagger}_{\alpha})(t)\\
&=\sum_{k=1}^{m}\delta(t-k\tau)\ln |\xi_{\alpha}^k\rangle\langle \tilde{\xi}_{\alpha}^{k-1}|~.
\end{split}
\end{equation}
where $W_{\tau}$ is the Dirac-Comb kernel. Addition of the operator $\hat{B}_{\alpha}$ to the Hamiltonian
$\hat{H}_{\alpha}(t)$ yields the non-Hermitian operator:
\begin{equation}
\hat{H}_{\mr{c},\alpha}[v_{\mr{p}}](t)=\hat{H}_{\alpha}[v_{\mr{p}}](t)+\mr{i}\hat{B}_{\alpha}[v_{\mr{p}}](t)~.
\end{equation}
The evolution of the system is now determined by $|\psi_{\alpha}[v_{\mr{p}}]\rangle$, which obeys
\begin{equation}
i\partial_t|\psi_{\alpha}[v_{\mr{p}}](t)\rangle=\hat{H}_{\mr{c},\alpha}[v_{\mr{p}}](t)|\psi_{\alpha}[v_{\mr{p}}](t)\rangle~.
\end{equation}
The total density is $n(\mb{r}t)=\sum_{\alpha}\langle \psi_{\alpha}(t)|\hat{n}(\mb{r})|\psi_{\alpha}(t)\rangle$ and the
number of particles in fragment $\alpha$ is $N_{\alpha}(t)=\langle \psi_{\alpha}(t)|\hat{N}|\psi_{\alpha}(t)\rangle$.
In general, any observable, $\hat{O}(t)$, is expressed as a functional of the partition potential,
$\langle \psi_{\alpha}[v_{\mr{p}}](t)|\hat{O}(t)|\psi_{\alpha}[v_{\mr{p}}]\rangle$.

Given the partition potential and occupation numbers as density-functionals, the
scheme to determine the evolution of the molecule is the following:
First, propagate the kets $\{\rket{\psi_{\alpha}}\}$ in the interval $[0,\tau)$ with fixed populations on each
fragment. Then, at $t=\tau$, obtain new occupation numbers according to Eq. (\ref{pdftmin}) as well
as new states to propagate; continue
the propagation in the block $[\tau,2\tau)$. The procedure continues similarly
for the rest of the propagation. The density of the system is then obtained
as $n({\bf r},t)=\sum_{\alpha}\lket{\psi_{\alpha}(t)}\hat{n}({\bf r})\rket{\psi_{\alpha}(t)}$.
The theorem discussed in section 2 also applies in this case. Therefore, the partition potential
for this scheme is uniquely determined by the TD electronic density, up to an arbitrary constant.
\begin{figure}[htb!]
\centering
\subfigure[$N_L$]{
\includegraphics[scale=0.25]{surf_occ.eps}}\\
\subfigure[$v_{\mr{p}}(t=40)$]{
\includegraphics[scale=0.20]{vp_t_40.eps}}
\subfigure[$n(t=40)$]{
\includegraphics[scale=0.20]{fragden_t_40.eps}}
\subfigure[$v_{\mr{p}}(t=70)$]{
\includegraphics[scale=0.20]{vp_t_70.eps}}
\subfigure[$n(t=70)$]{
\includegraphics[scale=0.20]{fragden_t_70.eps}}
\caption{Evolution of the fragments with TD electron populations. In a) the solid line
is the result from the inversion, the dashed line is the result from the two-state approximation, and
the dashed-dotted line is obtained by omitting the gluing potential: $v_{\mr{G}}=0$. In c) and e),
solid line: $n_{\mr{L}}$, dashed line: $n_{\mr{R}}$, dashed-dotted: $n$.}
\end{figure}

The partition potential is discontinuous at the relaxation nodes (points where
$t$ is an integer multiple of $\tau$). Discontinuities in time can be eliminated by using an integral transformation
that smooths the observable at the relaxation nodes. In practice, however,
it is convenient to propagate the occupation numbers and gluing potential assuming that they are
continuously differentiable functions of time. It can be shown, assuming that
the dynamics of the occupation numbers is much slower than that of the partition potential,
that the 1-1 map between the former and the density still holds. This follows
from the scheme we have shown here because the electronic populations are fixed in the first
block, allowing us to apply the Runge-Gross theorem in such block.

A density-functional approximation to the occupation numbers is needed
to apply the theory.
The dynamics of the occupation numbers can be investigated using master
equations, where the rate coefficients are determined by Fermi's golden rule,
or transition elements that couple the fragments. Here, we illustrate
a simple approach:
A trial wave function to investigate the evolution of the
occupation numbers is $\rket{\eta(t)}=a_{\mr{L}}(t)\rket{\varphi_{\mr{L}}}
+a_{\mr{R}}(t)\rket{\varphi_{\mr{R}}}$, where $\rket{\varphi_{\alpha}}$
is the ground-state of the electron described only by $\hat{H}_{\alpha}^0$
(This hamiltonian in coordinate representation is $-1/2\partial^2_x+v_{\alpha}(x)$).
The dynamics of electron transfer is governed
by a two-component wave-function $\mb{a}=(a_{\mr{L}},a_{\mr{R}})^{\mr{T}}$.
We assume that the Hamiltonian coupling that relates the two fragments is of the form:
\begin{equation}
\hat{\mc{H}}(t)=\hat{H}_{\mr{f}}^0+\int \mr{d} x~(v_{\mr{G}}(x,0)+v_{\mr{d}}(x,t))\hat{n}(x)
\end{equation}
where $\hat{H}_{\mr{f}}^0=\hat{H}_{\mr{L}}^0\oplus\hat{H}_{\mr{R}}^0$, is the uncoupled
Hamiltonian; $\hat{H}_{\alpha}^0\rket{\varphi_{\beta}}=0$ if $\alpha\neq\beta$.
For the sake of the illustration, the gluing field is frozen,
serving as a ``bridge'' for the charge to be transferred from one well into the other.

From the evolution equation: $\mr{i} \partial_t\rket{\eta(t)}=\hat{\mc{H}}(t)\rket{\eta(t)}$ we
infer that the state vector, $\mb{a}$, satisfies:
\begin{equation}\mr{i}\partial_t \mb{a}(t)=\mb{S}^{-1}(\bm{\epsilon}_0+\bm{\Delta}(t))\mb{a}(t)
\end{equation}
where $S_{\alpha\beta}=\int \mr{d} x~\varphi_{\alpha}^*(x)\varphi_{\beta}(x)$, $\bm{\epsilon}_0=\mr{diag}(\epsilon_0,\epsilon_0)$,
and
\begin{equation}
\Delta_{\alpha\beta}(t)=\int \mr{d} x~ \varphi_{\alpha}^*(x)(v_{\mr{G}}(x,0)+v_{\mr{d}}(x,t))\varphi_{\beta}(x)~.\end{equation}
The occupation numbers are obtained from the ``density'' of $\mb{a}$:
$N_{\alpha}(t)=|\mb{a}_{\alpha}|^2(t)+\mr{Re}(a^*_{\mr{L}}(t)a_{\mr{R}}(t)S_{\mr{LR}})$.
The last term arises from the overlap of the functions $\varphi_{\mr{L}}$ and $\varphi_{\mr{R}}$,
guaranteeing that $N_{\mr{L}}+N_{\mr{R}}=1$.

The example of the previous section, summarized in Fig. 1., illustrates how the partition potential
can be estimated even under
a strong laser field.
We now return to the example of section \ref{texample} and allow for variable occupations. The parameters for the propagation now
are $\tau=2$, $\Delta t=1$, $\omega=0.3$, $E_0=0.02$.
Comparison of the exact time-dependency of the average number of electrons of the left fragment with
the approximation described above is shown in Figure 3.a. The initial gluing potential
used here is that shown in Fig. 1.a.
The two-state approximation works well at short times, and displays deviations after $t=20$.
Capturing the results of the two-state approximation would be quite challenging by fixing
the occupation numbers and finding the corresponding partition potential.
Improvements over the two-state approximation can proceed by either refining the gluing potential (going beyond
the frozen approximation) or increasing the number of states
considered to couple the fragments. The first alternative has the advantage that the
equations can be solved very fast. Nonetheless, for functional development, the gluing potential is also a determinant
factor for the evolution of the shape of the electronic fragment-density ($\int n_{\alpha}/N_{\alpha}$).

Figure 3.b shows a snapshot of the ``exact'' partition potential at $t=40$. In contrast with
the results of section 3.2, the partition potential is now well localized (Figures 3.d and 3.f).
This suggests that the standard methods of ground-state PDFT can be used to estimate
the partition potential through the use of the adiabatic approximation.
The fragment densities also remain localized (Figure 3.c, 3.e, and 3.g). Qualitatively,
the partition potential accounts for the shape of the electronic densities
of the fragments, while the occupation numbers are responsible for their height.

\section{Concluding Remarks}
We formulated a TDDFT for treating molecules as composed of smaller composite
units. The successful application of this formulation requires
approximations to the partition potential and the occupation numbers. This can
be accomplished by a proper approximation to the Hamiltonians $\{\hat{H}_{{\mr{c}},{\alpha}}(t)\}$,
or the auxiliary evolution equations of the electron populations in the fragments.
The error analysis was also presented. It leads to a simple form of estimating the
errors in the potentials. The partition potential, obtained
by numerical inversion, can be uncertain in spatio-temporal regions where the density is small.
However, as time increases, the error
can propagate from the boundary areas into regions were the density is high.

\section{Acknowledgments}
We acknowledge support from the National Science Foundation CAREER program under grant No.CHE-1149968.
AW also acknowledges support from the Alfred P. Sloan Foundation and the Camille Dreyfus Teacher-Scholar
Awards Programs.


\section{Introduction}
\subsection{Background}
Hard X-ray Modulation Telescope (HXMT) is a planned Chinese space X-ray
telescope.
The telescope will perform an all-sky survey in hard X-ray band
($1$ -- $250\;\mathrm{keV}$), a series of deep imaging observations at small sky
regions and pointed observations.

We expect a large number of X-ray sources, e.g., AGNs, to be detected in its
all-sky survey.
We also expect through a series of deep imaging observations at the Galactic
plane X-ray transients can be detected\citep{li2007,lu2012}.
Therefore the point source detection performance is one of our concerns on
HXMT data analysis.

Methods and corresponding sensitivities of pointed observation have been
discussed by \citet{jin2010}.
In imaging observations such as all-sky survey and deep imaging at small
sky regions, a variety of data analysis aspects and methods are involved.
\textbf{First, images are computed instead of recorded directly by the optical
instrument.}
Mapping as well as image reconstruction methods are useful.
The direct scientific data products from imaging observations of HXMT are still scientific events, more specifically, X-ray photon arrival events. The attitude control system (ACS) of the spacecraft
reports the attitudes periodically. We take these reported attitudes as nodes to perform certain interpolations to determine the instantaneous attitude for each scientific event. In this way, a set of
parameters are assigned to each event, including the coordinates on the celestial sphere where the telescope is pointing at. Hence we call this process \textit{events mapping}, where scientific events
are mapped from time domain to the celestial sphere. The product of this process is refered to as the observed image, which implies the dimensionality of the data.
\textbf{Second, the exposure to a specific source is limited more strictly}
thus the signal-to-noise ratio (SNR) is tightly restricted.
Hence the importance of regularization methods becomes prominent.
\textbf{Finally, \textit{a picture is worth a thousand words}}.
Various information can be extracted from an image by numerical methods.

In this paper we investigate the point source detection performance of
the imaging and detecting system synthesised from the telescope as well
as diverse combinations of data analysis methods, especially the regularization
methods.

\subsection{Modulation in HXMT imaging observation}
The PSF of HXMT HE, which is a composite telescope consisting of 18 collimated
detectors, describes the response of the telescope to a point source when the
telescope is pointing at the source as well as different locations around it.
In other word, the PSF is a density function of the distribution of responses
occur in different observation states, which is denoted by the instantaneous
attitude of the telescope as a spacecraft.
So the PSF takes the attitude of the telescope as its input.

We use the proper Euler's angles to describe the attitude of the telescope,
i.e., $\phi$ and $\theta$ the longitude and latitude of the pointing, as well
as $\psi$ denoting the rotation angle of the telescope around its own pointing
axis, namely, the position angle.
The modulation equation that corresponds to the imaging observation over the
sphere surface is
\begin{equation}
d(\phi,\theta,\psi) = \iint\limits_{\Omega} p(\phi,\theta,\psi,\phi',\theta')
  f(\phi',\theta')\cos\theta'\mathrm{d}\phi'\mathrm{d}\theta' \text{,}
\end{equation}
where $f(\phi',\theta')$ is the object function (i.e., the image) defined in
a compact subset of the sphere surface $\Omega$,
$p(\phi,\theta,\psi,\phi',\theta')$ is the modulation kernel function that
relates the value of object function $f(\phi',\theta')$ defined on a
neighbourhood of the point $(\phi',\theta')$ of the subset $\Omega$
to the instantaneous response of the telescope, while its status is
$(\phi,\theta,\psi)$.

The modulation kernel function is determined by the point spread function (PSF)
of the telescope.
Suppose a unit point source is located at the zenith of the
celestial sphere, i.e., the point $(0,0,1)$ in the corresponding Cartesian
coordinate system.
Fix the position angle $\psi$, while slew the telescope across the polar cap,
and assign responses of the telescope to the unit source to pixels on the
celestial sphere.
Then we obtain a map $P(\phi,\theta)$ represents the PSF on the celestial
sphere with fixed rotation angle $\psi$.

The map is then projected to a tangent plane of the celestial sphere $z=1$ by
gnomonic mapping, i.e.,
\begin{equation}
\left\{
    \begin{aligned}
    u &= \cot\theta\cos\phi \\
    v &= \cot\theta\sin\phi
    \end{aligned}\right.
\label{eq-gnomonic}\text{,}
\end{equation}
where $u$ and $v$ are local Cartesian coordinates on the tangent plane.
Now we have
\begin{equation}
P_{\tan}(u,v) = P_{\tan}(\cot\theta\cos\phi,\cot\theta\sin\phi) =
  P(\phi,\theta)\text{,}
\end{equation}
i.e., the PSF defined on the tangent plane.

To provide for data analysis we have two discrete forms of $P_{\tan}$,
\begin{equation}
P_{i,j} = \frac{\iint_{\alpha_{i,j}} P_{\tan}(u,v)\mathrm{d}u\mathrm{d}v}
  {\iint_{\alpha_{i,j}} \mathrm{d}u\mathrm{d}v}\text{,}
\end{equation}
where $\alpha_{i,j}$ is the neighbourhood of the pixel $(i,j)$ of the
2-D pixel grid, and the normalized form
\begin{equation}
H_{i,j} = \frac{\iint_{\alpha_{i,j}} P_{\tan}(u,v)\mathrm{d}u\mathrm{d}v}
    {\iint P_{\tan}(u,v)\mathrm{d}u\mathrm{d}v}\text{.}
\label{eq-pixelization}
\end{equation}

Given the discrete image $\vect{F}$, the detection area $A$ and the
duration of exposure on each pixel $\tau$, the observed data is
\begin{equation}
\vect{D} = \tau \cdot \left[\frac{A}{\max\limits_{i,j}H_{i,j}}
    \left(\vect{F} \ast \vect{H}\right) + r_b\right] \text{,}
\label{eq-modulation}
\end{equation}
where $r_b$ is the constant background count rate,
$\ast$ denotes the convolution,
and $\vect{H}$ is the normalized PSF on the tangent plane.
Eq. \ref{eq-modulation} approximates modulations in HXMT imaging observations.

Distortion occurs when projecting a set of points from (or to) a sphere
surface to (or from) a plane. For example, distance between any two
points, area of any continuous subset, angle between any two crossing lines
(or tangents of curves) are altered non-uniformly.
On the other hand, the rotation angle $\psi$ is not always fixed during HXMT
imaging observations.
Both distortions in image reconstruction from HXMT observed data and rotations
during imaging observations have been discussed by \citet{huo2013raa}.
Here for the sake of simplicity, we ignore them in this article.

\section{Numerical methods}
\label{sect-methods}
\subsection{Single point source detection performance estimation}
\label{sect-estimation}
We estimate the single point source performance in terms of sensitivity and
position accuracy through the following procedures.
\begin{enumerate}
\item Determine flux threshold for point source detection.
    \label{item-threshold}
    \begin{enumerate}
    \item Simulate a frame of observed data contains only background counts.
        \label{item-simulate-data}

    \item Run denoise program on the simulated data to try to increase the
    signal-to-noise ratio.
        \label{item-denoise}

    \item Demodulate the denoised data.
        \label{item-demodulate}

    \item Run SExtractor, a source detection program, by
    \citet{bertin1996sextractor} on the demodulated image to detect point
    sources and extract their intensities, coordinates and other parameters.
    At this point a catalog of point sources is compiled from the simulated
    data.
    Point sources detected here, i.e., from images demodulated from background
    data are false detections.
        \label{item-extract-source}

    \item Repeat previous steps (from \ref{item-simulate-data} to
    \ref{item-extract-source}) and a series of catalogs are compiled.
    Draw a histogram of flux of false detections that could possibly been
    detected from background counts given a specific condition of both
    observation and detection.

    \item Choose a cut from the histogram as the flux threshold so that
    a certain percent of the false detections are rejected and the rejection
    percentage is precise enough.
    The rejection percentage, e.g., $95\%$ or $99.7\%$ etc., reflects the
    significance of detections above the corresponding threshold.
    \end{enumerate}

\item Estimate detection efficiency and position accuracy of a point source
of specific flux intensity.
    \label{item-efficiency}
    \begin{enumerate}
    \item Simulate observed data contains a single point source of given
    flux intensity $f_m$ located at $(x_m, y_m)$ in the model image.
        \label{item-simulate-data-pts}

    \item Perform steps from \ref{item-simulate-data} to
    \ref{item-extract-source}.
    A catalog is compiled.
        \label{item-extract-source-pts}

    \item Examine each detection is the catalog.
    Provided the $i-$th source in the catalog is detected at $(x_i, y_i)$
    in the demodulated image, the distance between the extracted source and
    the true point source
    \begin{equation}
      \delta_i = \sqrt{\left(x_m - x_i\right)^2 + \left(y_m - y_i\right)^2}
    \end{equation}
    as well as the flux of the extracted source $f_i$ are investigated to
    determine whether the $i-$th source is a true source or not.
    We define the score of the current catalog in detecting the single point
    source as
    \begin{equation}
        d_k =
        \begin{cases}
            1 & \exists \; i : \left(\delta_i \leq \Delta\right)\land
                \left(f_i \geq F_{thres}\right) \\
            0 & \text{otherwise}
        \end{cases}\text{,}
    \end{equation}
    where $k$ is the index of the current catalog, $\Delta$ and $F_{thres}$
    are position accuracy and flux thresholds,
    therefore the score $d_k$ reveals whether the $k-$th catalog contains
    the true source or not, in other word, through the previous steps
    (simulated observation, denoising, demodulation, source extraction and
    thresholding) whether we have detected the true source effectively.
    If we have, the outcome of these steps is counted as an effective
    detection of the true source, otherwise it is ineffective.
        \label{item-examine-source-pts}

    \item Repeat the previous steps (from \ref{item-simulate-data-pts} to
    \ref{item-examine-source-pts}) $N$ times and calculate the percentage of
    effective detections amongst all detections, namely,
    \begin{equation}
        \eta = \frac{1}{N}\sum_k^N d_k\text{,}
    \end{equation}
    which is defined as the detection efficiency.

    Let $(x_k, y_k)$ be the position of the brightest source in the $k-$th
    catalog, the position accuracy is calculated as
    \begin{equation}
    \rho = \frac{1}{\eta N}\sum_k^N d_k \sqrt{\left( x_k - x_m \right)^2 +
      \left( y_k - y_m \right)^2} \text{.}
    \end{equation}
    \end{enumerate}

\item Find an flux intensity $F_{0.5}$ so that the corresponding detection
effeciency $\eta=50\%$. The intensity $F_{0.5}$ marks the point source
sensitivity of the detecting system synthesised from both the telescope
in specific status and the data analysis program chain.
\end{enumerate}

\subsection{Imaging and mapping}

\subsubsection{Demodulation}
Direct demodulation (DD) method \citep{li1994dd} is used to estimate the
true image from observed data.
Residual map calculated with CLEAN algorithm \citep{hoegbom1974} is used
as lower limit constraint in DD method.
The skewness of the residual map is calculated in each CLEAN iteration and
its minimum absolute value serves as the main stopping criterion of iterations.

\subsubsection{Cross-correlation}
In contrary to detecting and extracting point sources from demodulated images,
it is also feasible to do so from cross-correlated maps as long as the point
sources are isolated with each other, compared to the FWHM of the PSF,
since the position of a peak of the expected value of such a map coincides
with the position of a source regardless of the PSF.
Cross-correlating the observed data and the PSF yields the correlated map
\begin{equation}
\vect{C} = \vect{D} \star \vect{H}\text{,}
\end{equation}
where $star$ denotes cross correlation, $\vect{D}$ and $\vect{H}$ denote the
observed data and the PSF respectively.
The peak of $\vect{C}$ that coincides with a point source of flux intensity
$F$ is
\begin{equation}
C = \frac{\tau \cdot A \cdot F}{\max\limits_{i,j}H_{i,j}}\sum_{i,j}H_{i,j}^2
  \text{,}
\end{equation}
according to Eq. \ref{eq-modulation}.
On the other hand, the background variance of the correlated map is
\begin{equation}
\sigma^2\left(\vect{C}\right) =
  \sum_{i,j}\sigma^2\left(\tau \cdot r_b \cdot H_{i,j}\right) =
  \tau \cdot r_b \cdot \sum_{i,j}H_{i,j}^2\text{,}
\end{equation}
since $\tau \cdot r_b$ follows Poisson distribution.
Hence the significance of the peak in term of numbers of $\sigma$ is
\begin{equation}
\mathrm{SI} =
  \frac{F \cdot A \cdot \sqrt{\tau \cdot \sum_{i,j}H_{i,j}^2}}
  {\max\limits_{i,j} H_{i,j}\sqrt{r_b}} =
  \frac{F \cdot A \cdot \sqrt{T}}{\sqrt{r_b}}
    \sqrt{\langle P_{\tan}^2 \rangle}\text{,}
\label{eq-significance}
\end{equation}
where $T$ is total duration of exposure on the 2-D pixel grid,
$\sqrt{\langle P_{\tan}^2 \rangle}$ is the square root of the arithme mean
of $P_{\tan}$ over the 2-D pixel grid, which is determined by the PSF as
well as the range of the pixel grid only, provided the pixel grid is fine
enough (see Eq. \ref{eq-pixelization}).

Cross-correlation significance of isolated point source defined Eq.
\ref{eq-significance} can be evaluated directly, given only the flux intensity
of the source, the background count rate, the duration of exposure,
the detection area and the PSF.
Hence it is determined by the object (i.e., the point source), the telescope
and the status of observation thus effects from data analysis programs are
minimized.

We use the cross-correlation significance as a reference.
For example, in our simulations an isolated point source of $1\;\mathrm{mCrab}$
flux has $2.42\sigma$ significance.

\subsection{Denoising}
\label{sect-denoise-methods}
\subsubsection{Linear methods}
Gaussian smoothing is often used in digital image processing to suppress the
noise at the cost of reduction in resolution.
Trade-off between noise suppression and resolution conservation is adjusted
through the standard deviation $\sigma$ of the Gaussian distribution serving
as the smoothing kernel function.
The best resolution of HXMT HE observed data is about $1.1^\circ$, limited by
the FWHM of its narrow-field collimator.
We set $\sigma$ to $28\mathrm{arcmin}$ so that the FWHM of the Gaussian kernel is
also $1.1^\circ$.

$N-$fold cross correlation transform ($N \ge 1$) can be used in DD method to
regularize the ill-posed problem, more specifically speaking, to ensure the
convergence as well as stability of the solution\citep{li2003dd}.
Here we put this technique in the denoising category.
We have tested $1-$fold and $2-$fold cross correlated DD methods respectively
in this article.

\subsubsection{Non-linear methods}
Non-local means denoising \citep{buades2005non} is an edge-preserving non-linear
denoising method.
To increase its performance we have implemented this method with
fast fourier transforms (FFTs).
The pixel-wisely evaluation of the general Euclidean distance between the $i-$th
pixel and other pixels of an image is replaced by
\begin{equation}
\vect{D}_i = \left\{\sum_k N_{i,k}\cdot W_k + \left(I^2 \ast W\right) -
  2\left[I \ast \left(N_i \cdot W\right)\right] \right\}^\frac{1}{2}\text{,}
  \label{eq-nlmeans}
\end{equation}
where $N_i$ is neighbourhood of the $i-$th pixel, $N_{i,k}$ is the $k-$ pixel
in the neighbourhood, $I$ is the image, $W$ is weight coefficient of the distance
function.
We use $7 \times 7$ pixels Gaussian kernel with standard deviation $\sigma=2$ (in
pixels) as the weight coefficient $W$ in our simulations.
We reduce the complexity of NLMeans method by computing convolutions in Eq.
\ref{eq-nlmeans} with FFTs instead of using searching windows.
 \citep{buades2008non}.

Median filter is another non-linear edge-preserving denoising method.
This method is effective in removing salt-and-pepper noise in digital images.
In HXMT observed data such noise is incurred typically by missing data or
charged particles of cosmic rays.
We fixed the size of the filter at $2^\circ \times 2^\circ$
(about $50 \times 50$ pixels) in HXMT HE data denoising.

The last non-linear denoising method included in this article is the
adaptive wavelet thresholding with multiresolution support \citep{starck2006}.
The multiresolution support of a noisy image is a subset that contains
significant coefficients only.
So wavelet coefficients that dominated by noise are discarded.
In this article we implemented the non-iterative algorithm.
The $5 \times 5$ $B_3$ spline wavelet is used for multiresolution
decomposition.

\section{Simulation and results}
\label{sect-simulations}
\subsection{In-orbit background simulation}
HXMT HE in-orbit background count rate $r_b$ ranges from $147.6\;\mathrm{counts/s}$ to
$210.7\;\mathrm{counts/s}$ \citep{li2009bgrd}.
We use a constant count rate $180\;\mathrm{counts/s}$ to simulate the average in-orbit
background of HXMT HE in this article.

\subsection{Source energy spectrum and telescope detection efficiency}
We use the formula proposed by \citet{massaro2000crab}  together with
parameters fitted by \citet{jourdain2009crab}
\begin{equation}
F(E) = 3.87 \times E^{-1.79-0.134\log_{10}\left(\sfrac{E}{20}\right)}
\end{equation}
to model energy spectra of Crab-like sources from $20\;\mathrm{keV}$ to $250\;\mathrm{keV}$,
where $E$ is in $\mathrm{keV}$ and the flux $F(E)$ is in $\mathrm{photons/s/cm^2}/\mathrm{keV}$.

The  detector efficiency of HXMT HE is derived from its simulated energy
response, as shown in Fig. \ref{fig-hxmt-he-eff}.
\begin{figure}[htbp]
\centering
\includegraphics[width=0.6\linewidth]{figures/hxmt-he-eff.eps}
\caption{HXMT HE detection efficiency}
\label{fig-hxmt-he-eff}
\end{figure}
Detection efficiency of HXMT HE to a Crab-like source is $67\%$.
Count rate of HXMT HE corresponds to the $1\;\mathrm{Crab}$ intensity is $1\,112\;\mathrm{counts/s}$,
given the detection area of HXMT HE is about $5\,100\;\mathrm{cm^2}$.

\subsection{PSF and modulation}
We use the diagram shown in Fig. \ref{fig-hxmt-he-psf} to simulate the PSF
$P_{\tan}(u,v)$ on the tangent plane.
\begin{figure}[htbp]
\centering
\includegraphics[width=0.6\linewidth]{figures/hxmt-he-psf.png}
\caption{Simulated HXMT HE PSF}
\label{fig-hxmt-he-psf}
\end{figure}

We use the concentric average of $P(\phi,\theta)$, namely,
\begin{equation}
S(\theta) = \frac{\int_{-\pi}^{\pi} P(\phi,\theta)\mathrm{d}\phi}{2\pi}\text{,}
\label{eq-psf-slope}
\end{equation}
and the cumulative sum of $S(\theta)$,
\begin{equation}
C(\theta) = \int_0^{\theta} S\left(\theta'\right)\mathrm{d}\theta'
\end{equation}
to characterize the radial fade-out of the PSF and the concentration of the PSF
respectively, as plotted in Fig.
\ref{fig-hxmt-he-psf-sl}.
\begin{figure}[htbp]
\centering
\includegraphics[width=0.6\linewidth]{figures/hxmt-he-psf-sl.eps}
\caption{Radial distribution of the simulated HXMT HE PSF}
\label{fig-hxmt-he-psf-sl}
\end{figure}
From Fig. \ref{fig-hxmt-he-psf-sl} we see that the FWHM of the simulated PSF is
about $1.7^\circ$ while $99.7\%$ responses occur in a diameter of $11^\circ$.

Despite the fact that the direct observed data is scientific events instead of 2-dimensional images, we start our simulation from \emph{simulated} observed data in the form of images defined on
2-dimensional cartesian pixel grid.
We use a $22^\circ \times 22^\circ$ model image for simulations. Given the
diameter of the PSF, the central $11^\circ \times 11^\circ$ region is fully
modulated, that is, all contributions to observed data in this region are from
the model image only. The surrounding $33^\circ \times 33^\circ$ region is
partially modulated, i.e., part of the contributions to observed data is this
region are from the model image.

The average exposure per unit solid angle is $382\;\mathrm{s/deg^2}$ in HXMT
half-year all-sky survey.
The partially modulated region is discretized by $N \times N$ pixel grid,
so $\tau \approx \sfrac{382 \times 33^2}{N^2}$.
The detection area of each HXMT HE detector is approximately
$300\;\mathrm{cm^2}$ so the total area of all the 17 HXMT HE detectors is
$5\,100\;\mathrm{cm^2}$.

\subsection{Results}
We have implemented several denoising methods (see Section \ref{sect-denoise-methods}).
Fig. \ref{fig-denoise-methods} shows denoising results by these methods.
\begin{figure}[htbp]
\centering
\includegraphics[width=0.6\linewidth]{figures/denoise-methods.png}
\caption{Denoising methods. From top to bottom, left to right:\
1. Model image;
2. Observed data of $10\;\mathrm{mCrab}$ point source.
3. Gaussian smoothed data, $\sigma = 28\mathrm{arcmin}$.
4. 1-fold cross correlated data.
5. NLMeans filter denoised data.
6. $2^\circ \times 2^\circ$ median filter denoised data.
7. $4^\circ \times 4^\circ$ median filter denoised data.
8. $3\sigma$ wavelet thresholding denoised data, with $B_3$ spline wavelet transform.}
\label{fig-denoise-methods}
\end{figure}

We have simulated $5\,000$ frames of observed data that contains the in-orbit
background counts only for each methods to estimate the corresponding flux
thresholds
by the method specified in Procedure \ref{item-threshold} of Section \ref{sect-estimation}.
From false detections we have obtained $2\sigma$ and $3\sigma$ flux thresholds,
see Table \ref{tab-threshold} for  results.

\begin{table}[htbp]
\centering
\begin{tabular}{c|r|r}
Denoising & $2\sigma$ thres., in $\mathrm{mCrab}$ & $3\sigma$ thres., in $\mathrm{mCrab}$ \\ \hline
w/o denoise & $0.520\pm0.001$ & $0.963\pm0.003$\\
1-fold CCT & $0.667\pm0.003$ & $1.192\pm0.010$\\
2-fold CCT & $0.919\pm0.014$ & $1.656\pm0.058$\\
Gaussian, $\sigma = 28\mathrm{arcmin}$ & $0.669\pm0.003$ & $1.133\pm0.008$\\
NLMeans & $0.845\pm0.004$ & $1.340\pm0.011$\\
Median filter, $2^\circ$ & $0.709\pm0.004$ & $1.185\pm0.015$\\
Median filter, $4^\circ$ & $0.665\pm0.005$ & $1.114\pm0.012$\\
Wavelet thres. & $0.0518\pm0.0002$ & $0.216\pm0.002$\\
\end{tabular}
\caption{$2\sigma$ and $3\sigma$ thresholds of point source detection}
\label{tab-threshold}
\end{table}

We have simulated $1\,000$ frames of observd data that contains a Crab-like
point source of $1\;\mathrm{mCrab}$, $1.25\;\mathrm{mCrab}$, $1.5\;\mathrm{mCrab}$, $1.75\;\mathrm{mCrab}$,
$2\;\mathrm{mCrab}$, $2.5\;\mathrm{mCrab}$, $3\;\mathrm{mCrab}$, $4\;\mathrm{mCrab}$, $5\;\mathrm{mCrab}$ and
$10\;\mathrm{mCrab}$ respectively, i.e., $10\,000$ frames of observed data in total.

With these simulated data we have estimated location accuracies as well as
detection efficiencies
by the method described in Procedure \ref{item-efficiency}
of the following methods:
\begin{enumerate}
\item DD without denoising,
\item DD with $1$-fold cross correlation,
\item DD with $2$-fold cross correlation,
\item DD with Gaussian smoothing, where $\sigma=28'$,
\item DD with NLMeans filtering, the size of the filter is
$7 \times 7$ and $\sigma=2$ (both parameters are in pixels),
\item DD with $2^\circ \times 2^\circ$ median filter,
\item DD with $4^\circ \times 4^\circ$ median filter, and
\item DD with adaptive wavelet thresholding.
\end{enumerate}
Implementation details of the above methods are in Sect. \ref{sect-denoise-methods}.

Table \ref{tab-accuracy} shows the location accuracies on simulated data
of $1\;\mathrm{mCrab}$, $2\;\mathrm{mCrab}$, $5\;\mathrm{mCrab}$ and $10\;\mathrm{mCrab}$ point sources.

\begin{table}[htbp]
\centering
\begin{tabular}{c|rrrrrrrrrr}
 & $1\;\mathrm{mCrab}$ & $2\;\mathrm{mCrab}$ & $5\;\mathrm{mCrab}$ & $10\;\mathrm{mCrab}$\\ \hline
w/o denoise & $95 \pm 6$ & $39 \pm 2$ & $9.6 \pm 0.2$ & $5.2 \pm 0.1$ \\
1-fold CCT & $104 \pm 3$ & $53 \pm 1$ & $20.5 \pm 0.4$ & $10.6 \pm 0.2$ \\
2-fold CCT & $149 \pm 5$ & $90 \pm 2$ & $36.9 \pm 0.8$ & $17.2 \pm 0.4$ \\
Gaussian, $\sigma = 28\mathrm{arcmin}$ & $90 \pm 4$ & $38 \pm 1$ & $11.1 \pm 0.2$ & $6.0 \pm 0.1$ \\
NLMeans & $107 \pm 4$ & $52 \pm 1$ & $21.0 \pm 0.3$ & $12.2 \pm 0.2$ \\
Median filter, $2^\circ$ & $91 \pm 4$ & $43 \pm 1$ & $14.3 \pm 0.3$ & $7.7 \pm 0.1$ \\
Median filter, $4^\circ$ & $100 \pm 3$ & $57 \pm 1$ & $23.3 \pm 0.4$ & $14.8 \pm 0.3$ \\
Wavelet thres. & $108 \pm 5$ & $48 \pm 2$ & $11.6 \pm 0.2$ & $5.9 \pm 0.1$ \\
\end{tabular}
\caption{Location accuracies, in arc minutes.}
\label{tab-accuracy}
\end{table}

Table \ref{tab-efficiency} shows the detection efficiencies on simulated data
of $1\;\mathrm{mCrab}$, $1.5\;\mathrm{mCrab}$, $2\;\mathrm{mCrab}$, $2.5\;\mathrm{mCrab}$ and $3\;\mathrm{mCrab}$
point sources.

Although the RL iteration itself we employed in DD method is total-counts-conservative, i.e., the sum of counts of each pixel is conserved after the iteration \citep{richardson1972},
the regularizations including the background constrains as well as various denoising techniques are not necessarily counts-conservative.
As a result, the absolute flux threshold
(Table \ref{tab-threshold})
for rejecting false detections doesn't reflect the sensitivity directly but the detection efficiency
(Table \ref{tab-efficiency})
does.

Errors of flux thresholds, location accuracies and detection efficiencies in Table \ref{tab-threshold},  Table \ref{tab-accuracy} and Table \ref{tab-efficiency} are calculated by bootstrapping. For
example, following Procedure \ref{item-threshold} a set of $5\,000$ frames of demodulated images are obtained, from which false detections are calculated and a histogram is plotted, where both the
$2\sigma$ and $3\sigma$ thresholds are determined. Now let's generate a new set with the same volume by resampling from the original set with replacement in order to calculate both the thresholds again.
We repeat this resampling process until we get enough thresholds calculated to estimate their standard deviations. In Table \ref{tab-threshold}, Table \ref{tab-accuracy} and Table \ref{tab-efficiency}
each of the errors is a standard deviation that calculated from $1\,000$ resampled sets.

\begin{table}[htbp]
\centering
\begin{tabular}{c|rrrrrrrrrr}
 & $1\;\mathrm{mCrab}$ & $1.5\;\mathrm{mCrab}$ & $2\;\mathrm{mCrab}$ & $2.5\;\mathrm{mCrab}$ & $3\;\mathrm{mCrab}$\\ \hline
w/o denoise & $29 \pm 2$ & $53 \pm 2$ & $77 \pm 1$ & $94 \pm 1$ & $98 \pm 0$ \\
1-fold CCT & $41 \pm 2$ & $71 \pm 1$ & $92 \pm 1$ & $99 \pm 0$ & $100 \pm 0$ \\
2-fold CCT & $27 \pm 2$ & $56 \pm 2$ & $79 \pm 2$ & $95 \pm 1$ & $99 \pm 0$ \\
Gaussian, $\sigma = 28\mathrm{arcmin}$ & $35 \pm 2$ & $66 \pm 2$ & $87 \pm 1$ & $98 \pm 0$ & $100 \pm 0$ \\
NLMeans & $36 \pm 2$ & $68 \pm 1$ & $91 \pm 1$ & $100 \pm 0$ & $100 \pm 0$ \\
Median filter, $2^\circ$ & $37 \pm 2$ & $64 \pm 2$ & $86 \pm 1$ & $97 \pm 1$ & $100 \pm 0$ \\
Median filter, $4^\circ$ & $39 \pm 2$ & $69 \pm 1$ & $89 \pm 1$ & $98 \pm 0$ & $100 \pm 0$ \\
Wavelet thres. & $27 \pm 1$ & $51 \pm 2$ & $79 \pm 1$ & $95 \pm 1$ & $100 \pm 0$ \\
\end{tabular}
\caption{Detection efficiencies, in percentages.}
\label{tab-efficiency}
\end{table}

A comprehensive summary of all the tested methods on all simulated data is
shown in Fig. \ref{fig-single-source}.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.995\linewidth]{figures/single-source-faint.eps}
\caption{Single source location accuracy and detection efficiency.
MF stands for median filter while WT for wavelet thresholding.}
\label{fig-single-source}
\end{figure}

\section{Conclusion}
According to the results from our tests no denoising method shows significant
advantage over $1$-fold cross correlated DD method in single point source
detection efficiency.
Therefore it's suggested that $1$-fold cross correlation should be the default
regularization method for single point source detection in HXMT imaging data
analysis.

In the other hand, the location accuracy can be improved with alternative
denoising methods, such as median filter, wavelet thresholding, Gaussian
smoothing with little kernel, or without denoising, according to the results
in this work.

This article is focused on the single point source detection of HXMT imaging
data analysis, where other interesting topics can not be all covered.
Although the alternative denoising methods have been out-performed more or less
by the default $1$-fold cross correlation in their contributions to the detection
efficiency, they have shown certain advantages in location accuracy, and these
features are promising for locating bright transients, multiple sources
resolving, and so on.

\section*{Acknowledgements}
In this work we made use of SciPy\citep{scipy}, PyFITS and AIRE.
PyFITS is a product of the Space Telescope Science Institute, which is operated
by AURA for NASA.
AIRE is a set of computing facilities initiated by Tsinghua Centre for
Astrophysics.
This work was supported by \emph{National Natrual Science Foundation of China}
(NSFC) under grants No. 11373025, No. 11173038 and No. 11403014 as well as
\emph{Tsinghua University Initiative Scientific Research Program} under grant
No. 20111081102.

\bibliographystyle{raa}