text:
\section{Introduction}
In recent years, fractional differential equations (FDEs) have attracted increasing interest due to
the fact that many mathematical problems in science and engineering can be modeled by fractional differential equations, see e.g.,\ \cite{Podlubny1999, Kilbas2006, Diethelm2010}.

In this paper, we consider a $d$-dimensional fractional differential equation involving the \emph{Caputo derivative of order $\alpha \in (0,1)$} on an finite interval $J:= [0,T] $ or on the half real line $J:=[0,\infty)$:
\begin{equation}\label{mainEq}
^{C\!}D_{0+}^\alpha x(t)=f(t,x(t)),
\end{equation}
where $x: J \rightarrow \mathbb{R}^d$ is a vector function,  $f: J\times \mathbb{R}^d\rightarrow \mathbb{R}^d$ is a continuous vector-valued function. A continuous function $x: J \rightarrow \mathbb{R}^d$ is called a {\em solution} of \eqref{mainEq} if it satisfied \eqref{mainEq} for all $0<t\in J$. The value $x(0)$  is called the initial value of the solution $x(\cdot)$.

We are interested in the dynamical properties of \eqref{mainEq}: we want to know whether \eqref{mainEq} generates a dynamical system so that the tools and methods of the classical theory of dynamical systems are applicable in the investigation of FDEs. Another important problem in the theory of FDEs, which is closely related to the problem of generation of dynamical systems by FDEs is the question of
 whether two different trajectories of a FDE can intersect.
We solved both these problems. Namely we show that
one-dimensional FDEs and triangular FDEs generate nonlocal fractional dynamical systems, whereas a higher dimensional FDE does, in general, not generate a nonlocal dynamical system. Correspondingly, two different trajectories of a one-dimensional or a triangular FDE cannot meet, whereas different trajectories of a high dimensional FDE may intersect each other. As a byproduct of our investigation we get  lower bounds for the solutions of one-dimensional FDEs and triangular linear FDEs.

The question of whether solutions of \eqref{mainEq} can intersect are treated in Diethelm~\cite{Diethelm2008,Diethelm2010}, Diethelm and Ford~\cite{DiethelmFord2012}, Agarwal \textit{et al.}~\cite{Agarwal2010}, Hayek \textit{et al.}~\cite{Hayek1999}, and Bonilla, Rivero and Trujillo~\cite{Bonilla2007}. Note that in the case of ordinary differential equations (ODEs) it is well known that two trajectories of an ODE either coincide or they do not intersect; the authors mentioned above proved that similar results hold for fractional differential equations of order $0<\alpha<1$. Note that the main difficulty of the problem for FDEs is the nonlocal nature (or history memory) of solutions of FDEs. The above authors used various tools to deal with the FDE case. However we notice flaws in their proofs which make their proofs incomplete. We will present some discussion about this matter in Section~\ref{sec.separation1d} of the paper.

The paper is organized as follows:
Section \ref{sec.preliminaries} is a preparatory section where we present some basic notions from fractional calculus and the theory of  FDEs.
Section \ref{sec.separation1d} is devoted to results on separation of solutions of one-dimensional FDEs; in this section we also give discussion about flaws in the proofs of separation of solutions of FDEs in the above mentioned papers.
In Section \ref{sec.1dimNonlocalDS} we present results on generation of nonlocal fractional dynamical systems by one-dimensional FDEs.
Section~\ref{sec.triangular} is devoted to high dimensional triangular systems of FDEs, where based on the results in Section~\ref{sec.1dimNonlocalDS} we show that a triangular system of FDEs does generate a nonlocal fractional dynamical system.
 In Section~\ref{sec.d-dim} we show that a higher dimensional FDE does, in general, not generate a nonlocal dynamical system; two different trajectories of a high dimensional FDE may intersect each other in finite time.

\section{Preliminaries}\label{sec.preliminaries}

We start this section by briefly recalling a framework of fractional calculus and fractional differential equations. We refer the reader to the books \cite{Diethelm2010,Kilbas2006} for more details.

Let $\mathbb{R}^d$ be the standard $d$-dimensional Euclidean space equipped with usual Euclidean norm.
We denote by $\mathbb{R}_+$ the set of all nonnegative real numbers.
Denote by $C([0,\infty);\mathbb{R}^d)$ the space of continuous functions from $[0,\infty)$ to $\mathbb{R}^d$,
 and by $\left(C_\infty(\mathbb{R}^d),\|\cdot\|_\infty\right)\subset C([0,\infty);\mathbb{R}^d)$ the space of all continuous functions $\xi:\mathbb{R}_+\rightarrow \mathbb{R}^d$
which are uniformly bounded on $\mathbb{R}_+$, i.e.\
\[
\|\xi\|_\infty:=\sup_{t\in \mathbb{R}_+}\|\xi(t)\|<\infty.
\]
It is well known that $\left(C_\infty(\mathbb{R}^d),\|\cdot\|_\infty\right)$ is a Banach space.

Let $\alpha>0$ and $[a,b]\subset \mathbb{R}$. Let $x:[a,b]\rightarrow \mathbb{R}$ be a measurable function such that $x\in L^1([a,b])$, i.e.,\ $\int_a^b|x(\tau)|\;d\tau<\infty$. Then, the Riemann-Liouville integral of order $\alpha$ is defined by
\[
I_{a+}^{\alpha}x(t):=\frac{1}{\Gamma(\alpha)}\int_a^t(t-\tau)^{\alpha-1}x(\tau)\;d\tau\quad \hbox{ for } t\in (a,b],
\]
where the Gamma function $\Gamma:(0,\infty)\rightarrow \mathbb{R}$ is defined as
\[
\Gamma(\alpha):=\int_0^\infty \tau^{\alpha-1}\exp(-\tau)\;d\tau,
\]
see e.g., Diethelm~\cite{Diethelm2010}. The corresponding Riemann--Liouville fractional derivative of order $\alpha $ is given by
\[
^{R\!}D_{a+}^\alpha x(t):=(D^m I_{a+}^{m-\alpha}x)(t) \qquad\hbox{ for } t\in (a,b],
\]
where $D=\dfrac{d}{dt}$ is the usual derivative and $m:=\lceil\alpha\rceil$ is the smallest integer bigger or equal to $\alpha$. On the other hand, the \emph{Caputo fractional derivative } $^{C\!}D_{a+}^\alpha x$ of a function $x\in C^m([a,b])$, which was introduced by Caputo (see e.g., Diethelm~\cite{Diethelm2010}), is defined by
\[
^{C\!}D_{a+}^\alpha x(t):=(I_{a+}^{m-\alpha}D^m x)(t),\qquad \hbox{ for } t\in (a,b].
\]
The Caputo fractional derivative of a $d$--dimensional vector function $x(t)=(x_1(t),\dots,x_d(t))^{\mathrm {T}}$ is defined component-wise as $$^{C\!}D^\alpha_{0+}x(t)=(^{C\!}D^\alpha_{0+}x_1(t),\dots,^{C\!}D^\alpha_{0+}x_d(t))^{\mathrm {T}}.$$

It is well-known that the initial value problem of the FDE \eqref{mainEq} is equivalent to a Volterra integral equation of the second kind. Namely we have the following

\begin{lemma}\label{lem.integral.form}
A continuous function $x: J\rightarrow \mathbb{R}$ is a solution of the FDE \eqref{mainEq} with the initial value condition $x(0)=x_0$ if and only if it is a solution of the Volterra integral equation of the second kind
\begin{equation}\label{eqn.volterra}
x(t) = x_0 +\frac{1}{\Gamma(\alpha)}\int_0^t (t-\tau)^{\alpha-1}f(\tau,x(\tau)) d\tau.
\end{equation}
\end{lemma}
For a proof we refer to Diethelm~\cite[Lemma 6.2, p. 86]{Diethelm2010} and Kilbas \textit{et al.}~\cite[Theorem 3.24, p. 199]{Kilbas2006} for the one-dimensional case; the multidimensional case is just componentwise application of the one-dimensional case.

\section{Separation of trajectories of solutions of one-dimensional FDE}\label{sec.separation1d}

\subsection{Two different trajectories of a one-dimensional FDE do not meet}
In this section, we consider the one-dimensional case of the system \eqref{mainEq}, i.e. the following FDE
\begin{equation}\label{1dimEq}
^{C\!}D_{0+}^{\alpha}x(t)=f(t,x(t)),
\end{equation}
where $f:J\times \mathbb{R} \to \mathbb{R}$ is a continuous function. Assume that $f$ satisfies the following Lipschitz condition on the second variable: there exists a nonnegative continuous function $L: J \rightarrow \mathbb{R}_+$ such that
\begin{equation}\label{eqn.Lip}
|f(t,x)-f(t,y)| \leq L(t)|x-y|\qquad \hbox{for all}\quad t\in J\quad\hbox{and all}\quad x,y\in\mathbb{R}.
\end{equation}
It is well known that under the Lipschitz condition \eqref{eqn.Lip} the initial value problem for \eqref{1dimEq} has unique solution defined on the whole interval $J$ for any given initial value (see, e.g. Baleanu and Mustafa~\cite[Theorem 2]{Baleanu2010}, Tisdell~\cite[Theorem 6.4]{Tisdell2012}) and Diethelm~\cite[Theorem 6.8]{Diethelm2010}).
We will show that any two solutions of \eqref{1dimEq} either coincide or do not intersect on $J$. For this we need the variation of constant formula for FDEs and the comparision principle, which are useful technical tools for investigation of FDEs.

\begin{lemma}[Variation of constants formula for FDEs]\label{lem.variationofconstant}
If the function $f$ in the FDE \eqref{1dimEq}
is of the form
$$
f(t,x) = Mx + g(t,x)
$$
for some fixed $M\in \mathbb{R}$ and all $t\in J$, $x\in \mathbb{R}$, then the solution $x(\cdot)$ of \eqref{1dimEq} with the initial value $x(0)=x_0$ satisfies for all $t\in J$ the formula
$$
x(t) = E_\alpha(Mt^\alpha)x_0 + \int_0^t (t-\tau)^{\alpha-1} E_{\alpha,\alpha}(M(t-\tau)^{\alpha})g(\tau,x(\tau)) d\tau,
$$
where
$$
E_{\alpha,\beta}(z):=\sum_{k=0}^\infty \frac{z^k}{\Gamma(\alpha k+\beta)},\qquad E_{\alpha}(z):=E_{\alpha,1}(z),\qquad z\in\mathbb{C},
$$
are Mittag-Leffler functions and $\Gamma(\cdot)$ is the Gamma function.
\end{lemma}
For a proof of this theorem we refer the reader to Kilbas \textit{et al.}~\cite[Theorem 5.15, p.~323]{Kilbas2006} and Diethelm~\cite[Theorem 7.2, p.~135]{Diethelm2010}.

\begin{theorem}[Comparision Principle~\cite{Lakshmikantham2008}]\label{comparision}
Let $0<q<1$,  $v,w\in C(J,\mathbb{R})$ are continuous function on $J$, $g\in C(J\times \mathbb{R},\mathbb{R})$ is a continuous two-variable funtion on $J\times \mathbb{R}$, and \\
(i) $v(t) \leq v(0) + \dfrac{1}{\Gamma(q)}\int_0^t (t-s)^{q-1}g(s,v(s)) ds$ for all  $t\in J$,\\
(ii) $w(t) \geq w(0) + \dfrac{1}{\Gamma(q)}\int_0^t (t-s)^{q-1}g(s,w(s)) ds$ for all $t\in J$.\\
Suppose further that $g(t,x)$ is nondeacreasing in $x$ for each $t\in J$, and
\begin{equation}\label{eqn.compar1}
v(0) < w(0).
\end{equation}
Then we have
\begin{equation}\label{eqn.compar2}
v(t) < w(t) \qquad\hbox{for all}\quad t\in J.
\end{equation}
\end{theorem}
For convenience of the reader we present here the proof of this theorem by
 Lakshmikantham and Vatsala~\cite[Theorem 2.1]{Lakshmikantham2008} with a small modification.
\begin{proof}
Suppose that \eqref{eqn.compar2} is not true. Then, because of the continuity of $v(\cdot)$ and $w(\cdot)$
there exists $0< t_1\in J$ such that
\begin{equation}\label{eqn.compar3}
v(t_1) = w(t_1), \qquad\hbox{and}\qquad v(t) <w(t)\quad\hbox{for all}\quad  0\leq t <t_1.
\end{equation}
Using the nondecreasing nature of $g$ and \eqref{eqn.compar3}, taking into account  \eqref{eqn.compar1}, we get
\begin{eqnarray*}
v(t_1) &\leq& v(0) + \dfrac{1}{\Gamma(q)}\int_0^{t_1} (t_1-s)^{q-1}g(s,v(s)) ds\\
&\leq& v(0) + \dfrac{1}{\Gamma(q)}\int_0^{t_1} (t_1-s)^{q-1}g(s,w(s)) ds\\
&<& w(0) + \dfrac{1}{\Gamma(q)}\int_0^{t_1} (t_1-s)^{q-1}g(s,w(s)) ds\\
&\leq& w(t_1),
\end{eqnarray*}
which contradict the condition $v(t_1) = w(t_1)$ in \eqref{eqn.compar3}. Hence \eqref{eqn.compar2} is valid and the proof is complete.
\end{proof}

Now we are in a position to prove our main theorem of this section.
\begin{theorem}[Different trajectories do not meet]\label{separation1dim}
 Assume that $f$ satisfies the Lipschitz condition \eqref{eqn.Lip}. Then for any two different initial values $x_{10}\not= x_{20}$ in $\mathbb{R}$ the trajectories of solutions of the FDE \eqref{1dimEq} do not meet on $J$, i.e., the solutions $x_1(\cdot),x_2(\cdot)$ of
\eqref{1dimEq} starting from $x_{10}=x_1(0)$ and $x_{20}=x_2(0)$ verify $x_1(t)\not= x_2(t)$ for all $t\in J$.
\end{theorem}
\begin{proof}
For definiteness we assume that
\begin{equation}\label{eqn.1dim1}
x_1(0)=x_{10} < x_{20}=x_2(0).
\end{equation}
To prove the theorem we show that $x_1(t)< x_2(t)$ for all $t\in J$.
Suppose that this is not true, then there is $0<T_1\in J$ such that $x_1(T_1)\geq x_2(T_1)$. By the continuity property of $x_1(\cdot), x_2(\cdot)$ and \eqref{eqn.1dim1} there is $T_2>0$ such that
\begin{eqnarray}
x_1(T_2)&=& x_2(T_2), \quad\hbox{and} \label{eqn.1dim2}\\
x_1(t) &<& x_2(t)\qquad\hbox{for all}\quad 0\leq t < T_2. \label{eqn.1dim3}
\end{eqnarray}
Put
\begin{eqnarray}
M &:=& \max_{0\leq t\leq T_2} L(t),  \label{eqn.1dim4}\\
g(t,x) &:=& f(t,x) + Mx\qquad \hbox{for all}\quad 0\leq t\leq T_2, x\in\mathbb{R}.\label{eqn.1dim5}
\end{eqnarray}
Then $g: [0,T_2]\times \mathbb{R} \rightarrow \mathbb{R}$ is nondecreasing in $x$ for each $t\in [0,T_2]$ since by \eqref{eqn.Lip}, \eqref{eqn.1dim4} and \eqref{eqn.1dim5}, for any $t\in [0,T_2]$ and $x\leq y$ we have
$$
g(t,y)-g(t,x) = M(y-x) +f(t,y)-f(t,x)\geq (M-L(t))(y-x)\geq 0.
$$
By virtue of Lemma~\ref{lem.variationofconstant}, since $x_1(\cdot)$ and $x_2(\cdot)$ are solutions of \eqref{1dimEq}, on the interval $[0,T_2]$ we have the following equations
\begin{eqnarray}
x_1(t) &=& E_\alpha(-Mt^\alpha)x_{10} \nonumber\\
&&\hspace*{0.1cm}+\;\int_0^t (t-\tau)^{\alpha-1}E_{\alpha,\alpha}(-M(t-\tau)^\alpha)g(\tau,x_1(\tau)) d\tau,\label{eqn.1dim6}\\[6pt]
x_2(t) &=& E_\alpha(-Mt^\alpha)x_{20} \nonumber\\
&&\hspace*{0.1cm}+\;\int_0^t (t-\tau)^{\alpha-1}E_{\alpha,\alpha}(-M(t-\tau)^\alpha)g(\tau,x_2(\tau)) d\tau.\label{eqn.1dim7}
\end{eqnarray}
Since $E_{\alpha,\alpha}(s)>0$ for all $s\in\mathbb{R}$ (see, e.g., Cong \textit{et al.}~\cite[Lemma 2]{Cong}), the function
$E_{\alpha,\alpha}(-M(t-\tau)^\alpha)g(\tau,x)$
is nondecreasing in the variable $x$ for all $0\leq\tau\leq t\leq T_2$. Therefore, a small modification of
Theorem~\ref{comparision} is applicable to the pair of integral equations \eqref{eqn.1dim6}-\eqref{eqn.1dim7} on $[0,T_2]$ and gives $x_1(T_2) < x_2(T_2)$. Thus we arrive at a contradiction. Consequently, conclusion of the theorem is true and the proof is completed.
\end{proof}

 \begin{remark}
 Theorem~\ref{separation1dim} provides a full solution to Conjecture 1.2 of
 Diethelm~\cite{Diethelm2008}. Note that in \cite{Diethelm2008} Diethelm provides a partial solution to the problem of separation of trajectories: he proved Conjecture 1.2 in \cite{Diethelm2008} under some restrictive conditions, see \cite[Theorem 2.2]{Diethelm2008}.
 \end{remark}

\subsection{Discussion on the problem of separation of solutions of a general FDE}\label{subsec.discussion}

Diethelm~\cite[Theorem 6.12]{Diethelm2010}
 formulated and proved a theorem on separation of solutions of a FDE, what is the same as our Theorem~\ref{separation1dim}. For the proof he used a fixed point theorem on a short interval of time, then used induction to extend the result to the whole (long) interval of time under consideration. However, his proof contains a flaw in the induction part: the passage to the next step from $N=1$ to $N=2$ does not work since the FDE is history dependent and on the next subinterval the argument leading to contractive mapping do not work like on the first subinterval of time.

Diethelm, in his joint work with Ford \cite[Theorems 3.1 and 4.1]{DiethelmFord2012}, gave and alternative (to \cite[Theorem 6.12]{Diethelm2010}) proof of the theorem on separation of solution of FDE. Namely, instead of induction forward in time from one subinterval to the next as in \cite{Diethelm2010} they used induction backward in time from one subinterval to the foregoing interval. However, their argument in the first step of induction does not work:
in the notation of \cite{DiethelmFord2012}, the equation (13) is equivalent to equation (8) only if we consider (8) and (13) in the whole interval $[0,b]$; actually the solution $y(t)$, considered as solution of the terminal problem (13), depends on the future value $y(b)=c$, hence we cannot say that the function $g(t)$ on pages 29--30 of \cite{DiethelmFord2012} are independent of the restriction of the function $y$ to $[T_{N-1},T_N]$, therefore the contraction property of the map $\mathcal F$  in the formula (14) of \cite{DiethelmFord2012} does not lead to existence and uniqueness of solution of (14) as claimed in \cite{DiethelmFord2012}. Therefore, the proof of Theorem 3.1 in \cite{DiethelmFord2012} is incomplete, hence so is Theorem 4.1 of \cite{DiethelmFord2012}.

Note that the proof of Theorem 6.12 of Diethelm~\cite{Diethelm2010}  is correct for a first "short" interval of time where the smallness of time combined with bounded Lipschitz condition of $f$ make certain operator contractive, hence two different solutions can not meet.
By the way, we notice that the continuity alone can also assure the non-intersection of different trajectories in a short time interval.

Hayek, Trujillo, Rivero, Bonilla and
Boreno \cite[Theorem 3.1]{Hayek1999} have proved the separation theorem on a "short" interval of time. We also note that Section 4 of \cite{Hayek1999} is false because of the history dependence of solutions of FDE, which prevents us from usage of "usual method of prolongation" (like in ODE) as claimed by the authors of \cite{Hayek1999}.

Bonilla, Rivero and Trujillo~\cite{Bonilla2007} treated higher dimensional linear systems of FDEs and in Section 3 of \cite{Bonilla2007} they relied on the above mention result by Hayek \textit{et al.}~\cite{Hayek1999}, hence there are gaps in their proofs of some results, namely the proofs of Theorem 1 on page 71 and of Propositions 1 and 2 on page 72 of \cite{Bonilla2007} are incomplete. For a counter example see Section~\ref{sec.d-dim}.

\section{One-dimensional FDEs generate nonlocal dynamical systems}\label{sec.1dimNonlocalDS}

In this section, based on the results on separation of trajectories presented in Section~\ref{sec.separation1d} we show that one-dimensional FDEs generate nonlocal dynamical systems, hence tools and methods of the classical theory of dynamical systems can be applied.

\subsection{Bounds for solutions of FDEs}
First we formulate and prove a lower bound for solutions of FDE, which provides us with a better understanding of geometry of solutions of one-dimensional fractional differential equations.
\begin{theorem}[Convergence rate for solutions of 1-dim FDEs]\label{thm.convergence.rate}
Assume that $f$ satisfies the Lipschitz condition \eqref{eqn.Lip}.
For any two solutions $x_1(\cdot), x_2(\cdot)$ of the FDE \eqref{1dimEq} and any $t\in J$ the following estimate holds
$$
|x_2(t)-x_1(t)| \geq |x_2(0) - x_1(0)| E_\alpha\big(- (\max_{0\leq \tau\leq t} L(\tau)) t^\alpha\big).
$$
\end{theorem}
\begin{proof}
 For definiteness we assume $x_2(0) \geq x_1(0)$. Then by Theorem \ref{separation1dim}, we have $x_2(t)\geq x_1(t)$ for any $t\in J$. For an arbitrary $t\geq 0$ but fixed, put
$M_t:=\sup_{0\leq s\leq t}L(s)$. On the interval $[0,t]$, repeat the arguments of the proof of Theorem~\ref{separation1dim}, we get
\begin{eqnarray*}
x_2(s) -x_1(s) &=& E_\alpha(-M_ts^\alpha)(x_2(0)-x_1(0)) +\\
&&\hspace*{-2.1cm}+\;\int_0^ s (s-\tau)^{\alpha-1}E_{\alpha,\alpha}(-M_t(s-\tau)^\alpha)( g(\tau,x_2(\tau))- g(\tau,x_1(\tau))) d\tau\\[6pt]
&\geq&\hspace*{0.1cm}  E_\alpha(-M_ts^\alpha)(x_2(0)-x_1(0)).
\end{eqnarray*}
Take $s=t$, and the proof is complete.
\end{proof}

\begin{corollary}[Lower bound for solutions of 1-dim FDEs]\label{cor.lower.bound}
Assume that $f$ satisfies the Lipschitz condition \eqref{eqn.Lip}. Assume additionally that $f(t,0)=0$ for all $t\in J$, $x\in\mathbb{R}$. Then for
 any solution $x(\cdot)$ of the FDE \eqref{1dimEq} and any $t\in J$ the following estimate holds
$$
|x(t)| \geq |x(0)| E_\alpha\big(- (\max_{0\leq \tau\leq t} L(\tau)) t^\alpha\big).
$$
\end{corollary}
\begin{proof}
Since $f(t,0)=0$ the FDE \eqref{1dimEq} has trivial solution. Apply Theorem~\ref{thm.convergence.rate} for the pair of $x(\cdot)$ and the trivial solution of \eqref{1dimEq}.
\end{proof}

For the divergence rate and upper bound for solutions of the FDEs the following statements are easy modifications of well known results.

\begin{theorem}[Divergence rate for solutions of 1-dim FDEs]\label{thm.divergence.rate}
Assume that $f$ satisfies the Lipschitz condition \eqref{eqn.Lip}.
For any two solutions $x_1(\cdot), x_2(\cdot)$ of the FDE \eqref{1dimEq} and for any $t\in J$ the following estimate holds
$$
|x_2(t)-x_1(t)| \leq |x_2(0) - x_1(0)| E_\alpha\big( (\max_{0\leq \tau\leq t} L(\tau)) t^\alpha\big).
$$
\end{theorem}
\begin{proof}
Let us denote $M_t = \max_{0\leq \tau\leq t} L(\tau)$. By Lemma~\ref{lem.integral.form}, for all $t\in J$ we have
\begin{eqnarray*}
x_1(t) &=& x_1(0) +\frac{1}{\Gamma(\alpha)}\int_0^t (t-\tau)^{\alpha-1}f(\tau,x_1(\tau)) d\tau,\\
x_2(t) &=& x_2(0) +\frac{1}{\Gamma(\alpha)}\int_0^t (t-\tau)^{\alpha-1}f(\tau,x_2(\tau)) d\tau.
\end{eqnarray*}
Therefore, for any $t\in J$ and  any $0\leq s\leq t$ we have
\begin{eqnarray*}
|x_2(s)-x_1(s)| &\leq& |x_2(0)-x_1(0)| + \\
&&\hspace*{0.5cm} + \frac{1}{\Gamma(\alpha)}\int_0^s (s-\tau)^{\alpha-1}(f(\tau,x_2(\tau)) - f(\tau,x_1(\tau)) ) d\tau\\[5pt]
&\leq& |x_2(0)-x_1(0)| + \frac{1}{\Gamma(\alpha)}\int_0^s (s-\tau)^{\alpha-1} M_t|x_2(\tau)-x_1(\tau)| d\tau.
\end{eqnarray*}
By virtue of the Gronwall-type inequality for FDE (see Diethelm~\cite[Theorem 6.19, p. 111]{Diethelm2010} and Tisdell~\cite[Lemma 3.1, p. 288]{Tisdell2012}), this implies
$$
|x_2(s)-x_1(s)| \leq |x_2(0) - x_1(0)| E_\alpha\big( M_t s^\alpha\big).
$$
Put $s=t$ and since $t\in J$ is arbitrary we have the theorem proved.
\end{proof}

\begin{corollary}[Upper bound for solutions of 1-dim FDEs]\label{cor.upper.bound}
 Assume that $f$ satisfies the Lipschitz condition \eqref{eqn.Lip}. Assume additionally that $f(t,0)=0$ for all $t\in J$, $x\in\mathbb{R}$. Then for
 any solution $x(\cdot)$ of the FDE \eqref{1dimEq} and any $t\in J$ the following estimate holds
$$
|x(t)| \leq |x(0)| E_\alpha\big((\max_{0\leq \tau\leq t} L(\tau)) t^\alpha\big).
$$
\end{corollary}
\begin{proof}
Since $f(t,0)=0$ the FDE \eqref{1dimEq} has trivial solution. Apply Theorem~\ref{thm.divergence.rate} for the pair of $x(\cdot)$ and the trivial solution of \eqref{1dimEq}.
\end{proof}

\begin{remark}\label{rem.multidim}
It is easily seen that Theorem~\ref{thm.divergence.rate} and Corollary~\ref{cor.upper.bound} hold true also for the case of a high dimensional system of FDEs.
\end{remark}

\subsection{One-dimensional FDEs generate two-parameter flows}\label{subsec.1dim}

Now we are in a position to show that one-dimensional FDEs generate two-parameter flows. First we define the evolution mappings of  \eqref{1dimEq}.
\begin{definition}\label{dfn.0evolution-mapping}
The mapping
\begin{equation}
\Phi_{0,T_1} :\mathbb{R}\rightarrow \mathbb{R}, \quad x_0 \mapsto x(T_1),
\end{equation}
where $x_0\in\mathbb{R}$ is an arbitrary initial value of \eqref{1dimEq}, $x(\cdot)$ is the solution of \eqref{1dimEq} starting from $x(0)=x_0$ and $x(T_1)$ is the evaluation of $x(\cdot)$ at $T_1$, is called {\em the evolution mapping of \eqref{1dimEq}}.
\end{definition}

\begin{definition}\label{2parameterflow}
A two-parameter family of mappings
$$
\phi_{s,t}(\cdot) : \mathbb{R}\rightarrow \mathbb{R}, \qquad s,t\in J,
$$
 is called  {\em a two-parameter flow in $\mathbb{R}$} if $\phi_{s,t}(x)$ is continuous as a function of three variables $s,t\in J, x\in\mathbb{R}$, for any fixed $s,t\in J$ the mapping $\phi_{s,t}$ is a homeomorphism of $\mathbb{R}$, and  this family satisfies the following flow property
$$
\phi_{s,t}\circ \phi_{u,s} = \phi_{u,t}\qquad\hbox{for all}\quad u,s,t\in J.
$$
\end{definition}

\begin{theorem}[One-dimensional FDEs generate two-parameter flows in $\mathbb{R}$]\label{thm.2param.flow}
The following statements hold for the one-dimensional FDE \eqref{1dimEq}.
\begin{enumerate}
\renewcommand{\labelenumi}{(\roman{enumi})}
\item The evolution mapping $\Phi_{0,t}$ generated by \eqref{1dimEq} is a bijection for any $t\in J$.
\item The FDE \eqref{1dimEq} generates a two-parameter family of bijections on $J$ by its evolution mappings as follows
\begin{equation}\label{eqn.2flow}
\Phi_{s,t} := \Phi_{0,t}\circ \Phi^{-1}_{0,s}\qquad\hbox{for all}\quad s,t\in J,
\end{equation}
where $\Phi_{0,\cdot}$ is the evolution mapping of \eqref{1dimEq} defined in Definition~\ref{dfn.0evolution-mapping}.
\item The family $\Phi_{s,t}$, $s,t\in J$, generated by the FDE \eqref{1dimEq} is a two-parameter flow in $\mathbb{R}$.
\item If $f$ is linear in $x$ then the two-parameter flow generated by the FDE \eqref{1dimEq} is a flow of linear operators.
\end{enumerate}
\end{theorem}
\begin{proof}
(i) Fix $T_1\in J$ arbitrary. By Theorem~\ref{separation1dim}, $\Phi_{0,T_1}$ is injective.

To show that
$\Phi_{0,T_1}$ is surjective, it suffices to show that for an arbitrary $x^*\in\mathbb{R}$ the boundary value problem
\begin{eqnarray}
&^{C\!}D_{0+}^{\alpha}x(t)=f(t,x(t)),\label{eqn.surj1}\\
&x(T_1) = x^*,\label{eqn.surj2}
\end{eqnarray}
where $t\in [0,T_1]$, $f$ is continuous and satisfies the Lipschitz condition \eqref{eqn.Lip}, has a continuous solution. Put
$$
M_1 := \max_{0\leq t\leq T_1} L(t),
$$
where $L(t)$ is determined from \eqref{eqn.Lip}. Let us denote by $\hat x(\cdot)$ the solution of the FDE \eqref{eqn.surj1} satisfying the initial value condition $\hat x(0)=0$. Put
$$
M_2 :=  \max_{0\leq t\leq T_1} |\hat x(t)|, \qquad M_3:= |x^*-\hat{x}(T_1)| +M_2,
$$
and
$$
M_4 := \frac{M_3}{E_\alpha(-M_1T_1^\alpha)}, \qquad M_5 := M_2 + M_4 E_\alpha(M_1T_1^\alpha)+1.
$$
Clearly $M_2<M_3<M_4<M_5$. Define a function on $[0,T_1]\times \mathbb{R}$ as follows
\begin{equation}\label{eqn.surj3}
\hat{f}(t,x) =
\begin{cases} f(t,x),&\quad\hbox{if}\quad |x| \leq M_5\\
f(t,M_5\frac{x}{|x|}),&\quad \hbox{if}\quad |x| > M_5,
\end{cases}
\end{equation}
and consider the boundary problem of the FDE
\begin{eqnarray}
&^{C\!}D_{0+}^{\alpha}x(t)=\hat{f}(t,x(t)),\label{eqn.surj4}\\
&x(T_1) = x^*.\label{eqn.surj5}
\end{eqnarray}
Clearly $\hat{f}$ is Lipschitz continuous with Lipschitz constant $L(t)$ and bounded on $[0,T_1]\times \mathbb{R}$,
hence the boundary problem \eqref{eqn.surj4}--\eqref{eqn.surj5} has at least one solution, say $x_1(\cdot)$ (see \cite[Theorem 8]{Benchohra2008}).
We show that $x_1(\cdot)$ is the required solution of \eqref{eqn.surj1}--\eqref{eqn.surj2} and we will be done. To this end notice that for all $t\in [0,T_1]$ we have
$$
|\hat{x}(t)| \leq M_2 < M_5 \qquad\hbox{hence}\quad f(t,\hat{x}(t)) = \hat{f}(t,\hat{x}(t)).
$$
Therefore, $\hat{x}(\cdot)$ is the solution of the FDE \eqref{eqn.surj4} satisfying the initial value condition $\hat{x}(0)=0$.
Apply Theorem~\ref{thm.convergence.rate} to the solutions $\hat{x}(\cdot), x_1(\cdot)$ of the FDE \eqref{eqn.surj4} we get for any $t\in [0,T_1]$ the inequality
$$
|\hat{x}(t)-x_1(t)| \geq |\hat{x}(0) - x_1(0)| E_\alpha\big(- (\max_{0\leq \tau\leq t} L(\tau)) t^\alpha\big) \geq |x_1(0)| E_\alpha(-M_1 T_1^\alpha).
$$
Substituting $t=T_1$ we get
$$
|\hat{x}(T_1) - x_1(T_1)| \geq |x_1(0)|E_\alpha(-M_1 T_1^\alpha),
$$
hence
$$
|x_1(0)| \leq \frac{|\hat{x}(T_1) - x^*|}{E_\alpha(-M_1 T_1^\alpha)}\leq M_4.
$$
Apply Theorem~\ref{thm.divergence.rate} to the solutions $\hat{x}(\cdot), x_1(\cdot)$ of the FDE \eqref{eqn.surj4}  we get  for any $t\in [0,T_1]$ the inequality
$$
|x_1(t)| \leq |\hat{x}(t)| + |x_1(0)| E_\alpha(M_1T_1^\alpha) \leq M_5,\qquad\hbox{hence}\quad f(t,x_1(t)) = \hat{f}(t,x_1(t)).
$$
Therefore $x_1(\cdot)$ is a solution of the FDE  \eqref{eqn.surj1}, and (i) is proved.

(ii) By (i), the evolution mappings of \eqref{1dimEq} are bijective, hence  $\Phi_{s,t}$ is well defined by \eqref{eqn.2flow}. The flow property is easily verified.

(iii) By (ii), the FDE \eqref{1dimEq} generates a two-parameter family of bijections $\Phi_{s,t}$ of $\mathbb{R}$ for all $s,t\in J$.
From Theorems \ref{thm.convergence.rate} and \ref{thm.divergence.rate} it follows that the bijections $\Phi_{s,t}$ are homeomorphisms and $\Phi$ depends continuously on three variables $s,t,x$.

(iv) Obvious.
\end{proof}

\begin{definition}
The two-parameter flow $\Phi_{s,t}$, specified in Theorem~\ref{thm.2param.flow}, generated by the FDE \eqref{1dimEq} is called the {\em nonlocal dynamical system generated by \eqref{1dimEq}}.
\end{definition}

\begin{remark}
Some distinguished features of the two-parameter flow generated by the FDE \eqref{1dimEq}:\\
(i) The flow has history memory. Though the past has impact on the behavior of the solutions, the solutions form a two-parameter flow of homeomorphisms.\\
(ii) The flow is in general $\alpha$-H\"older, but not $C^1$.
\end{remark}

\begin{remark}
Li and Ma~\cite[Theorem 2]{LiMa13} claimed that they constructed a dynamical system from a FDE. However their construction is false. For a counter example, see
Cong, Son and Tuan~\cite[Remark 12]{Cong_3}.
\end{remark}

\section{Triangular systems of FDEs generate nonlocal dynamical systems}\label{sec.triangular}
In this section, using the results of Section~\ref{sec.1dimNonlocalDS} we show that a higher dimensional triangular system of FDEs also generates a nonlocal dynamical system.

Let us consider a $d$-dimensional triangular system of (not necessarily linear) FDEs
\begin{equation}\label{eqn.triangular.nonlin}
\begin{cases}
^{C\!}D_{0+}^{\alpha}x_1(t)= &f_1(t,x_1(t)),\\
^{C\!}D_{0+}^{\alpha}x_2(t)= &f_2(t,x_1(t),x_2(t)),\\
\cdots&\cdots\\
^{C\!}D_{0+}^{\alpha}x_d(t)= &f_d(t,x_1(t),x_2(t),\cdots,x_d(t)),
\end{cases}
\end{equation}
where $t\in J$, $x(\cdot)=\big( x_1(\cdot),\ldots,x_d(\cdot)\big)^\mathrm {T} \in \mathbb{R}^d$, $f=(f_1,\ldots,f_d)^\mathrm {T}$ are Lipschitz functions in $x$ variables, i.e., there exists a continuous function $L: J\rightarrow [0,\infty)$ such that for all $i=1,\ldots, d$ and all $t\in J$ we have
\begin{equation}\label{eqn.Lip.tri.nonlin}
|f_i(t,x_1,\ldots,x_i)- f_i(t,y_1,\ldots,y_i)| \leq L(t) \sqrt{(x_1-y_1)^2 + \cdots (x_i-y_i)^2}.
\end{equation}
This triangular system has a distinguished property that it can be solved successively coordinate-wise and each time we have to solve only a one-dimensional FDE, hence the triangular system inherits many  features of the one-dimensional FDEs.

\begin{proposition}[Convergence rate for solutions of a triangular system of FDEs]\label{tri.nonlin.lower.bound}
Assume that the Lipschitz condition  \eqref{eqn.Lip.tri.nonlin} is satisfied.  Then for
 any two solutions $x(\cdot), y(\cdot)$ of the triangular FDE \eqref{eqn.triangular.nonlin} and any $t\in J$ the following estimate holds
 $$
\|x(t)-y(t)\| \geq \|x(0) - y(0)\| E_\alpha\big(- (\max_{0\leq \tau\leq t} L(\tau)) t^\alpha\big).
$$
\end{proposition}
\begin{proof}
Let $x(\cdot)=(x_1(\cdot),\ldots,x_d(\cdot))^\mathrm {T}$ and $y(\cdot)=(y_1(\cdot),\ldots,y_d(\cdot))^\mathrm {T}$ be arbitrary two solutions of  the triangular FDE \eqref{eqn.triangular.nonlin}.
Consider the first equation in \eqref{eqn.triangular.nonlin}: it is a 1-dimensional equation for the first coordinate. Applying Theorem~\ref{thm.convergence.rate} to this equation we get
$$
|x_1(t)-y_1(t)| \geq |x_1(0) - y_1(0)| E_\alpha\big(- (\max_{0\leq \tau\leq t} L(\tau)) t^\alpha\big).
$$
Since the first coordinate is solvable from the first equation we can substitute it into the second equation of the system \eqref{eqn.triangular.nonlin} and get a 1-dimensional FDE for the second coordinate
$$
^{C\!}D_{0+}^{\alpha} u(t)= f_2(t,x_1(t),u(t)) =: \hat{f}_2(t,u(t)),
$$
where due to \eqref{eqn.Lip.tri.nonlin} the function $\hat{f}_2(\cdot,\cdot): J\times \mathbb{R} \rightarrow \mathbb{R}$ is $L(t)$-Lipschitz continuous with respect to the second variable.
Applying Theorem~\ref{thm.convergence.rate} to the solutions $x_2(\cdot)$ and $y_2(\cdot)$ of  this 1-dimensional FDE we get
$$
|x_2(t)-y_2(t)| \geq |x_2(0) - y_2(0)| E_\alpha\big(- (\max_{0\leq \tau\leq t} L(\tau)) t^\alpha\big).
$$
Continue this process we get for any $i=1,\ldots,d$ and any $t\in J$ the inequality
\begin{equation}\label{eqn.tri.nonlin1}
|x_i(t)-y_i(t)| \geq |x_i(0) - y_i(0)| E_\alpha\big(- (\max_{0\leq \tau\leq t} L(\tau)) t^\alpha\big).
\end{equation}
This implies immediately the conclusion of the proposition.
\end{proof}

An important particular case of the triangular system of  FDEs \eqref{eqn.triangular.nonlin} is a linear triangular system of FDEs
\begin{equation}\label{triangularEq}
^{C\!}D_{0+}^{\alpha}x(t)=A(t)x(t),
\end{equation}
where $t\in J, x\in\mathbb{R}^d$,  $A:J\times \mathbb{R}^{d\times d}$ is a bounded continuous triangular $(d\times d)$-matrix function, i.e.,
$A(\cdot) = \Big( a_{ij}(\cdot)\Big)_{1\leq i,j\leq d}$
where either $a_{ij}=0$ for all $i>j$ (upper triangular) or $a_{ij}=0$ for all $i<j$ (lower triangular), and there exists a continuous function $L: J \rightarrow [0,\infty)$ such that
\begin{equation}\label{eqn.tri1}
\|A(t)\| \leq L(t)\qquad \hbox{for all}\quad t\in J.
\end{equation}
Clearly, Proposition~\ref{tri.nonlin.lower.bound} is applicable to the linear triangular system
\eqref{triangularEq}. Moreover, we also have a lower bound for solutions of \eqref{triangularEq}.

\begin{proposition}[Lower bound for solutions of a linear triangular system of FDEs]\label{tri.lower.bound}
Assume that the triangular matrix function $A$ satisfies the conditions \eqref{eqn.tri1}.  Then for
 any solution $x(\cdot)$ of the FDE \eqref{triangularEq} and for any $t\in J$ the following estimate holds
$$
\|x(t)\| \geq \|x(0)\| E_\alpha(-  (\max_{0\leq \tau\leq t} L(\tau))t^\alpha).
$$
\end{proposition}
\begin{proof}
Since the system \eqref{triangularEq} has trivial solution $0$, we can apply
Proposition~\ref{tri.nonlin.lower.bound} to the two solution $x(\cdot)$ and 0 of \eqref{triangularEq} and get the conclusion of the proposition.
\end{proof}

Similar to the one-dimensional case, for any $T_1\in J$ the mapping
\begin{equation}
\phi_{0,T_1} :\mathbb{R}^d\rightarrow \mathbb{R}^d, \quad x_0 \mapsto x(T_1),
\end{equation}
where $x_0\in\mathbb{R}^d$ is an arbitrary initial value of \eqref{eqn.triangular.nonlin}, $x(\cdot)$ is the solution of \eqref{eqn.triangular.nonlin} starting from $x(0)=x_0$ and $x(T_1)$ is the evaluation of $x(\cdot)$ at $T_1$, is called {\em the evolution mapping of \eqref{eqn.triangular.nonlin}}. By the same arguments as in Section~\ref{sec.1dimNonlocalDS} we can show that the evolution mapping $\phi_{0,t}$ of the triangular FDE \eqref{eqn.triangular.nonlin} is a bijection for any $t\in J$.

\begin{definition}\label{2parameterflow.tri.nonlin}
A two-parameter family of mappings
$$
\phi_{s,t}(\cdot) : \mathbb{R}^d\rightarrow \mathbb{R}^d, \qquad s,t\in J,
$$
 is called  {\em a two-parameter flow in $\mathbb{R}^d$} if the following three conditions are satisfied\\
(i) $\phi_{s,t}(x)$ is continuous as a function of three variables $s,t\in J, x\in\mathbb{R}^d$;\\
 (ii) For any fixed $s,t\in J$ the mapping $\varphi_{s,t}(\cdot)$ is a homeomorphism of $\mathbb{R}^d$;\\
(iii) This family satisfies the following flow property
$$
\phi_{s,t}\circ \phi_{u,s} = \phi_{u,t}\qquad\hbox{for all}\quad u,s,t\in J.
$$
\end{definition}

\begin{theorem}[Triangular systems of FDEs generate nonlocal dynamical systems]\label{thm.2param.flow.triangular}
 (i) The triangular FDE \eqref{eqn.triangular.nonlin} generates a two-parameter flow  in $\mathbb{R}^d$; namely it generates the two-parameter flow in $\mathbb{R}^d$
 $$
\phi_{s,t} : \mathbb{R}^d\rightarrow \mathbb{R}^d, \qquad s,t\in J,
$$
where $\phi_{0,t}$, $t\in J$, is the evolution mapping of \eqref{eqn.triangular.nonlin} and
$$
\phi_{s,t} := \phi_{0,t}\circ \phi_{0,s}^{-1}\qquad\hbox{for all}\quad s,t\in J.
$$
(ii) The linear triangular system of FDEs \eqref{triangularEq} generates a two-parameter flow of $d$-dimensional linear nonsingular operators.
\end{theorem}
\begin{proof}
Similar to Theorem~\ref{thm.2param.flow}.
\end{proof}

\begin{definition}
The two-parameter flow $\phi_{s,t}$ generated by the triangular FDE \eqref{eqn.triangular.nonlin} is called the {\em nonlocal dynamical system generated by \eqref{eqn.triangular.nonlin}}.
\end{definition}

\section{A general high dimensional system of FDEs does not generate a dynamical system}\label{sec.d-dim}

In this section we show that, in the high dimensional case different trajectories of a FDE can intersect each other. Thus a high dimensional system of FDEs does, in general, not generate a dynamical system. In order for a high dimensional FDE to generate a nonlocal dynamical system we need an additional property of the FDE, for example triangularity.
The results of this sections also provide a counter example to the assertions of Bonilla, Rivero and Trujillo~\cite[Theorem 1, Propositions 1 and 2]{Bonilla2007}.

\begin{theorem}[Different trajectories of a high dimensional system of FDEs can meet]\label{thm.intersection}
For any $d\geq 2$ there exists a system of type \eqref{mainEq} with a property that it has two different solutions $x_1(\cdot), x_2(\cdot)$ with $x_1(0)\not= x_2(0)$ which intersect each other at some finite time moment $0<T<\infty$, i.e., $x_1(T)=x_2(T)$.
\end{theorem}
\begin{proof}
It suffices to construct a two-dimensional system of type \eqref{mainEq} with the desired property.
Below we shall construct a two-dimensional linear autonomous system of FDEs with such a property.

Since $\alpha\in (0,1)$, the complex-valued Mittag-Leffler function $E_\alpha(\cdot)$ has infinitely many zeros in $\mathbb{C}$ (see Gorenflo \textit{et al.}~\cite[Corollary~3.10, p.~30]{GoKiMaRo14}). Take and fix $z^*\in \mathbb{C}$ such that $E_\alpha(z^*)=0$. Let $\phi:= \arg(z^*)\in (-\pi, \pi]$, put $\lambda := \cos\phi + i\sin\phi$, where $i = \sqrt{-1}\in\mathbb{C}$.
Note that since $\alpha\in \mathbb{R}$ we have  $E_\alpha(\overline{z^*}) = \overline{E_\alpha(z^*)} =0$, where $\overline{w}$ denotes the complex conjugate of the complex number $w$. Since $\alpha\in (0,1)$ we have $z^*\notin\mathbb{R}$, hence $\lambda\notin \mathbb{R}$. Let
$$
A:=
\left(
\begin{array}{cc}
 \cos\phi    &     \sin\phi \\
  -\sin\phi    & \cos\phi \\
   \end{array}
   \right).
$$
Then $A$ has two (complex) eigenvalues $\lambda, \overline\lambda$. We show that the linear FDE
\begin{equation}\label{eqn.2-dim}
    ^{C\!}D_{0+}^\alpha x(t)= Ax(t), \quad t\in\mathbb{R}_+,
\end{equation}
is the desired FDE. Indeed, it is known that this linear FDE is solvable (see Diethelm~\cite[Theorem~7.15, p.~152]{Diethelm2010}) and has two following independent real solutions
$$
 x_1(t):=  \left(
\begin{array}{c}
 \big(E_\alpha(\lambda t^\alpha) + E_\alpha(\overline\lambda t^\alpha)\big)
   \\ [6pt]
i\big(E_\alpha(\lambda t^\alpha) - E_\alpha(\overline\lambda t^\alpha)\big)  \\
   \end{array}
   \right),\;
  x_2(t):=   \left(
\begin{array}{c}
 -i\big(E_\alpha(\lambda t^\alpha) - E_\alpha(\overline\lambda t^\alpha)\big)  \\ [6pt]
 \big(E_\alpha(\lambda t^\alpha) + E_\alpha(\overline\lambda t^\alpha)\big)       \\
   \end{array}
   \right).
   $$
Put
$$
u(t) := E_\alpha(\lambda t^\alpha) + E_\alpha(\overline\lambda t^\alpha)
  \qquad
   \hbox{and}\quad
v(t) := i \big(E_\alpha(\lambda t^\alpha) - E_\alpha(\overline\lambda t^\alpha)\big).
$$
Then $u, v : \mathbb{R}_+ \rightarrow \mathbb{R}$, $u(0)=2, v(0)=0$, and
$$
 x_1(t) =  \left(
\begin{array}{c}
 u(t)
   \\
v(t) \\
   \end{array}
   \right),\;\;
  x_2(t) =   \left(
\begin{array}{c}
 -v(t) \\
 u(t)       \\
   \end{array}
   \right),
\quad
   x_1(0) =  \left(
\begin{array}{c}
 2
   \\
0\\
   \end{array}
   \right),\;\;
  x_2(0) =   \left(
\begin{array}{c}
 0\\
 2     \\
   \end{array}
   \right)
   $$
 The general solution of
\eqref{eqn.2-dim} is
\begin{equation}\label{eqn.2dim2}
x(t)= a x_1(t)  + b x_2(t)
\end{equation}
where $a,b\in\mathbb{R}$ are arbitrary real constants. Let $T>0$ be the unique finite positive number satisfying
$$
\lambda T^\alpha = z^*;
$$
such a $T$ exists uniquely due to the definition of $z^*$ and $\lambda$.
Clearly $u(T)=v(T)=0$, hence $x_1(T)=x_2(T)=(0,0)^\mathrm {T}$. From \eqref{eqn.2dim2} it follows that for any solution $x(t)$ of \eqref{eqn.2-dim} we have $x(T) =0$.
\end{proof}

From Theorem~\ref{thm.intersection} we obtain immediately the following corollary.
\begin{corollary}[A higher dimensional FDE does not generate a nonlocal dynamical system]
For $d\geq 2$ the FDE \eqref{mainEq} does, in general, not generate a two-parameter flow in $\mathbb{R}^d$. Hence it does, in general, not generate a nonlocal dynamical system.
\end{corollary}

\begin{remark}
(i) Actually for the FDE \eqref{eqn.2-dim} {\em all the solutions} are equal $(0,0)^\mathrm {T}$ at time $T$; hence all they meet each other.\\
(ii) Theorem~\ref{thm.intersection} shows that in contrary to the initial value problem, the terminal value problem for FDEs is not always solvable.\\
(iii) Theorem~\ref{thm.intersection} allows us to better understand the dynamics of the FDEs. It reveals a distinguished feature of FDEs in comparision with ODEs: different trajectories of an ODE cannot meet whereas for FDEs they may meet.\\
(iv) By small modification of $A$ in the proof of Theorem~\ref{thm.intersection} we can make the time of intersection $T$ small.\\
(v) By a small modification of the proof one can show that Theorem~\ref{thm.intersection} also holds for any positive real $\alpha \not=1$.
\end{remark}

\section*{Acknowledgement}
This research of the authors is funded by the Vietnam National Foundation for
Science and Technology Development (NAFOSTED) under Grant Number 101.03-2014.42.

redpajama_set_name: RedPajamaArXiv


text:
\section{INTRODUCTION}
\label{sec:intro}
The sky background plays a significant role in astronomy as it affects all forms of astronomical observations across the entire electromagnetic spectrum. The sky background in space and on the ground from the far ultraviolet to the far infrared wavelengths is well studied, and it is now known that a multitude of emission processes dominate different wavelength regimes\cite{leinert98}. The characterization of the sky background is crucial in astronomy because the sky brightness often determines the length of the time required to reach a given signal-to-noise for an astronomical observation. This is especially the case in the near-infrared (NIR) where the sky background is significantly brighter than in the visible. The NIR bands are mostly dominated by non-thermal emission, namely very narrow airglow lines that have been identified to be hydroxyl radical (OH) lines, Meinel bands \cite{meinel50,meinel50b}, and to a lesser extent O$_2$ ($^1 \Delta _g$ band) emission at 1.27 $\mu$m\cite{maihara93}. The airglow lines are generated by the vibrational de-excitation of excited OH molecules created in the mesosphere at an altitude of $80-105$km\cite{ramsay92}. The main reaction that produces this molecule is the following:
\begin{equation}
\textrm{H}+\textrm{O}_3\rightarrow\textrm{OH}^*+\textrm{O}_2,
\end{equation}
where OH is formed in a vibrationally excited ($\nu \le 9$) state\cite{bates50}. The excited OH radical cascades radiatively to lower energy states and produces the characteristic airglow spectrum observed in the NIR (see Figure \ref{fig:nearirsky}). This emission is highly variable in short (minutes) and longer (hours to years) timescales. The short time scale variations of tens of minutes have been identified to be caused by the passage of gravity (buoyancy) waves high in the atmosphere that modify the densities of reactants that produce the OH molecule\cite{ramsay92,swenson94}. Longer time scale diurnal variations along with correlations with the solar cycle have been observed. There have been studies that show that the OH airglow is brighter during periods of intense solar activity\cite{pertsev08}.
   \begin{figure}
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[height=8cm]{nearirsky.eps}
   \end{tabular}
   \end{center}
   \caption[nearirsky]
   { \label{fig:nearirsky}  The NIR atmospheric transmission for an excellent mid-latitude site, shown here for Mauna Kea (black curve) and over plotted with a corresponding high-resolution spectrum of the sky background (red). The purple, green, and orange curves are the Mauna Kea Observatories (MKO) $J$, $H$, and $K^\prime$ filter bandpasses. Most of the sky flux in this wavelength regime is completely dominated by OH airglow lines. Only at the very end of the NIR band does the thermal emission from the sky start to dominate. (Colour version available online)
}
   \end{figure}

Significant work has been done to find ways to suppress this OH background and improve NIR sensitivity. These techniques require complex hardware to accomplish this task\cite{herbst94,hawthorn11}. The goal of these methods is to reach the IR continuum in between the OH lines, which is significantly fainter and dominated by the zodiacal background though there may be other emission mechanisms at work\cite{maihara93}. Another alternative to reduce the OH background component is to find a site that has intrinsically lower OH background. A winter campaign carried out at the South Pole to measure the IR sky background from the near to mid-IR revealed that the sky background was significantly fainter than a mid-latitude site\cite{phillips99}. In the middle of winter at the South Pole, there should be continuous darkness resulting in little to no diurnal variations in the sky brightness. It is not particularly surprising that the thermal sky emission was measured to be an order of magnitude or more fainter than a comparable mid-latitude site because of the frigid temperatures ($\sim -60$C) and very low precipitable water vapour at the South Pole. But the $J$ and $H$-band sky background values were also fainter by factors of 2 to 3 ($\sim1$ mag fainter). This is particularly intriguing because it suggests the chemical process that produces the Meinel bands in the South Pole region is somehow less effective.
\par
Recent investigation of a potential Canadian astronomical site near Eureka, Nunavut at 80$^\circ$ N\cite{steinbring12}, and others on Ellesmere Island at higher latitude\cite{steinbring10}, prompted us to test the hypothesis that Arctic skies may also experience suppressed OH airglow in winter. There are a number of virtues of Eureka. Being located very close to the North magnetic pole\footnote[1]{The magnetic pole is not stationary and moves a few 10s of km per year.}, Eureka is within the auroral hole. It also remains accessible throughout the winter through chartered flights. Finally, the environment is not as harsh as the South Pole as average winter temperatures hover around $-40$C, the minimum rated temperature for most electronic equipment. Our primary objective was to carry out a long term monitoring of the NIR sky brightness of the site. We chose the Polar Environment Atmospheric Research Laboratory (PEARL), at an elevation of 610-meters, located approximately 15 km from the Eureka station, as the location of our site testing campaign due to both the availability of electrical power and internet access. We focussed primarily on the $J$-band as that was where the greatest improvement in the sky background was observed in the $1-2 \mu$m wavelength range in Antarctica. For this purpose, we constructed a ruggedized $J$-band photometer that monitored the sky brightness continuously. There was a short site testing campaign that took place in February 2012 when the photometer was set up for 11 days. During this time, we obtained four nights of good sky brightness data. In this paper, we present the results of the first sky brightness measurements in the High Arctic, along with the design, construction, and calibration of the photometer. We conclude with the issues we experienced and our future plans for NIR sky brightness characterization.

   \begin{figure}
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[height=6cm]{photometer.pdf}
   \end{tabular}
   \end{center}
   \caption[photometer]
   { \label{fig:photometer} A photograph of the $J$-band photometer. The black lens tube contains all of the optical components of the photometer. The electronics and the photodiode are housed in a sealed, protective box. Both power and communication are via the cable at left. There is a MIL-SPEC connector for the cable on one end of the photometer that provides a good environmental seal.
}
   \end{figure}

\section{Instrument Design}
\label{sec:design}
The goal of the design was to construct a simple photometer that can withstand temperatures reaching $\sim -40$C while being able to measure the zenith $J$-band sky brightness at a high signal-to-noise ratio (SNR). Eureka is a somewhat more forgiving environment than the South Pole because most common electronic components are rated to operate at temperatures down to $-40$C. Nevertheless, one needs to be mindful of using movable components. We already have significant experience constructing precise current measuring circuits that are central to commercial single pixel, visible-light sky brightness monitors constructed by our co-authors. These devices have been shown to meet the necessary operating temperature requirement. We used this experience to construct an infrared sky brightness photometer sensitive enough to accurately monitor the $J$-band sky with minimal moving parts. In the interest of cost, reliability, and ease of construction, we chose a zenith-pointing, single pixel photometer design with no moving parts.
\par
There were a number of design considerations that came into play in the construction of the photometer. The $J$-band passband was chosen for a number of reasons: (1) The $J$-band sky background is entirely dominated by OH airglow; (2) The greatest reduction in OH airglow-dominated sky brightness in Antarctica was seen in $J$-band\cite{phillips99}; and (3) The typical InGaAs diode sensitivity cuts off at $1.7\mu$m, which means part of $H$ band will not be accessible. Extended wavelength diodes can be purchased, but they are significantly more expensive, and have lower spectral detectivity, D$^*,$ values. The next design consideration is the accuracy to which we could measure the photocurrent generated by sky photons. From previous experience we knew that our current sensing circuits could measure currents with reasonable accuracy of about 0.25 pico amp (pA). We set a design requirement of sky generated photocurrent of $\sim5$ pA. This requirement drives the field-of-view (FOV) and aperture size of the photometer. Because we intend to measure the flux from an extended source, we need to maximize the photometer's etendue, i.e. $A\Omega$ where $A$ is the collecting area, and $\Omega$ is the solid angle of our photometer.
\par
For our first design we considered a simple radiometer, which consisted of a baffle tube with a photodiode at one end. In this case, $A$ is set by the size of the photodiode, which was 2mm in diameter, and $\Omega$ is defined by the baffle tube, which was approximately 15 degrees. We constructed a prototype of this design and found that one needed a larger FOV to reach a sufficient signal-to-noise ratio. A larger FOV would mean greater sensitivity to variations in sky brightness and atmospheric extinction across the sky. The bigger issue was that the angular response of the radiometer varied strongly with of field position and required fairly accurate lab characterization to determine the solid angle of the device. Our final revised design made use of a lens, which improved the sensitivity over the original baffle-tube prototype by a factor of 10. The lens radiometer provides a well-defined field-of-view with a uniform angular response across its FOV. In order to maximize the etendue, the lens needs to be fairly fast, i.e. low f/\#. This drove us to an aspheric lens design that used an InGaAs photodiode with a relatively large active area. For our lens, we used an off-the-shelf Edmund Optics precision asphere lens (part \#49115), which is 1-inch in diameter and has a focal length of $\sim30$mm. The lens is antireflection-coated to be less than 1.5\% reflective in the $0.6-1.05$ $\mu$m range. The diode was also an off-the-shelf part from Hamamatsu Corp (part \#G8605-22). The diode has a two stage thermoelectric cooler and an active area 2 mm in diameter. The $J$-band filter was obtained from Optec, which makes these filters for IR photometers. The $J$-band filter closely resembles the Mauna Kea Obsevatories (MKO) $J$-band filter\cite{tokunaga02}. A 30-mm uncoated Sapphire entrance window was used to protect the optical components from the elements. The optical layout of the system is shown in Figure \ref{fig:layout}. The optical parameters of the photometer are given in Table \ref{tab:parameters}.

  \begin{figure}
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[height=5cm]{layout.pdf}
   \end{tabular}
   \end{center}
   \caption[layout]
   { \label{fig:layout}
   The optical layout of the $J$-band photometer. Light enters from the left and the different coloured rays show the rays from different field angles. The red and green rays define the edge of the photometer's FOV. The flat sapphire window is shown to the left, followed by the aspheric lens, the $J$-band filter, and the focal plane. The stop of the system is located at the rear of the aspheric lens, which is defined by the inner diameter of the lens tube. The lens is mounted on a focusing stage, which is used to set the proper focus during assembly. The distance from the sapphire window to the lens is 35 mm. (Colour version available online)
   }
   \end{figure}

 \begin{table}[h]
\caption{System parameters for $J$-band photometer}
\label{tab:parameters}
\begin{center}
\small
\begin{tabular}{|l|l|}
\hline
\rule[-1ex]{0pt}{3.5ex}  {\bf Aperture Diameter} & 20.8 mm \\
\hline
\rule[-1ex]{0pt}{3.5ex}  {\bf f/\#} & 1.5\\
\hline
\rule[-1ex]{0pt}{3.5ex}  {\bf Field-of-view } &  3.5 degree diameter\\
\hline
\rule[-1ex]{0pt}{3.5ex}  {\bf Sky Area (Solid Angle)} & $1.32\times10^8$ arcsec$^2$   \\
\hline
\rule[-1ex]{0pt}{3.5ex}  {\bf Filter Center Wavelength} & 1.260 $\mu$m  \\
\hline
\rule[-1ex]{0pt}{3.5ex}  {\bf Filter Bandwidth (Full-Width-Half-Max)} & 0.1292 $\mu$m   \\
\hline
\rule[-1ex]{0pt}{3.5ex}  {\bf Filter Throughput} & $\sim 0.8$ \\
\hline
\rule[-1ex]{0pt}{3.5ex}  {\bf End-to-end Throughput (predicted/measured)} & $\sim 0.67/\sim 0.7$ \\
\hline
\rule[-1ex]{0pt}{3.5ex}  {\bf Average Responsivity} & $\sim0.84$ A/W   \\
\hline
\rule[-1ex]{0pt}{3.5ex}  {\bf Predicted Sky Photocurrent} & $\sim 4$ pA @ 16.0 mag arcsec$^{-2}$ \\
\hline
\end{tabular}
\end{center}
\end{table}

\par
The goal of the optical design was to obtain relatively uniform angular response, while maximizing throughput, over the full FOV of the photometer. This requires fairly small root-mean-squared (RMS) spot sizes that are much smaller than the size of the detector. This ensures all of the focussed light falls on the detector over the majority of the photometer's FOV. Conventional fast spherical lenses exhibit significant spherical aberration and are unable to satisfy this requirement. The aspheric design works remarkably well, and its spot diagram over the full field is shown in Figure \ref{fig:spotdiagram}. Typical RMS spot sizes are within $25-30$ $\mu$m across the entire field over the wavelengths of interest. We simulated the angular response function for and solid angle subtended by the photometer using ZEMAX. The results are shown in Figure \ref{fig:fov}. For the ideal case when the photometer is properly focussed at infinity, the angular response function (ARF) has a fairly sharp transition at 1.75 degrees as expected for our design. This is because the RMS spot sizes are much smaller than the detector size. We also show the effects of an improperly set focus and how that affects the ARF in the left panel of Figure \ref{fig:fov}. Small changes of $\pm 0.1$mm, which is much coarser than the our ability to set the lens separation, has a minor effect on the ARF. The solid angle subtended by the photometer for the slightly defocussed case does not change appreciably from the perfectly focussed case, see Figure \ref{fig:fov} right panel. However, we were unable to focus the photometer at infinity. Instead, the focussing was carried out by trying to maximize the observed signal of a fairly collimated ($\pm 5$ degrees), IR ($\lambda = 950$ nm) LED light source located 3-meters away. This means that the focus is not properly set to infinity, and the effect on the ARF is shown in Figure \ref{fig:fov} as the red curve. However, the solid angle of this as-built setup only differs by 2\% from the nominal infinity focus case, which means the focal shift does not significantly impact our knowledge of the photometer's solid angle (given in Table \ref{tab:parameters}). From this point on, we will use this solid angle to be the nominal solid angle of the photometer.
\par
We estimated the expected photocurrent for our $J$-band sky brightness measurements. Assuming the optimistic case of a sky brightness of 16.0 mag arcsec$^{-2},$ which was on the upper end of the typical sky brightness values measured in the Antarctic \cite{phillips99}, we estimate a photocurrent of 4 pA. This assumes a total photometer throughput of $\sim 0.67$ where we have a sapphire window transmission of 0.86, a lens transmission of 0.97, and filter transmission of 0.8. The diode responsivity is assumed to be 0.84 A/W, which was calculated by taking a weighted average value of a typical Hamamatsu InGaAs photodiode responsivity curve across the $J$-band bandpass. For the most optimistic case of a sky brightness of 16.8 mag arcsec$^{-2},$ the lower end of the Antarctic measurements, we would measure a photocurrent of 2 pA. Our photometer should detect this signal at the 5$\sigma$ level.

   \begin{figure}
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[height=8cm]{spotdiagram.pdf}
   \end{tabular}
   \end{center}
   \caption[spotdiagram]
   { \label{fig:spotdiagram} The spot diagram for the aspheric photometer design over multiple fields. The blue, green, and red colours represent the spots generated by rays at wavelengths of 1.1, 1.2, and 1.3 $\mu$m, respectively. These wavelengths span the $J$-band bandpass. The spot diagram shows that the performance is fairly uniform across the full field at $J$-band wavelengths. RMS spot radii are typically in the $25-30$ $\mu$m range. (Colour version available online) }
   \end{figure}

    \begin{figure}
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[height=6cm]{fov_v2.eps} \includegraphics[height=6cm]{solidangle.eps}
   \end{tabular}
   \end{center}
   \caption[fov]
   { \label{fig:fov} {\it Left:} Simulated angular response function for the photometer for different focus settings of the photometer. The throughput at angles smaller than 1 degree is 1.0. This region has not been plotted to accentuate the changes in the response function at larger angles. The nominal case is shown by the black curve where the lens is properly set to infinity focus. The blue and green curves show the effect of a slight focal shift on the angular response function. The red curve shows the as-built simulated angular response function. {\it Right:} Solid angle of photometer as a function of focal shift of lens. 0 mm is defined to be nominal infinity focus case. The as-built focal position of the photometer is expected to be +0.2 mm from nominal. In all cases, the solid angle of the photometer does not change more than a few percent from the infinity focus case. (Colour version available online) }
   \end{figure}

 \par
The final component of the photometer is the photocurrent measurement device. The photocurrent is measured by a custom-built electronic circuit that accomplishes the following tasks: (1) Measures and digitizes the photocurrent; (2) Allows one to change the gain and integration time for these measurements; and (3) Communicates with a computer over a serial line. To communicate with the photometer, we use a netbook computer that executes a PERL script to setup the photometer's operating parameters and retrieve data. A protocol that uses a suite of text commands is used to communicate with the photometer. The netbook then creates a text file of the measurements that can be retrieved at later time via an Internet connection.
\par
The heart of the custom circuit consists of a current input 20-bit analog-to-digital converter (ADC) - Texas Instruments DDC112 - and a PIC microprocessor, which provides a bidirectional link to a computer over a serial line and sets up the DDC112 to carry out the current measurements. The photodiode is configured in zero bias mode, which means there is no dark current. We do not require high speed measurements, which makes the lack of biasing a worthwhile compromise because at zero bias the diode response to light is very linear over a wide intensity range. The DDC112 is able to make unipolar current measurements to a high degree of accuracy over a wide temperature range. Schematically, the DDC 112 is an integrator with a known capacitance (5\% tolerance) and its output is fed into a built-in ADC.   There are several integration capacitors that can be chosen depending on the amount of charge that needs to be measured for a given integration time, which essentially changes the gain of the measurement. After one integration time, the digitized photocurrent value is retrieved by the microprocessor and can be read by the control computer. In our case, due to the exceedingly small photocurrent to be measured, we use the smallest integrating capacitor (12.5 pF) and the longest possible integration time of one second. The full well charge level at this setting is 50pC, which means we can measure a current up to 50 pA at 1s integration time without saturating the integrator. At this setting, we estimate a 1$\sigma$ current measurement of 0.25 pA. This error estimate is obtained by characterizing the stability of measurements and shifts in the bias values as a function of temperature. This means that the predicted photocurrent of $\sim 4$ pA is detected at a high level ($> 10\sigma$) of significance.
\par
Another function that we implemented in the circuit is the temperature control of the photodiode. The photodiode in the photometer has a thermoelectric cooler (TEC) that can be used to cool the diode 40C below ambient. This feature improves the stability of the photodiode and reduces the leakage current, see Section \ref{sec:issues} for further discussion. The circuit also consists of a temperature sensing circuit that determines the diode temperature by measuring the resistance of its internal thermistor. The TEC driver is a constant current source that is turned on and off by a thermostat control implemented in the photometer firmware. This implementation is able to maintain the temperature of the diode to within 1C.

\section{Calibration and Testing}
\label{sec:calibration}
In order to know the conversion between counts measured by the photometer and the true sky brightness, one needs to know the flux calibration of and the solid angle subtended by the photometer. The ZEMAX simulations show that the solid angle of the photometer is insensitive to slight errors in the focus adjustment. We did carry out crude measurements of the angular response that were not at the level necessary to confirm the solid angle of the photometer to high accuracy, but we found our measurements to be consistent with simulations. It is possible that there might be additional sensitivity at larger angles due to scattering inside the lens barrel. This will only increase the solid angle of the photometer, and therefore reduce the total sky brightness reported by the photometer. Given that our goal was to determine if the High Arctic was as dark as the South Pole, our measurements would then set an upper limit on the measured sky brightness. We estimate a 10\% uncertainty in our solid angle determination. Our calibration procedure focused entirely on the flux calibration of the photometer.
\par
The flux calibration of the photometer was done in a lab environment using a calibrated light source. It is not possible to carry out a flux calibration by observing an astrophysical source since even the brightest stars are too faint, and the Sun is too bright. The solution then was to use a calibrated lab source to carry out the flux calibration. The lab characterization setup used a 1310 nm pigtailed laser with a collimator at its fibre end. The 1310 nm laser is ideal because it is within the $J$-band bandpass. The output of the collimator produced a Gaussian beam with a measured beam waist of 1.5 mm and a divergence angle of $\sim0.3$ degrees. The laser was set to operate in constant power mode, and an InGaAs photodiode 2mm in diameter with a calibrated response curve that was NIST traceable was used to measure the laser's power output by placing the photodiode very close to the collimator output. An ND filter was inserted at the output of the laser collimator to attenuate the light to a level that did not cause the photodiode readout to saturate at the lowest gain and shortest integration time. To determine what fraction of the light from the laser fell within the diode active area, we scanned the diode across the laser beam and recorded its readout as a function of scan distance. The best fit to the scan curve indicated that $95-97$\% of the laser light fell within the diode active area. To test the photometer, we kept the ND filter in the beam and set up the photometer such that the optical axis of the collimated laser beam and the photometer optics were co-aligned.  Coalignment was thought to have been achieved when the highest readout was obtained from the photometer while the beam was centred on its lens. Note that the width of the laser beam was only 5 mm at the photometer entrance aperture and we did not fill the full clear aperture. Nevertheless, all of the light incident on the photometer should be focussed on the photodiode, and one should be able to obtain a flux calibration of the system. This calibration obtained a throughput of $\sim 0.7,$ which assumes the responsivity of the calibration diode and the photometer diode are the same at 1.31 $\mu$m. This measurement should be to good to within a few percent. This is also roughly consistent with our predicted throughput for the system. Based on the stability of our measurements, we expect this flux calibration to be good to 10\%. The weakness of this calibration technique is that it does not test the photometer at the same flux levels and spectrum as the sky background, and the same environmental conditions as in the Arctic. Furthermore, the calibrated light source does not fill the full aperture of the photometer, which means spatial variations in transmission of the optical system and the responsivity of the diode can affect our calibration. More careful calibrations need to be done in the future to verify the accuracy of this method. The total uncertainty in the sky brightness measurement, dominated entirely by systematic errors in the photometer's flux calibration and solid angle determination, is therefore $\sim 15\%$ or 0.15 mag.
\par
The system was tested in the dark to determine its bias value. We discovered that its bias shifted significantly as a function of temperature. This was already discovered when the first prototype was constructed, which did not have a thermoelectric cooler. The photometer's onboard temperature measurement was used to correlate the changes in the bias value. While a general trend was discovered, it was difficult to accurately predict the expected bias value for a given temperature. Further testing discovered that two main components affected the bias value of the diode: (1) A leakage current flowing in the opposite direction of the photocurrent that was highly sensitive to temperature; (2) Thermal flux detected by the photodiode from emission from the optical and mechanical components near it; and (3) The operation of the thermoelectric cooler. The first issue can be mitigated by cooling the diode, but the second two are intrinsic to the photometer design. We discuss these issues further in Section \ref{sec:issues}. Fortunately, none of these issues are of concern at the very low temperatures the photometer operates at in the High Arctic where active cooling is not necessary. However, they are of considerable concern for operation in warmer temperatures if, for example, the photometer used to test a warmer mid-latitude site.

\section{High Arctic Sky Background Measurement}
\label{sec:irskyback}
We measured the $J$-band sky brightness at PEARL over a period of 11 days during which we continuously operated the photometer. During this period, we obtained four good days of data. Some days were lost due to poor weather depositing snow on the entrance window while others were clouded out. We also installed an off-the-shelf Unihedron visible light sky quality meter (Model: SQM-LU), which measured the $V$-band sky brightness at the site and served as an additional source of information. The $V$-band sky monitor had a larger field that was 20 degrees (half-power beam width). The IR measurements were taken in 10 second intervals, while the $V$-band measurements were taken in one minute intervals. We show our results in Figure \ref{fig:skyback}. During this period of four days, both the Moon and Sun remained below the horizon producing the darkest skies when the Sun was below $-18^\circ$ elevation. The $J$-band sky brightness varied from 15.5 to 15.9 mag arcsec$^{-2}$ during this time. On the best night, which stretched from the 17th to the 18th, the $J$-band sky brightness remained consistently at 15.8 mag arcsec$^{-2}.$ This is not corrected for atmospheric extinction, which is typically $< 0.1$ mag in $J$-band at a good astronomical site\cite{lombardi11}. Also, the $V$-band sky brightness values ranged within 21.0 to 21.3 mag arcsec$^{-2}$ during this period. If we account for a typical $V$-band atmospheric extinction of 0.2 mag, these values agree with previous measurements\cite{steinbring12}. Both of these measurements suggest that the site is of excellent quality in terms of sky brightness. The $J$-band sky brightness is consistent with typical values at the best mid-latitude astronomical sites, but is still brighter than the Antarctic measurements. A compilation of NIR sky-brightness values for mid-latitude astronomical sites is available elsewhere\cite{sanchez08}.

  \begin{figure}
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[height=12cm]{irsky_v4.eps}
   \end{tabular}
   \end{center}
   \caption[skyback]
  {The measured sky background over a four day period at the PEARL station near Eureka, NU. Time is shown in Eastern Standard Time (EST). The black curve represents the $J$-band sky brightness, while the red curve represents the $V$-band sky brightness. The gaps in the data streams are regions with unreliable data, i.e. snow covered windows or spikes in brightness due to personnel moving near the photometers. The orange and light blue curves represent the Sun and Moon elevations during that time. Neither the Sun nor the Moon were above the horizon during the period, but the Sun did spend a significant fraction of the time above $-12^\circ$ twilight. The green hatched region represents typical Mauna Kea values while the blue hatched region cover the typical range of values observed at the South Pole\cite{phillips99}. (Colour version available online) \label{fig:skyback} }
   \end{figure}

\par
There could be a number of reasons why our measurements are higher than at the South Pole; one being that the OH airglow is not suppressed in the High Arctic. More data are needed to be convinced of that though, as our statistics are not yet robust. The Antarctic measurements were carried out continuously over a course of four winter months; however, we only obtained approximately four days of good data. While the majority of the Antarctic $J$-band sky brightness measurements fall within the $16.0-16.8$ mag arcsec$^{-2}$ range, there is a long tail on the bright end that stretches out to 15.5 mag arcsec$^{-2}$\cite{phillips99}. Our current sky brightness measurements do not have enough statistical significance to suggest that we are comparable or worse than the Antarctic measurements. A longer term monitoring program is required. Another explanation may be one related to solar activity. It has been shown that the OH airglow brightness is correlated with solar activity measured in the form of 10.7-centimeter (cm) solar microwave flux\cite{pertsev08}. It was found that the OH airglow varied by $30-40$\% per 100 solar flux units (sfu) measured at 10.7cm\cite{pertsev08}. The original South Pole measurements were taken close to the solar minimum during 1995. However, our measurements were taken in 2012 as the Sun was approaching its maximum period of activity. We obtained the sfu measurements for both periods from Space Weather Canada database\footnote[2]{http:$//$spaceweather.gc.ca$/$data-donnee$/$sol\_flux$/$sx-eng.php}. During the period of Antarctic measurements, the solar microwave flux varied from $65-70$ sfu's, while during our measurement period the solar flux varied from $91-93$ sfu's. If we were to correct our sky brightness measurements to match the solar activity in 1995, we would only see a reduction of 0.1 mag. This cannot be the sole explanation. Another possibility mentioned in Section \ref{sec:calibration} is the accuracy of the solid angle measurement for the photometer. It is possible that we underestimate the solid angle of the photometer, which will result in an overestimate of the sky brightness. Careful characterization in the lab can resolve this issue. It is also worth noting that the Antarctic measurements were subject to a $30-40$\% uncertainty in the J and $H$-bands, which could lead to an underestimate of the sky brightness by up to $0.3-0.4$ mags\cite{phillips99}. Although this could bring the South Pole measurements into closer agreement with ours, we still have not seen sky as dark as from the Antarctic glacial plateau.

\section{Issues Encountered and Lessons Learnt}
\label{sec:issues}
A number of issues were identified in the design of the photometer, which were mainly due to the non-ideality of the components. To make accurate radiometric measurements, the dark value of the photometer needs to be known to a high degree of accuracy. It was discovered that this dark value changed significantly with temperature. One puzzling aspect was that the photodiode should not produce any dark current because it was run in zero bias mode. Whatever bias was in the system was inherent to the readout circuit, which was designed to have a well-defined fat zero value. For sufficiently long integration times of 1s, the bias value would become negative and be clipped because it had fallen below the low end of the ADC's range. We discovered that this was because there was a current flowing through the photodiode that was in the opposite direction of the photocurrent. This current was also very sensitive to the diode temperature. It is worth noting that photodiodes do not have infinite resistances and they have a shunt resistance. The shunt resistance is inversely proportional to the diode active area and material properties. In our case, we were using a fairly large InGaAs diode, which meant its shunt resistance was only 20 M$\Omega$ at 20C. This value increases significantly to 200 M$\Omega$ when the diode temperature is $-10$C. There is an exponential relation between shunt resistance and temperature. Of course, a current would not flow through a resistor unless there was potential difference. It turns out the diode was not exactly zero-biased due to a small input offset voltage introduced by the DDC 112. This voltage is at the level of a few 10s of $\mu$V. This is enough to generate a leakage current of $\sim0.5$ pA at room temperature. However, this diminishes significantly at temperatures lower than $-10$C where it has a negligible effect.
\par
To mitigate this issue, we used a TEC-cooled diode in order to maintain its temperature at $<-10$C at all times. Unfortunately, when we turned on the TEC we observed a rather large dark current. We managed to isolate this problem to two separate issues: (1) Leakage of the TEC driver current into the photocurrent sensing circuit; and (2) Thermal light sensitivity of the photodiode. The TEC requires a current of a few hundred mAs to adequately cool the diode. We attempted to isolate the TEC from the photo sensing circuit by driving it using a battery supply. The problem remained. This meant that there were stray impedances, most likely between the terminal posts of the photodiode, that were coupling a small portion of the TEC current into the photo sensing circuit. The second issue was originally an oversight. A typical InGaAs diode's wavelength sensitivity cuts off at 1.7$\mu$m. Since the $J$-band filter blocks all light long ward of $\sim1.3$ $\mu$m, the filter is close to 100\% emissive at wavelengths outside of its bandpass. The diode then sees blackbody emission across $2\pi$ sr at room temperature outside of the filter bandpass, and the light long ward of $1.3$ $\mu$m makes a significant contribution to the signal. Our calculations show that at an ambient temperature of 25C, we expect to see approximately 1 pA from thermal light. This was very close to the measured value in our tests. At $-10$ and $-40$C, one would expect a photocurrent of 0.02 pA and 0.001 pA from thermal photons, respectively. Therefore, the leakage current and thermal photons make negligible contributions to a photometer operating at Arctic temperatures. However, for a general purpose photometer that operates under all environmental conditions, these effects need to be carefully characterized or mitigated.
\par
The last issue was an operational one, which was the result of ice crystals and snow collecting on the entrance window of the photometer. Since the photometer is not heated and was designed to operate at ambient temperatures, the accumulation remained on the entrance window. Fortunately, during this observing campaign, we were able to periodically check up on the photometers and physically remove any accumulation. The ultimate objective is to leave the photometer unattended over an entire winter campaign, which requires a way to automatically remove accumulation or prevent it from forming.

\section{Future Plans}
\label{sec:futureplans}
Our goal is to obtain a full winter of measurements. This would provide a sufficient baseline for solid statistics of the NIR sky brightness near PEARL. An interesting comparison might be to the Summit Station at the apex of the Greenland ice sheet. It may be that the darkest sky brightness will occur during the middle of winter when there is little to no sunlight incident on the upper atmosphere for long periods. In addition to carrying out these measurements in the Arctic, we are planning to install another version of the photometer at a mid-latitude astronomical site to obtain a long baseline comparison dataset. However, to accomplish this task, we need to resolve the thermal issues that we discussed earlier. As it stands, more work needs to be done for the photometer to function in all environmental conditions.
\par
One other area we need to improve upon is the calibration of the photometer. We have not yet accurately measured the angular response function or carried out a flux calibration across the full bandpass of the photometer. There are two possible avenues to improve the calibration. We can use an integrating sphere to generate a photometric standard and measure the response of the photometer over all angles and illuminate the entire entrance aperture. The second method requires the simultaneous measurement of the $J$-band sky background with both the photometer and a NIR camera with a $J$-band filter. The NIR camera can be first used to observe at a stellar photometric standard to carry out a flux calibration of the camera system. By observing blank sky and then a known flux calibrator, one can carry out an absolute measurement of the $J$-band sky background with the camera system.

\section{Conclusions}

We set out to test the hypothesis that the NIR sky brightness is lower in polar regions due to suppressed OH airglow emission. Measurements in the high Antarctic plateau revealed that the $J$-band sky brightness was $2-3$ times (1 mag) fainter than the best mid-latitude sites like Mauna Kea. For this reason, we designed a rugged photometer that was suitable for operation in the High Arctic to measure the $J$-band sky background. We installed this device at PEARL near Eureka, Nunavut, Canada and carried out a short sky brightness measurement campaign in February 2012. Our measurements indicate that the $J$-band sky brightness varies between $15.5-15.9$ mag arcsec$^2,$ with a measurement uncertainty of 0.15 mag, which is dominated by systematic errors in our radiometric calibration. On our best night, we measured a fairly consistent sky brightness of $15.8\pm0.15$ mag arcsec$^2.$ This is not corrected for atmospheric extinction, which is usually $<0.1$ mag in the $J$-band at a good site. This sky brightness is comparable to the best mid-latitude sites, but is not as dark as the Antarctic measurements. We tested the idea that the increased solar activity at present times may explain the increased sky brightness. However, we discovered that using known correlations between OH line flux and 10.7cm solar flux to correct for this effect can only account for a 0.1 mag decrease in the sky brightness. One other possible explanation is that we simply do not have sufficient statistics on the variation of the sky brightness in the Arctic due to the short span of our measurements. Finally, we may underestimate the solid angle of the photometer, which means our measured sky brightness value is actually somewhat higher than the true value. Future improvements in our calibration techniques and a winter-long sky brightness measurement campaign should resolve these issues.

\label{sec:conclusions}

\acknowledgments

S.S. was supported by the Dunlap Institute through the Dunlap fellowship program. S.S. would like to thank Bill Reeve and Mark Anderson of Lockheed Martin Corporation for our long discussions about making precise pico amp level current measurements and the intricacies of op amps. S.S. would also like to thank Adam Palmentieri of Hamamatsu Corporation for his assistance in troubleshooting spurious dark current in the readout circuit. We are grateful to the Canadian Network for the Detection of Atmospheric Change, particularly Pierre Fogal and James Drummond for their assistance in deploying the photometer on the PEARL roof. The Eureka weather station is operated by Environment Canada, and E.S. and W.N. greatly appreciate the warm hospitality of their staff during our \"observing run." This research uses data from Space Weather Canada solar flux database.