\section{Introduction}\label{sec:intro}
Kenzo~\cite{Kenzo} is a Computer Algebra System devoted to Algebraic Topology which was developed by F. Sergeraert. This
system has computed some homology and homotopy groups which cannot be easily obtained by theoretical or computational
means; some examples can be seen in~\cite{RR2012}. Therefore, in this situation, it makes sense to analyze the Kenzo
programs in order to ensure the correctness of the mathematical results which are obtained thanks to it. To this aim,
two different research lines were launched some years ago to apply formal methods in the study of Kenzo.

On the one hand, the ACL2 theorem prover has been used to verify the correctness of actual Kenzo \emph{programs},
see~\cite{LOPSTR10,LMRR2}. ACL2 fits perfectly to this task since Kenzo is implemented in Common Lisp~\cite{Gra96}, the same
language in which ACL2 is built on. Nevertheless, since the ACL2 logic is first-order, the full verification of Kenzo
is not possible, because it uses intensively higher order functional programming. On the other hand, some instrumental
Kenzo \emph{algorithms}, involving higher-order logic, have been formalized in the proof assistants Isabelle/HOL and
Coq. Namely, we can highlight the formalizations of the Basic Perturbation Lemma in Isabelle/HOL, see~\cite{ABR08},
and the Effective Homology of Bicomplexes in {\sc Coq}, published in~\cite{DR2011}.

The work presented in this paper goes in the same direction that the latter approach, formalizing Kenzo \emph{algorithms}.
In particular, we have focused on the formalization of \emph{Discrete Vector Fields}, a powerful notion which will play
a key role in the new version of Kenzo; see the Kenzo web page~\cite{Kenzo}. To carry out this task, we will use the {\sc Coq }
proof assistant~\cite{Coq} and its {\sc SSReflect } library~\cite{SSReflect}.

The importance of Discrete Vector Fields, which were first introduced in~\cite{Forman}, stems from the fact that they
can be used to considerably reduce the amount of information of a discrete object but preserving homological properties.
In particular, we can use discrete vector fields to deal with biomedical images inside {\sc Coq } in a reasonable amount of time.

The rest of this paper is organized as follows. In the next section, we introduce some mathematical preliminaries,
which are encoded abstractly in {\sc Coq } in Section~\ref{sec:foandaics}. Such an abstract version is \emph{refined} to an effective one
in Section~\ref{sec:aeifhtc}; namely, the implementation and formal verification of the main algorithm involved in our
developments are presented there. In order to ensure the feasibility of our programs, a major issue when applying formal
methods, we use them to study a biomedical problem in Section~\ref{sec:atbi}. The paper ends with a section of
Conclusions and Further work, and the Bibliography.

The interested reader can consult the complete development in
\url{http://wiki.portal.chalmers.se/cse/pmwiki.php/ForMath/ProofExamples#wp3ex5}.

\section{Mathematics to formalize}\label{sec:mtf}
In this section, we briefly provide the minimal mathematical background needed to understand
the rest of the paper. We mainly focus on definitions which, mainly, come from the algebraic setting of
discrete Morse theory presented in~\cite{cvd} and the Effective Homology theory~\cite{RS02}.
We assume as known the notions of \emph{ring}, \emph{module} over a ring and \emph{module morphism}
(see, for instance,~\cite{Basic-Algebra}).

First of all, let us introduce one of the main notions in the context of Algebraic Topology: \emph{chain complexes}.

\begin{Defn}
A \emph{chain complex} $C_\ast$ is a pair of sequences $(C_n,d_n)_{n\in\mathbb{Z}}$ where for every $n\in\mathbb{Z}$,
$C_n$ is an ${\mathcal R}$-module and $d_n:C_n\rightarrow C_{n-1}$ is a module morphism, called the \emph{differential map},
such that the composition $d_nd_{n+1}$ is null. In many situations the ring ${\mathcal R}$ is either the integer ring,
${\mathcal R} =\mathbb{Z}$, or the field $\mathbb{Z}_2$. In the rest of this section, we will work with $\mathbb{Z}$ as
ground ring; later on, we will change to $\mathbb{Z}_2$.

\noindent The module $C_n$ is called the module of \emph{$n$-chains}. The image
$B_n = im~d_{n+1} \subseteq C_n$ is the (sub)module of $n$-boundaries. The kernel $Z_n = \ker d_n \subseteq C_n$ is the
(sub)module of \emph{$n$-cycles}.
\end{Defn}

Given a chain complex $C_\ast=(C_n,d_n)_{n \in \mathbb{Z}}$, the identities $d_{n-1}\circ d_n= 0$ mean the inclusion
relations \(B_n \subseteq Z_n\): every boundary is a cycle (the converse in general is not true). Thus the next definition makes sense.

\begin{Defn}
The \emph{$n$-homology group} of $C_\ast$, denoted by $H_n(C_\ast)$, is defined as the quotient $H_n(C_\ast)=Z_n/B_n$
\end{Defn}

Chain complexes have a corresponding notion of morphism.

\begin{Defn}
Let $C_\ast=(C_n,d_n)_{n\in\mathbb{Z}}$ and $D_\ast=(D_n,\widehat{d}_n)_{n\in\mathbb{Z}}$ be two chain complexes. A \emph{chain
complex morphism} $f:C_\ast \rightarrow D_\ast$ is a family of module morphisms, $f=\{f_n:C_n \rightarrow D_n\}_{n\in\mathbb{Z}}$,
satisfying for every $n\in\mathbb{Z}$ the relation $f_{n-1}d_n = \widehat{d}_nf_n$. Usually, the sub-indexes are skipped, and we just write
$fd_C = d_Df$.
\end{Defn}

Now, we can introduce one of the fundamental notions in the effective homology theory.

\begin{Defn}
\label{defn:rdct}
A \emph{reduction} $\rho$ between two chain complexes $C_\ast$ and $D_\ast$, denoted in this paper by $\rho : C_\ast \mbox{\,\(\Rightarrow\hspace{-9pt}\Rightarrow\)\,} D_\ast$, is a triple $\rho=(f,g,h)$
where $f:C_\ast\rightarrow D_\ast$ and $g:D_\ast\rightarrow C_\ast$ are chain complex morphisms, $h = \{h_n:C_n \rightarrow C_{n+1}\}_{n\in\mathbb{Z}}$
is a family of module morphism, and the following relations are satisfied:
\begin{enumerate}
\item [1)]$f \circ g = Id_{D_\ast}$;
\item [2)]$d_{C}\circ h + h\circ  d_{C} = Id_{C_\ast}- g\circ f$;
\item [3)]$f \circ h = 0$; \quad $h \circ g = 0$; \quad $h \circ h = 0$.
\end{enumerate}
\end{Defn}

\noindent The importance of reductions lies in the fact that given a reduction $\rho~:~C_\ast \mbox{\,\(\Rightarrow\hspace{-9pt}\Rightarrow\)\,} D_\ast$, then
$H_n(C_\ast)$ is isomorphic to $H_n(D_\ast)$ for every $n\in \mathbb{Z}$. Very frequently, $D_\ast$ is a much smaller
chain complex than $C_\ast$, so we can compute the homology groups of $C_\ast$ much faster by means of those of $D_\ast$.

Let us state now the main notions coming from the algebraic setting of Discrete Morse Theory~\cite{cvd}.

\begin{Defn}
Let $C_\ast=(C_n,d_n)_{n \in \mathbb{Z}}$ be a free chain complex with distinguished $\mathbb{Z}$-basis $\beta_n \subset C_n$.
A \emph {discrete vector field} $V$ on $C_\ast$ is a collection of pairs $V = \{(\sigma_i; \tau_i)\}_{i\in I}$ satisfying the conditions:
\begin{itemize}
\item Every $\sigma_i$ is some element of $\beta_n$, in which case $\tau_i\in \beta_{n+1}$.
The degree $n$ depends on $i$ and in general is not constant.
\item Every component $\sigma_i$ is a \emph{regular face} of the corresponding $\tau_i$ (regular face means that the coefficient of
$\sigma_i$ in $d_{n+1}\tau_i$ is $1$ or $-1$).
\item Each generator (\emph{cell}) of $C_\ast$ appears at most one time in $V$.
\end{itemize}
\end{Defn}

It is not compulsory all the cells of $C_\ast$ appear in the vector field $V$.

\begin{Defn}
A cell $\chi$ which does not appear in a discrete vector field $V=\{(\sigma_i; \tau_i)\}_{i\in I}$ is called a \emph{critical cell}.
\end{Defn}

From a discrete vector field on a chain complex, we can introduce $V$-paths.

\begin{Defn}
A $V$-path of degree $n$ and length $m$ is a sequence $((\sigma_{i_k}, \tau_{i_k}))_{0 \leq k < m}$ satisfying:
\begin{itemize}
\item Every pair $((\sigma_{i_k}, \tau_{i_k}))$ is a component of $V$ and $\tau_{i_k}$ is a $n$-cell.
\item For every $0 < k < m$, the component $\sigma_{i_k}$ is a face of $\tau_{i_{k-1}}$ (the coefficient of $\sigma_{i_k}$ in $d_{n}\tau_{i_{k-1}}$ is non-null)
different from $\sigma_{i_{k-1}}$.
\end{itemize}
\end{Defn}

Now we can present the notion of \emph{admissible} discrete vector field on a chain complex, a concept which can be understood as
a \emph{recipe} indicating both the ``useless'' elements of the chain complex (in the sense, that they can be removed without
changing its homology) and the \emph{critical} ones (those whose removal modifies the homology).

\begin{Defn}
A discrete vector field $V$ is \emph{admissible} if for every $n \in \mathbb{Z}$, a function $\lambda_n: \beta_n \rightarrow \mathbb{N}$ is
provided satisfying the following property: every $V$-path starting from $\sigma \in \beta_n$ has a length bounded by $\lambda_n(\sigma)$.
\end{Defn}

Finally, we can state the theorem where Discrete Morse Theory and Effective Homology converge.

\begin{Thm}~\cite[Theorem~19]{cvd}\label{thm:advred}
Let $C_\ast=(C_n,d_n)_{n \in \mathbb{Z}}$ be a free chain complex and $V$ be an admissible discrete vector
field on $C_\ast$. Then the vector field $V$ defines a canonical reduction $\rho : (C_n, d_n)\mbox{\,\(\Rightarrow\hspace{-9pt}\Rightarrow\)\,} (C^c_n, d^c_n)$ where $C^c_n = \mathbb{Z}[\beta^c_n]$
is the free $\mathbb{Z}$-module generated by $\beta^c_n$, the critical $n$-cells.
\end{Thm}

\noindent Therefore, as the bigger the admissible discrete vector field $V$ the smaller the chain complex $C^c_\ast$, we need algorithms
which produce admissible discrete vector fields as large as possible.

If we consider the case of \emph{finite type} chain complexes, where there is a finite number of generators in each dimension of the chain complex,
the differential maps can be represented as matrices. In that case, the problem of finding an admissible discrete vector field on the chain complex
can be solved through the computation of an admissible vector field for those matrices.

\begin{Defn}\label{def:dvfM}
Let $M$ be a matrix with coefficients in $\mathbb{Z}$, and with $m$ rows and $n$ columns. A discrete vector field $V$
for this matrix is a set of natural pairs $\{(a_i, b_i)\}$ satisfying these conditions:
\begin{enumerate}
 \item $1\leq a_i \leq m$ and $1\leq b_i \leq n$.
 \item The entry $M[a_i,b_i]$ of the matrix is $\pm 1$.
 \item The indexes $a_i$ (resp. $b_i$) are pairwise different.
\end{enumerate}
\end{Defn}

Given $V$ be a vector field for our matrix $M$, we need to know if $V$ is admissible. If $1\leq a,a'\leq m$, with $a\neq a'$, we can
decide $a > a'$ if there is an elementary $V$-path from $a$ to $a'$, that is, if a vector $(a,b)$ is present in $V$ and the entry $M[a',b]$
is non-null. In this way, a binary relation is obtained. Then the vector field $V$ is admissible if and only if this binary relation
generates a partial order, that is, if there is no loop $a_1 > a_2 > \ldots > a_k = a_1$.

Eventually, given a matrix and an admissible discrete vector field on it, we can construct a new matrix, smaller than the original one,
preserving the homological properties. This is the equivalent version of Theorem~\ref{thm:advred} for matrices, a detailed description
of the process can be seen in~\cite[Proposition~14]{cvd}.

In the rest of the paper, we will focus on the formally certified construction of an admissible discrete vector field from a matrix.
The task of verifying the reduction process remains as further work.

\section{A non deterministic algorithm in {\sc SSReflect }}\label{sec:foandaics}
First of all, we have provided in {\sc Coq}/{\sc SSReflect } an abstract formalization of admissible discrete vector fields on matrices
and a non deterministic algorithm to construct an admissible discrete vector field from a matrix\footnote{Thanks
are due to Maxime D\'en\`es and Anders M\"ortberg which guided us in this development.}. {\sc SSReflect}
is an extension for the {\sc Coq} proof assistant, which was developed by G. Gonthier while formalizing
the Four Color Theorem~\cite{FCT}. Nowadays, it is used in the formal proof of the Feit-Thompson theorem~\cite{MathComp}.

{\sc SSReflect} provides all the necessary tools to achieve our goal. In particular, we take
advantage of the \verb"matrix", \verb"ssralg" and \verb"fingraph" libraries, which formalize,
respectively, matrix theory, the main algebraic structures and the theory of finite graphs.

First of all, we are going to define admissible discrete vector field on a matrix $M$ with coefficients in
a ring $\mathcal{R}$, and with $m$ rows and $n$ columns. It is worth noting that our matrices
are defined over a generic ring instead of working with coefficients in $\mathbb{Z}$ since
the {\sc SSReflect} implementation of $\mathbb{Z}$, see~\cite{CM11}, is not yet included in the {\sc SSReflect } distributed
version. The vector fields are represented by a sequence of pairs where the first component is an ordinal $m$ and the second one an
ordinal $n$.

\begin{verbatim}
Variable R : ringType.
Variables m n : nat.
Definition vectorfield := seq ('I_m * 'I_n).
\end{verbatim}

Now, we can define in a straightforward manner a function, called \verb"dvf", which given a matrix \verb"M" (with coefficients in
a ring $\mathcal{R}$, and with $m$ rows and $n$ columns, \verb"'M[R]_(m,n)") and an object \verb"V" of type \verb"vectorfield"
checks whether \verb"V" satisfies the properties of a discrete vector field on \verb"M" (Definition~\ref{def:dvfM}).

\begin{verbatim}
Definition dvf (M : 'M[R]_(m,n)) (V : vectorfield) :=
  all [pred p | (M p.1 p.2 == 1) || (M p.1 p.2 == -1)] V &&
  (uniq (map (@fst _ _) V) && uniq (map (@snd _ _) V)).
\end{verbatim}

\noindent It is worth noting that the first condition of Definition~\ref{def:dvfM} is implicit in the \verb"vectorfield" type. Now, as we have
explained at the end of the previous section, from a discrete vector field \verb"V" a binary relation is obtained between the first elements of each pair
of \verb"V". Such a binary relation will be encoded by means of an object of the following type.

\begin{verbatim}
Definition orders := (simpl_rel 'I_m).
\end{verbatim}

Finally, we can define a function, which is called \verb"advf", that given a matrix \verb"'M[R]_(m,n)", \verb"M", a \verb"vectorfield",
\verb"V" and an \verb"orders", \verb"ords", as input, tests whether both \verb"V" satisfies the properties of a discrete vector field
on \verb"M" and the admissibility property for the relations, \verb"ords", associated with the vector field, \verb"V". In order to test the admissibility
property we generate the transitive closure of \verb"ords", using the \verb"connect" operator of the \verb"fingraph" library,
and subsequently check that there is not any path between the first element of a pair of \verb"V" and itself.

\begin{verbatim}
Definition advf (M:'M[R]_(m,n)) (V:vectorfield) (ords:orders) :=
 dvf M V && all [pred i|~~(connect ords i i)] (map (@fst _ _) V).
\end{verbatim}

Now, let us define a non deterministic algorithm which construct an admissible discrete vector field from a matrix. Firstly,
we define a function, \verb"gen_orders", which generates the relations between the elements of the discrete vector
field as we have explained at the end of the previous section.

\begin{verbatim}
Definition gen_orders (M0 : 'M[R]_(m,n)) (i:'I_m) j :=
  [rel i x | (x != i) && (M0 x j != 0)].
\end{verbatim}

Subsequently, the function, \verb"gen_adm_dvf", which generates an admissible discrete vector field from a matrix is introduced.
This function invokes a recursive function, \verb"genDvfOrders", which in each step adds a new component to the vector field
in such a way that the admissibility property is fulfilled. The recursive algorithm stops when either there is not any new element
whose inclusion in the vector field preserves the admissibility property or the maximum number of elements of the discrete vector field
(which is the minimum between the number of columns and the number of rows of the matrix) is reached.

\begin{verbatim}
Fixpoint genDvfOrders M V (ords : simpl_rel _) k :=
  if k is l.+1 then
    let P := [pred ij | admissible (ij::V) M
                        (relU ords (gen_orders M ij.1 ij.2))] in
    if pick P is Some (i,j)
       then genDvfOrders M ((i,j)::V)
                         (relU ords (gen_orders M i j)) l
    else (V, ords)
  else (V, ords).

Definition gen_adm_dvf M :=
  genDvfOrders M [::] [rel x y | false] (minn m n).
\end{verbatim}

Eventually, we can certify in a straightforward manner (just 4 lines) the correctness of the function \verb"gen_adm_dvf".

\begin{verbatim}
Lemma admissible_gen_adm_dvf m n (M : 'M[R]_(m,n)) :
 let (vf,ords) := gen_adm_dvf M in admissible vf M ords.
\end{verbatim}

As a final remark, it is worth noting that the function \verb"gen_adm_dvf" is not executable. On the one
hand, {\sc SSReflect } matrices are locked in a way that do not allow direct computations since they may trigger
heavy computations during deduction steps. On the other hand, we are using the \verb"pick" instruction, in
the definition of \verb"genDvfOrders", to choose the elements which are added to the vector field; however,
this operator does not provide an actual method to select those elements.

\section{An effective implementation: from Haskell to {\sc Coq }}\label{sec:aeifhtc}

In the previous section, we have presented a non deterministic algorithm to construct an admissible discrete vector field from a matrix.
Such an abstract version has been described on high-level datastructures; now, we are going to obtain from it a \emph{refined} version,
based on datastructures closer to machine representation, which will be executable.

The necessity of an executable algorithm which construct an admissible discrete vector field stems from the fact that they will
will play a key role to study biomedical images. There are several algorithms to construct an admissible discrete vector field;
the one that we will use is explained in~\cite{cvd} (from now on, called RS's algorithm; RS stands for Romero-Sergeraert). The
implementation of this algorithm will be executable but the proof of its correctness will be much more difficult than the one
presented in the previous section.

\subsection{The Romero-Sergeraert algorithm}\label{subsec:trsa}

\begin{comment}
{This algorithm returns an admissible discrete vector field built in the following way. It consists of running the rows of the matrix $M$ in the usual reading order looking for the entries whose value is $1$. This entry will be in a position of $M$, for instance $(i,j)$ which will be exactly the vector to add to the admissible discrete vector field. Let us note that every the entries can not be selected because of the properties of an admissible discrete vector field depicted in Section \ref{sec:mtf}. Moreover, to avoid that cycles are not created, that is to say, satisfy the admissibility property, we have to pay attention to the entries of the column \emph{j} different from $0$ which are in a different position from \emph{i}.}
\end{comment}

The underlying idea of the RS algorithm is that given an admissible discrete vector field, we try to enlarge it adding new vectors which
preserve the admissibility property. We can define algorithmically the RS algorithm as follows.

\begin{Alg}[The RS Algorithm]\label{alg:rs}
\textcolor{white}{.}\\
{\itshape Input: } a matrix $M$ with coefficients in $\mathbb{Z}$.\\
{\itshape Output: } an admissible discrete vector field for $M$.\\
{\itshape Description: }
  \begin{enumerate}
    \item Initialize the vector field, $V$, to the void vector field.
    \item Initialize the relations, $ords$, to nil.
    \item For every row, $i$, of $M$:
      \begin{enumerate}
        \item Search the first entry of the row equal to $1$ or $-1$, $j$.
        \item \verb"If" $(i,j)$ can be added to the vector field; that is, if we add it to $V$ and generate all
              the relations, the properties of an admissible discrete vector field are preserved.
        \begin{itemize}
          \item[] \verb"then":
                     \begin{itemize}
                        \item[-] Add $(i,j)$ to $V$.
                        \item[-] Add to $ords$ the corresponding relations generated from $(i,j)$.
                        \item[-] Go to the next row and repeat from Step 3.
                      \end{itemize}
          \item[] \verb"else": look for the next entry of the row whose value is $1$ or $-1$.
                     \begin{itemize}
                        \item[-] \verb"If" there is not any.
                         \begin{itemize}
                           \item[] \verb"then": go to the next row and repeat from Step 3.
                           \item[] \verb"else": go to Step 3.2 with $j$ the column of the entry whose value is $1$ or $-1$.
			\end{itemize}
                      \end{itemize}
        \end{itemize}
      \end{enumerate}
  \end{enumerate}
\end{Alg}

In order to clarify how this algorithm works, let us construct an admissible discrete vector field from the following matrix.

 \[ \left( \begin{array}{cccc}
1 & 1 & 0 & 0  \\
1 & 1 & 1 & 0 \\
0 & 0 & 1 & 1 \\
0 & 1 & 1 & 0
 \end{array} \right)\]

We start with the void vector field $V= \left\{\right\}$. Running the successive rows, we find $M\left[1,1\right] = 1$, and we include the
vector $(1,1)$ to $V$, obtaining $V = \left\{(1,1)\right\}$. Then, let us add the relations that, in this case is $1 > 2$ because
$M\left[2,1\right]\neq 0$. So, it will be forbidden to incorporate the relation $2>1$ as the cycle $1 > 2 > 1$ would appear.
Besides, the row 1 and the column 1 are now used and cannot be used anymore. So, we go on with the second row and find $M\left[2,1\right] = 1$,
but we cannot add $(2,1)$ as we have just said. Moreover, the element $(2,2)$ can not be incorporated because the cycle $1 > 2 > 1$ would be created.
So, we continue and find the next element, $M\left[2,3\right] = 1$. This does not create any cycle and satisfies the properties of a discrete vector field.
Then, we obtain $V=\left\{(2,3),(1,1)\right\}$ and the relations $2 > 3$ and $2 > 4$. Running the next row, the first element equal to $1$ is in the position $(3,3)$,
but we cannot include it due to the admissibility property. Therefore, we try with the last element of this row $(3,4)$. No relation is generated in this case
because in this column the only non null element is in the chosen position. So, $V=\left\{(3,4),(2,3),(1,1)\right\}$. Finally, we run the last row.
The elements that could be added are $(4,2),(4,3)$, but in both cases we would have to append the relation $4>2$. This would generate a cycle with one
of the previous restrictions, $2>4$. So, we obtain $V= \left\{(3,4),(2,3),(1,1)\right\}$ and the relations are: $1>2$, $2>3$ and $2>4$.

In general, Algorithm~\ref{alg:rs} can be applied over matrices with coefficients in a general ring ${\cal R}$. From now on, we will work with
${\cal R}=\mathbb{Z}_2$, since this is the usual ring when working with monochromatic images in the context of Digital Algebraic Topology.

The development of a formally certified implementation of the RS algorithm has followed the methodology presented in~\cite{and10}. Firstly,
we implement a version of our programs in \emph{Haskell}~\cite{haskell98}, a \emph{lazy} functional programming language. Subsequently, we
intensively test our implementation using \emph{QuickCheck}~\cite{Quickcheck}, a tool which allows one to automatically test properties
about programs implemented in Haskell. Finally, we verify the correctness of our programs using
the {\sc Coq } \emph{interactive} proof assistant and its {\sc SSReflect } library.

\subsection{A Haskell program}

The choice of Haskell to implement our programs was because both the code and the programming style is similar to the
ones of {\sc Coq }. In this programming language, we have defined the programs which implement the RS algorithm. The
description of the main function is shown here:

\begin{description}
\item[{\small gen\_admdvf\_ord $M$:}]
From a matrix $M$ with coefficients in $\mathbb{Z}_2$, represented as a list of lists, this function generates
an admissible discrete vector field for $M$, encoded by a list of natural pairs, and the relations, a list of lists of
natural numbers.
\end{description}

Let us emphasize that the function \verb"gen_admdvf_ord" returns a pair of elements. The former one, \verb"(gen_admdvf_ord M).1",
is a discrete vector field and the latter one, \verb"(gen_admdvf_ord M).2", corresponds to the relations associated with the vector field.
To provide a better understanding of these tools, let us apply them in the example presented in Subsection~\ref{subsec:trsa}.

\begin{verbatim}
> gen_admdvf_ord [[1,1,0,0],[1,1,1,0],[0,0,1,1],[0,1,1,0]]
[([(3,4),(2,3),(1,1)], [[2,4],[2,3],[1,2],[1,2,4],[1,2,3]])]
\end{verbatim}

Let us note that we return the transitive closure of the relations between the first components of the pairs of the discrete vector field. This will
make the proof of the correctness of our programs easier.

\subsection{Testing with QuickCheck}\label{subsec:twq}

Using QuickCheck can be considered as a good starting point towards the formal verification of our programs.
On the one hand, a specification of the properties which must be satisfied by our programs is given (a necessary
step in the formalization process). On the other hand, before trying a formal verification of our programs
(a quite difficult task) we are testing them, a process which can be useful in order to detect bugs.

In our case, we want to check that the output by \verb"gen_admdvf_ord" gives us an admissible discrete vector field.
Then, let $M$ be a matrix over $\mathbb{Z}_2$ with $m$ rows and $n$ columns, $V=(a_i, b_i)_i$ be a discrete
vector field from $M$ and $ords$ be the transitive closure of the relations associated with $V$, the properties to test are
the ones coming from Definition~\ref{def:dvfM} and the admissibility property adapted to the $\mathbb{Z}_2$ case.

\begin{enumerate}
  \item $1\leq a_i \leq m$ and $1\leq b_i \leq n$.
  \item $\forall i\;, M(a_i,b_i) = 1$.
  \item $(a_i)_i$ (resp. $(b_i)_i$) are pairwise different.
  \item $ords$ does not have any loop (admissibility property).
\end{enumerate}

These four properties has been encoded in Haskell by means of a function called \verb"isAdmVecfield". To test in QuickCheck
that our implementation of the RS algorithm fulfills the specification given in \verb"isAdmVecfield", the following
\emph{property definition}, using QuickCheck terminology, is defined.

\begin{verbatim}
condAdmVecfield M =
let advf = (gen_admdvf_ord M) in isAdmVecfield M (advf.1) (advf.2)
\end{verbatim}

The definition of \verb"condAdmVecfield" states that given a matrix \verb"M", the output
of \verb"gen_admdvf_ord", both the discrete vector field (first component) and the relations (second component) from $M$, fulfill
the specification of the property called \verb"isAdmVecfield". Now, we can test whether \verb"condAdmVecfield" satisfies such
a property.

\begin{verbatim}
> quickCheck condAdmVecfield
+ + + OK, passed 100 tests.
\end{verbatim}

The result produced by QuickCheck when evaluating this statement, means that QuickCheck has generated $100$ random values for \verb"M",
checking that the property was true for all these cases.

\subsection{Formalization in {\sc Coq }/{\sc SSReflect }}\label{subsec:fics}

After testing our programs, and as a final step to confirm their correctness, we can undertake the challenge of formally verify them.

First of all, we define the data types related to our programs, which are effective matrices, vector fields and relations.
We have tried to keep them as close as possible to the Haskell ones; therefore, a matrix is represented by means of a list
of lists over $\mathbb{Z}_2$, a vector field is a sequence of natural pairs and finally, the relations is a
list of lists of natural numbers.

\begin{verbatim}
Definition Z2 := Fp_fieldType 2.
Definition matZ2 := seq (seq Z2).
Definition vectorfield := seq (prod nat nat).
Definition orders := seq (seq nat).
\end{verbatim}

Afterwards, we translate both the programs and the properties, which were specified during the testing of the programs,
from Haskell to {\sc Coq}, a task which is quite direct since these two systems are close.

\begin{comment}Let us emphasize that any modification in the code {\sc Coq} regarding to Haskell version has to be modified in this one. Moreover, the algorithms will have to be tested again to check that these changes do not affect to the hopeful result.\end{comment}

Then, we have defined a function \verb"isAdmVecfield" which receives as input a matrix over $\mathbb{Z}_2$, a vector field and the relations and
checks if the properties, explained in Subsection~\ref{subsec:twq}, are satisfied.
\begin{verbatim}
Definition isAdmVecfield (M:matZ2)(vf:vectorfield)(ord:orders):=
  ((longmn (size M) (getfirstElemseq vf)) /\
  (longmn (size (nth [::] M 0)) (getsndElemseq vf))) /\
  (forall i j:nat, (i , j) \in vf -> compij i j M = 1) /\
  ((uniq (getfirstElemseq vf)) /\ (uniq (getsndElemseq vf))) /\
  (admissible ord).
\end{verbatim}

Finally, we have proved the theorem \verb"genDvfisVecfield" which says that given a matrix $M$, the output produced by
\verb"gen_admdvf_ord" satisfies the properties specified in \verb"isAdmVecfield".

\begin{verbatim}
Theorem genDvfisVecfield (M:matZ2):
  let advf := (gen_admdvf_ord M) in
  isAdmVecfield M (advf.1) (advf.2).
\end{verbatim}

\noindent We have split the proof of the above theorem into $4$ lemmas which correspond with each one of the properties that should
be fulfilled to have an admissible discrete vector field. For instance, the lemma associated with the first property of the definition of a
discrete vector field is the following one.

\begin{verbatim}
Lemma propSizef (M:matZ2):
  let advf := (gen_admdvf_ord M).1 in
  (longmn (size M) (getfirstElemseq advf) /\
  (longmn (size (nth nil M 0))(getsndElemseq advf)).
\end{verbatim}

Both the functions which implement the RS algorithm and the ones which specify the definitional properties of admissible discrete vector fields
are defined using a \emph{functional style}; that is, our programs are defined using \emph{pattern-matching} and \emph{recursion}. Therefore,
in order to reason about our recursive functions, we need elimination principles which are fitted for them. To this aim,
we use the tool presented in~\cite{BC02} which allows one to reason about complex recursive definitions since {\sc Coq } does not directly generate
elimination principles for complex recursive functions. Let us see how the tool presented in~\cite{BC02} works.

In our development of the implementation of the RS algorithm, we have defined a function, called \verb"subm", which takes as arguments
a natural number, \verb"n", and a matrix, \verb"M", and removes the first \verb"n" rows of \verb"M". The inductive
scheme associated with \verb"subm" is set as follows.

\begin{verbatim}
Functional Scheme subm_ind := Induction for subm Sort Prop.
\end{verbatim}

Then, when we need to reason about \verb"subm", we can apply this scheme with the corresponding parameters using the instruction
\verb"functional induction". However, as we have previously said both our programs to define the RS algorithm and the ones which
specify the properties to prove are recursive. Then, in several cases, it is necessary to merge several inductive schemes to induction
simultaneously on several variables.  For instance, let $M$ be a matrix and $M'$ be a submatrix of $M$ where we have removed the
$(k-1)$ first rows of $M$; then, we want to prove that $\forall j,\; M (i, j) = M' (i - k + 1, j)$. This can be stated in {\sc Coq } as follows.

\begin{verbatim}
Lemma Mij_subM (i k: nat) (M: matZ2):
  k <= i -> k != 0 -> let M' := (subm k M) in
  M i j ==  M' (i - k + 1) j.
\end{verbatim}

To prove this lemma it is necessary to induct simultaneously on the parameters \verb"i", \verb"k" and \verb"M", but
the inductive scheme generated from \verb"subm" only applies induction on \verb"k" and \verb"M". Therefore,
we have to define a new recursive function, called \verb"Mij_subM_rec", to provide a proper inductive scheme to prove this theorem.

\begin{verbatim}
Fixpoint Mij_subM_rec (i k: nat) (M: matZ2) :=
match k with
|0 => M
|S p => match M with
        |nil => nil
        |hM::tM => if (k == 1)
                    then a::b
                  else (Mij_subM_rec p (i- 1) tM)
        end
end.
\end{verbatim}

This style of proving functional programs in {\sc Coq } is the one followed in the development of the proof of Theorem \verb"genDvfisVecfield".

\subsection{Experimental results}

Using the same methodology presented throughout this section, we are working in the formalization of the algorithm which, from a matrix and an
admissible discrete vector field on it, produces a reduced matrix preserving the homological properties of the original one. Up to now,
we have achieved a Haskell implementation which has been both tested with QuickCheck and translated into {\sc Coq}; however, the proof of
its correctness remains as further work.

Anyway, as we have a {\sc Coq } implementation of that procedure, we can execute some examples inside this proof assistant. Namely,
we have integrated the programs presented in this paper with the ones devoted to the computation of homology groups of digital images
introduced in~\cite{HDMMPS}; and we have considered matrices coming from $500$ randomly generated images.

The size of the matrices associated with those images was initially around $100\times 300$, and after the reduction process the average
size was $5\times 50$. Using the original matrices {\sc Coq } takes around 12 seconds to compute the homology from the matrices; on the contrary,
using the reduced matrices {\sc Coq } only needs milliseconds. Furthermore, as we will see in the following section, we have studied some images
which are originated from a real biomedical problem.

\section{Application to biomedical images}\label{sec:atbi}
Biomedical images are a suitable benchmark for testing our programs. On the one hand, the amount of information
included in this kind of images is usually quite big; then, a process able to reduce those images but keeping
the homological properties can be really useful. On the other hand, software systems dealing with biomedical
images must be trustworthy; this is our case since we have formally verified the correctness of our programs.

As an example, we can consider the problem of counting the number of \emph{synapses} in a neuron. Synapses~\cite{Neuroscience}
are the points of connection between neurons and are related to the computational capabilities of the brain. Therefore,
the treatment of neurological diseases, such as Alzheimer, may take advantage of procedures modifying the number
of synapses~\cite{C11}.

Up to now, the study of the \emph{synaptic density evolution} of neurons was a time-consuming task since it was performed, mainly,
manually. To overcome this issue, an automatic method was presented in~\cite{HMPR11-ACA}. Briefly speaking, such a process can be
split into two parts. Firstly, from three images of a neuron (the neuron with two antibody markers and the structure of the neuron),
a monochromatic image is obtained, see Figure~\ref{fig:neuron}\footnote{The same images with higher resolution can be seen in \url{http://www.unirioja.es/cu/joheras/synapses/}}.
In such an image, each connected component represents a synapse. So, the problem of measuring the number of synapses
is translated into a question of counting the connected components of a monochromatic image.

\begin{figure}
\centering
 \includegraphics[scale=0.14]{combinacion.png}
\caption{Synapses extraction from three images of a neuron}\label{fig:neuron}
\end{figure}

In the context of Algebraic Digital Topology, this issue can be tackled by means of the computation of the homology group $H_0$
of the monochromatic image. This task can be performed in {\sc Coq } through the formally verified programs
which were presented in~\cite{HDMMPS}. Nevertheless, such programs are not able to handle images like the one of the right side
of Figure~\ref{fig:neuron} due to its size (let us remark that {\sc Coq } is a proof assistant tool and not a computer algebra system). In order
to overcome this drawback, as we have explained at the end of the previous section, we have integrate our reduction programs with
the ones presented in~\cite{HDMMPS}. Using this approach, we can successfully compute the homology of the biomedical images in just $25$ seconds,
a remarkable time for an execution inside {\sc Coq}.

\section{Conclusions and Further work}
In this paper, we have given the first step towards the formal verification of a procedure
which allows one to study homological properties of big digital images inside {\sc Coq}. The underlying idea consists
in building an admissible discrete vector field on the matrices associated with an image and,
subsequently reduce those matrices but preserving the homology.

Up to now, we have certified the former step of this procedure, the construction of an admissible discrete
vector field from a matrix, both in an abstract and a concrete way. The reason because the abstract formalization is
useful is twofold: on the one hand, it provides a high-level theory close to usual mathematics, and, on the other hand,
it has been refined to obtain the effective construction of admissible discrete vector fields. As we have explained, there are several
heuristics to construct an admissible discrete vector field from a matrix, the one that we have
chosen is the RS algorithm~\cite{cvd} which produces, as we have experimentally seen, quite large discrete
vector fields, a desirable property for these objects. The latter step, the process to reduce the matrices,
is already specified in {\sc Coq}, but the proof of its correctness is still an ongoing work.

The suitability of our approach has been tested with several examples coming from randomly generated images
and also real images associated with a biomedical problem, the study of synaptic density evolution. The results
which have been obtained are remarkable since the amount of time necessary to compute homology groups
of such images inside {\sc Coq } is considerably reduced (in fact, it was impossible in the case of biomedical images).

As further work, we have to deal with some formalization issues. Namely, we have to verify that a reduction
can be constructed from a matrix and an admissible discrete vector field to a reduced matrix. Moreover,
we have hitherto worked with matrices over the ring $\mathbb{Z}_2$; the more general case
of matrices with coefficients in a ring ${\cal R}$ (with convenient constructive properties)
should be studied.

As we have seen in Subsection~\ref{subsec:fics}, it is necessary the definition of inductive schema
which fits to our complex recursive programs. Then, this opens the door to an integration between
{\sc Coq } and the ACL2 Theorem Prover~\cite{ACL2}. ACL2 has good heuristics to generate inductive schemes
from recursive functions; so, we could translate our functional programs from {\sc Coq } to ACL2, generate
the inductive schemes in ACL2; and finally return such inductive schemes to {\sc Coq}. Some preliminary
experiments have been performed to automate that process, obtaining encouraging results.

In a different research line, we can consider the study of more complex biomedical problems
using our certified programs. As an example, the recognition of the structure of neurons seems to
involve the computation of homology groups in higher dimensions; a question which could be tackled
with our tools.

\bibliographystyle{abbrv}


\section{Introduction}

Open clusters (OCs) have been widely used to investigate the existence of trends
in the Galactic disk such as radial and vertical metallicity
gradients or an age-metallicity relationship
\citepads[e.g.][]{1979ApJS...39..135J,1980A&A....81..375P,1997AJ....114.2556T,
2002AJ....124.2693F,2003AJ....125.1397C,2010A&A...511A..56P,2011A&A...535A..30C}
.
However, our knowledge of OC properties
is still far from complete. Age and distance estimations, mainly obtained
from isochrone fitting, are available for about $\sim$70\% of the more than
2100 OCs known in our Galaxy \citepads{2002A&A...389..871D}\footnote{The updated
version of this catalog can be found at {\tt
http://www.astro.iag.usp.br/~wilton/}.}. However, radial velocities have been
determined only for a 24\% of them. The picture is even worse in the case of
chemical compositions, which have been obtained for only 9\%, mainly
by means of different studies in the Washington
\citepads[e.g.][]{2005MNRAS.363.1247P}, DDO
\citepads[e.g.][]{1977AJ.....82...35J,1999A&AS..134..301C}, Str\"omgren
\citepads[e.g.][]{2003A&A...403..937P}, UBV
\citepads[e.g.][]{1998A&AS..128..131C}, IR
\citepads[e.g.][]{1999A&A...343..825C,2000A&A...353..147V}, and Vilnius
\citepads[e.g.][]{2011BaltA..20...27B,2011BaltA..20....1Z} photometric systems,
often giving rise to
considerable differences among them and to those obtained from spectroscopy.

Reliable information about the chemical composition can be retrieved only
from spectroscopy. However, the acquisition of high-resolution spectra
(R$\geq$20000), which is the only
way to derive detailed abundances, needs large amounts of telescope time. This
together with the complexity of the associated analysis explain why this kind of
study has been performed for only 4\% of known OCs
\citepads{2011A&A...535A..30C}. The alternative is to perform low- and
medium-resolution spectroscopy, although these data
provide information only information about the metallicities and radial
velocities of
stars
\citepads[e.g.][]{1993A&A...267...75F,2002AJ....124.2693F,2007AJ....134.1298C,
2009MNRAS.393..272W}.

\begin{table*}[htb]
\caption{Observing logs and program star information.}
\label{obsstars}
\centering
\renewcommand{\footnoterule}{}
\begin{scriptsize}
\begin{tabular}{l c c c c c c c c c c c c c c}
\hline\hline
	Cluster & Star\footnote{Identification taken from WEBDA database.} & $\alpha_{2000}$ & $\delta_{2000}$ & V & V-I & B-V & V$_{r}$ & EW$_{8498}$ & EW$_{8542}$ & EW$_{8662}$ & $\Sigma$ Ca & t$_{exp}$ & S/N$^{tot}$ & Note\\
	 &  & (hrs) & (deg) & (mag) & (mag) &(mag) & (km s$^{-1}$) & (\AA) & (\AA) & (\AA) & (\AA) & (sec) &  &\\
            \hline
Be~26    & 1038 & 06:50:02.9 & +05:44:59 & 18.845 & 1.283 &  ---  &  66.1$\pm$11.9 &      ---      &      ---      &      ---      &       ---      &  2$\times$950 &  25 & 2 \\
         & 1155 & 06:50:09.5 & +05:44:28 & 15.754 & 1.729 &  ---  &  71.0$\pm$7.9  & 1.17$\pm$0.13 & 3.39$\pm$0.12 & 2.74$\pm$0.12 &  7.30$\pm$0.21 &  2$\times$950 &  29 & 1 \\
         & 1231 & 06:50:08.6 & +05:44:07 & 16.702 & 1.819 &  ---  &  66.1$\pm$7.7  & 1.13$\pm$0.14 & 3.21$\pm$0.21 & 2.77$\pm$0.10 &  7.11$\pm$0.27 & 2$\times$1500 &  25 & 1 \\
         & 1288 & 06:50:05.7 & +05:43:53 & 15.500 & 1.735 &  ---  &  68.3$\pm$7.2  & 1.27$\pm$0.11 & 3.51$\pm$0.12 & 2.93$\pm$0.12 &  7.71$\pm$0.20 &  2$\times$850 &  29 & 1 \\
         & 1421 & 06:50:18.4 & +05:43:22 & 15.556 & 2.063 &  ---  &  68.3$\pm$6.7  & 1.45$\pm$0.08 & 3.94$\pm$0.09 & 2.62$\pm$0.07 &  8.01$\pm$0.14 &  2$\times$900 &  40 & 1 \\
         & 1650 & 06:50:14.9 & +05:42:17 & 14.974 & 2.092 &  ---  &  97.1$\pm$7.8  &      ---      &      ---      &      ---      &       ---      &  2$\times$850 &  45 & 9 \\
Be~70    & 1105 & 05:25:44.8 & +41:56:44 & 12.651 & 2.035 & 2.014 & -13.8$\pm$7.6  & 1.94$\pm$0.05 & 4.66$\pm$0.05 & 3.50$\pm$0.06 & 10.10$\pm$0.09 &  2$\times$700 & 186 & 1 \\
         & 1088 & 05:25:45.0 & +41:55:55 & 14.097 & 1.994 & 1.849 &  -4.3$\pm$7.5  & 1.61$\pm$0.02 & 4.32$\pm$0.02 & 3.20$\pm$0.03 &  9.13$\pm$0.04 &  2$\times$800 &  98 & 1 \\
         & 0820 & 05:25:49.9 & +41:56:51 & 14.482 & 1.944 & 1.816 & -17.3$\pm$7.9  & 1.53$\pm$0.04 & 4.18$\pm$0.03 & 2.99$\pm$0.04 &  8.70$\pm$0.06 &  2$\times$800 &  69 & 1 \\
         & 0609 & 05:25:54.6 & +41:58:22 & 14.608 & 1.734 & 1.585 & -11.7$\pm$7.2  & 1.56$\pm$0.04 & 3.69$\pm$0.05 & 2.80$\pm$0.06 &  8.05$\pm$0.09 &  2$\times$850 &  44 & 1 \\
         & 0811 & 05:25:50.1 & +41:57:46 & 14.774 & 1.796 & 1.645 & -19.7$\pm$7.7  & 1.46$\pm$0.08 & 3.87$\pm$0.07 & 2.94$\pm$0.08 &  8.27$\pm$0.12 &  2$\times$850 &  58 & 1 \\
         & 0690 & 05:25:52.7 & +41:57:13 & 14.923 & 1.617 & 1.457 & -20.4$\pm$7.3  & 1.49$\pm$0.05 & 3.71$\pm$0.07 & 2.71$\pm$0.06 &  7.91$\pm$0.11 &  2$\times$950 &  58 & 1 \\
         & 0735 & 05:25:51.9 & +41:57:35 & 15.142 & 1.762 & 1.624 & -19.7$\pm$7.9  & 1.46$\pm$0.08 & 3.89$\pm$0.08 & 2.80$\pm$0.10 &  8.15$\pm$0.15 &  2$\times$950 &  40 & 1 \\
         & 0171 & 05:26:03.3 & +41:58:54 & 15.304 & 1.716 & 1.547 & -12.9$\pm$8.9  & 1.35$\pm$0.12 & 4.01$\pm$0.13 & 2.35$\pm$0.13 &  7.71$\pm$0.22 &  2$\times$900 &  30 & 1 \\
NGC~1798 & 0005 & 05:11:36.7 & +47 41 49 & 14.380 & 1.633 & 1.483 &   4.3$\pm$7.3  & 1.54$\pm$0.06 & 3.79$\pm$0.07 & 2.91$\pm$0.07 &  8.24$\pm$0.12 &  2$\times$800 &  56 & 1 \\
         & 0009 & 05:11:42.8 & +47:41:32 & 14.777 & 1.653 & 1.486 &   0.2$\pm$7.3  & 1.31$\pm$0.10 & 3.45$\pm$0.10 & 2.60$\pm$0.08 &  7.36$\pm$0.16 &  2$\times$750 &  40 & 1 \\
         & 0011 & 05:11:37.0 & +47:40:15 & 15.305 & 1.666 & 1.529 &  -3.4$\pm$7.0  & 1.40$\pm$0.09 & 3.85$\pm$0.11 & 2.71$\pm$0.08 &  7.96$\pm$0.16 &  2$\times$800 &  41 & 1 \\
         & 0013 & 05:11:41.2 & +47:40:41 & 15.325 & 1.660 & 1.524 &  -5.6$\pm$7.2  & 1.38$\pm$0.07 & 3.52$\pm$0.10 & 2.64$\pm$0.13 &  7.54$\pm$0.18 &  2$\times$850 &  40 & 1 \\
         & 0043 & 05:11:47.9 & +47:40:26 & 12.875 & 2.252 & 1.970 &  11.2$\pm$7.8  & 1.89$\pm$0.02 & 4.68$\pm$0.03 & 3.44$\pm$0.03 & 10.01$\pm$0.05 &  2$\times$700 & 110 & 1 \\
         & 0608 & 05:11:41.0 & +47:40:43 & 16.653 & 0.909 & 0.816 &  12.8$\pm$11.0 &      ---      &      ---      &      ---      &       ---      &  2$\times$850 &  25 & 2 \\
NGC~2266 & 0034 & 06:43:16.0 & +26 57 45 & 12.454 &  ---  & 0.957 &  21.0$\pm$7.7  &      ---      &      ---      &      ---      &       ---      &  2$\times$350 &  55 & 9 \\
         & 0067 & 06:43:20.2 & +26:57:32 & 10.538 &  ---  & 1.288 &  30.2$\pm$7.2  &      ---      &      ---      &      ---      &       ---      &   2$\times$50 &  55 & 9 \\
         & 0075 & 06:43:12.9 & +26:57:07 & 15.242 &  ---  & 0.341 & -18.0$\pm$9.2  &      ---      &      ---      &      ---      &       ---      &  2$\times$350 &  24 & 2\\
         & 0096 & 06:43:28.0 & +26:58:10 & 12.287 &  ---  & 1.253 &  -8.9$\pm$7.5  & 1.47$\pm$0.05 & 3.84$\pm$0.05 & 2.65$\pm$0.06 &  7.96$\pm$0.09 &  2$\times$250 &  58 & 1 \\
         & 0101 & 06:43:25.3 & +26:57:49 & 12.978 &  ---  & 1.187 & -19.1$\pm$7.8  & 1.46$\pm$0.06 & 3.71$\pm$0.06 & 2.25$\pm$0.10 &  7.42$\pm$0.13 &  2$\times$500 &  50 & 1 \\
         & 1011 & 06:43:25.3 & +26:57:50 & 12.117 &  ---  & 1.273 & -18.5$\pm$7.1  & 1.43$\pm$0.06 & 4.13$\pm$0.05 & 2.97$\pm$0.06 &  8.53$\pm$0.10 &  2$\times$500 &  58 & 1\\       \hline\hline
\end{tabular}
\tablefoot{ (1)  member star used both for radial velocity and metallicity determination; (2) member star used only for radial velocity determination; (9) non-member star. }
\end{scriptsize}
\end{table*}

The goal of this paper is to increase the number of clusters with radial
velocities and metallicities determined from spectroscopy. For this purpose, we
have selected four OCs, Berkeley~ 26, Berkeley~70, NGC~1798, and NGC~2266, that
had not
been studied spectroscopically before. These clusters are located towards the
Galactic anticenter, at relatively  large distances (R$_{GC}\geq$10 kpc) where
the trend of the radial metallicity seems to flatten
\citepads[e.g.][]{2011A&A...535A..30C}. Medium-resolution spectra in the region
of the near-infrared \ion{Ca}{II} triplet (CaT) at $\sim$8500\AA~were obtained
for several red giant branch (RGB) stars in each cluster. The
observations and data reduction are described in Sect.~\ref{sec2}, radial
velocities and metallicities are obtained in Sects.~\ref{sec3}--\ref{sec4}, the
results for each cluster are discussed and compared with the literature in
Sect.~\ref{sec5}, and the comparison of the results obtained here with the
trends observed in the disk is performed in Sec.~\ref{sec6}. Finally, our main
conclusions are summarised in
Sect.~\ref{sec7}.

\section{Observations and data reduction}\label{sec2}

\begin{table}
\begin{minipage}[htb]{\columnwidth}
\caption{Adopted cluster parameters.}
\label{clusterspropierties}
\centering
\renewcommand{\footnoterule}{}
\begin{tabular}{l c c c c c}
\hline\hline
Cluster   & E(B--V) & (m-M)$_o$ & Age & R$_{gc}$ & z \\
          & (mag)   & (mag)     & (Gyr)  & (kpc) & (kpc) \\
\hline
Be~26\footnote{Average of the values by \citetads{2010MNRAS.402.2720P} and  \citetads{2008PASJ...60.1267H}. See text for details.}
          & 0.68$\pm$0.06 & 13.8$\pm$0.5 & 4.25$\pm$0.25 & 13.9 & 0.24\\
Be~70\footnote{Values taken from \citetads{2002AJ....123..905A} which are similar to those derived by \citetads{2004PASJ...56..295H} with the exception that they derived an age of 2.8 Gyr. }
          & 0.48$\pm$0.10 & 13.1$\pm$0.3 & 4.0$\pm$0.8 & 12.6 & 0.26 \\
NGC~1798\footnote{Values taken from \citetads{1999MNRAS.304..883P} which are similar to others available in the literature within error bars \citepads[e.g.][]{2000A&A...359..347D,2002NewA....7..553T,2004A&A...414..163S,2007A&A...467.1065M}.}
          & 0.51$\pm$0.04 & 13.1$\pm$0.2 & 1.4$\pm$0.3 & 12.5 & 0.35 \\
NGC~2266\footnote{Values taken from \citetads{2012A&A...539A.125D}. Similar values have been obtained by other investigations in this cluster \citepads{1991AcA....41..191K,2000A&A...359..347D,2001NewA....6..293T,2002NewA....7..553T,2002A&A...388..158L,2004A&A...414..163S,2006MNRAS.371.1641P,2007A&A...467.1065M,2008BaltA..17...51M,2010A&A...517A..32P,2010A&A...516A...2M}}
          & 0.15$\pm$0.01  & 12.6$\pm$0.2 & 0.7$\pm$0.1 & 11.8 & 0.59 \\
\hline
\end{tabular}
\end{minipage}
\end{table}

The target stars were selected from the color-magnitude diagram of each cluster
taken from \citetads{2010MNRAS.402.2720P}, \citetads{2002AJ....123..905A},
\citetads{1999MNRAS.304..883P}, and \citetads{2007A&A...467.1065M} for
Berkeley~26, Berkeley~70, NGC~1798, and NGC~2266, respectively
(Fig.~\ref{CMDs}). These data together with the coordinates of each target stars
were obtained from the WEBDA\footnote{{\tt http://www.univie.ac.at/webda}}
database \citepads{1995ASSL..203..127M}. In total, we observed 26 stars: 8 in
Berkeley~70, and 6 in the others. Table~\ref{obsstars} summarizes the coordinates,
magnitudes, exposure times, and signal--to--noise  ratios for each target star.
The global properties of each cluster are listed in
Table~\ref{clusterspropierties}.

Observations were carried out between the nights of 5 and 9 of February 2012
using the Intermediate Dispersion Spectrograph (IDS) mounted at the Cassegrain
focus of the 2.5 m Isaac Newton Telescope (INT) located at the Roque de los
Muchachos Observatory, Spain. We used the R1200R grism centered at 8500\AA~and
the RED+2 CCD, providing a spectral resolution of about 8000. For each target, we
obtained two exposures with the star shifted along the slit in each of them. In
some cases, another star was observed together with the main target because it
was aligned with the slit. Since stellar spectra have reasonable signal-to-noise ratios, we
also analyzed these stars. In all cases, these are main-sequence stars, which
therefore, were used only for the radial velocity analysis.

The data reduction was performed using \textit{imred} and \textit{specred}
packages within IRAF\footnote{Image Reduction and Analysis Facility, IRAF is
distributed by the National Optical Astronomy Observatories, which are operated
by the Association of Universities for Research in Astronomy, Inc., under
cooperative agreement with the National Science Foundation.}. Firstly, each
image was overscan-subtracted and flat-field corrected. Then, since we had acquired
two images of each target with the star shifted along the slit, we subtracted one
from the other, obtaining a positive and a negative spectrum in the same image.
With this procedure, the sky was subtracted from the same physical pixel in which
the star was observed, thus minimizing the effects of pixel-to-pixel sensitivity
variations. A time dependency of course remained, since the two spectra had not
been taken simultaneously. These sky residuals were eliminated in the following
step, in which the spectrum was extracted in the traditional way and the
remaining sky background subtracted from the information on both sides of the
stellar spectra. In the next step, the spectrum was wavelength-calibrated. We then
again subtracted the negative from the positive spectrum (so we effectively added both spectra
because one is negative) to obtain the final spectrum. Finally, each spectrum
was normalized by fitting a polynomial, excluding the strongest lines in the
wavelength range, such as those of the CaT.

\begin{figure}
\centering
\includegraphics[height=14cm,width=\columnwidth]{fig_dcm.eps}
\caption{Color-magnitude diagrams of clusters studied in this paper (small gray
dots). Filled circles are stars assumed to be cluster members based on their
radial velocities and used in the metallicity determination. Filled triangles
are main-sequence stars assumed to be members based on their radial velocity but not
used to determine the metallicity of these clusters. Open circles are
stars rejected as possible members based on their radial velocities.}
\label{CMDs}%
\end{figure}

\section{Radial velocities}\label{sec3}

The radial velocity of each star was calculated using the \textit{fxcor}
task in IRAF, which performs a cross-correlation between the target and template
spectra of known radial velocity. As templates, we used two bright stars,
HD15656 and HD74462, observed in the same nights and with a high
signal--to--noise ratio. The final radial velocity for each target star was
obtained as the average of the velocities obtained for each template, weighted
by the width of the correlation peaks. The measured  radial velocities for each star are
listed in column 8 of Table~\ref{obsstars}. The velocity distribution of each
cluster is shown in Fig.~\ref{Vrdist} and discussed in Sec.~\ref{sec5}.

\begin{figure}
\centering
\includegraphics[height=\columnwidth,width=\columnwidth]{fig_velo.eps}
\caption{Radial velocity distributions of the four clusters studied.}
\label{Vrdist}%
\end{figure}

\section{CaT equivalent widths and metallicity determination}\label{sec4}

The equivalent widths of the CaT lines of the observed RGB stars were used
to determine the metallicity of each cluster. This technique was originally
developed to determine the metallicities of old and metal-poor populations such as
those of globular clusters
\citepads[e.g.][]{1988AJ.....96...92A,1991AJ....101.1329A,1997PASP..109..907R}.
However, more recent investigations have extended its validity to much larger
ranges of ages and metallicities, including extremely metal-poor regimes
\citepads{2010A&A...513A..34S} and the young and metal-rich OC populations
\citepads{2004MNRAS.347..367C,2007AJ....134.1298C}.

Although the procedure used to determine metallicities from the CaT lines
strengths is described in depth by \citetads{2007AJ....134.1298C}, we summarize
here the main steps. Firstly, we determine the equivalent width of each line by
fitting its profile with a Gaussian plus a Lorentzian. This combination provides
the best fit to the line core and wings. Although in our case the spectra had already been
normalized, the position of the continuum was recalculated by performing a linear fitting
to the mean values of each continuum bandpass. The bandpasses used to fit the
line profile and determine the continuum position are those described by
\citetads{2001MNRAS.326..959C}. The CaT index, denoted $\Sigma$Ca, is defined as the sum of the equivalent widths of the three CaT
lines. We list the equivalent widths of the three CaT lines and the $\Sigma$Ca values obtained for each target star in
Table~\ref{obsstars}.

The strengths of the CaT lines do not only depend on the chemical abundance,
but also on the  temperature and gravity of the stellar atmosphere. Since we are
only interested in the abundance dependence, the temperature and gravity
dependences must be removed. To do this, we used the finding that, for a given
chemical abundance, the stars define a sequence in the luminosity-$\Sigma$Ca
plane when the temperature and/or gravity change. It is known that this sequence
is not linear for extremely metal-poor regimes
\citepads{2010A&A...513A..34S} and in the bottom of the RGB
\citepads{2007AJ....134.1298C}. However, we can neglect this trend since our
targets are located in the upper part of the RGB and because for OCs we expect
metallicities of around solar. Traditionally, the difference between the
$V$ magnitude of a star and that of the Horizontal Branch (HB), $V-V_{HB}$,
has been used as a luminosity indicator, which also removes any dependence on
distance and reddening
\citepads[e.g.][]{1991AJ....101.1329A,1997PASP..109..907R,2004MNRAS.347..367C}.
However, in many stellar populations, it is difficult to define accurately the
position of the HB, such as poorly populated OCs. For this reason, other authors
have used the absolute magnitude in either the $V$ or $I$ bandpasses
\citepads[e.g.][]{1991AJ....101.1329A,2004AJ....127..840P,2007AJ....134.1298C}.
We used this second approach in this paper. The CaT metallicity calibration
based on M$_I$ is less sensitive to age. Unfortunately, $I$ magnitudes are unavailable for the stars of one of the clusters, NGC~2266. Therefore, we derive the metallicities of each cluster using both M$_V$ and M$_I$ to minimize the effect of age on the metallicity determination of NGC~2266.

In the M$_{V,I}$-$\Sigma$Ca planes, the sequences can be parametrized in the form
$\Sigma Ca = W'_{V,I} + \beta_{V,I}\times M_{V,I}$. The slopes $\beta_{V,I}$ are
independent of metallicity as demonstrated by
\citetads{2007AJ....134.1298C}, whereas $W'_{V,I}$, the so-called reduced
equivalent width, varies as a function of the chemical abundance. We used the
slopes $\beta_V=-0.677\pm0.004$ \AA~mag$^{-1}$ and $\beta_I=-0.611\pm0.002$
\AA~mag$^{-1}$ calculated by \citetads{2007AJ....134.1298C}. The position of the
observed stars in each cluster in the M$_V$-$\Sigma$Ca and M$_I$-$\Sigma$Ca
planes are shown in Fig.~\ref{MvSigmaCa} and \ref{MiSigmaCa}, respectively. The
$W'_{V,I}$ values obtained for each cluster are listed in
Table~\ref{clustersresults}.

\begin{figure}
\centering
\includegraphics[height=7cm,width=\columnwidth]{fig_mv_sigmaca.eps}
\caption{Observed stars in each cluster in the M$_V$-$\Sigma$Ca plane. Dashed lines are the linear fit to the stars in each cluster.}
\label{MvSigmaCa}%
\end{figure}

\begin{figure}
\centering
\includegraphics[height=7cm,width=\columnwidth]{fig_mi_sigmaca.eps}
\caption{Same as Fig.~\ref{MvSigmaCa} but for M$_I$-$\Sigma$Ca plane. }
\label{MiSigmaCa}%
\end{figure}

Relationships between $W'_{V,I}$ and metallicity on different metallicity
scales, including the widely used one from \citetads{1997A&AS..121...95C} were
derived by \citetads{2007AJ....134.1298C}. \citetads{2009A&A...508..695C} updated their metallicity scale based on more
homogeneous and higher resolution spectra of a larger sample of bright giants in
24 Galactic globular clusters. They provided a detailed comparison between their
old and new scales. The existence of this new metallicity scale for globular
clusters and the new determinations of metallicities available for some open clusters used
by \citetads{2007AJ....134.1298C} to derive their calibrations, motivated
\citetads{2011AJ....142...61C} to obtain new relationships between $W'_{V,I}$
and metallicity.
We used these updated relations to obtain the metallicity, on the
\citetads{2009A&A...508..695C} scale, of the clusters studied in this paper. The
values obtained for each cluster are listed in Table~\ref{clustersresults} and
discussed below. We have to note that the metallicities derived are
related to the reddenings and distance moduli used as references in each case.
For these reason, the metallicity error bars were estimated by taken into
account the uncertainties in these quantities.

\section{Cluster-by-cluster discussion}\label{sec5}

\begin{table*}
\begin{minipage}[t]{17.5cm}
\caption{Cluster metallicities and radial velocities}
\label{clustersresults}
\centering
\renewcommand{\footnoterule}{}
\begin{tabular}{l c c c c c}
\hline\hline
Cluster   & W'$_V$ & [Fe/H]$_V$ & W'$_I$ & [Fe/H]$_I$ & $\langle V_r \rangle$ \\
          & (mag)   & (dex)     & (mag)  & (dex) & (km s$^{-1}$) \\
\hline
Be~26 & 7.42$\pm$0.20 & -0.43$\pm$0.20 & 6.88$\pm$0.23 & -0.35$\pm$0.17 & 68$\pm$12\\
Be~70 & 8.37$\pm$0.01 & -0.02$\pm$0.18 & 7.67$\pm$0.01 & -0.01$\pm$0.14 & -15$\pm$7\\
NGC~1798 & 8.05$\pm$0.04 & -0.16$\pm$0.10 & 7.41$\pm$0.03 & -0.12$\pm$0.07 & 2$\pm$10\\
NGC~2266 & 7.54$\pm$0.02  & -0.38$\pm$0.07 & --- & --- & -16$\pm$15\\
\hline
\end{tabular}
\end{minipage}
\end{table*}

\subsection{Berkeley~26}

Also known as Biurakan~12, Berkeley~26 is a poorly studied cluster in Monoceros.
To our knowledge, there are only three studies of this cluster in the
literature. \citetads{2008PASJ...60.1267H} and \citetads{2010MNRAS.402.2720P}
obtained the color-magnitude diagram of Berkeley~26 in the BVI photometric
system. \citetads{2008MNRAS.389..285T} analyzed near-infrared JHK photometry
obtained from 2MASS\footnote{2MASS (Two Mi\-cron All Sky Survey) is a
joint project of the University of Massachusetts and the Infrared Processing and
Analysis Center/California Institute of Technology, funded by the National
Aeronautics and Space Administration and the National  Science Foundation.}. All
of these authors derived the age, reddening, distance, and metallicity by means of
isochrone fitting. There are some discrepancies between the values obtained by each of
them. There is good agreement, within the uncertainties, in the values derived
by \citetads{2008PASJ...60.1267H} and \citetads{2010MNRAS.402.2720P}, but not
between them and the \citetads{2008MNRAS.389..285T} measurements. For example,
\citetads{2008PASJ...60.1267H} and \citetads{2010MNRAS.402.2720P} derived an age
of 4.5 Gyr and  4 Gyr, respectively. In contrast, \citetads{2008MNRAS.389..285T}
obtained a much younger age of 0.6 Gyr. A reddening of E(B-V)=0.61 was derived
by \citetads{2008PASJ...60.1267H}, which is between the values derived by
\citetads{2008MNRAS.389..285T} of E(B-V)=0.5, and \citetads{2010MNRAS.402.2720P}
of E(B-V)=0.75. In the case of the distance modulus, the obtained values range
from (m-M)$_0$=12.6 \citepads{2008MNRAS.389..285T} to 14.45
\citepads{2008PASJ...60.1267H}. \citetads{2010MNRAS.402.2720P} derived a value
between both of (m-M)$_0$=13.17. As was noted by
\citetads{2010MNRAS.402.2720P},  the poor quality of the data used by \citetads{2008MNRAS.389..285T}, mainly owing to the techniques that they used to decontaminate the foreground
contribution and select the cluster members, may explain the discrepancies with other investigations in the literature. For this reason, we used in our
work the distance modulus and reddening obtained by averaging the values derived
by \citetads{2008PASJ...60.1267H} and \citetads{2010MNRAS.402.2720P} listed in
Table~\ref{clusterspropierties}.

A total of six stars were observed along the Berkeley~26 line of sight. The
spectrum of the star 1650 contains clear molecular bands. Moreover, this star has a larger
radial velocity than the other observed stars. For these reasons, we use
this star in neither the radial velocity nor the metallicity determinations.
The star 1038 is not a RGB and was used only for the radial velocity
determination. Therefore, excluding star 1650, we derived a mean radial velocity
for this cluster of $\langle V_r \rangle$=68$\pm$12  km s$^{-1}$

Only \citetads{2010MNRAS.402.2720P} provided an estimation of the Berkeley~26
metallicity of [Fe/H]=-0.70 dex. \citetads{2008PASJ...60.1267H} assumed a
metallicity of [Fe/H]$\sim$-0.35 for this cluster. We derived
[Fe/H]=-0.43$\pm$0.20 and  -0.35$\pm$0.17 dex using M$_V$ and M$_I$,
respectively. These values are in relatively good agreement with the photometric
estimations taken into account the uncertainties. As explained above,
\citetads{2008MNRAS.389..285T} derived a different distance modulus and reddening
for Berkeley~26 than those used here. When these values were used in the
metallicity determination, we obtained [Fe/H]=+0.09 and +0.06 in the $V$ and $I$
bandpasses, respectively.

\subsection{Berkeley~70}

Berkeley~70 is a rich OC, which probably belongs to the Perseus arm
\citepads{2008A&A...482..777C}. It shows a highly populated main-sequence and a
well-developed RGB. A significant amount of blue stranglers can also be
observed. To our knowledge, only two studies have determined the properties of
this cluster, using both $UBV$ photometry
\citepads{2002AJ....123..905A,2004PASJ...56..295H}. The reddenings and distances
determined by them are the same within the uncertainties. There is a slight
difference in the ages obtained, although they are consistent with the uncertainties.
While \citetads{2002AJ....123..905A} determined an age of 4 Gyr,
\citetads{2004PASJ...56..295H} obtained 2.8 Gyr.

The eight observed stars have  similar radial velocities so we conclude that all
of them are members of Berkeley~70. They define a narrow velocity
distribution (top-right panel of Fig.~\ref{Vrdist}). We derive a mean radial
velocity for this clusters of $\langle V_r \rangle$=-15$\pm$7  km s$^{-1}$.

For Berkeley~70, we derived a metallicity of [Fe/H]=-0.02$\pm$0.18 and
-0.01$\pm$0.14 dex using $V$ and $I$ magnitudes, respectively. To our knowledge,
the only determinations of metallicity for this cluster available in the
literature come from isochrone fitting, and they are all lower than obtained here. \citetads{2002AJ....123..905A} determined a metallicity
of -0.32 dex, while \citetads{2004PASJ...56..295H} obtained a slightly more
metal-poor value of -0.48 dex. This discrepancy may be explained by the
large uncertainties involved in the metallicity estimation from isochrone fitting owing to the limitation of the method itself and to the sparse color magnitude-diagram used as a reference.

\subsection{NGC~1798}

NGC~1798, also denoted as Berkeley~16, is another poorly studied OC located in
the Auriga constellation. Its color-magnitude diagram shows a well-populated
main-sequence and red clump, and a poorly defined RGB. The reddening towards
the NGC~1798 is relatively high, E(B-V)=0.51$\pm$0.04
\citepads{1999MNRAS.304..883P,2000A&A...359..347D,2002NewA....7..553T,
2004A&A...414..163S,2007A&A...467.1065M}, which is in part explained by its
proximity to the Galactic plane. An age of 1.4$\pm$0.2 Gyr has been estimated by
most of the previous investigations in this cluster
\citepads[e.g.][]{1994AJ....108.1773J,1999MNRAS.304..883P,2002A&A...388..158L,
2004A&A...414..163S,2007A&A...467.1065M}. These works have also obtained a similar
distance modulus for this cluster of (m-M)$_0$=13.1$\pm$0.2.

The six stars observed in the NGC~1798 field seem to be members of this cluster
based on their radial velocity. The main-sequence star 608 was observed because it
was aligned with the slit in the observations of star 13 but it was only used for the radial velocity analysis. We determined a mean radial velocity of
$\langle V_r \rangle$=2$\pm$10 km s$^{-1}$ for this cluster.

For the five RGB stars observed in NGC~1798, we inferred metallicities of [Fe/H]=-0.16$\pm$0.10 and
-0.12$\pm$0.07 dex using $V$ and $I$ magnitudes, respectively. There are a
handful of metallicity estimations for NGC~1798 in the literature. All of them
assign a metallicity lower than that obtained here, [Fe/H]=-0.46$\pm$0.10 dex
\citepads{1999MNRAS.304..883P,2001NewA....6..293T,2003NewA....8..737T,
2010A&A...517A..32P}.  As in the case of Berkeley~70, this disagreement may be explained by the
large uncertainties in the metallicities derived from isochrone fitting.

\subsection{NGC~2266}

NGC~2266, also known as Melotte~50, is a rich and well-condensed system. It is
the youngest cluster of our sample with an age similar to that of the Hyades,
\citepads[0.7$\pm$0.1 Gyr;
][]{1991AcA....41..191K,2004A&A...414..163S,2010A&A...516A...2M,
2012A&A...539A.125D}, although some authors have obtained a slightly older age
\citepads[1.2$\pm$0.1 Gyr; ][]{2007A&A...467.1065M,2008BaltA..17...51M}.  Its
color-magnitude diagram has a populated main-sequence, especially in the turn-off
area, and a well-defined red-clump located at $V\sim$13.7 with a few stars above
it (Fig.~\ref{CMDs}). The contamination of the color-magnitude diagram seems
to be almost negligible. We used the reddening and distance derived
by \citetads{2012A&A...539A.125D} (see Table~\ref{clusterspropierties}),
which are similar, within the uncertainties to others available in the
literature
\citepads[e.g.][]{2002A&A...388..158L,2006MNRAS.371.1641P,2008BaltA..17...51M,
2010A&A...516A...2M}.

Six stars were observed in this cluster. However, it seems that only four of
them, including the main-sequence object 75, are true members of the cluster.
The other two, stars 34 and 67, have a radial velocity differing from the average of the other stars much more than the expected uncertainties, hence may not be members of this cluster. Excluding
these two stars, we obtained a mean radial velocity of $\langle V_r \rangle$=-16$\pm$15  km
s$^{-1}$ for NGC~2266.

Since one of the four stars confirmed as a NGC~2266 member is a main-sequence
object, its metallicity determination is based on only three stars. Moreover,
$I$ magnitudes are unavailable for these stars. We therefore, obtained a
metallicity of [Fe/H]=-0.38$\pm$0.07 dex derived only from $V$ magnitudes. This
value is in very good agreement with the following estimates available in the literature:
-0.26$\pm$0.20 dex \citepads{2010A&A...517A..32P}; -0.68
\citepads{2008BaltA..17...51M}; -0.26$\pm$-0.18 \citepads{1991AcA....41..191K};
and -0.35 \citepads{2010A&A...516A...2M,2012A&A...539A.125D}. We note, however, that
\citetads{2004A&A...414..163S} assigned a more metal-rich metallicity of
0.00$\pm$0.20 dex.
\begin{figure}
\centering
\includegraphics[height=10cm,width=\columnwidth]{fig_trends.eps}
\caption{Run of [Fe/H] with R$_{gc}$ (top), age (middle), and $|z|$ (bottom) of
the OCs in the \citetads{2011A&A...535A..30C} compilation. Dark and light gray
points are clusters in the inner 10 kpc and from there off, respectively. Black
points are the clusters studied in this paper.}
\label{trends}%
\end{figure}

\section{Comparison with trends in the Galactic disk}\label{sec6}

As mentioned earlier, OCs are fundamental test particles in the study of the
formation and evolution of the Galactic disk. In particular, they are used to
trace the existence of trends in metallicity with Galactocentric distance,
Rgc, the vertical distance to the Galactic plane, z, and age, providing
fundamental constraints on chemical evolution models. All the investigations
based on OCs agree that the iron content decreases almost linearly
with increasing radius
\citepads[e.g.][]{2002AJ....124.2693F,2010A&A...511A..56P}. However,
investigations based on samples containing clusters at larger distances
 found that the metallicity decreases as a function of
increasing radius to R$_{gc}\sim$12.5 kpc and appears to flatten from there
outwards \citepads[e.g.][]{2005AJ....130..597Y,2008A&A...488..943S,2011MNRAS.412.1265A,
2011A&A...535A..30C}. Moreover, it seems that this gradient has evolved with time being, having been
steeper in the past
\citepads[e.g.][]{2009A&A...494...95M,2011MNRAS.412.1265A,2011A&A...535A..30C}.
However, there are no hints of the existence of an age-metallicity relationship
\citepads[e.g.][]{1995ARA&A..33..381F}. This indicates that the chemical
enrichment of OCs is modulated by their location in the Galaxy and not by the
moment at which they formed \citepads{2011A&A...535A..30C}.

It is helpful to check how the clusters studied in this paper fit the general trends
observed in the Galactic disk. This comparison, of course, should be interpreted with caution
because the results obtained here should be confirmed by analyses based on
high-resolution spectra. In all cases, this comparison is
interesting because of the location of these clusters at relatively large
distances from the Galactic anticenter (R$_{gc}<$10 kpc). To perform this comparison, we used the OC compilation
obtained by \citetads{2011A&A...535A..30C}, which contains data for 89 clusters whose
chemical abundances were derived by investigations based on high-resolution spectra
(R$\ge$20000) available in the literature. We refer the reader to this paper for
details. The run of metallicity with R$_{gc}$, age, and $|z|$ for this sample
have been plotted in the top, middle, and bottom panels of Fig.~\ref{trends},
respectively. The solid lines in the top panel are the linear fits to the sample
assuming a change in the slope at R$_{gc}\sim$12.5 kpc, while the long-dashed
line is the linear fit to the whole sample. Clusters in the inner 10 kpc and
from there outwards have been plotted with dark and light gray circles, respectively, in
the three panels.

The four clusters studied here have been overplotted with
black symbols. Berkeley~26, in spite of the uncertainty in its distance, and NGC~1798 closely follow the general trends
observed in other coeval clusters located at similar Galactocentric distances.
NGC~2266 seems to be slightly younger than other metal-poor systems of similar metallicity situated at the
same distance. Nevertheless, its metal content is akin to other OCs located away from the Galactic plane, such as NGC~2266. In all cases, this comparison must be considered be taken with caution because its
metallicity has been obtained from only three stars. In contrast, Berkeley~70,
which has a more robust metallicity determination based on six stars, appears more metal-rich than
other systems located at $\sim$12.5 kpc. There are a few clusters with similar ages
and metallicities to Berkeley~70, but all of them are situated in the inner 10
kpc at the same distance. Although it is beyond the scope of this paper, one explanation could be that this cluster was formed at a smaller
Galactocentric distance and then migrated to its current location. However, more information than is currently available is needed to confirm or disprove this hypothesis.

\section{Conclusions}\label{sec7}

We have analyzed medium-resolution spectra (R$\sim$8000) in the infrared CaT region ($\sim$8500 \AA) of several stars in four OCs: Berkeley~26, Berkeley~70, NGC~1798, and NGC~2266. To our knowledge, this is the first time that these clusters have been studied spectroscopically. Our
main results can be summarised as follows:

\begin{itemize}
 \item For Berkeley~26, we derived a mean radial velocity of $\langle V_r \rangle$=68$\pm$12 km s$^{-1}$
based on six stars. One of them is a main-sequence star that was not used in the
metallicity determinations. Using $V$ and $I$ magnitudes, we determined a
metallicities of [Fe/H]=-0.43$\pm$0.20 and -0.35$\pm$0.17 dex, respectively.
\item The eight observed stars in Berkeley~70 were confirmed as members based on
their radial velocity. For them, we derived a mean radial velocity of $\langle V_r \rangle$=-15$\pm$7
km s$^{-1}$ and metallicities of [Fe/H]=-0.02$\pm$0.18 and -0.01$\pm$0.14 dex based on $V$ and
$I$ magnitude data, respectively.
\item We derived a mean radial velocity for NGC~1798 of $\langle V_r \rangle$=2$\pm$10 km s$^{-1}$
from six stars, although one object is on the main sequence. For the
RGB stars, we obtained a metallicity of [Fe/H]=-0.16$\pm$0.10 and -0.12$\pm$0.07
dex in $V$ and $I$ bandpasses, respectively.
\item In the case of NGC~2266, its mean radial velocity, $\langle V_r \rangle$=-16$\pm$15 km s$^{-1}$
was obtained for four objects. A metallicity of [Fe/H]=-0.38$\pm$0.07 dex was
derived from the $V$ magnitudes of the three RGB stars observed in this cluster.
\end{itemize}

Finally, we investigated how the four analyzed clusters fit the trends
defined by other well-studied OCs. This comparison is motivated by our
clusters being situated at distances where other investigations have observed a
change in the slope of the metallicity gradient. In general, Berkeley~26 and
NGC~1798 follow the trends described by other coeval systems situated at the
same distance. In the case of NGC~2266, its height above the Galactic plane can
explain its low metallicity compared to other clusters of the same age.
In contrast, Berkeley~70 seems to be more metal-rich than other coeval clusters
situated at similar distances. We suggest that this cluster may have formed
at relatively small Galactocentric distances and has migrated outwards with time.
Nevertheless, more information is needed to confirm or discard this hypothesis.

\begin{acknowledgements}
R. C acknowledge the support of the Spanish Ministry of Economy and Competitiveness
(Plan Nacional de Investigaci\'on Cient\'{\i}fica, Desarrollo, e Investigaci\'on
Tecnol\'ogica, AYA2010-16717). R. C. also acknowledges the
funds of the Spanish Ministry of Science and Innovation under the Juan de la
Cierva fellowship. This research has made use of the WEBDA database, operated at the Institute for Astronomy of the University of Vienna, and the SIMBAD database,
operated at CDS, Strasbourg, France
\end{acknowledgements}


\section{Introduction}
The Andreev reflection in normal metal-superconductor(NS) junctions has  attracted great interests of researchers in the field of superconductivity in recent years\cite{Blonder,Satoshi, Yang, wliming}. The tunneling spectrum of electrons on a NS junction is sensitive to the superconducting gap of the superconductor(SC), thus providing an important technique to measure the gaps and gap properties of SCs. When an electron tunnels into a SC from a metal there appears a hole reflection to conserve the electric charge, which is called the Andreev reflection, since the incident electron may form a Cooper pair in the SC. A commonly used theoretical method to investigate the Andreev reflection is the Blonder-Tinkham-Klapwijk(BTK) approach\cite{Blonder}, which takes the interface of a NS junction as a $\delta(x)$ potential barrier. This theory has been widely and successfully applied to systems like NS junctions\cite{Tanaka}, ferromagnet-SC junctions\cite{Dong2}, etc.  Wei, Dong and Xing {\it et al} studied the tunneling spectrum in metal-SC-metal(NSN) junctions using the BTK approach\cite{Wei,Dong}. They found that the incident and outgoing currents in both sides do not balance each other\cite{Wei}. Thus they claimed earlier that the BTK approach were not suitable for the Andreev reflection in NSN junctions and treated them in the Landauer-Buttiker formalism\cite{Buttiker}. Later they proposed a mechanism to balance the currents by adjusting the chemical potential of the SC inside the NSN junctions\cite{Dong}. This mechanism seems reasonable since the unbalance of currents will change the charge density and thus the Fermi surface of the SC. In a similar idea we proposed recently another mechanism that the interface of the NS junctions is charged by the unbalanced currents\cite{wliming}.

Although these mechanisms successfully balanced the incident current and the outgoing current, however, they mis-considered the phase between the incident electron and hole. We found that different phases just balance the incident current and the ougoing current. After averaging over the phases the incident and outgoing currents coincide each other through the whole range of bias voltage. The charge density, the Fermi surface, and the interface charge state won't change at all. This provides us a more natural and physical mechanism for the theory of multiple NS junctions. In the present paper we study the superlattice of multiple NS junctions using this new mechanism.

\section{formalism}
 A NS junction superlattice is shown in Fig.1, where metals and SCs with width $L_N$ and $L_S$ are sandwiched periodically but the outmost layers are metals working as two leads.  A thin insulating interface exists in every junction to be treated as a $\delta(x)$ potential barrier in the BTK approach. An electron and a hole are incident into the superlattice under a bias voltage on both sides. There is no incident electron in the right metal, thus $b_{n+1}=0$ and in the left metal no incident hole $d_1=0$. Due to the potential barriers the electron and the hole in one layer are partially reflected and partially transmitted into the next layer. In the SC only a pseudo-particle can propagate above the energy gap $\Delta$. When the bias voltage is smaller than $\Delta$ only Cooper pairs can do in the SC. This causes a hole reflected from the interface, which is called the Andreev reflection in literatures.

\begin{figure}
\includegraphics[width=6cm,height=4cm]{nsnsn.eps}
\includegraphics[width=8cm,height=5cm]{bands.eps}
\caption{(a) a NS junction superlattice with an incident electron and an incident hole under a bias voltage. The metal  and the SC are colored by blue and yellow, respectively. An electron is incident from the left metal and partially tunnels into the adjacent SC and is partially reflected. In the right side there is no incident electron, $b_{n+1}=0$ and in the left side $d_1=0$. (b) The fermi surfaces and energy bands of metals and superconductors. Except for the two metal leads the metals and SCs have 0 bias voltage. }
\end{figure}

The tunneling process of an electron and a hole in the NS junction superlattice is described by the Bogoliubov-de Gennes (BdG) equation\cite{BdG}
\begin{align}
i\hbar {\partial \over \partial t} \begin{pmatrix} u({\bf r},t)\\v({\bf r},t)\end{pmatrix} = \begin{pmatrix} H_0 &\Delta \\ \Delta^* & -H_0\end{pmatrix} \begin{pmatrix} u({\bf r},t)\\v({\bf r},t)\end{pmatrix}\\
H_0 = -{\hbar^2\over 2m} \nabla^2 + V({\bf r}) -\mu
\end{align}
where $V({\bf r})= U\delta(x)$ on the interfaces.

When the bias voltage satisfies $eV > \Delta$ an incident particle in metals tunnels into the SC through the interfaces and propagate as a quasi-particle inside the SC. It is then partially reflected by the following interface and partially tunnels into the following metal.  The wave functions of the electron and the hole inside the $i$th metal layer and the $i$th SC layer are written as
\begin{align}
\psi^i_{N}= \left(
  \begin{array}{c}  a_i e^{ik_+ x}+b_i e^{-ik_+ x}\\   c_i e^{i k_- x}+ d_i e^{-i k_- x}\\
  \end{array}\right) \label{wf1}
\\
\psi^i_{SC}=\biggl[(e_i e^{iq_+ x}+ f_i e^{-iq_+ x})
\left(
  \begin{array}{c}   u_{+}\\  v_{+}\\  \end{array} \right)   \nonumber\\
  + (g_i e^{iq_- x} + h_i e^{-iq_- x})
\left(   \begin{array}{c}     u_{-}\\    v_{-}\\   \end{array} \right) \biggr]
\label{wf2}
\end{align}
where $\left(
  \begin{array}{c}   u_{+}\\  v_{+}\\  \end{array} \right)$ and $\left(
  \begin{array}{c}   u_{-}\\  v_{-}\\  \end{array} \right)$ are the quasi-eletron ($\epsilon_q>0$) and quasi-hole ($\epsilon_q<0$) wave functions with the same energy $E$ in the SC layers, respectively. They are given by
\begin{align}
u = \sqrt{{1\over 2} + {\epsilon_{q} \over 2E}}, v = {|\Delta|\over \Delta}\sqrt{{1\over 2} - {\epsilon_{q}\over 2E}}\\
\epsilon_q ={\hbar^2 q^2\over {2m}}-\mu=\pm \sqrt{E^2 - |\Delta|^2}
\end{align}
The wave vectors are determined by
 \begin{align}
 E\equiv eV=eV+{\hbar^2 k_{+,1}^2\over {2m}} - \mu =  \mu - {\hbar^2 k_{-,1}^2 \over {2m}}-eV\\
 \quad k_{+,1}=k_F = \sqrt{2m\mu/\hbar^2}\\ k_{-,1}/k_F=\sqrt{1- {2E/\mu}}\\
\quad k_{\pm, i}/k_F = \sqrt{1\pm {E/\mu}}, 1<i<n+1\\
\quad q_\pm/k_F = \sqrt{1\pm \sqrt{E^2-\Delta_q^2}/\mu}\\
2eV={\hbar^2 k_{+,n+1}^2\over {2m}} - \mu ,\quad 0 =  \mu - {\hbar^2 k_{-,n+1}^2 \over {2m}}\\
\quad k_{+,n+1}/k_F=\sqrt{1+{2E/\mu}}, \quad k_{-,n+1}= k_F
\end{align}
When $eV<\Delta$, $ q_\pm$ become complex numbers, so that Cooper pairs appear and the wave functions in the SC become damping traveling waves.
In this case one has
 \begin{align}
\epsilon_q  = \pm i\sqrt{ |\Delta|^2- E^2}\\
\begin{pmatrix}u\\v\end{pmatrix}_\pm ={1\over \sqrt 2}\begin{pmatrix}1 \\{E- \epsilon_q\over \Delta}\end{pmatrix}
\end{align}
The coefficients $a,b,c,d,e,f,g,h$ in (\ref{wf1}) to (\ref{wf2}) are obtained from the boundary conditions of wave functions on the interfaces, where the wave functions are continuous but due to the $\delta(x)$ barriers
the first derivatives of them are not, i.e.,
\begin{align}
\psi_{1}(X_i)=\psi_{2}(X_i)\\
\psi'_{1}(X_i)-\psi'_{2}(X_i)+Zk_F\psi_{1}(X_i)=0
\end{align}
where 1 and 2 denote the adjacent metal and SC layers, and $Z=\frac{2mU}{\hbar^{2}k_F}$ and $X_i$ is the position of the interface. These boundary conditions give the following matrix equations
\begin{align}
(a_i, b_i, c_i, d_i)^T = M_1(X_i)(e_i,f_i,g_i,h_i)^T\\
(e_i,f_i,g_i,h_i)^T=M_2(Y_{i}) (a_{i+1}, b_{i+1}, c_{i+1}, d_{i+1})^T
\end{align}
for $i = 1,2, ..., n$, where $M_1(X_i)$ and $M_2(Y_{i})$ are called the transfer matrices at interface $X_i$ between the metal and SC layers and interface $Y_i$ between the SC and metal layers. Using this matrices we find the relation between the coefficients of the incident and outgoing particles in the following form
\begin{align}
\begin{pmatrix}a_1\\ b_1\\c_1\\ d_1\end{pmatrix} = \prod_i M_1(X_i) M_2(Y_i) \begin{pmatrix}a_{n+1}\\ b_{n+1}\\ c_{n+1}\\ d_{n+1}\end{pmatrix}
\end{align}
Assigning $a_1=1, d_1=0, b_{n+1}=0$, and $c_{n+1}=exp(i\phi)$ for an incident electron and an incident hole with a phase difference $\phi$,  other coefficients can be obtained relative to $\phi$ from this matrix equation. These coefficients determine the the incident and outgoing current densities and the differential electric conductance\cite{wliming}.
\begin{align}
J_1 \sim[k_+(1-|b_1(eV)|^2)+k_- |c_1(eV)|^2]\\
J_{n+1} \sim[k_-(1-|d_{n+1}(eV)|^2)+k_+|a_{n+1}(eV)|^2]\\
{dI_1(eV)\over dV}\sim 1-|b_1(eV)|^2+|c_1(eV)|^2\\
{dI_2(eV)\over dV}\sim 1-|d_{n+1}(eV)|^2+|a_{n+1}(eV)|^2
\end{align}

In general the incident current and the outgoing current do not match each other in a NS junction superlattice even in NSN junctions as pointed by Xing et al\cite{Dong}. The key idea of this work is that the currents depend on the phase difference $\phi$. After averaging over the phase the incident current just balances the outgoing current and so does the differential conductance.

\section{Results of computation}
The phase difference $\phi$ has strong effect on the currents on the NS superlattice, as shown in Fig.2 for a 2-NS-junction lattice. It is found that the probabilities of the reflected electron $|b_1|^2$, the outgoing electron $|a_{n+1}|^2$, the Andreev reflection $|c_1|^2$ and the hole reflection $|d_{n+1}|^2$ vary periodically with the the phase difference $\phi$, and so do the incident and outgoing currents. The total probability  is always conserved at different phase differences $|a_{n+1}|^2+|b_1|^2+|c_1|^2+|d_{n+1}|^2=2$. However, the incident and outgoing currents do not match each other at different phase differences(except two points). This result has never been reported in literatures and no one noticed that these two currents oscillate about the phase difference. Xing et al simply averaged the two currents from only in-phase incidences of an electron and a hole to balance the two currents\cite{Dong}. One of the present authors assumed in an earlier paper a charge accumulation on the interfaces to conserve the currents\cite{wliming}. They all neglected the effect of the phase difference on the currents. An important result of the present computation is that the incident and outgoing currents have almost the same average value for different $\phi$. This indicates that a great number of electrons and holes with random phase differences incident on the two leads of the NS junction superlattice will balance the two currents naturally. This provides a new but final mechanism for the balance of currents on the NS junction superlattice.

Using the above mechanism the bias voltage is scanned for the currents as shown in Fig.3. It is seen that the currents from in-phase incidences of an electron and a hole obviously deviate from each other above the superconducting  gap as reported by earlier works\cite{Wei, wliming}. After averaging over the phase difference the currents meet each other perfectly except small deviations at the three peaks of the currents. Therefore, this averaging process successfully recovers the balance of currents in a NS superlattice. The differential conductance for different interface barriers under this phase-averaging mechanism is shown in Fig.4(a). For smaller and smaller barrier potentials the differential conductance approaches the saturating value, i.e. fully conducting NS junctions. As the barrier increases the differential conductance shows stronger and stronger tunneling effect.  The first resonance peak takes place slightly higher than the superconducting gap. This provides a technique for measuring superconducting gaps. The other two resonance peaks occur at about $eV = 1.35\Delta, 1.8\Delta$. These resonances need experimental verification. For a comparison to Dong's work\cite{Dong}, where only in-phase incidence of electrons and holes is considered, Fig.4(b) shows the differential conductance using the phase averaging method and Dong's parameters.  For weak barriers the present result is quite different from Dong's result, but for strong barriers the present result approaches to Dong's result.

\begin{figure}
\includegraphics[width=4.cm,height=3.8cm]{scanphase1.eps}
\includegraphics[width=4.cm,height=3.8cm]{scanphase2.eps}
\caption{(a) Probabilities of the reflected electron, the outgoing electron, the Andreev reflection, and the hole reflection in NSN junctions. The total probability of these terms is exactly equal to 2. (b) The incident current and the outgoing current $J_1$, $J_{n+1}$ at different phase  $\phi$. Parameters in the computation are set to be $\mu=0.5, \Delta = 0.01\mu, Z = 1, k_F L=1000, eV = 2\Delta$. }
\end{figure}

\begin{figure}
\includegraphics[width=4.cm,height=3.8cm]{j12.eps}
\includegraphics[width=4.cm,height=3.8cm]{j12a.eps}
\caption{(a) Incident and outgoing currents in NSN junctions with an in-phase incidence of an electron and a hole. (b) The currents averaged over the phase $\phi$. Parameters see Fig.2.}
\end{figure}

\begin{figure}
\includegraphics[width=4.cm,height=3.8cm]{di12z.eps}
\includegraphics[width=4.cm,height=3.8cm]{xing.eps}
\caption{Differential conductance for different interface barriers. (a) parameters see Fig.2; (b) A comparison to Fig.6 of Dong's work\cite{Dong} with parameters $\mu=0.5, \Delta = 0.001\mu, k_F L=5000$. }
\end{figure}

Next superlattices of 2,4,8 and 16 NS junctions are studied in the above phase averaging mechanism. The differential conductances of them with barriers $Z=0.5$ and $Z=4$ are shown in Fig.5(a,b). It is seen that in these two cases the NS junction superlattices have some sharp zero conductance resonances. Especially those of the superlattice of 16 junctions has a full width at half maximum of $0.012\Delta$, as seen in the inset of (b), which is equivalent to a voltage of about $0.5mV$. This effect may find applications such as  sensitive millivolt electronic switches, where a weak voltage increase induces a strong current response.

\begin{figure}\includegraphics[width=4.cm,height=3.8cm]{248junctionsZ05.eps}\includegraphics[width=4.cm,height=3.8cm]{248junctionsZ40.eps}
\caption{Differential conductances of superlattices of 2,4,8, and 16 NS junctions with (a) Z = 0.5 and (b) Z=4.0. Other parameters see Fig.2.}
\end{figure}

Finally it should be pointed out that the first resonance peak of the differential conductance does not occur exictly at voltages corresponding to the energy gap of the SC at small values of $k_fL$. As shown in Fig.6 until $k_fL = 1200 $ the first resonance peak moves to the  gap position. This indicates that the measuring density of states using the technique of tunneling spectrum of metal-SC junctions requires a thickness of the SC layer above the minimum value.
\begin{figure}
\includegraphics[width=4.cm,height=3.8cm]{kL_2J.eps}
\caption{Differential conductances of multiple NS junction lattices with different width $k_f L$ of superconducting layers. Other parameters see Fig.2.}
\end{figure}

In summary, the tunneling spectrum of an electron and a hole on a superlattice of NS junctions is computed using the BTK approach. It shows a more abundant structure compared to that on a single NS junction, such as the sharp resonance above the gap. In particular, the sharper the resonance is the more layers the superlattice has. This effect may find applications in the electronic industry in the future. We find  for the first time a  mechanism to balance the incident and outgoing currents on the superlattice  by averaging over the phase between the incident electron and the incident hole. Compared to other mechanisms such as the adjusting of the Fermi surface and the charge accumulation the present mechanism is much more natural and physical.

\section{acknowledgment}
 This work was supported by the National Natural Science Foundation of
China (Grant No. 10874049),  the State Key Program for Basic Research of China (No.
2007CB925204).