text:
\section{Introduction}

Spacelike hypersurfaces with constant mean curvature in Lorentz
manifolds have been object of great interest in recent years,
both from physical and mathematical points of view. In~\cite{ABC:03},
the authors studied the uniqueness of spacelike hypersurfaces  with CMC in
generalized Robertson-Walker (GRW) spacetimes, namely, Lorentz warped
products with 1-dimensional negative  definite base and Riemannian fiber.
They proved that in a GRW spacetime obeying the timelike convergence condition
(i.e, the Ricci curvature is non-negative on timelike directions), every compact
spacelike hypersurface with CMC must be umbilical. Recently, Al\'{i}as and  Montiel obtained,
in~\cite{AM:01}, a more general condition on the warping function $f$ that is
sufficient in order to guarantee uniqueness. More precisely, they proved the following

\begin{theorem} Let $f:I\rightarrow\mathbb R$ be a positive smooth function defined on
an open interval, such that $ff^{''}-(f^{'})^{2}\leq 0$, that is,
such that $-\log f$ is convex. Then, the only compact spacelike
hypersurfaces immersed into a generalized Robertson-Walker spacetime
$I\times_fF^{n}$ and having constant mean curvature are the slices
$\{t\}\times F$, for a {\rm(}necessarily compact{\rm)} Riemannian
manifold $F$.
\end{theorem}

Stability questions concerning CMC, compact hypersurfaces in
Riemannian space forms began with Barbosa and do Carmo
in~\cite{BdC:84}, and Barbosa, Do Carmo and Eschenburg
in~\cite{BdCE:88}. In the former paper, they introduced the notion
of stability and proved that spheres are the only stable critical
points for the area functional, for volume-preserving variations. In
the setting of spacelike hypersurfaces in Lorentz manifolds, Barbosa
and Oliker proved in~\cite{Barbosa:93} that CMC spacelike
hypersurfaces are critical points of volume-preserving variations.
Moreover, by computing the second variation formula they showed that
CMC embedded spheres in the de Sitter space $S_1^{n+1}$ maximize the
area functional for such variations. In this paper, we give a
characterization of {\em strongly stable}, CMC spacelike
hypersurfaces in GRW spacetimes, the essential tool for the proof
being a formula for the Laplacian of a new support function. More
precisely, it is our purpose to show the following

\begin{theorem}
Let $\overline M^{n+1}=I\times_{\phi}F^n$ be a generalized
Robertson-Walker spacetime, and $x:M^n\rightarrow\overline M^{n+1}$
be a closed spacelike hypersurface of $\overline M^{n+1}$, having
constant mean curvature $H$. If the warping function $\phi$
satisfies $\phi''\geq\max\{H\phi',0\}$ and $M^n$ is strongly
stable, then $M^{n}$ is either minimal or a spacelike slice
$M_{t_0}=\{t_0\}\times F$, for some $t_0\in I$.
\end{theorem}

\section{Stable spacelike hypersurfaces}

In what follows, $\overline M^{n+1}$ denotes an orientable, time-oriented Lorentz manifold
with Lorentz metric $\overline g=\langle\,\,,\,\,\rangle$ and semi-Riemannian connection
$\overline\nabla$. If $x:M^n\rightarrow\overline M^{n+1}$ is a spacelike hypersurface of
$\overline M^{n+1}$, then $M^n$ is automatically orientable (\cite{O'Neill:83}, p. 189),
and one can choose a globally defined unit normal vector field $N$ on $M^n$ having the
same time-orientation of $V$, that is, such that
$$\langle V,N\rangle<0$$
on $M$. One says that such an $N$ {\em points to the future}.

A {\em variation} of $x$ is a smooth map
$$X:M^n\times(-\epsilon,\epsilon)\rightarrow\overline M^{n+1}$$
satisfying the following conditions:
\begin{enumerate}
\item[(1)] For $t\in(-\epsilon,\epsilon)$, the map $X_t:M^n\rightarrow\overline M^{n+1}$
given by $X_t(p)=X(t,p)$ is a spacelike immersion such that $X_0=x$.
\item[(2)] $X_t\big|_{\partial M}=x\big|_{\partial M}$, for all $t\in(-\epsilon,\epsilon)$.
\end{enumerate}

The {\em variational field} associated to the variation $X$ is the vector field
$\frac{\partial X}{\partial t}$. Letting
$f=-\langle\frac{\partial X}{\partial t},N\rangle$, we get
$$\frac{\partial X}{\partial t}\Big|_M=fN+\left(\frac{\partial X}{\partial t}\right)^T,$$
where $T$ stands for tangential components. The {\em balance of volume} of the variation $X$
is the function $\mathcal V:(-\epsilon,\epsilon)\rightarrow\mathbb R$ given by
$$\mathcal V(t)=\int_{M\times[0,t]}X^*(d\overline M),$$
where $d\overline M$ denotes the volume element of $\overline M$.

The {\em area functional} $\mathcal A:(-\epsilon,\epsilon)\rightarrow\mathbb R$ associated
to the variation $X$ is given by
$$\mathcal A(t)=\int_MdM_t,$$
where $dM_t$ denotes the volume element of the metric induced in $M$ by $X_t$. Note that
$dM_0=dM$ and $\mathcal A(0)=\mathcal A$, the volume of $M$.
The following lemma is classical:

\begin{lemma}\label{lemma:first variation}
Let $\overline M^{n+1}$ be a time-oriented Lorentz manifold and
$x:M^n\rightarrow\overline M^{n+1}$ a spacelike closed hypersurface having mean curvature $H$.
If $X:M^n\times(-\epsilon,\epsilon)\rightarrow\overline M^{n+1}$ is a variation of $x$, then
$$\frac{d\mathcal V}{dt}\Big|_{t=0}=\int_MfdM,\ \   \ \ \frac{d\mathcal A}{dt}\Big|_{t=0}=\int_MnHfdM.$$
\end{lemma}

Set $H_0=\frac{1}{\mathcal A}\int_MdM$ and
$\mathcal J:(-\epsilon,\epsilon)\rightarrow\mathbb R$ given by
$$\mathcal J(t)=\mathcal A(t)-nH_0\mathcal V(t).$$
$\mathcal J$ is called the {\em Jacobi functional} associated to the variation,
and it is a well known result~\cite{BdCE:88} that $x$ has constant mean curvature $H_0$
if and only if $\mathcal J'(0)=0$ for all variations $X$ of $x$.

We wish to study here immersions $x:M^n\rightarrow\overline M^{n+1}$ that
maximize $\mathcal J$ for all variations $X$. Since $x$ must be a critical point of
$\mathcal J$, it thus follows from the above discussion that $x$ must have constant mean
curvature. Therefore, in order to examine whether or not some critical immersion $x$ is
actually a maximum for $\mathcal J$, one certainly needs to study the second variation
$\mathcal J''(0)$. We start with the following

\begin{proposition} Let $x:M^n\rightarrow\overline M^{n+1}$ be a closed spacelike
hypersurface of the time-oriented Lorentz manifold $\overline M^{n+1}$, and
$X:M^n\times(-\epsilon,\epsilon)\rightarrow\overline M^{n+1}$ be a variation of $x$. Then,
\begin{equation}\label{eq:fundamental relation}
n\frac{\partial H}{\partial t}=\Delta f-\left\{\overline{Ric}(N,N)+|A|^2\right\}f-n\langle\left(\frac{\partial X}{\partial t}\right)^T,\nabla H\rangle.
\end{equation}
\end{proposition}

Although the above proposition is known to be true, we believe there is a lack,
in the literature, of a clear proof of it in this degree of generality, so we present
a simple proof here.

\begin{proof} Let $p\in M$ and $\{e_k\}$ be a moving frame on a neighborhood $U\subset M$ of
$p$, geodesic at $p$ and diagonalizing $A$ at $p$, with $Ae_k=\lambda_ke_k$ for
$1\leq k\leq n$. Extend $N$ and the $e_k's$ to a neighborhood of $p$ in $\overline M$, so
that $\langle N,e_k\rangle=0$ and $(\overline\nabla_Ne_k)(p)=0$. Then
\begin{eqnarray*}
n\frac{\partial H}{\partial t}&=&-{\rm tr}\left(\frac{\partial A}{\partial t}\right)=-\sum_k\langle\frac{\partial A}{\partial t}e_k,e_k\rangle=-\sum_k\langle\left(\overline\nabla_{\frac{\partial X}{\partial t}}A\right)e_k,e_k\rangle\\
&=&-\sum_k\left\{\langle\overline\nabla_{\frac{\partial X}{\partial t}}Ae_k,e_k\rangle-\langle A\overline\nabla_{\frac{\partial X}{\partial t}}e_k,e_k\rangle\right\}\\
&=&\sum_k\langle\overline\nabla_{\frac{\partial X}{\partial t}}\overline\nabla_{e_k}N,e_k\rangle+\sum_k\langle A\overline\nabla_{e_k}\frac{\partial X}{\partial t},e_k\rangle,
\end{eqnarray*}
where in the last equality we used the fact that $[\frac{\partial X}{\partial t},e_k]=0$.
Letting
$$I=\sum_k\langle\overline\nabla_{\frac{\partial X}{\partial t}}\overline\nabla_{e_k}N,e_k\rangle\ \ \text{and}\ \ II=\sum_k\langle A\overline\nabla_{e_k}\frac{\partial X}{\partial t},e_k\rangle,$$
we have
\begin{eqnarray*}
I&=&\sum_k\left\{\langle\overline\nabla_{\frac{\partial X}{\partial t}}\overline\nabla_{e_k}N-\overline\nabla_{e_k}\overline\nabla_{\frac{\partial X}{\partial t}}N+\overline\nabla_{[e_k,\frac{\partial X}{\partial t}]}N,e_k\rangle+\overline\nabla_{e_k}\overline\nabla_{\frac{\partial X}{\partial t}}N,e_k\rangle\right\}\\
&=&\sum_k\left\{\langle\overline R\left(e_k,\frac{\partial X}{\partial t}\right)N,e_k\rangle+\langle\overline\nabla_{e_k}\overline\nabla_{\frac{\partial X}{\partial t}}N,e_k\rangle\right\}\\
&=&-\overline{Ric}\left(\frac{\partial X}{\partial t},N\right)+\sum_k\langle\overline\nabla_{e_k}\overline\nabla_{\frac{\partial X}{\partial t}}N,e_k\rangle.
\end{eqnarray*}

Since the frame $\{e_k\}$ is geodesic at $p$, it follows that
$$\langle\overline\nabla_{\frac{\partial X}{\partial t}}N,\overline\nabla_{e_k}e_k\rangle=\langle\overline\nabla_{\frac{\partial X}{\partial t}}N,N\rangle\langle\overline\nabla_{e_k}e_k,N\rangle=0$$
at $p$, and hence
\begin{eqnarray*}
\langle\overline\nabla_{e_k}\overline\nabla_{\frac{\partial X}{\partial t}}N,e_k\rangle&=&e_k\langle\overline\nabla_{\frac{\partial X}{\partial t}}N,e_k\rangle=-e_k\langle N,\overline\nabla_{\frac{\partial X}{\partial t}}e_k\rangle=-e_k\langle N,\overline\nabla_{e_k}\frac{\partial X}{\partial t}\rangle\\
&=&-e_ke_k\langle N,\frac{\partial X}{\partial t}\rangle+e_k\langle\overline\nabla_{e_k}N,\frac{\partial X}{\partial t}\rangle\\
&=&e_ke_k(f)+e_k\langle\overline\nabla_{e_k}N,\left(\frac{\partial X}{\partial t}\right)^T\rangle\\
&=&e_ke_k(f)+\langle\overline\nabla_{e_k}\overline\nabla_{e_k}N,\left(\frac{\partial X}{\partial t}\right)^T\rangle-\langle Ae_k,\overline\nabla_{e_k}\left(\frac{\partial X}{\partial t}\right)^T\rangle.
\end{eqnarray*}

For $II$, we have
\begin{eqnarray*}
II&=&\sum_k\langle Ae_k,\overline\nabla_{e_k}\frac{\partial X}{\partial t}\rangle=\sum_k\langle Ae_k,\overline\nabla_{e_k}(fN+\left(\frac{\partial X}{\partial t}\right)^T)\rangle\\
&=&\sum_k\langle Ae_k,f\overline\nabla_{e_k}N\rangle+\sum_k\langle Ae_k,\overline\nabla_{e_k}\left(\frac{\partial X}{\partial t}\right)^T\rangle\\
&=&-f|A|^2+\sum_k\langle Ae_k,\overline\nabla_{e_k}\left(\frac{\partial X}{\partial t}\right)^T\rangle\\
\end{eqnarray*}

Therefore,
\begin{equation}\label{eq:aux_I}
n\frac{\partial H}{\partial t}=-\overline{Ric}\left(\frac{\partial X}{\partial t},N\right)+\Delta f-f|A|^2+\sum_k\langle\overline\nabla_{e_k}\overline\nabla_{e_k}N,\left(\frac{\partial X}{\partial t}\right)^T\rangle.
\end{equation}

Now, letting
$$\frac{\partial X}{\partial t}=\sum_l^n\alpha_le_l+fN$$
and $Ae_k=\sum_jh_{jk}e_j$, one successively gets
\begin{eqnarray*}
\overline{Ric}\left(\frac{\partial X}{\partial t},N\right)&=&\sum_l\alpha_l\overline{Ric}(N,e_l)+f\overline{Ric}(N,N)\\
&=&\sum_{k,l}\alpha_l\langle\overline R(e_k,e_l)e_k,N\rangle+f\overline{Ric}(N,N)\\
\end{eqnarray*}
and, since $(\overline\nabla_Ne_k)(p)=0$,
\begin{eqnarray*}
\langle\overline R(e_k,e_l)e_k,N\rangle_p&=&\langle\overline\nabla_{e_l}\overline\nabla_{e_k}e_k-\overline\nabla_{e_k}\overline\nabla_{e_l}e_k,N\rangle_p\\
&=&e_l\langle\overline\nabla_{e_k}e_k,N\rangle_p-\langle\overline\nabla_{e_k}e_k,\overline\nabla_{e_l}N\rangle_p-e_k\langle\overline\nabla_{e_l}e_k,N\rangle_p\\
&=&-e_l\langle e_k,\overline\nabla_{e_k}N\rangle_p+e_k\langle e_k,\overline\nabla_{e_l}N\rangle_p\\
&=&e_l(h_{kk})-e_k(h_{kl}),
\end{eqnarray*}
so that
\begin{equation}\label{eq:aux_II}
\overline{Ric}\left(\frac{\partial X}{\partial t},N\right)_p=\sum_{k,l}\alpha_le_l(h_{kk})-\sum_{k,l}\alpha_le_k(h_{kl})+f\overline{Ric}(N,N)_p.
\end{equation}

Also,
\begin{eqnarray*}
\alpha_l\langle\overline\nabla_{e_k}\overline\nabla_{e_k}N,e_l\rangle&=&\alpha_l\langle\nabla_{e_k}\overline\nabla_{e_k}N,e_l\rangle=-\alpha_l\sum_j\langle\nabla_{e_k}h_{kj}e_j,e_l\rangle\\
&=&-\alpha_l\sum_j\left\{e_k(h_{kj})\delta_{lj}+h_{kj}\langle\nabla_{e_k}e_j,e_l\rangle\right\}\\
&=&-\alpha_le_k(h_{kl}),
\end{eqnarray*}
and hence
\begin{equation}\label{eq:aux_III}
\sum_k\langle\overline\nabla_{e_k}\overline\nabla_{e_k}N,\left(\frac{\partial X}{\partial t}\right)^T\rangle=-\sum_{k,l}\alpha_le_k(h_{kl}).
\end{equation}

Substituting (\ref{eq:aux_II}) and (\ref{eq:aux_III}) into (\ref{eq:aux_I}), we finally arrive
at
\begin{eqnarray*}
n\frac{\partial H}{\partial t}&=&-\sum_{k,l}\alpha_le_l(h_{kk})-f\overline{Ric}(N,N)_p+\Delta f-f|A|^2\\
&=&-\left(\frac{\partial X}{\partial t}\right)^T(nH)-f\overline{Ric}(N,N)_p+\Delta f-f|A|^2.
\end{eqnarray*}
\end{proof}

\begin{proposition} Let $\overline M^{n+1}$ be a Lorentz manifold and
$x:M^n\rightarrow\overline M^{n+1}$ be a closed spacelike hypersurface having constant
mean curvature $H$. If $X:M^n\times(-\epsilon,\epsilon)\rightarrow\overline M^{n+1}$
is a variation of $x$, then
\begin{equation}\label{eq:second variation of J}
\mathcal J''(0)(f)=\int_Mf\left\{\Delta f-\left(\overline{Ric}(N,N)+|A|^2\right)f\right\}dM.
\end{equation}
\end{proposition}

\begin{proof} In the notations of the above discussion, set $f=f(0)$ and note that $H(0)=H$.
It follows from lemma~\ref{lemma:first variation} that
$$\mathcal J'(t)=\int_Mn\left\{H(t)-H\right\}f(t)dM_t.$$
Therefore, differentiating with respect to $t$ once more
\begin{eqnarray*}
\mathcal J''(0)&=&\int_MnH'(0)f(0)dM_0+\int_Mn\left\{H(0)-H\right\}\frac{d}{dt}f(t)dM_t\Big|_{t=0}\\
&=&\int_MnH'(0)fdM.
\end{eqnarray*}

Taking into account that $H$ is constant, relation (\ref{eq:fundamental relation}) finally
gives formula~\ref{eq:second variation of J}
\end{proof}

It follows from the previous result that $\mathcal J''(0)=\mathcal J''(0)(f)$ depends only
on $f\in C^{\infty}(M)$, for which there exists a variation $X$ of $M^n$ such that
$\left(\frac{\partial X}{\partial t}\right)^{\bot}=fN$. Therefore, the following definition
makes sense:

\begin{definition} Let $\overline M^{n+1}$ be a Lorentz manifold and
$x:M^n\rightarrow\overline M^{n+1}$ be a closed spacelike hypersurface having constant
mean curvature $H$. We say that $x$ is strongly stable if, for every function
$f\in C^{\infty}(M)$ for which there exists a variation $X$ of $M^n$ such that
$\left(\frac{\partial X}{\partial t}\right)^{\bot}=fN$, one has $\mathcal J''(0)(f)\leq 0$.
\end{definition}

\section{Conformal vector fields}

As in the previous section, let $\overline M^{n+1}$ be a Lorentz manifold.
A vector field $V$ on $\overline M^{n+1}$ is said to be {\em conformal} if
\begin{equation}
\mathcal L_V\langle\,\,,\,\,\rangle=2\psi\langle\,\,,\,\,\rangle
\end{equation}
for some function $\psi\in C^{\infty}(\overline M)$, where $\mathcal L$ stands for the Lie
derivative of the Lorentz metric of $\overline M$. The function $\psi$ is called the
{\em conformal factor} of $V$.

Since $\mathcal L_V(X)=[V,X]$ for all $X\in\mathcal X(\overline M)$, it follows from the
tensorial character of $\mathcal L_V$ that $V\in\mathcal X(\overline M)$ is conformal if
and only if
\begin{equation}\label{eq:1.1}
\langle\overline\nabla_XV,Y \rangle+\langle X,\overline\nabla_YV\rangle=2\psi\langle X,Y\rangle,
\end{equation}
for all $X,Y\in\mathcal X(\overline M)$. In particular, $V$ is a Killing vector field
relatively to $\overline g$ if and only if $\psi\equiv 0$.

Any Lorentz manifold $\overline M^{n+1}$, possessing a globally defined, timelike conformal
vector field is said to be a {\em conformally stationary spacetime}.

\begin{proposition}\label{prop:Laplacian of conformal vector field}
Let $\overline M^{n+1}$ be a conformally stationary Lorentz manifold, with conformal vector
field $V$ having conformal factor $\psi:\overline M^{n+1}\rightarrow\mathbb R$. Let also
$x:M^n\rightarrow\overline M^{n+1}$ be a spacelike hypersurface of $\overline M^{n+1}$,
and $N$ a future-pointing, unit normal vector field globally defined on $M^n$.
If $f=\langle V,N\rangle$, then
\begin{equation}\label{eq:Laplacian formula_I}
\Delta f=n\langle V,\nabla H\rangle+f\left\{\overline{Ric}(N,N)+|A|^2\right\}+n\left\{H\psi-N(\psi)\right\},
\end{equation}
where $\overline{Ric}$ denotes the Ricci tensor of $\overline M$, $A$ is the second fundamental
form of $x$ with respect to $N$, $H=-\frac{1}{n}{\rm tr}(A)$ is the mean curvature of $x$ and
$\nabla H$ denotes the gradient of $H$ in the metric of $M$.
\end{proposition}

\begin{proof} Fix $p\in M$ and let $\{e_k\}$ be an orthonormal moving frame on $M$, geodesic
at $p$. Extend the $e_k$ to a neighborhood of $p$ in $\overline M$, so that
$(\overline\nabla_Ne_k)(p)=0$, and let
$$V=\sum_l^n\alpha_le_l-fN.$$
Then
\begin{eqnarray*}
f=\langle N,V\rangle\Rightarrow e_k(f)&=&\langle\overline\nabla_{e_k}N,V\rangle+\langle N,\overline\nabla_{e_k}V\rangle\\
&=&-\langle Ae_k,V\rangle+\langle N,\overline\nabla_{e_k}V\rangle,
\end{eqnarray*}
so that
\begin{eqnarray}\label{eq:I}
\Delta f&=&\sum_ke_k(e_k(f))=-\sum_ke_k\langle Ae_k,V\rangle+\sum_ke_k\langle N,\overline\nabla_{e_k}V\rangle\nonumber\\
&=&-\sum_k\langle\overline\nabla_{e_k}Ae_k,V\rangle-2\sum_k\langle Ae_k,\overline\nabla_{e_k}V\rangle+\sum_k\langle N,\overline\nabla_{e_k}\overline\nabla_{e_k}V\rangle.
\end{eqnarray}

Now, differentiating $Ae_k=\sum_lh_{kl}e_l$ with respect to $e_k$, one gets at $p$
\begin{eqnarray}\label{eq:II}
\sum_k\langle\overline\nabla_{e_k}Ae_k,V\rangle&=&\sum_{k,l}e_k(h_{kl})\langle e_l,V\rangle+\sum_{k,l}h_{kl}\langle\overline\nabla_{e_k}e_l,V\rangle\nonumber\\
&=&\sum_{k,l}\alpha_le_k(h_{kl})-\sum_{k,l}h_{kl}\langle\overline\nabla_{e_k}e_l,N\rangle\langle V,N\rangle\nonumber\\
&=&\sum_{k,l}\alpha_le_k(h_{kl})-\sum_{k,l}h_{kl}^2f\nonumber\\
&=&\sum_{k,l}\alpha_le_k(h_{kl})-f|A|^2.
\end{eqnarray}

Asking further that $Ae_k=\lambda_ke_k$ at $p$ (which is always possible), we have at $p$
\begin{equation}\label{eq:III}
\sum_k\langle Ae_k,\overline\nabla_{e_k}V\rangle=\sum_k\lambda_k\langle e_k,\overline\nabla_{e_k}V\rangle=\sum_k\lambda_k\psi=-nH\psi.
\end{equation}

In order to compute the last summand of (\ref{eq:I}), note that the conformality of $V$ gives
$$\langle\overline\nabla_NV,e_k\rangle+\langle N,\overline\nabla_{e_k}V\rangle=0$$
for all $k$. Hence, differentiating the above relation in the direction of $e_k$, we get
$$\langle\overline\nabla_{e_k}\overline\nabla_NV,e_k\rangle+\langle\overline\nabla_NV,\overline\nabla_{e_k}e_k\rangle+\langle\overline\nabla_{e_k}N,\overline\nabla_{e_k}V\rangle+\langle N,\overline\nabla_{e_k}\overline\nabla_{e_k}V\rangle=0.$$
However, at $p$ one has
\begin{eqnarray*}
\langle\overline\nabla_NV,\overline\nabla_{e_k}e_k\rangle&=&-\langle\overline\nabla_NV,\langle\overline\nabla_{e_k}e_k,N\rangle N\rangle=-\langle\overline\nabla_NV,\lambda_kN\rangle\\
&=&-\lambda_k\psi\langle N,N\rangle=\lambda_k\psi
\end{eqnarray*}
and
$$\langle\overline\nabla_{e_k}N,\overline\nabla_{e_k}V\rangle=-\lambda_k\langle e_k,\overline\nabla_{e_k}V\rangle=-\lambda_k\psi,$$
so that
\begin{equation}\label{eq:IV}
\langle\overline\nabla_{e_k}\overline\nabla_NV,e_k\rangle+\langle N,\overline\nabla_{e_k}\overline\nabla_{e_k}V\rangle=0
\end{equation}
at $p$. On the other hand, since
$$[N,e_k](p)=(\overline\nabla_Ne_k)(p)-(\overline\nabla_{e_k}N)(p)=\lambda_ke_k(p),$$
it follows from (\ref{eq:IV}) that
\begin{eqnarray*}
\langle\overline R(N,e_k)V,e_k\rangle_p&=&\langle\overline\nabla_{e_k}\overline\nabla_NV-\overline\nabla_N\overline\nabla_{e_k}V+\overline\nabla_{[N,e_k]}V,e_k\rangle_p\\
&=&-\langle N,\overline\nabla_{e_k}\overline\nabla_{e_k}V\rangle_p-N\langle\overline\nabla_{e_k}V,e_k\rangle_p+\langle\overline\nabla_{\lambda_ke_k}V,e_k\rangle_p\\
&=&-\langle N,\overline\nabla_{e_k}\overline\nabla_{e_k}V\rangle_p-N(\psi)+\lambda_k\psi,
\end{eqnarray*}
and hence
\begin{equation}\label{eq:V}
\sum_k\langle N,\overline\nabla_{e_k}\overline\nabla_{e_k}V\rangle_p=-nN(\psi)-nH\psi-\overline{Ric}(N,V)_p
\end{equation}
Finally,
\begin{eqnarray*}
\overline{Ric}(N,V)&=&\sum_l\alpha_l\overline{Ric}(N,e_l)-f\overline{Ric}(N,N)\\
&=&\sum_{k,l}\alpha_l\langle\overline R(e_k,e_l)e_k,N\rangle-f\overline{Ric}(N,N),
\end{eqnarray*}
and
\begin{eqnarray*}
\langle\overline R(e_k,e_l)e_k,N\rangle_p&=&\langle\overline\nabla_{e_l}\overline\nabla_{e_k}e_k-\overline\nabla_{e_k}\overline\nabla_{e_l}e_k,N\rangle_p\\
&=&e_l\langle\overline\nabla_{e_k}e_k,N\rangle_p-\langle\overline\nabla_{e_k}e_k,\overline\nabla_{e_l}N\rangle_p-e_k\langle\overline\nabla_{e_l}e_k,N\rangle_p\\
&&+\langle\overline\nabla_{e_l}e_k,\overline\nabla_{e_k}N\rangle_p\\
&=&-e_l\langle e_k,\overline\nabla_{e_k}N\rangle_p+e_k\langle e_k,\overline\nabla_{e_l}N\rangle_p\\
&=&e_l(h_{kk})-e_k(h_{kl}),
\end{eqnarray*}
so that
$$\overline{Ric}(N,V)_p=\sum_{k,l}\alpha_le_l(h_{kk})-\sum_{k,l}\alpha_le_k(h_{kl})-f\overline{Ric}(N,N)_p,$$
and it follows from (\ref{eq:V}) that
\begin{eqnarray}\label{eq:VI}
\sum_k\langle N,\overline\nabla_{e_k}\overline\nabla_{e_k}V\rangle_p&=&-nN(\psi)-nH\psi+V^T(nH)\nonumber\\
&&+\sum_{k,l}\alpha_le_k(h_{kl})+f\overline{Ric}(N,N).
\end{eqnarray}

Substituting (\ref{eq:II}), (\ref{eq:III}) and (\ref{eq:VI}) into (\ref{eq:I}), one gets
the desired formula (\ref{eq:Laplacian formula_I}).
\end{proof}

\section{Applications}

A particular class of conformally stationary spacetimes is that of
{\em generalized Robertson-Walker} spacetimes~\cite{ABC:03}, namely, warped products
$\overline M^{n+1}=I\times_{\phi}F^n$, where $I\subseteq\mathbb R$ is an interval with the
metric $-dt^2$, $F^n$ is an $n$-dimensional Riemannian manifold and
$\phi:I\rightarrow\mathbb R$ is positive and smooth. For such a space, let
$\pi_I:\overline M^{n+1}\rightarrow I$ denote the canonical projection onto the $I-$factor.
Then the vector field
$$V=(\phi\circ\pi_I)\frac{\partial}{\partial t}$$
is conformal, timelike and closed (in the sense that its dual
$1-$form is closed), with conformal factor $\psi=\phi'$, where the
prime denotes differentiation with respect to $t$.
Moreover, according to~\cite{Montiel:99}, for $t_0\in I$, orienting the
(spacelike) leaf $M_{t_0}^n=\{t_0\}\times F^n$ by using the
future-pointing unit normal vector field $N$, it follows that
$M_{t_0}$ has constant mean curvature
$$H=\frac{\phi'(t_0)}{\phi(t_0)}.$$

If $\overline M^{n+1}=I\times_{\phi}F^n$ is a generalized Robertson-Walker spacetime and
$x:M^n\rightarrow\overline M^{n+1}$ is a complete spacelike hypersurface of
$\overline M^{n+1}$, such that $\phi\circ\pi_I$ is limited on $M$, then
$\pi_F\big|_M:M^n\rightarrow F^n$ is necessarily a covering map (\cite{ABC:03}).
In particular, if $M^n$ is closed, then $F^n$ is automatically closed.

One has the following corollary of proposition~\ref{prop:Laplacian of conformal vector field}:

\begin{corollary}
Let $\overline M^{n+1}=I\times_{\phi}F^n$ be a generalized Robertson-Walker spacetime,
and $x:M^n\rightarrow\overline M^{n+1}$ a spacelike hypersurface of $\overline M^{n+1}$,
having constant mean curvature $H$. Let also $N$ be a future-pointing unit normal vector
field globally defined on $M^n$. If $V=(\phi\circ\pi_I)\frac{\partial}{\partial t}$ and
$f=\langle V,N\rangle$, then
\begin{equation}\label{eq:Laplacian formula_II}
\Delta f=\left\{\overline{Ric}(N,N)+|A|^2\right\}f+n\left\{H\phi'+\phi''\langle N,\frac{\partial}{\partial t}\rangle\right\}.
\end{equation}
where $\overline{Ric}$ denotes the Ricci tensor of $\overline M$, $A$ is the second fundamental
form of $x$ with respect to $N$, and $H=-\frac{1}{n}{\rm tr}(A)$ is the mean curvature of $x$.
\end{corollary}

\begin{proof} First of all,
$f=\langle V,N\rangle=\phi\langle N,\frac{\partial}{\partial t}\rangle$, and
it thus follows from (\ref{eq:Laplacian formula_I}) that
$$\Delta f=\left\{\overline{Ric}(N,N)+|A|^2\right\}f+n\left\{H\phi'-N(\phi')\right\}.$$
However,
$$\overline\nabla\phi'=-\langle\overline\nabla\phi',\frac{\partial}{\partial t}\rangle\frac{\partial}{\partial t}=-\phi''\frac{\partial}{\partial t},$$
so that
$$N(\phi')=\langle N,\overline\nabla\phi'\rangle=-\phi''\langle N,\frac{\partial}{\partial t}\rangle$$
\end{proof}

We can now state and prove our main result:

\begin{theorem}
Let $\overline M^{n+1}=I\times_{\phi}F^n$ be a generalized
Robertson-Walker spacetime, and $x:M^n\rightarrow\overline M^{n+1}$
be a closed spacelike hypersurface of $\overline M^{n+1}$, having
constant mean curvature $H$. If the warping function $\phi$
satisfies $\phi''\geq\max\{H\phi',0\}$ and $M^n$ is strongly
stable, then $M^{n}$ is either minimal or a spacelike slice
$M_{t_0}=\{t_0\}\times F$, for some $t_0\in I$.
\end{theorem}

\begin{proof} Since $M^n$ is strongly stable, we have
$$0\geq\mathcal J''(0)(g)=\int_Mg\left\{\Delta g-\left(\overline{Ric}(N,N)+|A|^2\right)g\right\}dM$$
for all $g\in C^{\infty}(M)$ for which $gN$ is the normal component of the variational field of
some variation of $M^n$. In particular, if
$f=\langle V,N\rangle=\phi\langle N,\frac{\partial}{\partial t}\rangle$, where
$V=(\phi\circ\pi_I)\frac{\partial}{\partial t}$, and $g=-f=-\langle V,N\rangle$, then
$$\Delta g=\left\{\overline{Ric}(N,N)+|A|^2\right\}g-n\left\{H\phi'+\phi''\langle N,\frac{\partial}{\partial t}\rangle\right\}.$$
Therefore, $M^n$ stable implies
$$0\geq\int_M\phi\langle N,\frac{\partial}{\partial t}\rangle\left\{H\phi'+\phi''\langle N,\frac{\partial}{\partial t}\rangle\right\}dM$$
Letting $\theta$ be the hyperbolic angle between $N$ and $\frac{\partial}{\partial t}$, it
follows from the reversed Cauchy-Schwarz inequality that
$\cosh\theta=-\langle N,\frac{\partial}{\partial t}\rangle$, with $\cosh\theta\equiv 1$ if and only
if $N$ and $\frac{\partial}{\partial t}$ are collinear at every point, that is, if and only if
$M^n$ is a spacelike leaf $M_{t_0}$ for some $t_0\in I$. Hence,
$$0\geq\int_M\phi\cosh\theta\left\{-H\phi'+\phi''\cosh\theta\right\}dM.$$

Now, notice that $-H\phi'+\phi''\cosh\theta\geq-\phi''+\phi''\cosh\theta$, which gives
$$\phi\cosh\theta(-H\phi'+\phi''\cosh\theta)\geq\phi\phi''\cosh\theta(\cosh\theta-1).$$
Therefore,
$$0\geq\int_M\phi\cosh\theta(-H\phi'+\phi''\cosh\theta)dM\geq\int_M\phi\phi''\cosh\theta(\cosh\theta-1)\geq 0,$$
and hence
$$\phi''(\cosh\theta-1)=0\ \ \text{and}\ \ \phi''=H\phi'$$
on $M$. If, for some $p\in M$, one has $\phi''(p)=0$, then $\phi'H=0$ at $p$. If $H\neq 0$,
then $\phi'(p)=0$. But if this is the case, then
proposition 7.35 of~\cite{O'Neill:83} gives that
$$\overline\nabla_V\frac{\partial}{\partial t}=\frac{\phi'}{\phi}V=0$$
at $p$ for any $V$, and $M$ is totally geodesic at $p$. In particular, $H=0$, a
contradiction. Therefore, either $\phi''(p)=0$ for some $p\in M$, and $M$ is minimal,
or $\phi''\neq 0$ on all of $M$, whence $\cosh\theta=1$ always, and $M$ is an umbilical
leaf such that $\phi''=H\phi'$.
\end{proof}

\remark Note that $\frac{\phi''}{\phi'}=H=\frac{\phi'}{\phi}$, i.e.,
$\phi''\phi-(\phi')^2=0$, which is a limit case of Al\'{i}as and Montiel's timelike convergent condition.


text:
\section{Introduction}
 \label{intro}
 Nitrogen, the fifth most abundant element in the Universe and the fourth most important biogenic element, exists in the form of two stable isotopes: $^{14}$N (the main one) and $^{15}$N. The $^{14}$N/$^{15}$N ratio is thought to be an important indicator of the chemical evolution of the Galaxy, although its value across the Galaxy is still uncertain. The ratio measured for the Proto-Solar Nebula (PSN), from the Solar wind, is 441$\pm$6 (Marty et al.~\citeyear{marty10}), and this value is about two times larger than that measured in the terrestrial atmosphere (TA), derived from N$_{2}$, $\sim$272 (Marty et al.~\citeyear{marty09}). The PSN value is also larger than that measured in some comets (a factor two) and in carbonaceous chondrites (until a factor five). A $^{14}$N/$^{15}$N ratio of 139$\pm$26, from HCN, was estimated in the comet 17P/Holmes after an outburst (Bockelée-Morvan et al.~\citeyear{bockelee-morvan08}). Manfroid et al.~(\citeyear{manfroid09}) estimated a mean value of $^{14}$N/$^{15}$N=148$\pm$6 from CN in 18 comets observed with VLT+UVES. Different values of $^{14}$N/$^{15}$N were estimated in the carbonaceous chondrite Isheyevo: $\sim$424 in a osbornite-bearing calcium-aluminium-rich inclusion (CAI) (Meibom et al.~\citeyear{meibom07}), and from 44 up to 264 (van Kooten et al.~\citeyear{vankooten17}) in lethic clasts.
Measurements in carbonaceous chondrites are made in specific parts of the objects, as described above, while measurements in comets refers to values integrated over a whole part of comets's comas. For these reasons values obtained from extraterrestrial material studied on Earth could be different from observations of comets. At the same time, different values from different part of the same chondrite are measured, and this indicated that different part of a chondrite are formed during different time of the star formation. Meibom et al.~(\citeyear{meibom07}) infer that the carbon-bearing titanium-nitride (osbornite) in the Isheyevo CAI formed by gas-solid condensation, started in a high-temperature ($\sim$2000 K) inner region (<0.1 AU) of the PSN, where all solids were initially evaporated and the gas homogenized. Under such high temperature conditions, little isotopic fractionation is expected between gas and solids, and the nitrogen isotopic composition of osbornite in the Isheyevo CAI must be representative of the Solar nebula. In fact, winds associated with bipolar outflows, or turbulent transport of the hot inner nebula (silicate) dust, may have carried small refractory particles out to colder zones where CAI formed during protoplanet formation. van Kooten et al.~(\citeyear{vankooten17}) infer, from isotope data of Isheyevo lethic clasts, that variations in N isotopes are consistent with the accretion of multiple organic precursors and subsequent alteration by fluids with different (isotopic) compositions, that is in stages more evolved than the enrichments in CAI.

The relation between the $^{15}$N-enrichments in pristine Solar System materials and the natal core is still uncertain. The number of observational works is increasing, but observations towards both low- and high-mass star-forming cores can not reproduce the low $^{14}$N/$^{15}$N values found in comets or meteorites (e.g. Fontani et al.~\citeyear{fontani15}, Colzi et al.~\citeyear{colzi18}, Zeng et al.~\citeyear{zeng17}, Daniel et al.~\citeyear{daniel16}, Kahane et al.~\citeyear{kahane18}). Recent observations of PSN analogs, demonstrated that multiple isotopic reservoirs were present in the early phases of the formation of the Solar System (Hily-Blant et al.~\citeyear{hily-blant17}). In low-mass pre-stellar cores or protostellar objects, values of $^{14}$N/$^{15}$N were found comparable with that of the PSN: 330$\pm$150 from N$_{2}$H$^{+}$ (Daniel et al.~\citeyear{daniel16}), 350 -- 850 from NH$_{3}$ (Gerin et al.~\citeyear{gerin09}), 334$\pm$50 (Lis et al.~\citeyear{lis10}), and an average value of 270$\pm$20 from different molecules, CN, HCN, HNC, HC$_{3}$N and N$_{2}$H$^{+}$ (Kahane et al.,~\citeyear{kahane18}) or even larger, 1000$\pm$200, in N$_{2}$H$^{+}$ (Bizzocchi et al.~\citeyear{bizzocchi13}). On the contrary, observations of nitrile-bearing species such as HCN and HNC in low-mass sources show lower values of the $^{14}$N/$^{15}$N ratio (140--360, Hily-Blant et al.~\citeyear{hily-blant13}; 160--290, Wampfler et al.~\citeyear{wampfler14}), but still not as low as in comets or protoplanetary discs ( 80--160, Guzmán et al.~\citeyear{guzman17}). Daniel et al.~(\citeyear{daniel13}) measured the nitrogen isotopic ratio in several nitrogenated species towards B1b. They found 300$^{+55}_{-40}$ from NH$_{3}$, 400$^{+100}_{-65}$ from N$_{2}$H$^{+}$, 330$^{+60}_{-50}$ from HCN, 225$^{+75}_{-45}$ from HNC and 290$^{+160}_{-80}$ from CN.

The $^{15}$N-enrichment is still not well understood, even from a theoretical point of view. Wirström et al.~(\citeyear{wirstrom12}) suggested that the isotopic enrichments measured in the primitive organic matter probably had their chemical origin in a low-temperature environment in which ion-molecule isotopic exchange reactions are favoured. Roueff et al.~(\citeyear{roueff15}) discusses some possible chemical pathways in cold ({\it T}$=$10 K) and dense gas ({\it n({\rm H}$_{2}$)}=2$\times$10$^{4-5}$cm$^{-3}$) that bring $^{15}$N-enrichments in molecules such as N$_{2}$H$^{+}$, HCN and HNC. However, the $^{14}$N/$^{15}$N found in the model after 10$^{7}$ yr is not different from that of the PSN value, and chemical reactions that bring $^{15}$N in molecules are uncertain or incomplete. Therefore, it seems that time does not play a role in N-fractionation, as found by Colzi et al.~(\citeyear{colzi18}) in a sample of massive cores in different evolutionary stages.
In fact, comparing the $^{14}$N/$^{15}$N ratios derived from HCN and HNC for these different high-mass star-forming cores, they do not find any trend with the evolutionary stages. This result suggests that chemical reactions may still be missing in existing models, or that the enrichment is a local phenomenon occurring in a spatial region much smaller than the beam size of the observations.

The abundance ratio $^{14}$N/$^{15}$N is also considered a good indicator of stellar nucleosynthesis since the two elemental isotopes are not originated in the same way. Both isotopes have indeed an important secondary production in the CNO cycles. There are two types of CNO cycles: a cold cycle and a hot cycle. The cold CNO cycle takes place in main-sequence stars and in the H-burning shells of red giants: $^{14}$N is created from $^{13}$C or $^{17}$O and brought in stellar surfaces through dredge-up on the red giant branch. The hot CNO cycle occurs instead in novae outbursts (Clayton~\citeyear{clayton03}), and is the main way to produce $^{15}$N.
However, there is also a strong primary component of $^{14}$N created in the so-called Hot Bottom Burning (HBB) of asymptotic giant branch (AGB) stars (e.g. Schmitt \& Ness~\citeyear{schmitt02}, Izzard et al.~\citeyear{izzard04}), and an (over-)production of $^{15}$N with respect to $^{14}$N in the relative role played by massive stars and novae (Romano \& Matteucci \citeyear{romano03}, Romano et al.~\citeyear{romano17}). Therefore, $^{15}$N is principally a secondary element.

Different origins of the two nitrogen isotopes lead to an increase of $^{14}$N/$^{15}$N ratio with the Galactocentric distance, up to a distance of $\sim$8 kpc, as predicted by models of Galactic chemical evolution (GCE, Romano \& Matteucci \citeyear{romano03}, Romano et al.~\citeyear{romano17}). However, the relative importance of these processes, or the eventual existence of additional processes, are still unclear. The only way to test this is to provide more observational contraints.
The GCE model made by Romano et al.~(\citeyear{romano17}) was compared  with measurements of the $^{14}$N/$^{15}$N ratios in a sample of warm molecular clouds observed by Adande \& Ziurys (\citeyear{adande12}). The model seems able to reproduce the trend found by Adande \& Ziurys (\citeyear{adande12}), but this latter was obtained from a small sample (22 sources), and mixing up data from different molecules and different instruments (i.e. different observational parameters).

In this work we report the measurements of the $^{14}$N/$^{15}$N ratio derived in a sample of 66 dense cores that span Galactocentric distances in the range $\sim$ 2--12 kpc. Moreover, together with the sources of Colzi et al.~(\citeyear{colzi18}), we present a new Galactocentric behavior of $^{14}$N/$^{15}$N, based now on 87 sources, which makes this new Galactocentric dependence statistically more robust than that of past works. We compare also this trend to GCE models. The results are presented in Sect.~\ref{res} and a discussion of the results is presented in Sect.~\ref{disc}.
\begin{figure*}
\centering
\includegraphics[width=30pc]{spectra_18151.eps}
\caption{Spectra of the source 18151-1208M1 in all the transitions detected: HN$^{13}$C(1--0) (top-left), H$^{15}$NC(1--0) (top-right), H$^{13}$CN(1--0) (bottom-left) and HC$^{15}$N(1--0) (bottom-right). For each spectrum the x-axis represents the systemic velocity listed in Tables \ref{fit1} and \ref{fit2}. The y-axis shows the intensity in main beam temperature units. The red curves are the best fits obtained with MADCUBA (from the method described in Sect.~\ref{res}).}
\centering
\label{text-spectra}
\end{figure*}

 \section{Source selection and observations}
 \label{observations}
 \subsection{Source sample}
The source list of the 66 sources is reported in Table \ref{tab-coo}, where the source coordinates and the distances from the Sun are shown.
A part of the source sample was compiled from a list of high-mass young stellar objects from northern millimeter continuum surveys (Sridharan et al.~\citeyear{sridharan02}, Beuther et al.~\citeyear{beuther02}). We selected 35 sources with (i) distances < 5 kpc, to sample comparable spatial scales, (ii) strong (peak intensity > 0.5 Jy/beam) and compact 1.2 mm continuum emission, suggestive of high-mass young objects embedded in dust, and (iii) observable with the IRAM-30m telescope. Note that we do not consider any criteria related to the line intensity or molecular richness of the source, nor do we consider the centimeter continuum emission and luminosity, in order to avoid possible biases. The sources cover a broad range of luminosities and evolutionary stages (traced with different star formation signposts, e.g. presence of masers, infrared sources, HII regions) which will permit to study the dependence of chemical properties with luminosity and evolutionary stage. This sample is not classified in evolutionary stages (like that of Colzi et al.~\citeyear{colzi18}), yet. This classifation will be done in a forthcoming paper (Mininni et al., in prep). The whole sample and their properties will be presented in the forthcoming paper, while in this paper we focus on the lines of the isotopologues of HCN and HNC. The other 31 sources, indicated with $^{**}$ in Table \ref{tab-coo}, were selected from the literature (Fontani et al.~\citeyear{fontani11}, Wouterloot et al.~\citeyear{wouterloot93}, Beltrán et al.~\citeyear{beltran16}, Cesaroni et al.~\citeyear{cesaroni17}, Tan et al.~\citeyear{tan14}).

The distances are kinematic distances based on the rotation curve of the Galaxy, and are taken from Sridharan et al.~(\citeyear{sridharan02}), Beuther et al.~(\citeyear{beuther02}) and from the other works cited above.
For the sources with the ambiguity of distance we have chosen the near one.

The sample of 66 sources was implemented with that observed by Colzi et al.~(\citeyear{colzi18}) of 27 sources, to obtain a statistically significant sample.   The sample of 27 sources was divided in the three evolutionary stages: 11 high-mass starless cores (HMSCs), 9 high-mass protostellar objects (HMPOs), and 7 ultra-compact H{\sc ii}\ regions (UC HIIs), respectively (see Fontani et al.~\citeyear{fontani11}). Because 6 sources of the two samples were in common: 00117+6412M1, 18089-1732M1, 18517+0437M1, G75, 20081+3122M2 and NGC7538IRS9, we have decided to keep the results found in Colzi et al.~(\citeyear{colzi18}) as the S/N of their observations is better than that of the observations presented here by a factor of $\sim$2 (compare spectra in Appendix \ref{ap-spectra} of this paper and Appendix B of Colzi et al.~\citeyear{colzi18}).
Despite the ratios found in Colzi et al.~(\citeyear{colzi18}) were derived using the kinetic temperatures calculated from ammonia (given in Fontani et al.~\citeyear{fontani11}), that are different from the excitation temperature ({\it T}$_{\rm ex}$) adopted in this work, the $^{14}$N/$^{15}$N ratios are almost independent (differences of 1--10\%) from the chosen temperature (see Sect.~\ref{res}), and then their results could be used to increment our sample.

We finally have a sample of 87 sources: the 27 of Colzi et al.~(\citeyear{colzi18}) and the 60 of this work not in common with the previous one.

 \subsection{Observations}
The 66 sources analyzed in this work were observed with the IRAM 30m telescope on August 2014. We have used the FTS spectrometer (Fast Fourier Transform Spectrometer; Klein et al.~\citeyear{klein12}) to cover the frequency range 85.7--93.5 GHz with a frequency resolution of 200 kHz, corresponding to 0.6 km s$^{-1}$ at the observed frequencies. The half power beam width (HPBW) is about 29 arcsec. The observations were made in position-switching mode. Pointing was checked every 1.5 hours, and focus was corrected at the beginning of the observations and every 4-6 hours. The integration time ranged from 10 to 30 minutes, except for some sources, as 04579+4703, for which it was less (about 5 minutes). The r.m.s. noise levels of the spectra are presented in Appendix \ref{fit1}, \ref{fit2}, \ref{fit3} and \ref{fit4}. The system temperature was in the range 100--150 K during the observations. The observed spectra were first calibrated and reduced using the CLASS/GILDAS package (Pety et al.~\citeyear{pety05}), the data were corrected for platforming and spike channels were removed.
The beam efficiency ($B_{\rm eff}$) for the IRAM 30m telescope at these frequencies is 0.86, while the forward efficiency ($F_{\rm eff}$) is 0.95. The antenna temperatures were converted to main beam temperatures by using the expression: $T_{\rm A}^{*}$ = $T_{\rm MB} \eta _{\rm MB}$, where $\eta_{\rm MB}$ = $B_{\rm eff}$/$F_{\rm eff}$.
 The spectra obtained were exported from the software package CLASS of GILDAS\footnote{The GILDAS
software is available at http://www.iram.fr/IRAMFR/GILDAS} to MADCUBA\footnote{Madrid Data Cube Analysis on ImageJ is a software developed in the
Center of Astrobiology (Madrid, INTA-CSIC) to visualize and analyze
single spectra and datacubes (Martín et al., in prep.).} (see e.g.  Rivilla et al.~\citeyear{rivilla16},  Rivilla et al.~\citeyear{rivilla17}). With MADCUBA we have selected the part of the spectra around the four targeted lines, and we have fitted the baselines with a maximum order of 1. MADCUBA was then used to identify the lines and analyze them. For the analysis we have used the spectroscopic parameters from the Jet Propulsion Laboratory (JPL) molecular catalog\footnote{http://spec.jpl.nasa.gov/}. Rest frequencies are taken from the laboratory works and quantum-chemical calculations cited in Table \ref{observation}.

\begin{figure*}
\centering
\includegraphics[width=40pc]{ColDen.eps}
\caption{Column densities of HCN against those of HC$^{15}$N (left panel) and of HNC against those of H$^{15}$NC (right panel). HCN and HNC column densities are derived from those of H$^{13}$CN and HN$^{13}$C, respectively (see Sect.~\ref{res}). In all panels the filled circles represent the detected sources and the open triangles are the upper limits for the corresponding column densities. In green the 60 sources of this work, in black, red and blue the 27 sources of Colzi et al.~(\citeyear{colzi18}) divided in the three evolutionary stages: high-mass starless cores (HMSCs), high-mass protostellar objects (HMPOs), and ultra-compact H{\sc ii}\ regions (UC HIIs), respectively. The dashed line indicates the mean atomic composition as measured in the terrestrial atmosphere (TA) and the solid line that measured in the Proto-Solar Nebula (PSN). In the lower right corner of each panel we show the distribution of $^{14}$N/$^{15}$N values for HCN (left) and HNC (right). Here the blue vertical line represents the TA value and the red vertical line the PSN value.}
\centering
\label{fig-colden}
\end{figure*}

  \section{Analysis and results}
 \label{res}

 \subsection{Line detection and fitting procedure}
 \label{fitting}
  We have considered detections of HN$^{13}$C(1--0), H$^{15}$NC(1--0), H$^{13}$CN(1--0), and HC$^{15}$N(1--0) if the spectral lines show a peak main beam temperature ({\it T}$_{\rm MB}^{\rm peak}$)$\ge$3$\sigma$, where $\sigma$ is the r.m.s noise of the spectrum; for the detections close to the 3$\sigma$ limit, we have distinguished between tentative detections ($2.5\sigma \leq T_{\rm MB}^{\rm peak}< 3\sigma$) and non detections ($T_{\rm MB}^{\rm peak}<2.5\sigma$). For the non detections we have computed the upper limits for the column densities as explained in Colzi et al.~(\citeyear{colzi18}), with average value of the full width at half maximum ($\Delta v_{\rm 1/2}$) of the lines: $\Delta v_{{\rm HN}^{13}{\rm C}}=2.6\pm 0.1$ km s$^{-1}$, $\Delta v_{{\rm H}^{13}{\rm CN}}=3.3\pm 0.2$ km s$^{-1}$, $\Delta v_{{\rm H}^{15}{\rm NC}}=2.3\pm 0.1$ km s$^{-1}$ and $\Delta v_{{\rm HC}^{15}{\rm N}}=3.1\pm 0.2$ km s$^{-1}$.
 The percentage of detected lines (together with tentative tetections) is the following: the HN$^{13}$C line is clearly detected in 58 cores (88\%), the H$^{13}$CN line is detected in 58 cores (88\%) which, however, are not exactly the same (see Tables \ref{fit1} and \ref{fit2}), the H$^{15}$NC line is detected in 45 cores (68\%) and the HC$^{15}$N line is detected in 49 cores (74\%).

Both HN$^{13}$C(1--0) and H$^{13}$CN(1--0) have hyperfine structure,
but we can resolve it only for H$^{13}$CN(1--0) because the line widths found ($\sim$3 km s$^{-1}$) are smaller than the separations in velocity of the hyperfine components (e.g. Fig.~\ref{text-spectra}). Conversely, H$^{15}$NC(1--0) and HC$^{15}$N(1--0) do not have hyperfine structure. All the lines were fitted with single Gaussians, except in the case of H$^{13}$CN(1--0) for which we take into account the hyperfine structure, using same {\it T}$_{\rm ex}$ and $\Delta v_{1/2}$ for each transition. Fitting the HN$^{13}$C(1--0) with a single Gaussian can overestimate the line width. To estimates how much the line widths are overestimated by this simplified approach, when possible, we fitted a line both with one gaussian, and considering the hyperfine structure. We found that with a single Gaussian, the line width is about 10\% larger than with the hyperfine structure method.

One of the sources, G31.41+0.31, requires some comments: it is a prototypical hot molecular core (HMC) which harbours deeply embedded young stellar objects (YSOs), and the HMC is separated by $\sim$5'' from an ultracompact H{\sc ii}\ region (UC HII). Observations with the IRAM Plateau de Bure demonstrated the existence of a velocity gradient across the core (Beltrán et al.~\citeyear{beltran04} and \citeyear{beltran05}). Subsequently, observations with the Submillimeter Array confirmed the velocity gradient (Cesaroni et al.~\citeyear{cesaroni11}), supporting the claim of the presence of a rotating toroid. Moreover, observations of $^{12}$CO and $^{13}$CO revealed self-absorption caused by infalling gas close to the HMC (Cesaroni et al.~\citeyear{cesaroni11}). Therefore, G31.41+0.31 is a very complicated source and, as can be noted from the spectra in Appendix \ref{ap-spectra}, there is an evidence of possible self-absorption at the center of the lines (for H$^{13}$CN and HN$^{13}$C in particular) and a distorted line profile due to the possible presence of the rotating toroid at larger scales, together with the possibility of outflows from the center of the HMC (Cesaroni et al.~\citeyear{cesaroni17}, Beltrán et al.,~\citeyear{beltran18}).
For these reasons we were not able to fit the lines neither with a single Gaussian nor with the three hyperfine components, because of the complex line shape. All of this makes it very difficult to fit the lines with the approximation used for the other sources; therefore, we decided to exclude this source from the analysis.

For G14.99-0.67, instead, we have identified two velocity components, and fit them simultaneously (G14.99-0.67 and G14.99-2).

All the spectra are shown in Appendix \ref{ap-spectra}.

\subsection{Column density calculation}
To compute the colum densities we have used MADCUBA assuming local thermal equilibrium (LTE) conditions. MADCUBA takes five parameters to model the LTE line profile into account: total column density of the molecule ({\it N}), {\it T}$_{\rm ex}$, $\Delta v_{\rm 1/2}$, peak velocity ($v$), and source angular dimension ($\theta$).
The total column densities of all species have been evaluated assuming that the emission fills the telescope beam (i.e. no beam dilution has been applied), since we do not have any measurement of the emitting region in our sources. Therefore, the derived column densities are beam-averaged values. Observations towards low-mass protostars at high-angular resolution with the NOEMA interferometer indicate that the H$^{13}$CN and HC$^{15}$N(1--0) line emission is indeed very similar (Zapata et al.~\citeyear{zapata13}, Wampfler et al., personal comunication). Therefore, the assumed filling factor should not influence the $^{14}$N/$^{15}$N ratio.
We have fitted all the lines fixing also the {\it T}$_{\rm ex}$. We have used the excitation temperature derived from the K-ladder of the $J$=5-4 transition of the CH$_{3}$CN molecule, which is a good thermometer of massive dense cores (e.g. Araya et al.~\citeyear{araya05}, Purcell et al.~\citeyear{purcell06}). The detailed analysis of this species will be presented in a forthcoming paper (Mininni et al., in prep.).
The critical densities of the analysed transitions are relatively high ($\sim$3$\times$10$^{5}$ cm$^{-3}$), thus the observed transitions could be in non-LTE conditions. As shown by Daniel et al.~(\citeyear{daniel16}), if the excitation temperatures are lower than the kinetic temperatures, the excitation temperature of the less abundant  isotopologue (with $^{15}$N) could be different from that of the more abundant one and the isotopic ratio could also be affected. Therefore, we have checked with RADEX\footnote{ http://var.sron.nl/radex/radex.php} how T$_{ex}$ would change in non-LTE conditions and we have found that the T$_{ex}$ are lower comparing with those obtained with CH$_{3}$CN. However, the T$_{ex}$ of both isotopologues are similar (within $\sim$16\%), being that of the more abundant isotopologue always higher.
Therefore, the hypothesis that the rotational levels of the two molecules are populated with the same T$_{ex}$ is valid, and the column densities are computed with the same formula used in LTE approximation (Caselli \& Ceccarelli \citeyear{caselli12}). Furthermore, we have tested how the ratios can change if we use a non-LTE analysis (RADEX), and we have found that the ratios would be, on average, lower than a factor from 1.1 up to 1.6. However we stress that a precise non-LTE analysis can not be made because we don't know H$_{2}$ column densities and kinetic temperatures of all the sources.

After that, leaving free the other three parameters ({\it N}, {\it v} and $\Delta v_{\rm 1/2}$), we have used the AUTOFIT tool, which compares automatically the synthetic LTE spectra with the observed spectra, and provides the best non-linear least-squared fit using the Levenberg-Marquardt algorithm. Column densities, line widths, peak velocities and opacity ($\tau$) derived from the fit are given in Tables \ref{fit1} and \ref{fit2}. It can be seen that all the transitions are optically thin as expected. It should also be noted that the opacities are lower limits of the real opacity of the lines, since we assume that the emission of the lines comes from the entire beam. Despite this, we do not expect much higher opacities, and the transitions can be considered still optically thin. To test this, we have assumed a source size of 10'' and found that the $\tau$ values increase at most by a factor 10, so that the values would be $\sim0.01$--$0.1$. The assumption of optically thin transitions is hence still valid.

 In Table \ref{colden} column densities of the main isotopologues and of the corresponding $^{15}$N-bearing species are shown.

\begin{figure*}
\centering
\includegraphics[width=40pc]{fig-ratio.eps}
\caption{$^{14}$N/$^{15}$N ratios for HCN (left panel) and HNC (right panel) as a function of Galactocentric distances of the sources listed in Table \ref{tab-coo}. Symbols are the same of Fig.~\ref{fig-colden} but the open triangles are lower limits for the corresponding $^{14}$N/$^{15}$N ratio. The cyan solid line is the linear regression fit computed for the two data sets and the pink solid line is the model of Romano et al.~(\citeyear{romano17}). The dashed orange line is the regression linear fit obtained by Adande \& Ziurys (\citeyear{adande12}). Note that we used the scale from 50 to 1000 of y-axes for visualization purposes, and that one source for HCN, with value 1305$\pm$650, falls outside.}
\centering
\label{fig-ratio}
\end{figure*}

\subsection{$^{14}$N/$^{15}$N ratios}
The reason why we have observed H$^{13}$CN and HN$^{13}$C is that HNC and HCN are usually optically thick in their lower energy rotational transitions (e.g. Padovani et al.~\citeyear{padovani11}). Hence, we decided to use the $^{13}$C-bearing species that have been found to be optically thin in high-mass star-forming cores with similar densities (Colzi et al.~\citeyear{colzi18}), and then, the $^{14}$N/$^{15}$N ratio was computed from the column density ratio of the two isotopologues and corrected by the $^{12}$C/$^{13}$C Galactic trend given in Milam et al.~(\citeyear{milam05}):
\begin{equation}
^{12}\textrm{C}/^{13}\textrm{C}=(6.01\pm1.19) \textrm{ kpc}^{-1}\times\textrm{D}_{\rm GC} +(12.28\pm9.33).
\label{milam}
\end{equation}
This trend is measured with observations of CN, CO and H$_{2}$CO in a sample of 18 molecular clouds that lie in the range of Galactocentric distances 0.09--16.41 kpc.
Note that $^{13}$C-fractionation may affect the abundances of $^{13}$C-bearing species. In fact nitriles and isonitriles are predicted to be significantly depleted in $^{13}$C (Roueff et al.~\citeyear{roueff15}). However, this depletion is at most of a factor 2 and is derived from a chemical model with a fixed kinetic temperature of 10 K, which probably is not the average kinetic temperature of our sources (see Table \ref{tab-coo}). The predictions of this model may therefore not be appropriate for our objects and observational tests to verify whether this theoretical effect is real have yet to be performed.

We have computed the $^{14}$N/$^{15}$N ratios for HNC and HCN, along with uncertainties derived propagating the error on total column densities (Table \ref{colden}). The uncertainties do not include the calibration errors, which cancels out in the ratio because the two lines were obtained in the same spectral setup. In Table \ref{colden} the $^{14}$N/$^{15}$N ratios are given.

As discussed in Sect.~\ref{observations}, we have implemented these observations with those made by Colzi et al.~(\citeyear{colzi18}). The observations of this sample were made with the same receiver and spectrometer of this work, giving us the possibility to combine the two samples with consistency. For further information about the observations and the results concerning the sample of 27 sources, see Colzi et al.~(\citeyear{colzi18}).
The isotopic ratios $^{14}$N/$^{15}$N, of the whole sample are in the range $\sim$115--1305 from HCN, and in the range $\sim$185--780 from HNC, as can be noted in Table \ref{colden}. All of the ratios found in this work are on average consistent with those computed by Colzi et al.~(\citeyear{colzi18}) (Fig.~\ref{fig-colden}). The most direct result of this work, based on an unprecedented statistics, is a firm confirmation of the finding suggested by Colzi et al.~(\citeyear{colzi18}): the lowest values found in pristine Solar System material (e.g. $\sim$150 in nitrogen hydrides of some comets, Manfroid et al.~\citeyear{manfroid09}) are not typical at the core scales.

In Fig.~\ref{fig-colden}, where the column densities of the two samples are shown, at the bottom right of each panel we have made a histogram, with bin width of 40, that shows the distribution on the ratios. We can note that for both molecules the distribution is asymmetric and centered in the bin 310$\leq^{14}$N/$^{15}$N$\leq$350, namely below the PSN value ($\sim$441) and just above the TA value ($\sim$272).

\begin{figure*}
\centering
\includegraphics[width=43pc]{fig-bayes.eps}
\caption{$^{14}$N/$^{15}$N ratios for HCN (left panel) and HNC (right panel) as a function of Galactocentric distances of the sources. In all panels only the detection (not the upper limits), are represented with blue filled circles. The red solid lines are all the linear fit found with the Bayesian method, described in text, and the green solid line is the best fit performed with \emph{linmix}. For comparison we have added with black solid lines the linear regression fits (in grey are represented the errors on this fit), computed as explained in Sect.~\ref{res}. Note that we used the scale from 50 to 1000 of y-axes for visualization purposes, and that one source for HCN, with value 1305$\pm$650, falls outside.}
\centering
\label{fig-bayes}
\end{figure*}

\begin{figure*}
\centering
\includegraphics[width=40pc]{fig-mappe.eps}
\caption{Plan view of the Milky Way. The Galactic center (red asterisk) is at (0,0) and the Sun (red Sun symbol) is at (0,8.34). The background grey discs correspond to the Galactic bar region ($\sim$4 kpc), the solar circle ($\sim$8 kpc), co-rotation of spiral pattern ($\sim$12 kpc) and the edge of major star-formation regions ($\sim$16 kpc). The solid black lines are the center of spiral arms traced by masers, and the dotted lines the 1$\sigma$ widths. For more details see Reid et al.~(\citeyear{reid14}). The filled circles represent sources of our sample implemented with those of Colzi et al.~(\citeyear{colzi18}), and the three colors are the $^{14}$N/$^{15}$N ratios measured for HCN (left panel) and HNC (right panel): in red $^{14}$N/$^{15}$N$\le$300, in green 300$\le^{14}$N/$^{15}$N$\le$390 and in blue $^{14}$N/$^{15}$N$\ge$390.}
\centering
\label{fig-mappe}
\end{figure*}

\subsection{Galactocentric behavior: linear analysis}
The merged sample of 87 objects, given by the sources studied in this work and those analysed by Colzi et al.~(\citeyear{colzi18}), are located at different distances from the Galactic Centre ({\it D}$_{\rm GC}$). This gives us the opportunity to study the Galactocentric trend of the $^{14}$N/$^{15}$N ratio.

Besides our study, one of the most recent works on the $^{14}$N/$^{15}$N Galactocentric gradient is that of Adande \& Ziurys (\citeyear{adande12}). These authors evaluated the ratio from rotational lines of CN, HCN and HNC towards a sample of mixed low- and high-mass star-forming regions observed with different telescopes. In total, their statistics is poorer, as the sample only contains 22 objects. They found an increasing trend of the $^{14}$N/$^{15}$N with {\it D}$_{\rm GC}$ and with an unweighted linear regression fit to the data set they obtained:
\begin{equation}
^{14}\textrm{N}/^{15}\textrm{N}=(21.1\pm5.2)\textrm{ kpc}^{-1}\times\textrm{D}_{\rm GC} +(123.8\pm37.1).
\end{equation}
Our total sample increases the statistics by a factor 4, and we can properly constrain for the first time the $^{14}$N/$^{15}$N Galactocentric trend with this robust statistics. We have chosen to separate the results of the two molecules because the distributions presented in Fig.~\ref{fig-colden} suggest a possible difference between the two isotopologues. We have first computed an unweighted linear regression fit to the data (excluding the lower limits), in order to compare it to the results of Adande \& Ziurys (\citeyear{adande12}), and found:
\begin{equation}
\textrm{HCN}/\textrm{HC}^{15}\textrm{N}=(21\pm9)\textrm{ kpc}^{-1}\times\textrm{D}_{\rm GC} +(250\pm67),
\end{equation}
\begin{equation}
\textrm{HNC}/\textrm{H}^{15}\textrm{NC}=(20\pm6)\textrm{ kpc}^{-1}\times\textrm{D}_{\rm GC} +(221\pm42).
\end{equation}
These trend are shown in Fig.~\ref{fig-ratio}.
Both slopes are consistent with that found by Adande \& Ziurys (\citeyear{adande12}), and are also consistent with each other. However, we found a general trend that shows an offset of about 100 toward larger values for the $^{14}$N/$^{15}$N, as can be noted also in both panels of Fig.~\ref{fig-ratio}.
We have also tried to fit the data with a single power-law relation using the maximum-likelihood Bayesian tool \emph{linmix} by Kelly et al.~(\citeyear{kelly07}), which derives the linear dependence considering measurement uncertainties (Fig.~\ref{fig-bayes}) and excluding the lower limits. The best fits we found are:
\begin{equation}
\label{bay1}
\textrm{HCN}/\textrm{HC}^{15}\textrm{N}=(17\pm6) \textrm{ kpc}^{-1}\times\textrm{D}_{\rm GC} +(223\pm12),
\end{equation}
\begin{equation}
\label{bay2}
\textrm{HNC}/\textrm{H}^{15}\textrm{NC}=(22\pm6) \textrm{ kpc}^{-1}\times\textrm{D}_{\rm GC} +(198\pm12).
\end{equation}
The linear regression fit is consistent with the best Bayesian fit within the uncertainties, as can be seen in Fig.~\ref{fig-bayes}, mostly for HNC. Note that in the case of the Bayesian analysis the errors are lower than in the linear regression fit because in the standard analysis the errors on the parameters are computed separately, as they are considered independent. On the contrary in the Bayesian analysis an optimised procedure, that takes into account the dependence between the two parameters, is performed and converges to the best parameters with lower errors than in the other analysis.

\subsubsection{Dependence with spatial distribution in the Galaxy}
In order to investigate a dependence with the Galactic longitude ($l$, listed in Table \ref{tab-coo}), we have plotted our results in the Galactic plan view of the Milky Way made by Reid et al.~(\citeyear{reid14}). The solid curved lines trace the center of spiral arms as measured with masers associated with young, high-mass stars. We have plotted in the Galactic plane the position of our sources and we have divided the ratio in three intervals: $^{14}$N/$^{15}$N$\le$300, 300$\le$ $^{14}$N/$^{15}$N$\le$390 and $^{14}$N/$^{15}$N$\ge$390 (Fig.~\ref{fig-mappe}). We have decided to use these ranges because each of them contains a similar number of objects. For both molecules no clear trend along spiral arms was found. We have also checked if there was any residual trend after substracting to all the $^{14}$N/$^{15}$N ratios the values obtained at the correponding {\it D}$_{\rm GC}$ of the source with the Bayesian fit, but we have not found trends either. Thus, we conclude that the Galactocentric gradient, probably, is related only to the source distances and not to other local processes.

\begin{figure*}
\centering
\includegraphics[width=40pc]{fig-par.eps}
\caption{Upper panels: $^{13}$C$^{14}$N/$^{12}$C$^{15}$N ratios for HCN (left panel) and HNC (right panel) as a function of Galactocentric distances of the sources. Symbols are the same of Fig.~\ref{fig-ratio} and the cyan solid line is the linear regression fit computed from the two data sets. Note that the negative slope is the consequence of the two combined ratios, $^{13}$C/$^{12}$C and $^{14}$N/$^{15}$N. Bottom panels: same of Fig.~\ref{fig-ratio} with all the sources in black, while the three parabolas are described in Sect.~\ref{res-par}. The yellow one is found using the minimum value in slope of Eq.~\eqref{milam} (Parabola 1), the red one is found using the central value of the slope (Parabola 2) and the pink one is found using the maximum value of the slope (Parabola 3). Note that we used the scale from 50 to 1000 of y-axes for visualization purposes, and that one source for HCN, with value 1305$\pm$650, falls outside.}
\centering
\label{fig-par}
\end{figure*}

\subsection{Galactocentric behavior: parabolic analysis}
\label{res-par}
The linear analysis relies on the $^{12}$C/$^{13}$C relation by Milam et al.~(\citeyear{milam05}) (Eq.~\eqref{milam}),
that we have assumed to compute $^{14}$N/$^{15}$N ratios.
Note that what can be measured directly from our observations are the H$^{13}$CN/HC$^{15}$N and HN$^{13}$C/H$^{15}$NC ratios. In the upper panels of Fig.~\ref{fig-par} these ratios have been plotted as a function of the Galactocentric distances, and, since the points suggest a negative linear trend, we have computed unweighted linear regression fits to these data, obtaining:
\begin{equation}
\textrm{H}^{13}\textrm{CN}/\textrm{HC}^{15}\textrm{N}=(-0.5\pm0.2) \textrm{ kpc}^{-1}\times\textrm{D}_{\rm GC} +(11.0\pm1.2),
\end{equation}
\begin{equation}
\textrm{HN}^{13}\textrm{C}/\textrm{H}^{15}\textrm{NC}=(-0.4\pm0.1) \textrm{ kpc}^{-1}\times\textrm{D}_{\rm GC} +(9.4\pm0.8).
\end{equation}
We have tested the linearity of this relation also using a non-parametric statistics that allows to explore the data sets without assuming an underlying model that describes the data, finding a linear trend of the ratios with the Galactocentric distances.
The $^{14}$N/$^{15}$N ratios, that have been derived in this work, are the product between $^{14}$N$^{13}$C/$^{15}$N$^{12}$C and $^{12}$C/$^{13}$C. Considering that $^{14}$N$^{13}$C/$^{15}$N$^{12}$C ratios are well fitted by a linear trend, and the assumption that $^{12}$C/$^{13}$C is also linear with {\it D}$_{\rm GC}$ (Milam et al.~\citeyear{milam05}), then the derived $^{14}$N/$^{15}$N must have a parabolic one.
The parabolas have equations, for HNC and HCN, respectively:
\begin{equation}
\begin{split}
&\textrm{HCN}/\textrm{HC}^{15}\textrm{N}=[(-0.5\pm0.2) \textrm{ kpc}^{-1}\times\textrm{D}_{\rm GC} +(11.0\pm1.2)]\times\\
&\times[^{12}\textrm{C}/^{13}\textrm{C}],
\end{split}
\end{equation}
\begin{equation}
\begin{split}
&\textrm{HNC}/\textrm{H}^{15}\textrm{NC}=[(-0.4\pm0.1) \textrm{ kpc}^{-1}\times\textrm{D}_{\rm GC} +(9.4\pm0.8)]\times\\
&\times[^{12}\textrm{C}/^{13}\textrm{C}].
\end{split}
\end{equation}
In the bottom panels of Fig.~\ref{fig-par} the three parabolas obtained choosing the slope of $^{12}$C/$^{13}$C as the minimum, the centered and the maximum in the error, are plotted.
Interestingly, with this analysis we are able to reproduce the flattening trend above $\sim$8 kpc predicted by the GCE model of Romano et al.~(\citeyear{romano17}). This is the distance where the parabolic trend start to disagree with the linear trend of the linear analysis.
This result remains the same whatever the assumption of the $^{12}$C/$^{13}$C.

\section{Discussion and Conclusions}
\label{disc}
We have observed the rotational transitions \emph{J}=1-0 of HN$^{13}$C, H$^{15}$NC, H$^{13}$CN and HC$^{15}$N towards 66 massive star-forming cores to derive the $^{14}$N/$^{15}$N ratios. We have complemented this sample with that observed by Colzi et al.~(\citeyear{colzi18}), obtaining a total sample of 87 sources. The isotopic ratios measured range from 115 up to 1305 for HCN and from 185 up to 780 for HNC, which are higher than the low values in the pristine Solar System material.

With the strong statistics of 87 sources a new Galactocentric trend has been derived. From the linear fits given in Eq.~\eqref{bay1} and \eqref{bay2}, we derive that the $^{14}$N/$^{15}$N ratios in the local interstellar medium (ISM), i.e. the values calculated at the Galactocentric distance of the Sun (8.4 kpc, Reid et al.~\citeyear{reid09}), are 383$\pm$51 for HCN and 366$\pm$51 for HNC. These values are both higher than 301$\pm$44 which was obtained from the gradient of Adande \& Ziurys (\citeyear{adande12}). Despite the large uncertainties, our new local interstellar values are more similar to the PSN value of $\sim$441, measured from the Solar wind adopting a totally different approach (Marty et al.~\citeyear{marty10}).

Let us now discuss briefly our findings in the general context of the origin of $^{14}$N and $^{15}$N in the Galaxy. As discussed in Sect.\ref{intro}, primary contributors to $^{14}$N are massive stars in He-shell burning, while the main contributor of $^{15}$N, which is a secondary element, is the hot CNO cycle in novae outbursts. In our Galaxy the abundance of heavy elements is found to decrease with the distance, as observed by Maciel \& Costa (\citeyear{maciel09}). Therefore the $^{14}$N/$^{15}$N ratio should increase with Galactocentric distance as the main contributors of $^{15}$N are novae and some supernovae events (for example, Romano et al.~\citeyear{romano17}). The positive trend found with the Galactocentric distance suggests that there is not a primary component production of $^{15}$N, that is instead important for $^{14}$N in the environment with an high rate of star formation, like at {\it D}$_{\rm GC}$<8 kpc.
The slope of the GCE model of Romano et al.~(\citeyear{romano17}), which takes novae as the main contributors of $^{15}$N into account, is in agreement with our observational results (Fig.~\ref{fig-ratio}), up to 8 kpc. Then, up to 12 kpc the curve in the model flattens. With our linear fit we can not say what is happening after 8 kpc as the linear trend can be biased by points at smaller distances.

A second analysis points out a parabolic trend of $^{14}$N/$^{15}$N ratios with the Galactocentric distance. This is based taking into account the linear assumption made for the $^{12}$C/$^{13}$C ratios (Milam et al.~\citeyear{milam05}). With this analysis we were able to reproduce the flattening trend above 8 kpc, as predicted by the GCE model of Romano et al.~(\citeyear{romano17}).
The flattening of the $^{14}$N/$^{15}$N gradient for {\it D}$_{\rm GC}>$ 8 kpc in the model of Romano et al.~(\citeyear{romano17}) is due to the absence of substantial $^{14}$N-enrichment from low-metallicity massive stars (Nomoto et al.~\citeyear{nomoto13}).

Note that the model is a lower limit of the $^{14}$N/$^{15}$N ratios found (Fig.~\ref{fig-ratio} and bottom panels of Fig.~\ref{fig-par}). An interpretation of this result could be that the first N-fractionation is regulated by the nucleosynthesis processes, and afterwards there could be local (at core level) enrichments of $^{14}$N/$^{15}$N.

Of course, to confirm our results, more observations of sources in the outer Galaxy are needed. New observations could place important constraints on the amount of $^{14}$N produced by massive stars in low-metallicity environments.
The trend we have found in our Milky Way can be relevant
also as "template" for extra-galactic sources, in which the $^{14}$N/$^{15}$N
has been measured in a few objects. For example, toward the massive star forming region N113 in the Large Magellanic Cloud a $^{14}$N/$^{15}$N ratio of 111$\pm$17 was found (Chin et al.~\citeyear{chin99}). Toward the nuclear region of nearby Seyfert 2 galaxy NGC 1068 a $^{14}$N/$^{15}$N ratio greater than 419 was found (Wang et al.~\citeyear{wang14}). Toward the z=0.89 spiral galaxy, located on the line of sight to the quasar PKS1830-211, a ratio of 152$\pm$27 was measured (Muller et al.~\citeyear{muller11}). Moreover, through millimeter molecular absorption lines arising in the z=0.68 absorber toward B0218+357, a ratio of $\sim$120 was measured (Wallström et al.~\citeyear{wallstrom16}). These measurements indicate that even in external galaxies the $^{14}$N/$^{15}$N is far from being homogeneous. We stress, however, that the comparison between local clouds in the Milky Way and external galaxies should be taken with caution. In this kind of sources the relative contribution of the processes responsible for changes in the $^{14}$N/$^{15}$N (nucleosynthesis and/or chemical fractionation), has not been explored. Moreover,
it is difficult to determine local trends because the single clouds cannot be resolved easily.

\section*{Acknowledgements}
We are grateful to the IRAM-30m telescope staff for their help during the observations.
Many thanks to the anonymous referee for the careful reading of the paper and the comments that improved the work.
We would like to thank Mark Reid and the BeSSel Project for the plot of the structure of our Galaxy.
Paola Caselli acknowledges support from the European Research Council (project PALs 320620).
V. M. Rivilla has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk\l{}odowska-Curie grant agreement No 664931. \'A.\, S.-M. acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) in the framework of the collaborative research grant SFB\,956. project A6.

redpajama_set_name: RedPajamaArXiv