\section{Introduction}
Our paper is devoted to spectral theory of quantum graphs -- ordinary differential operators on metric graphs. This is one of the most rapidly growing areas of modern mathematical physics
due to its important applications in physics and applied sciences as well as interesting mathematical problems that emerge \cite{kuchment,gutkin,wires,kurasovschrodinger,kurasovBirk}.
It appeared that symmetries of the underlying metric graphs play a very important role in constructing counterexamples and proving spectral estimates
\cite{parzan,boman,pavel2015,kurasov2014rayleigh,post2012}. If the theory of self-adjoint operators on metric graphs is
rather well-understood, the corresponding theory of non-self-adjoint operators is in its incubatory stage \cite{hussein}. The main subject of our studies is precisely non-self-adjoint quantum graphs.

The main goal of our current paper is to understand connections to the
theory of $ \mathcal P \mathcal T$-symmetric operators -- yet another area of mathematical physics that has got a lot of attention recently
\cite{MR1627442,MR1686605,MR3035410,MR1958251,MR3485039}.
 Standard quantum mechanics in one dimension is described
by self-adjoint differential operators leading to purely real spectrum. Extending the set of allowed operators by including  $ \mathcal P \mathcal T$-symmetric ones leads to the spectrum with
reflection symmetry with respect to the real axis, not only as a set but also including multiplicities. (In the classical studies $ \mathcal P $ is the reflection operator $ (\mathcal P f)(x) = f(-x) $ and
$ \mathcal T $ is the time-reversal operator of complex conjugation $ (\mathcal T f) (x) = \overline{f(x)}.$) If a metric graph possesses a certain automorphism
(symmetry) $ \mathcal P $, then the corresponding differential
operator can be chosen to have  $ \mathcal P \mathcal T$-symmetry, leading to reflection-symmetric spectrum. We would like to understand whether this mechanism is unavoidable for a
quantum graph to have a reflection-symmetric spectrum.

We consider the case of an equilateral star-graph $\Gamma$, as the one in Figure~\ref{fig:stargraph},  formed by  $ N $ identical edges joined  at the central vertex, together with the Laplace operator acting
on it. The metric graph $ \Gamma $ has a rich symmetry group generated by the permutations of the edges.
If, the so-called standard vertex conditions (continuity of the function and vanishing of the sum of normal derivatives), are introduced, then the corresponding operator is self-adjoint
and the spectrum is real (an infinite set of discrete eigenvalues tending to $+ \infty$). One may break the self-adjointness by introducing Robin conditions with non-real parameters at the degree-one
vertices. It is relatively easy to see that if the set of Robin parameters is invariant under conjugation, then the corresponding Laplace operator is  $ \mathcal P \mathcal T$-symmetric with respect
to a certain automorphism
 $ \mathcal P $ of the underlying metric graph $ \Gamma. $ Then the spectrum possesses reflection symmetry with respect to the real axis. Our main question is whether the opposite
statement holds, namely, whether the reflection symmetry of the spectrum implies  $ \mathcal P \mathcal T$-symmetry of the operator with respect to a certain
automorphism $ \mathcal P $
of the metric graph $ \Gamma $. For an operator with discrete spectrum and eigenfunctions of the operator and its adjoint building a biorthogonal basis, one may easily construct a symmetry
operator in the Hilbert space, provided the spectrum is reflection symmetric. Such symmetry operator is defined using the eigenfunctions and there is no guarantee that it comes from an
automorphism of the metric graph.
Our main goal is to prove the following theorem:

\begin{theorem}[Main Theorem]\label{thm:maintheorm}
Let  $ L_{\boldsymbol{h}} $ be the Laplace operator acting on the equilateral star-graph $ \Gamma $ with the domain given by standard vertex conditions at the internal vertex and
Robin conditions including complex valued parameters $ h_i $  at the degree-one vertices. For the spectrum of the operator $  L_{\boldsymbol{h}} $ to possess reflection symmetry
with respect to the real line it is necessary and sufficient that the quantum graph $ L_{\boldsymbol{h}} $ is $ \mathcal{PT}$-symmetric, where $ \mathcal P$ is a certain
automorphism of the
underlying metric star-graph $ \Gamma.$
\end{theorem}

We doubt that the theorem holds for arbitrary metric graphs, but it might be interesting to characterize the class of graphs and vertex conditions for which a similar
statement holds at least for the Laplace operator. In the case under consideration the metric graph $ \Gamma $ is an equilateral star-graph and the group of its automorphisms is
clear. Moreover, we do not consider all possible differential operators on $ \Gamma $ but just the family $ L_{\boldsymbol{h}} $ containing a finite
number of (Robin) parameters. Therefore, we are going to prove first a slightly different theorem and obtain the Main Theorem \ref{thm:maintheorm} as
its corollary:

\begin{theorem}\label{thm:2}
Let  $ L_{\boldsymbol{h}} $ be the Laplace operator acting on the equilateral star-graph $ \Gamma $ with the domain given by standard vertex conditions at the internal vertex and
Robin conditions including complex valued parameters $ h_i $  at the degree-one vertices. For the spectrum of the operator $  L_{\boldsymbol{h}} $ to possess reflection symmetry
with respect to the real line it is necessary and sufficient that the  set of Robin parameters $ \{ h_i \}_{i=1}^N $
 is invariant under conjugation.
\end{theorem}

Spectral theory of non-self-adjoint operators on the star-graph has already been discussed by one of the co-authors\cite{astudillo}. The main difference to the current paper is that
Neumann conditions were assumed at the degree-one vertices, while one considered rotation-invariant conditions at the central vertex. Operators with the Robin
conditions at degree-one vertices were considered by M.\.Znojil\cite{znojil1,znojil2}. All these three papers are devoted to the
 problem of proving that the spectrum of the operator
possesses reflection symmetry (with respect to the real or imaginary axis), provided the vertex conditions have a certain special form. Our interest is in proving the opposite
statement.

The paper is organized as follows. In Section \ref{sec:statement} we describe the model and discuss its elementary spectral properties.
Section \ref{sec:sufficientcondition} is devoted to an easier direction of Theorem \ref{thm:2}, proving if the set of Robin parameters $ \{ h_i \}_{i=1}^N $
is invariant under conjugation, then the spectrum has reflection symmetry with respect to the real axis.
The opposite (harder) statement is proven in the Section \ref{sec:necessarycondition}.
The last section is devoted to the proof of the Main Theorem \ref{thm:maintheorm}. The proofs of elementary Propositions \ref{prop:symmetry} and \ref{prop:symmetricpolynomial}
can be found  in the Appendix.

\section{Description of the Model and Definitions}\label{sec:statement}
In this section we define rigorously the Laplace operator on the equilateral star-graph $\Gamma$, depicted in Figure~\ref{fig:stargraph}, with complex Robin parameters at the degree one vertices.
The operator we define should be consistent with the geometric picture -- the vertex conditions should connect together values of the functions at each vertex separately.
\begin{figure}[H]
\begin{center}
\begin{tikzpicture}
	\draw [thick] (-2,0)--(2,0);
	\draw [thick] (1,-1.73205)--(-1,1.73205);
	\draw [thick] (-1,-1.73205)--(1,1.73205);
	\draw [fill] (1,-1.73205) circle [radius=0.08];
	\draw [fill] (1,1.73205) circle [radius=0.08];
	\draw [fill] (0,0) circle [radius=0.08];
	\draw [fill] (2,0) circle [radius=0.08];
	\draw [fill] (-1,1.73205) circle [radius=0.08];
	\draw [fill] (-1,-1.73205) circle [radius=0.08];
	\draw [fill] (-2,0) circle [radius=0.08];
	\draw [fill] (0,0) circle [radius=0.08];
	\node [below] at (0,-0.1) {$1$};
	\node [right] at (2.05,0) {$0$};
	\node [left] at (-2.05,0) {$0$};
	\node [above left] at (-1,1.73205) {$0$};
	\node [below right] at (1,-1.73205) {$0$};
	\node [above right] at (1,1.73205) {$0$};
	\node [below left] at (-1,-1.73205) {$0$};
	\node [above] at (1,0) {$e_1$};
	\node [below] at (-1,0) {$e_4$};
	\node [above left] at (-0.6,1.73205/2) {$e_3$};
	\node [below right] at (-.6,-1.73205/2) {$e_5$};
	\node [above left] at (.6,1.73205/2) {$e_2$};
	\node [below right] at (0.6,-1.73205/2) {$e_N$};
	\draw [thick, dotted] (-.5,-1.5)--(.5,-1.5);
\end{tikzpicture}
\end{center}\caption{Star-graph $\Gamma$ with $N$ edges.}\label{fig:stargraph}
\end{figure}
More explicitly let us assume that the graph $ \Gamma$ is built of $ N $ edges $e_1$ till $e_N$, each of unit length. We identify each edge with a  separate copy
of the interval $ [0,1] $,
where $ 1 $ corresponds to the unique internal vertex.
Consider the Hilbert space  $\mathscr{H}=\text{L}^2\big((0,1), \mathbb C^N \big) \ni \boldsymbol{u}=(u_1,\ldots,u_N) $
 equipped with the inner product:
\begin{align}
	\langle\,\boldsymbol{u}\,,\boldsymbol{v}\,\rangle:=\sum_{i=1}^{N}\int_{0}^{1}\overline{u}_i(x)v_i(x)\,\mathrm{d} x.
\end{align}
Let $h_1$ till $h_N$, known as \emph{Robin parameters}, be given complex numbers and also let $\boldsymbol{h}:=(h_1,\ldots,h_N)$.  Consider the \emph{Laplace} operator $L_{\boldsymbol{h}}=-\mathrm{d}^2/\mathrm{d} x^2$ with the domain $\dom(L_{\boldsymbol{h}})$, consisting of all  complex-valued functions in the Sobolev space $W_2^2 \big( (0,1), \mathbb C^N \big) $ fulfilling the following vertex conditions:
\begin{itemize}
	\item \emph{standard} vertex conditions at the internal vertex:
	\begin{align}
	u_1(1)=u_2(1)=\cdots=u_N(1),\label{eqn:continuity}\\
	u'_1(1)+u'_2(1)+\cdots+ u'_N(1)=0;\label{eqn:normalderivative}
\end{align}
    \item \emph{complex Robin} conditions at the external vertices:
    \begin{align}\label{eqn:robin}
	u'_j(0)=h_j\,u_j(0),\qquad 1\leqslant j\leqslant N.
\end{align}
\end{itemize}

The reader should note that the operator $L_{\boldsymbol{h}}$ is not necessarily self-adjoint. This stems from the fact that complex values for Robin parameters are allowed. To see this, let $\boldsymbol{u}$ and $\boldsymbol{v}$ be two elements of $\dom(L_{\boldsymbol{h}})$.  Using integration by parts, we have:
\begin{align*}
	\langle\,L_{\boldsymbol{h}}\boldsymbol{u}\,,\boldsymbol{v}\,\rangle-\langle\,\boldsymbol{u}\,,L_{\boldsymbol{h}}\boldsymbol{v}\,\rangle&=\sum_{i=1}^{N}\bigg(-\int_{0}^{1}\overline{u}''_i(x)v_i(x)\,\mathrm{d} x+\int_{0}^{1}\overline{u}_i(x)v''_i(x)\,\mathrm{d} x\bigg)\\&=
	\sum_{i=1}^{N}\Big(-\overline{u}'_i(x)v_i(x)\big|_{0}^{1}+\overline{u}'_i(x)v_i(x)\big|_{0}^{1}\Big)\\&=
	\sum_{i=1}^{N}(\overline{h}_i-h_i)\overline{u}(0)v_i(0)+\sum_{i=1}^{N}\big[\overline{u}_i(1)v'_i(1)-\overline{u}'_i(1)v_i(1)\big],
\end{align*}
where Equations~\eqref{eqn:robin} are used. From Equations~\eqref{eqn:continuity} and \eqref{eqn:normalderivative}, one can see that the second sum above vanishes and we get:
\begin{align}
	\langle\,L_{\boldsymbol{h}}\boldsymbol{u}\,,\boldsymbol{v}\,\rangle-\langle\,\boldsymbol{u}\,,L_{\boldsymbol{h}}\boldsymbol{v}\,\rangle=\sum_{i=1}^{N}(\overline{h}_i-h_i)\overline{u}_i(0)v_i(0),
\end{align}
which does not in general vanish, since we allow for complex Robin parameters. Therefore, the spectrum of $L_{\boldsymbol{h}}$ might contain complex values.

In what follows, we use a boldface letter to denote an $N$-tuple and  the corresponding Roman letter with subscripts, running from $1$ till $N$, to denote its components. For instance, $\boldsymbol{c}$ is used to denote the $N$-tuple $(c_1,\ldots,c_N)$. In this context, we use  $\overline{\boldsymbol{c}}$ to indicate $(\overline{c}_1,\ldots,\overline{c}_N)$ and use  $\boldsymbol{c}_i$, $1\leqslant i\leqslant N$, to indicate the $(N-1)$-tuple obtained from $\boldsymbol{c}$ by suppressing its $i$th component. \noindent For example,  if $\boldsymbol{c}=(-\mathrm{i},1,3\mathrm{i},-3\mathrm{i},\mathrm{i})$, then $\overline{\boldsymbol{c}}=(\mathrm{i},1,-3\mathrm{i},3\mathrm{i},-\mathrm{i})$ and $\boldsymbol{c}_3=(-\mathrm{i},1,-3\mathrm{i},\mathrm{i})$.
\begin{definition}
	An $N$-tuple $\boldsymbol{c}$ is called \emph{invariant under conjugation}, if non-real components of this tuple can be grouped in conjugate pairs.
\end{definition}
\noindent For example, $\boldsymbol{c}$ introduced above is invariant under conjugation but $\boldsymbol{d}=(1,\mathrm{i},\mathrm{i},-\mathrm{i})$ is not.

\begin{definition}\label{dfn:time-reversal}
	The \emph{time-reversal} operator $\mathcal{T}$ is an operator defined on  $\mathscr{H}$ by
\begin{align}\label{eqn:Treversal}
	\mathcal{T}\boldsymbol{u}=\overline{\boldsymbol{u}}.
\end{align}
This is an anti-linear operator.
\end{definition}
\begin{definition}\label{dfn:parity}
    Let $i$ and $j$ be integers and $1\leqslant i\leqslant j\leqslant N$. A linear operator  $\mathcal{P}_{i,j}$ acting on $\mathscr{H}$ as follows
\begin{align}
	(\mathcal{P}_{i,j}\boldsymbol{u})_k=
	\left\{\hspace{-.2cm}
	\begin{array}{rl}
		u_k,&k\not= i, j\\
		u_j,&k=i\\
		u_i,&k=j
	\end{array}
    \right.,  \; \; 1\leqslant k\leqslant N
\end{align}
is called a \emph{permutation} operator. Note that, for any $1\leqslant i\leqslant N$, $\mathcal{P}_{i,i}$ is the identity operator $\id_{\mathscr{H}}$.
\end{definition}
\noindent Permutations $\mathcal{P}_{i,j}$ are elements of the automorphism group for the metric star-graph $\Gamma$. This group naturally induces
 a group of unitary transformations (symmetries) in the Hilbert space $\mathscr{H}$:
\begin{align}
	\mathcal{S}f(x)=f(\mathcal{S}^{-1}x)
\end{align}
for any $ \mathcal S $ from the graph automorphism group.
If all $h_i$'s are equal then the quantum graph, namely, the operator $L_{\boldsymbol{h}}$ possesses the same symmetry group as $ \Gamma$. This is not the case if $ h_i $'s are different.

\begin{definition}\label{dfn:S-symmetric}
	A unitary operator $\mathcal{S}$ defined on $\mathscr{H}$ is called a \emph{symmetry} of the quantum graph $L_{\boldsymbol{h}}$, if $\mathcal{S}L_{\boldsymbol{h}}=L_{\boldsymbol{h}}\mathcal{S}$. In this case, the quantum graph  is called $\mathcal{S}$-\emph{symmetric}. For $\mathcal{S}$ to be a symmetry of $L_{\boldsymbol{h}}$, it is necessary that $\mathcal{S}\big(\dom(L_{\boldsymbol{h}})\big)\subseteq\dom(L_{\boldsymbol{h}})$.
\end{definition}

\section{Sufficient Condition for Theorem \ref{thm:2}}\label{sec:sufficientcondition}

In this section we prove that the invariance of the $N$-tuple  $\boldsymbol{h}=(h_1,\ldots,h_N)$ under conjugation implies that the spectrum of $L_{\boldsymbol{h}}$ possesses reflection symmetry with respect to the real axis. To
establish this, one needs the following proposition whose proof is given in Appendix \ref{app:proposition}.

\begin{prop}\label{prop:symmetry}
	Let $\mathcal{AT}$, where $\mathcal{A}$ is an invertible linear operator defined on $ {H}$ and $\mathcal{T}$ is the time-reversal operator, be a symmetry for a linear operator $ L $ acting in a Hilbert space $  H $, namely, $\mathcal{AT}L=L \mathcal{AT}$. Then if $\lambda$ is an eigenvalue of $L $ with degeneracy $d$, then $\overline{\lambda}$ is also an eigenvalue of $L$ with the same degeneracy $d$.
\end{prop}

\begin{lemma}\label{lma:PTsymmetric}
	If the $N$-tuple $\boldsymbol{h}=(h_1,\ldots,h_N)$, consisting of Robin parameters, is invariant under conjugation, then there exists an automorphism  $\mathcal{P}$ of the graph $ \Gamma $ and hence an invertible linear operator on $\mathscr{H}$ such that the mentioned quantum graph  $L_{\boldsymbol{h}}$ is $\mathcal{PT}$-symmetric, that is, $\mathcal{PT}L_{\boldsymbol{h}}=L_{\boldsymbol{h}}\mathcal{PT}$.
\end{lemma}

\begin{proof}
	If all components of $\boldsymbol{h}$ are real, then it is straightforward to see that, say, $\mathcal{P}_{1,1}\mathcal{T}=\id_{\mathscr{H}}\mathcal{T}$ is a symmetry of the quantum graph. Therefore, suppose that this is not the case and, since $\boldsymbol{h}$ is assumed to be invariant under conjugation, suppose that $\boldsymbol{h}$ has at least two non-real components.  Let $2m$, for some integer $m$ where $1\leqslant m\leqslant\lfloor N/2\rfloor$,  be the number of non-real components of $\boldsymbol{h}$.  Without loss of generality, one can assume that the components of $\boldsymbol{h}$ are coming in an order in which $h_{2j}=\overline{h}_{2j-1}$, for all integers $1\leqslant j\leqslant m$, and all other components (if any) are real.

	We now show that $\mathcal{PT}$, where  $\mathcal{P}:=\prod_{j=1}^{m}\mathcal{P}_{2j-1,2j}$, is a symmetry of $L_{\boldsymbol{h}}$. Expressed in words, $\mathcal{P}$ is the operator that interchanges the $(2j-1)$th and $(2j)$th ($1\leqslant j\leqslant m$) components of the $N$-tuple it is acting on, and leaves the rest (if any) of the components unchanged.

	First, we show that $\dom(L_{\boldsymbol{h}})$ is invariant under the action of $\mathcal{PT}$. Let $\boldsymbol{u}$ be in $\dom(L_{\boldsymbol{h}})$ and let $\boldsymbol{v}=\mathcal{PT}\boldsymbol{u}$. Then, for every integer $j$ in the interval $[1,m]$, $v_{2j-1}=\overline{u}_{2j}$ and $v_{2j}=\overline{u}_{2j-1}$ and, for every integer $j$ in the interval $[2m+1,N]$,  $v_j=\overline{u}_j$. Therefore, Equations~\eqref{eqn:continuity} and \eqref{eqn:normalderivative} are fulfilled for components of $\boldsymbol{v}$, since the corresponding equations for $\boldsymbol{v}$ components are complex conjugates of, at most, a rearrangement of their $\boldsymbol{u}$-component counterparts. To see that Robin vertex conditions are also satisfied for $\boldsymbol{v}$ components, we show that $v'_k(0)=h_k\,v_k(0)$, for every integer $k$ in the interval $[1,N]$.

	If $k=2j-1$, for some integer $j$ in the interval $[1,m]$, then $v_k=\overline{u}_{k+1}$ and $h_k=\overline{h}_{k+1}$. Hence $v'_k(0)=h_k\,v_k(0)$ if and only if $\overline{u}'_{k+1}(0)=\overline{h}_{k+1}\,\overline{u}_{k+1}(0)$ or, equivalently,  $u'_{k+1}(0)=h_{k+1}\,u_{k+1}(0)$. For even integers $k$ in this interval, the proof is similar. If $k$ is an integer in $[2m+1,N]$, provided that this interval is non-empty, then $v_k=\overline{u}_k$ and $h_k=\overline{h}_k$. Thus $v'_k(0)=h_k\,v_k(0)$ is equivalent to $u'_k(0)=h_k\,u_k(0)$ in this case.

	Checking that $L_{\boldsymbol{h}}\mathcal{PT}\boldsymbol{u}=\mathcal{PT}L_{\boldsymbol{h}}\boldsymbol{u}$, for any $\boldsymbol{u}$ in the domain of $L_{\boldsymbol{h}}$, is striaghtforward.
\end{proof}

The fact that the invariance of the $N$-tuple of Robin parameters under complex conjugation implies
reflection symmetry of the spectrum, follows simply as a combination of Proposition~\ref{prop:symmetry} with Lemma~\ref{lma:PTsymmetric}.

\section{Necessary Condition  for Theorem \ref{thm:2} }\label{sec:necessarycondition}
In this section we assume that the spectrum of $L_{\boldsymbol{h}}$ is invariant under complex conjugations
 and we prove that $\boldsymbol{h}$ is invariant under conjugation. To this end, we start by determining the secular equation  for $L_{\boldsymbol{h}}$.
\subsection{The Secular Equation}

Let $\lambda:=z^2$ be a \emph{non-zero} eigenvalue of $L_{\boldsymbol{h}}$.   Hence, there exists a non-zero $\boldsymbol{u}$ in $\dom(L_{\boldsymbol{h}})$ such that
\begin{align*}
	L_{\boldsymbol{h}}\boldsymbol{u}=z^2\,\boldsymbol{u}.
\end{align*}
In particular, after writing it in component form, we have:
\begin{align}\label{eqn:eigenvalue}
	-\frac{\mathrm{d}^2}{\mathrm{d} x^2}u_i(x)=z^2\,u_i(x),\qquad i=1,2,\ldots,N.
\end{align}
The general solutions of Equations~\eqref{eqn:eigenvalue}, for non-zero $z$, have the following form:
\begin{align}\label{eqn:generalsolution}
	u_i(x)=A_i\,\cos zx+B_i\,\sin zx,\qquad i=1,2,\ldots,N,
\end{align}
where $A_i$'s and $B_i$'s are complex constants.
Taking into account the vertex conditions \eqref{eqn:robin} on them, we obtain
$ B_i=(h_iA_i)/z$  and
\begin{equation} \label{alpha}
u_i(x)=\frac{A_i}{z}  \underbrace{ \left(z \cos zx+ h_i\,\sin zx \right) }_{\displaystyle := \alpha_{h_i} (z)} .
\end{equation}
Taking into account the vertex conditions \eqref{eqn:continuity} on these functions, one gets:
\begin{align}\label{eq:008}
	\alpha_{h_1}(z)\,A_1=\alpha_{h_2}(z)\,A_2=\cdots=\alpha_{h_N}(z)\,A_N,
\end{align}
where $\alpha_{h_i}(z) $ were introduced in (\ref{alpha}).
Vertex condition \eqref{eqn:normalderivative} give rise to the following equation:
\begin{align}\label{eq:010}
	\beta_{h_1}(z)\,A_1+\beta_{h_2}(z)\,A_2+\cdots+\beta_{h_N}(z)\,A_N=0,
\end{align}
where
\begin{align}
	\beta_{h_i}(z):=-z\sin z+h_i\,\cos z,\qquad(i=1,2,\ldots,N).
\end{align}

Equations~\eqref{eq:008} along with Equation~\eqref{eq:010} form a system of $N$ linear equations for $N$ unknowns $A_1$ till $A_N$. For this system to have a non-trivial solution, it is necessary and sufficient that the determinant of the coefficients vanishes, namely,
\begin{align}\label{eqn:secularequations}
	D_{\boldsymbol{h}}(z)=0,
\end{align}
where
\begin{align}
    D_{\boldsymbol{h}}(z):=\left|
	\begin{matrix}
		\alpha_{h_1}(z)&-\alpha_{h_2}(z)&0&0&\cdots&0&0\\
		0&\alpha_{h_2}(z)&-\alpha_{h_3}(z)&0&\cdots&0&0\\
		0&0&\alpha_{h_3}(z)&-\alpha_{h_4}(z)&\cdots&0&0\\
		0&0&0&\alpha_{h_4}(z)&\ddots&0&0\\
		\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\
		0&0&0&0&\cdots&\alpha_{h_{N-1}}(z)&-\alpha_{h_N}(z)\\
		\beta_{h_1}(z)&\beta_{h_2}(z)&\beta_{h_3}(z)&\beta_{h_4}(z)&\cdots&\beta_{h_{N-1}}(z)&\beta_{h_N}(z)
	\end{matrix}
	\right|.
\end{align}
Equation~\eqref{eqn:secularequations} is the \emph{secular} equation of $L_{\boldsymbol{h}}$. Expanding the determinant above along the last row and noting that each of the $N$ resultant determinants is in triangular form, one gets the following relation for $D_{\boldsymbol{h}}(z)$:

\begin{align}\label{eqn:D}
	D_{\boldsymbol{h}}(z)=\sum_{i=1}^N\Bigg(\beta_{h_i}(z)\prod_{\substack{j=1\\j\not=i}}^{N}\alpha_{h_j}(z)\Bigg).
\end{align}
The function $ D_{\boldsymbol{h}}(z) $ is an entire function, since the functions $\alpha$ and $\beta$  are entire. Secondly,
\begin{align}\label{eqn:Dconjugate}
	\overline{D_{\boldsymbol{h}}(\overline{z})}=D_{\overline{\boldsymbol{h}}}(z),
\end{align}
since $\overline{\alpha_{h_i}(\overline{z})}=\alpha_{\overline{h}_i}(z)$ and $\overline{\beta_{h_i}(\overline{z})}=\beta_{\overline{h}_i}(z)$.

\subsection{On the Roots of the Secular Equation}
When all the Robin parameters are zero, that is, when $\boldsymbol{h}=\boldsymbol{0}$, Equation~\eqref{eqn:D} reduces to
\begin{align}\label{eqn:D0}
	D_{\boldsymbol{0}}(z)=-N\,z^N&\sin z(\cos z)^{N-1}.
\end{align}
The non-zero roots of $D_{\boldsymbol{0}}(z)$ then are $n\,\uppi$ (for any non-zero integer $n$), each  of multiplicity one, and $n\,\uppi+\uppi/2$ (for any integer $n$), each of multiplicity $N-1$. We concentrate ourselves here on the former set of roots, which are easier to handle.  Now consider the case in which $\boldsymbol{h}\not=\boldsymbol{0}$, and let $\widetilde{z}_n(\boldsymbol{h})$, for some fixed positive integer $n$,  be a root of the secular equation $D_{\boldsymbol{h}}(z)=0$. Then we write $\widetilde{z}_n(\boldsymbol{h})$  in the following form:
\begin{align}\label{eqn:z_n}
	\widetilde{z}_n(\boldsymbol{h})=n\,\uppi+\mathit{\Delta}_n(\boldsymbol{h}),
\end{align}
where $\mathit{\Delta}_n(\boldsymbol{h})$ denotes the deviation of $\widetilde{z}_n(\boldsymbol{h})$, as a root  of equation $D_{\boldsymbol{h}}(z)=0$, from $n\,\uppi$.
 We show that, for sufficiently large $n$,  the correction term $\mathit{\Delta}_n(\boldsymbol{h})$ takes the following form:
\begin{align}\label{eqn:delt_n}
	\mathit{\Delta}_n(\boldsymbol{h})=  \frac{a_1(\boldsymbol{h})}{n}+\frac{a_3(\boldsymbol{h})}{n^3}+\frac{a_5(\boldsymbol{h})}{n^5}+\cdots\raisepunct{,}
\end{align}
with some coefficients $a_i(\boldsymbol{h})$. The general structure of these coefficients  will be determined later.

First we engineer a complex-valued function $K_\epsilon$, in which $\epsilon$ is a positive real parameter,  so that
\begin{align}\label{eqn:implication}
	D_{\boldsymbol{h}}(n\,\uppi+z)=0\Leftrightarrow K_{\frac{1}{n\,\uppi}}(z)=0,
\end{align}
for any complex number $z$ and for any positive integer $n$. Employing the elementary formulas:
\begin{align*}
	\sin(n\,\uppi+z)&=(-1)^n\,\sin z,\\
	\cos(n\,\uppi+z)&=(-1)^n\,\cos z,
\end{align*}
it is straightforward to see that
\begin{align}
	K_{\epsilon}(z):=\sum_{i=1}^{N}\Bigg\{\big[-(1+\epsilon\,z)\sin z+\epsilon h_i\,\cos z\big] \prod_{\substack{j=1\\j\not=i}}^{N}\big[(1+\epsilon\,z)\cos z+\epsilon h_j\,\sin z\big] \Bigg\}\raisepunct{,}
\end{align}
does the job.

We now show that there exists a neighborhood $\mathcal{N}$ of the origin and a positive number $\epsilon_0$ such that for every $\epsilon$ in the interval $(0,\epsilon_0)$, the function $K_\epsilon$ has a unique root $\mathit{\Delta}(\epsilon)$ in $\mathcal{N}$. More importantly, we show that the dependence of $\mathit{\Delta}$ on $\epsilon$ is analytic. Note that $\mathit{\Delta}$ is, in general, a complex-valued function of $\epsilon$.

The equation $K_{\epsilon}(z)=0$ can be considered as a polynomial equation in $\epsilon$:
\begin{align}\label{eqn:g-eqn}
	g_0(z)+g_1(z)\,\epsilon+g_2(z)\,\epsilon^2+\cdots+g_N(z)\,\epsilon^N=0,
\end{align}
in which the coefficient-functions $g_i(z)$ are entire trigonometric functions of $z$. Among all these coefficient-functions, what we actually need is the explicit form of just the first function $g_0(z)$. It is straightforward to see that
\begin{align*}
	g_0(z)=-N\,\sin z(\cos z)^{N-1}.
\end{align*}
For this function we have
\begin{align*}
	g_0(0)=0,\quad g_0'(0)=-N,
\end{align*}
and, consequently, the entire function $g_0$ is invertible in a neighborhood of the origin, where its first derivative is separated from zero.  Using the inverse of $g_0$, the Equation~\eqref{eqn:g-eqn} can be written in the following manner:
\begin{align}\label{eqn:Gepsilon}
	G_\epsilon(z)=z,
\end{align}
 where
 \begin{align}\label{eqn:Gepsilon-g0}
 	G_\epsilon(z):=-g_0^{-1}\Big(g_1(z)\,\epsilon+g_2(z)\,\epsilon^2+\cdots+g_N(z)\,\epsilon^N \Big).
 \end{align}
Of course, $\epsilon$ in ~\eqref{eqn:Gepsilon-g0} should be chosen sufficiently small.  The reader also notes that for any $z$ in the suitable domain,  $g_0^{-1}(z)=-g_0^{-1}(z)$, since $g_0$ itself is an odd function.

 We now have:
\begin{align}\label{eqn:Gprimeepsilon}
 	G_\epsilon'(z)=-\frac{1}{g'_0\Big(g_1(z)\,\epsilon+g_2(z)\,\epsilon^2+\cdots+g_N(z)\,\epsilon^N\Big)}\,\epsilon\big(g'_1(z)+g'_2(z)\,\epsilon+\cdots+g'_N(z)\,\epsilon^{N-1} \big),
\end{align}
 where prime everywhere denotes the derivative with respect to $z$. As Equation~\eqref{eqn:Gprimeepsilon} shows, $|G_\epsilon'(z)|$ is controlled by $\epsilon$. In particular, there exists  $\epsilon_0>0$ and there exists a neighborhood $\mathcal{N}$ such that  $|G_\epsilon'(z)|<1$, for all $z$ in $\mathcal{N}$ and all $\epsilon$ in the interval $(0,\epsilon_0)$. The fixed-point theorem then implies that the Equation~\eqref{eqn:Gepsilon}, or equivalently, the Equation $K_\epsilon(z)=0$ has a unique root $\mathit{\Delta}(\epsilon)$ in $\mathcal{N}$. To see the analytic dependence of $\mathit{\Delta}$ on $\epsilon$, note that the solution curve $\mathit{\Delta}=\mathit{\Delta}(\epsilon)$ can be viewed as the intersection of the following two analytic manifolds:
 \begin{align}
 	x=G_\epsilon(\mathit{\Delta}),\quad x=\mathit{\Delta},
 \end{align}
 in the three-dimensional space, and therefore it is analytic and one can write:
\begin{align} \label{expD}
	\mathit{\Delta}(\epsilon)=b_0(\boldsymbol{h})+b_1(\boldsymbol{h})\,\epsilon+b_2(\boldsymbol{h})\,\epsilon^2+\cdots\raisepunct{,}
\end{align}
with some coefficients $b_i(\boldsymbol{h})$. The function $ K_\epsilon (z) $ possesses the property $K_{-\epsilon}(-z)=-K_\epsilon(z)$
and, consequently, $\mathit{\Delta} (\epsilon) $ is an odd function of $\epsilon$. Therefore,  all even coefficients in the expansion (\ref{expD}) vanish:
\begin{align} \label{expD2}
	\mathit{\Delta}(\epsilon) = b_1(\boldsymbol{h})\,\epsilon+b_3(\boldsymbol{h})\,\epsilon^3+\cdots\,\raisepunct{.}
\end{align}

Now let $n$ be a positive integer such that $1/(n\,\uppi)$ lies in $(0,\epsilon_0)$. Therefore, the unique root $\mathit{\Delta}_{n}(\boldsymbol{h})$ of $K_{\frac{1}{n\,\uppi}}$ in $\mathcal{N}$ or, considering \eqref{eqn:implication},  equivalently the unique root of Equation $D_{\boldsymbol{h}}(n\,\uppi+z)=0$ in $\mathcal{N}$ is
$\mathit{\Delta}_{n}(\boldsymbol{h})=\mathit{\Delta}\big(\frac{1}{n\,\uppi} \big)$, and expansion (\ref{eqn:delt_n}) holds with
\begin{align}
	a_j(\boldsymbol{h}):=\frac{1}{\uppi^j}\,b_j(\boldsymbol{h}),
\end{align}
for all odd positive integers $j$.

Our next goal is to discuss formulas for the expansion coefficients, this can be done in terms of
symmetric polynomials in $ h_i $.

\subsection{Elementary Symmetric Polynomials}
In this subsection, we introduce a special kind of symmetric polynomials that allows us to write Equation~\eqref{eqn:D} in a more flexible form, which fits better to our purpose.

\begin{definition}\label{dfn:symmetricpolynomial}
	For a non-negative integer $m$, the $m$th \emph{elementary symmetric polynomial} in $N$ variables $x_1$ till $x_N$ is denoted by $s_m(x_1,\ldots,x_N)$ and it is defined by
\begin{align}\label{eqn:symmetricpolynomial}
	s_m(x_1,\ldots,x_N)=\left\{\hspace{-.2cm}
	\begin{array}{cc}
	1,&m=0\\[3mm]
	\sum x_{i_1}\!\!\cdots x_{i_m},&1\leqslant m\leqslant N\\[3mm]
	0,&m\geqslant N+1
	\end{array},
    \right.
\end{align}
 where the sum above is over all indices $i_1$ till $i_m$ such that $1\leqslant i_1<\cdots<i_m\leqslant N$.
\end{definition}
For instance, the elementary symmetric polynomials in three variables are:
\begin{align*}
    s_0(x_1,x_2,x_3)&=1,\\
	s_1(x_1,x_2,x_3)&=x_1+x_2+x_3,\\
	s_2(x_1,x_2,x_3)&=x_1x_2+x_2x_3+x_1x_3,\\
	s_3(x_1,x_2,x_3)&=x_1x_2x_3,\\
	s_m(x_1,x_2,x_3)&=0,\qquad(m\geqslant4).
\end{align*}

\noindent It is readily seen that the symmetric polynomial $s_m(x_1,\ldots,x_N)$,  $m\leqslant N$, is  homogeneous of degree $m$ in the sense that for any number $\lambda$,
\begin{align}\label{eqn:homogenity}
	s_m(\lambda\,x_1,\ldots,\lambda\,x_N)=\lambda^m\,s_m(x_1,\ldots,x_N).
\end{align}

In the coming sections, we make use of the properties of elementary symmetric polynomials mentioned in the following proposition, which proof
can be found in Appendix~\ref{app:proposition}.

\begin{prop}\label{prop:symmetricpolynomial}
	Let $\boldsymbol{c}=(c_1,\ldots,c_N)$, and let $i$ and $j$ be integers such that $1\leqslant i\leqslant N$ and   $0\leqslant j\leqslant N-1$. Then
\begin{gather}
    c_i\,s_j(\boldsymbol{c}_i)=s_{j+1}(\boldsymbol{c})-s_{j+1}(\boldsymbol{c}_i)\label{eqn:reccursive},\\
	\sum_{i=1}^{N}s_j(\boldsymbol{c}_i)=(N-j)\,s_{j}(\boldsymbol{c})\label{eqn:023}.
\end{gather}
\end{prop}
\noindent Moreover, we also make use of the following well-known identity:
\begin{align}\label{eqn:identity}
	\prod_{j=1}^{N}(x+c_j)=\sum_{k=0}^{N}s_k(\boldsymbol{c})\,x^{N-k}.
\end{align}

\noindent The following lemma is in some sense the key observation that makes the proof to work.

\begin{lemma}\label{lem:symmetricpolynomials}
An $N$-tuple $\boldsymbol{c}=(c_1,\ldots,c_N)$ of complex numbers is invariant under conjugation if and only if, for all $1\leqslant m\leqslant N$, $s_m(\boldsymbol{c})$ is real.
\end{lemma}

\begin{proof}
First, suppose that $\boldsymbol{c}$ is invariant under conjugation and let $m$, $1\leqslant m\leqslant N$, be a given integer. From Equation~\eqref{eqn:symmetricpolynomial}, we have: \label{page:A}
\begin{align*}
	\overline{s_m(c_1,\ldots,c_N)}=\sum \overline{c_{i_1}\!\!\cdots c}_{i_m}.
\end{align*}
Since $\boldsymbol{c}$ is invariant under conjugation, for each term in the sum above $\overline{c_{i_1}\!\!\cdots c}_{i_m}\!\!\!\!=c_{j_1}\!\!\cdots c_{j_m}$, where $1\leqslant j_1<\cdots<j_m\leqslant N$. Of course,  not all of these $j$ indices are necessarily distinct from their corresponding $i$ indices. Thus, each term in this sum is either real or its complex conjugate is also a term in this sum and, hence, $s_m(c_1,\ldots,c_N)$ is real.

Conversely, suppose that $s_m(\boldsymbol{c})$'s, $1\leqslant m\leqslant N$, are all real. Consider the following polynomial:\label{page:property of A}
\begin{align}
	P(z):=\prod_{i=1}^{N}(z-c_i).
\end{align}
Using Equations~\eqref{eqn:identity} and \eqref{eqn:homogenity}, $P(z)$ can be written as follows:
\begin{align}
	P(z)=z^n+\sum_{m=1}^{N}(-1)^ms_m(\boldsymbol{c})\,z^{N-m},
\end{align}
and, therefore, $P(z)$ is a polynomial with real coefficients. Consequently,  a complex number $z_0$ is a root of $P(z)$ if and only if $\overline{z}_0$ is a root and, since the roots of $P$ are $c_1$, $c_2$, \dots, and $c_N$. This would then mean that $\boldsymbol{c}$ is invariant under conjugation.
\end{proof}

\subsection{The Secular Equation, New Guise}
We are now ready to write Equation~\eqref{eqn:D} in terms of elementary symmetric polynomials. Plugging $z\cos z$ for $x$ and  $h_j\,\sin z$ for $c_j$ in  Equation~\eqref{eqn:identity}, and using Equation~\eqref{eqn:homogenity}, one can write:
\begin{align*}
	\smash[b]{\prod_{\substack{j=1\\j\not=i}}^{N}}\alpha_{h_j}(z)&=\smash[b]{\prod_{\substack{j=1\\j\not=i}}^{N}}(z\cos z+h_j\,\sin z)\\&\hspace{1.3cm}=
	\smash[b]{\sum_{k=0}^{N-1}}s_k(h_1\,\sin z,\ldots,h_{i-1}\,\sin z,h_{i+1}\,\sin z,\ldots,h_N\,\sin z)\,(z\cos z)^{N-1-k}\\&\hspace{6.8cm}=
	\sum_{k=0}^{N-1}s_k(\boldsymbol{h}_i)\,z^{N-1-k}(\sin z)^k(\cos z)^{N-1-k}.
\end{align*}
Hence,
\begin{align*}
	\beta_{h_i}(z)\prod_{\substack{j=1\\j\not=i}}^{N}\alpha_{h_j}(z)&=(-z\sin z+h_i\,\cos z)\sum_{k=0}^{N-1}s_k(\boldsymbol{h}_i)\,z^{N-1-k}(\sin z)^k(\cos z)^{N-1-k},
\end{align*}
and after multiplying the terms, we get:
\begin{align}\label{eqn:024}
	\beta_{h_i}(z)\smash[b]{\prod_{\substack{j=1\\j\not=i}}^{N}}\alpha_{h_j}(z)&=-\smash[b]{\sum_{k=0}^{N-1}}s_k(\boldsymbol{h}_i)\,z^{N-k}(\sin z)^{k+1}(\cos z)^{N-1-k}\notag\\&\hspace{2cm}+\sum_{k=0}^{N-1}h_is_k(\boldsymbol{h}_i)\,z^{N-k-1}(\sin z)^{k}(\cos z)^{N-k}.
\end{align}
The first sum, after the term corresponding to $k=0$ is separated,  can be written in the following form:
\begin{align}\label{eqn:025}
	z^N\sin z(\cos z)^{N-1}+\sum_{k=1}^{N-1}s_k(\boldsymbol{h}_i)\,z^{N-k}(\sin z)^{k+1}(\cos z)^{N-1-k}\raisepunct{.}
\end{align}
Using Equation~\eqref{eqn:reccursive}, the second sum in Equation~\eqref{eqn:024} can be rewritten as follows:
\begin{align}\label{eqn:026}
	\smash[b]{\sum_{k=0}^{N-1}}[s_{k+1}&(\boldsymbol{h})-s_{k+1}(\boldsymbol{h}_i)]\,z^{N-k-1}(\sin z)^{k}(\cos z)^{N-k}\notag\\&\hspace{-.5cm}=
	\smash[b]{\sum_{k=1}^{N}}[s_{k}(\boldsymbol{h})-s_{k}(\boldsymbol{h}_i)]\,z^{N-k}(\sin z)^{k-1}(\cos z)^{N-k+1}\notag\\&\hspace{.5cm}=
	\sum_{k=1}^{N-1}\big\{[s_{k}(\boldsymbol{h})-s_{k}(\boldsymbol{h}_i)]\,z^{N-k}(\sin z)^{k-1}(\cos z)^{N-k+1}\big\}+s_N(\boldsymbol{h})\,(\sin z)^{N-1}\cos z,
\end{align}
where in the last line, the $N$th term of the sum is written separately, and $s_N(\boldsymbol{h}_i)=0$ has been employed. Plugging Equations~\eqref{eqn:025} and \eqref{eqn:026} back into Equation~\eqref{eqn:024}, we obtain:
\begin{multline}\label{eqn:0028}
	\beta_{h_i}(z)\smash[b]{\prod_{\substack{j=1\\j\not=i}}^{N}}\alpha_{h_j}(z)=-z^N\sin z(\cos z)^{N-1}\\\hspace{-1cm}+\smash[b]{\sum_{k=1}^{N-1}}\big\{-s_k(\boldsymbol{h}_i)\,z^{N-k}(\sin z)^{k+1}(\cos z)^{N-k-1}\\\hspace{3cm}+[s_{k}(\boldsymbol{h})-s_{k}(\boldsymbol{h}_i)]\,z^{N-k}(\sin z)^{k-1}(\cos z)^{N-k+1}\big\}\\+s_N(\boldsymbol{h})\,(\sin z)^{N-1}\cos z.
\end{multline}
The summand in the middle term above can be written as:
\begin{align*}
	\big\{[s_k(\boldsymbol{h})-s_k(\boldsymbol{h}_i)]\,\cos^2z-s_k(\boldsymbol{h}_i)\,\sin^2z\big\}\,z^{N-k}(\sin z)^{k-1}(\cos z)^{N-k-1},
\end{align*}
or
\begin{align*}
	[s_k(\boldsymbol{h})\,\cos^2z-s_k(\boldsymbol{h}_i)]\,z^{N-k}(\sin z)^{k-1}(\cos z)^{N-k-1}.
\end{align*}
Plugging back the last expression into Equation~\eqref{eqn:0028}, yields:
\begin{align*}
	\beta_{h_i}(z)&\smash[b]{\prod_{\substack{j=1\\j\not=i}}^{N}}\alpha_{h_j}(z)=-z^N\sin z(\cos z)^{N-1}\\&\hspace{1cm}+\smash[b]{\sum_{k=1}^{N-1}}s_k(\boldsymbol{h})\,z^{N-k}(\sin z)^{k-1}(\cos z)^{N-k+1}\\&\hspace{3cm}-
	\smash[b]{\sum_{k=1}^{N-1}}s_k(\boldsymbol{h}_i)\,z^{N-k}(\sin z)^{k-1}(\cos z)^{N-k-1}\\&\hspace{7.5cm}
	+s_N(\boldsymbol{h})\,(\sin z)^{N-1}\cos z.
\end{align*}
Finally, summing over $i$ and using Equation~\eqref{eqn:023}, one gets:
\begin{align}
	D_{\boldsymbol{h}}(z)=-N\,z^N&\sin z(\cos z)^{N-1}\notag\\&+\,\smash[b]{\sum_{k=1}^{N-1}}\big\{s_k(\boldsymbol{h})\,(k-N\,\sin^2z)z^{N-k}(\sin z)^{k-1}(\cos z)^{N-k-1}\big\}\notag\\&\hspace{6cm}+N\,s_N(\boldsymbol{h})(\sin z)^{N-1}\cos z,
\end{align}
which can be written more compactly as:
\begin{align}\label{eqn:secularequation}
	D_{\boldsymbol{h}}(z)=-N\,z^N&\sin z(\cos z)^{N-1}+z^N\sum_{k=1}^N\frac{s_k(\boldsymbol{h})}{z^k}\,f_k(z),
\end{align}
where
\begin{align}
	f_k(z):=(k-N\,\sin^2z)(\sin z)^{k-1}(\cos z)^{N-k-1},\qquad(k=1,2,\ldots,N).
\end{align}

\subsubsection{A Short Detour}\label{subsub:detour}
Before continuing the main theme, let us have a short detour here to see how practical the formula  given in Equation~\eqref{eqn:secularequation} for the secular equation actually is. Here we use this formula to demonstrate directly that for large $ n$, eigenvalues of $ L_{\boldsymbol{h}} $ can be found close to $n\,\uppi$.
 Assuming that $z\not=0$, Equation~\eqref{eqn:secularequation} can be written as follows:
\begin{align}\label{eqn:perturbed}
	\frac{D_{\boldsymbol{h}}(z)}{z^N}=-N\,\sin z(\cos z)^{N-1}+\sum_{k=1}^N\frac{s_k(\boldsymbol{h})}{z^k}\,f_k(z).
\end{align}
For any positive integer $n$, $n\,\uppi$ is a root of the first term on the right hand side of Equation~\eqref{eqn:perturbed}. As Equation~\eqref{eqn:D0} shows, this term actually gives the secular equation when $\boldsymbol{h}=\boldsymbol{0}$ and $z\not=0$.

Let $\mathcal{C}_n$ denote a circle of some \emph{small}  but fixed radius $\rho_n$ centered at $n\,\uppi$. We have:
\begin{align}\label{eqn:estimate}
	\Big|\frac{s_k(\boldsymbol{h})}{z^k}\Big|\leqslant\binom{N}{k}\bigg(\frac{h}{|z|}\bigg)^k,
\end{align}
where $h:=\max\big\{\,|h_1|,\ldots,|h_N|\,\big\}$. Thus, for any $z$ chosen inside $\mathcal{C}_n$ with large $n$, the expression on the left in  Equation~\eqref{eqn:estimate} becomes small. On the other hand, $f_k$'s are all uniformly bounded inside every $\mathcal{C}_n$. To see this, note that  all $f_k$'s are uniformly bounded inside, say, $\mathcal{C}_1$ and they are periodic with period $2\uppi$. Hence,  the second term on the right of Equation~\eqref{eqn:perturbed} becomes negligible with respect to the first term, for all $z$ inside $\mathcal{C}_n$ with large $n$. Since the second term is continuous inside any of the circles $\mathcal{C}_n$,
 one then expects the secular equation, when $\boldsymbol{h}\not=\boldsymbol{0}$, has a root that is arbitrary close to $n\,\uppi$, if $n$ is sufficiently large.

\subsection{Structure of the  Coefficients $a_j(\boldsymbol{h})$}\label{subsec:coeff}

In this subsection we deal with the practical problem of determining the coefficients $a_j(\boldsymbol{h})$, which appear in the expansion of the correction term $\mathit{\Delta}_n(\boldsymbol{h})$ in Equation~\eqref{eqn:delt_n}. We remind that all even coefficients in the expansion are zero and we  employ a recursive calculation to determine the general structure of odd  coefficients.

We Taylor expand, around $z=0$,  the sine and cosine functions appearing in $D_{\boldsymbol{h}}(z)$ given by Equation~\eqref{eqn:secularequation}  and plug $\widetilde{z}_n(\boldsymbol{h})=n\,\uppi+\mathit{\Delta}_n(\boldsymbol{h})$, with $\mathit{\Delta}_n(\boldsymbol{h})$ given by Equation~\eqref{eqn:delt_n}, into this Taylor's expansion and then collect the resultant series in terms of different powers of $n$. The powers of $n$ that appear in this expansion are:
\begin{align*}
	n^{N-1},n^{N-3}, n^{N-5},n^{N-7},\ldots=\{n^{N-2i+1}\}_{i\in\mathbb{N}=\{0,1,2,3,\ldots\}}
\end{align*}
For $\widetilde{z}_n(\boldsymbol{h})$ to be a root of the secular equation, all the coefficients of the powers of $n$ introduced above, must vanish.
 Equating all coefficients to zero and solving the equations recursively, we get the following result for coefficients $a_i(\boldsymbol{h})$ with odd indices:
\begin{align}\label{eqn:acoefficents}
	a_{2i-1}(\boldsymbol{h})=s_1(\boldsymbol{h})\,F_i\big(s_1(\boldsymbol{h}),\ldots,s_{i-1}(\boldsymbol{h})\big)+\frac{i}{\uppi^{2i-1}N^i}\,\big(s_1(\boldsymbol{h})\big)^{i-1}s_i(\boldsymbol{h}),\qquad i=2,3,4,\ldots,
\end{align}
where $F_i$'s are some polynomial of degree $2i$ in $i-1$ variables. The coefficients of $F_i$'s are real-valued rational functions of $N$. For $i=1$, the corresponding function $F_1$ is zero. More explicitly, for instance:
\begin{align}\label{eqn:a1}
	a_1(\boldsymbol{h})&=\frac{1}{\uppi N}\,s_1(\boldsymbol{h}),
\end{align}
and
\begin{align*}
	a_3(\boldsymbol{h})&=s_1(\boldsymbol{h})F_1\big((s_1(\boldsymbol{h})\big)+\frac{2}{\uppi^3 N^2}\,s_1(\boldsymbol{h})s_2(\boldsymbol{h}),\\
	a_5(\boldsymbol{h})&=s_1(\boldsymbol{h})F_2\big(s_1(\boldsymbol{h}),s_2(\boldsymbol{h})\big)+\frac{3}{\uppi^5N^3}\,\big(s_1(\boldsymbol{h})\big)^2s_3(\boldsymbol{h}),
\end{align*}
where
\begin{align*}
	F_1(x_1)&=-\frac{3N-2}{3\uppi^3N^3}x_1^2-\frac{1}{\uppi^3N^2}x_1,\\
	F_2(x_1,x_2)&=\frac{10N^2-15N+6}{5\uppi^5N^5}x_1^4+\frac{12N-8}{3\uppi^5N^4}x_1^3-\frac{7N-6}{\uppi^5N^4}x_1^2x_2+\frac{2}{\uppi^5N^3}x_1^2-\frac{8}{\uppi^5N^3}x_1x_2+\frac{4}{\uppi^5N^3}x_2^2.
\end{align*}
A closer look at Equations~\eqref{eqn:acoefficents} reveals that $s_i(\boldsymbol{h})$ appears, for the first time, in  the equation for $a_{2i-1}(\boldsymbol{h})$ and, moreover, it appears there with power equal to one. This we use to prove the main result of this section.
\subsection{Necessary Condition, The Proof}
Finally, everything is in order and we finilize the proof for the necessary condition of the main theorem.  In fact, here we  assume that the eigenvalues of $L_{\boldsymbol{h}}$ possess reflection symmetry with respect to the real axis and we prove that the $N$-tuple $\boldsymbol{h}=(h_1,\ldots,h_N)$, consisting of Robin parameters, is invariant under conjugation.

Since the eigenvalues of $L_{\boldsymbol{h}}$ possesses reflection symmetry with respect to the real axis, then $D_{\boldsymbol{h}}(z)=0$ if and only if $D_{\boldsymbol{h}}(\overline{z})=0$ and, by Equation~\eqref{eqn:Dconjugate}, this is equivalent to $D_{\overline{\boldsymbol{h}}}(z)=0$. Therefore, $\widetilde{z}_n(\boldsymbol{h})$ introduced in Equation~\eqref{eqn:z_n}, is a root of $D_{\overline{\boldsymbol{h}}}(z)$ as well. If one uses the same procedure that gave rise to Equation~\eqref{eqn:acoefficents}, but this time for $D_{\overline{\boldsymbol{h}}}(z)$, one gets:
\begin{align}\label{eqn:a'1}
	a_1(\boldsymbol{h})&=\frac{1}{\uppi N}\,\overline{s_1(\boldsymbol{h})},
\end{align}
and
\begin{align}\label{eqn:a'coefficents}
	a_{2i-1}(\boldsymbol{h})=\overline{s_1(\boldsymbol{h})}\,F_i\big(\overline{s_1(\boldsymbol{h})},\ldots,\overline{s_{i-1}(\boldsymbol{h})}\big)+\frac{i}{\uppi^{2i-1}N^i}\,\big(\overline{s_1(\boldsymbol{h})}\big)^{i-1}\overline{s_i(\boldsymbol{h})},\qquad i=2,3,4,\ldots.
\end{align}
Here we used $s_m(\overline{\boldsymbol{h}})=\overline{s_m(\boldsymbol{h})}$ for $m=1,\ldots,N$.

Comparing Equations~\eqref{eqn:a1} and \eqref{eqn:a'1}, one finds that $s_1(\boldsymbol{h})$ is real. Now assume that $s_1(\boldsymbol{h})$ till $s_{k-1}(\boldsymbol{h})$, $1\leqslant k\leqslant N$, are real. From Equation~\eqref{eqn:a'coefficents}, we have:
\begin{align*}
	a_{2k-1}(\boldsymbol{h})=s_1(\boldsymbol{h})\,F_k\big(s_1(\boldsymbol{h}),\ldots,s_{k-1}(\boldsymbol{h})\big)+\frac{k}{\uppi^{2k-1}N^k}\,\big(s_1(\boldsymbol{h})\big)^{k-1}\overline{s_k(\boldsymbol{h})},
\end{align*}
and from Equation~\eqref{eqn:acoefficents}, we have
\begin{align*}
	a_{2k-1}(\boldsymbol{h})=s_1(\boldsymbol{h})\,F_k\big(s_1(\boldsymbol{h}),\ldots,s_{k-1}(\boldsymbol{h})\big)+\frac{k}{\uppi^{2k-1}N^k}\,\big(s_1(\boldsymbol{h})\big)^{k-1}s_k(\boldsymbol{h}).
\end{align*}
Thus $s_k(\boldsymbol{h})$ is also real. Consequently, all polynomials $s_1(\boldsymbol{h})$ till $s_N(\boldsymbol{h})$ are real and by Lemma~\ref{lem:symmetricpolynomials}, $\boldsymbol{h}$ is invariant under conjugation.

\section{Proof of the Main Theorem \ref{thm:maintheorm} and perspectives}

We need to prove that existence of a certain $ \mathcal{PT}$-symmetry is necessary and sufficient for the spectrum of the operator $ L_{\boldsymbol{h}} $
to possess reflection symmetry with respect to the real axis.
Sufficiency (it is a general property of $ \mathcal{PT}$-symmetric operators) is already proven -- it is given by Proposition \ref{prop:symmetry}. Hence it remains to prove the necessity.
Theorem \ref{thm:2} implies that the spectrum possesses reflection symmetry only if the set of Robin parameters $ \{ h_i \} $ is
invariant under complex conjugation. Then  Lemma \ref{lma:PTsymmetric} implies that there exists an automorphism $ \mathcal P$ of $ \Gamma $
(given by a product of edge permutations) such that the operator $ L_{\boldsymbol{h}} $ is $ \mathcal{PT}$-symmetric. It follows that every operator $ L_{\boldsymbol{h}} $
with reflection-symmetric spectrum possesses $ \mathcal{PT}$-symmetry with a certain automorphism $ \mathcal P $
given as a product of permutations. The Main Theorem is proven.

We have already mentioned in the Introduction that we doubt that the statement that reflection symmetry of
the spectrum implies $ \mathcal{PT}$-symmetry of the quantum graph operator $ L $ with respect to a certain automorphism $ \mathcal P $
of the underlying metric graph holds in full generality. This would imply in particular that all quantum graphs with reflection symmetric
spectrum necessarily are defined on symmetric metric graphs. It might happen that this statement holds if one restricts consideration to
vertex conditions given only via complex delta interactions. We are planning to return back to this question in one of our future publications.

Note that in the proof of Theorem \ref{thm:2} and therefore in the proof of the Main Theorem we did not use that the multiplicities of
the conjugated eigenvalues coincide. The reason is that only a small part of the eigenvalues was used. One may show that the
eigenvalues of $ L_{\boldsymbol{h}} $ are situated close to the eigenvalues of the standard Laplace operator $ L_{\boldsymbol{0}}.$ More precisely,
there is exactly one eigenvalue close to $ \uppi n $ and $ N-1 $ eigenvalues close to $ \uppi (n+1/2). $ In our proof we used
just the series of eigenvalues close to $ \uppi n $ and such a series is invariant under conjugation only if it is real. Such eigenvalues cannot be
degenerate. Generalising our results for
arbitrary graphs would require to take into account that the multiplicities in the conjugated pairs of eigenvalues coincide.

Researchers working with $\mathcal{PT}$-symmetric operators are often interested when the operator is not self-adjoint, but the
spectrum is real anyway. One may study this question for quantum graphs getting explicit examples.

\section*{Acknowledgements}
	B.M.G. would like to thank E. Ardonne, I. Mahyaeh,  and F. \v{S}tampach for discussions. He was sponsored, in part, by the Swedish Research Council. P.K. was partially supported by the Center for
  Interdisciplinary Research (ZiF) in Bielefeld in the framework of
  the cooperation group on ``Discrete and continuous models in the
  theory of networks'' and by the Swedish Research Council grant D0497301.\nopagebreak[4]


\section{Introduction}

Massive star formation is a challenging field of observational astronomy, with massive young stellar objects (MYSOs) residing in deeply embedded environments, often at great distances from Earth.
From a theoretical perspective, high- and low-mass star formation were once considered distinctly different, with MYSOs initiating nuclear fusion while still deeply embedded in their parent cores, accretion persisting into the main sequence, in addition to producing strong radiation capable of limiting accretion onto the central star (see reviews by \citealt{Zinnecker07,Tan14}).

However, with recent advances in instrument sensitivity, resolution and data reduction techniques used in astronomy there has emerged a gradual convergence in the observed features present in low- and high-mass star formation.
MYSOs are increasingly found to harbour disk-jet systems - a feature ubiquitous to low mass star formation. Phenomena such as circumstellar disks \citep{Beltran04,Hirota14,Chen16,Ilee16}, episodic ejection \citep{Burns16b} and even jet rotation \citep{Burns15a} are now becoming common targets of investigation in MYSOs.

In disk-jet systems, ejection rates correlate with accretion rates \citep{Corc98,Garatti15} - thus the ejection history of an MYSO allows inference of its accretion history - which is itself unobservable in a practical sense due to the timescales involved. Accretion history, based on ejection history, can therefore be used to compare low- and high-mass star formation.

In addition to accretion mechanisms, outflows serve as another point of comparison - with the outflows of massive stars considered to be less collimated than those of low mass stars \citep{Wu04}.
However, such conclusions are often based on surveys of large scale outflows which could have formed via entrainment by a collimated, precessing jet as was observed in IRAS 20126+4104 \citep{Shep00}. Evidently, all physical scales must be considered when comparing the outflows of low- and high-mass YSOs.

AFGL 5142 is a massive star forming region which has multiple outflows on multiple scales.
The region contains nine millimeter cores, by far the most predominant of which being MM1 and MM2 which exhibit hot core chemistry \citep{Zhang07,Palau11,Palau13}. MM1 has a mass of about 6.5 M$_{\odot}$ \citep{Liu16} and houses an embedded massive star, indicated by the presence of 6.7 GHz methanol masers which trace an infalling disk \citep{Goddi11}. An ionized bipolar jet extends to the NW-SE - at an angle near perpendicular to the disk. Water masers tracing the leading tip of the ionised jet reveal expanding motions in the line of sight and sky plane \citep{Goddi06,Goddi11}, thus giving AFGL 5142 MM1 the image of a prototypical disk-jet system in an MYSO. Furthermore, the combination of features associated with both low- and high-mass star formation make AFGL 5142 MM1 a good target for a comparative investigation.

About 1$^{\prime \prime}$ to the south is MM2, with a mass of about 6.2 M$_{\odot}$ \citep{Liu16} but no centimeter emission or 6.7 GHz methanol masers. Both MM1 and MM2 exhibit similar hot molecular core chemical compositions \citep{Zhang07,Palau11}, suggesting similar evolutionary stages. Water masers have also been detected near AFGL 5142 MM2 \citep{Hunter95,Goddi06}.

Interferometric observations of the AFGL 5142 region in CO ($2-1$) \citep{Zhang07,Palau11}, SiO ($2-1$) and HCO$^+$ ($1-0$) \citep{Hunter99} reveal at least four distinct collimated molecular outflows on $10^{ \prime \prime}$ scales; outflows A, B, C \citep{Zhang07} and outflow D \citep{Palau11}. Several of these outflows intersect both MM1 and MM2, barring a simple assignment of progenitors. Single dish observations reveal the presence of larger, arc-minute scale outflows traced in CO ($2-1$) which also intersect the MM1 and MM2 region \citep{Hunter95} these propagate primarily in the skyplane. Recently \citet{Liu16} also find a $\sim 9^{ \prime \prime}$ length, extremely wide-angle bipolar outflow (EWBO), with an opening angle of $\sim$180$^{\circ}$, driven from the MM1 - MM2 region, leading the authors to interpret its formation as being driven by a precessing jet rather than entrainment at wide angles from a narrow jet.

In this work we investigate the outflows seen in AFGL 5142 at the smallest scales by performing new multi-epoch very long baseline interferometry (VLBI) observations of water masers at 22 GHz. VLBI is an ideal observational approach to investigating protostellar ejections as multi-epoch observations provide sky-plane and line of sight information, thus revealing the 3D kinematics of gas in the vicinity of MYSOs.
Our investigation aims to provide observational characterisation of outflow behaviour (launching mechanisms and episodic ejection) for a prototypical MYSO, to provide comparison between low- and high-mass formation mechanisms.
Our report concentrates on the kinematics and energetics of the jets associated with AFLG 5142 MM1. We use our findings to offer a reinterpretation of the recent history of jets and outflows in this region and attempt to assign progenitors to each of the known outflows from the literature. Via annual parallax we also provide the first precise measurement of the distance to AFGL 5142; a quantity essential for converting observations to physical quantities such as ejection velocities, ages and momentum rates.

\section{Observations and data reduction}

VLBI observations of AFGL 5142 were carried out with VERA (VLBI exploration of radio astrometry). All observations were conducted in dual-beam mode \citep{dual} with beams centered on the maser target, AFGL 5142, and a reference quasar, J0533+3451. For both beams, left-hand circular polarisation signals were recorded at each of the 4 VERA station. Data were collected and correlated using the Mitaka FX correlator \citep{Chikada91}, adopting a rest frequency of 22.235080 GHz. For AFGL 5142 data we used a phase tracking center of
$(\alpha, \delta)_{\mathrm{J}2000.0}=(05^{\mathrm{h}}30^{\mathrm{m}}48^{\mathrm{s}}.01733742$,
+33$^{\circ}$47'54".56750)
and for J0533+3451 data we set a phase tracking center at
$(\alpha, \delta)_{\mathrm{J}2000.0}=(05^{\mathrm{h}}33^{\mathrm{m}}12^{\mathrm{s}}.76510600$,
+34$^{\circ}$51'30".336995). Further phase corrections, including more accurate atmospheric models and antennae positions than those used at correlation, were made and applied post-calibration.

\begin{table}
\scriptsize
\caption{\small Summary of observations\label{obs}}
\begin{center}
\begin{tabular}{cccc}
\hline
Epoch& Observation & Modified& Number of \\
number& date& Julian date &features\\ \hline
1&  21 Apr 2014 & 56768 & 12\\
2&  20 May 2014 & 56797 & 9\\
3& 2 Oct 2014 &   56932 & 17\\
4& 25 Nov 2014 &  56986 & 22\\
5& 31 Jan 2015 &  57053 & 24\\
6& 29 Mar 2015 &  57110 & 25\\
7& 29 May 2015 &  57171 & 19\\
\hline
\end{tabular}
\end{center}
\end{table}

\begin{figure}
\begin{center}
\includegraphics[scale=0.88]{SpectraFull444.eps}
\caption{\small Scalar averaged cross-power spectra of the maser emission in AFGL 5142 as a function of time, coloured arbitrarily.
\label{spectra}}
\end{center}
\end{figure}

The total correlator bandwidth of 240 MHz was shared into 16 intermediate frequencies (IFs). One IF was allocated to the maser data with a bandwidth of 8 MHz and 15.63 kHz channel spacing, providing a 0.12 km s$^{-1}$ velocity spacing. The remaining 15 IFs were allocated to the data of the quasar reference source - one IF had similar properties to the maser data while the other 14 IFs had 16 MHz bandwidth and spacings of 125 kHz. The 15 quasar data IFs were manipulated into common form and merged - resulting in almost continuous frequency coverage spanning 232 MHz.

Data reduction was performed using AIPS (Astronomical Image Processing System), developed by National Radio Astronomy Observatory. VLBI phase referencing data reduction made use of the \emph{inverse phase referencing} technique, customised for the dual-beam data of VERA. The technique was first introduced in \citet{Imai12} and further development, in addition to a guide to its implementation, is given in \citet{Burns15a}.

\begin{table*}
\vspace{-0.2cm}
\begin{center}
\caption{\small The general properties of H$_{2}$O masers in AFGL 5142, detected with VERA. \label{TAB}}
\begin{tabular}{cclcccccccc}
\hline
Maser&$V_{\rm LSR}$&Detected &$\Delta \alpha \cos \delta$&$\Delta \delta$ &$\mu_{\alpha}\cos\delta$&$\mu_{\delta}$& F$_{\rm int}$&$\pi$\\
ID &(km s$^{-1}$)&epochs&(mas)&(mas)&(mas yr$^{-1}$)&(mas yr$^{-1}$)&(Jy)& (mas)\\
\hline
N.W. 1	&	$	2.47	$	&	1234567 &	-35.216	&	35.133	&	$	-0.55	\pm	0.03	$	&	$	0.89	\pm	0.26	$	&23.67&	$	0.460 \pm 0.014	$	\\
N.W. 2	&	$	0.42	$	&	1234567 &	-39.637	&	29.447	&	$	-1.03	\pm	0.02	$	&	$	0.23	\pm	0.11	$	&11.09&	$	0.479 \pm 0.021	$	\\
N.W. 3	&	$	3.49	$	&	****567 &	-28.365	&	38.437	&	$	0.07	\pm	0.27	$	&	$	0.16	\pm	0.05	$	&1.27&	$		$	\\
N.W. 4	&	$	3.07	$	&	1*34*67 &	-37.371	&	32.979	&	$	-0.58	\pm	0.08	$	&	$	1.07	\pm	0.08	$	&2.99&	$		$	\\
N.W. 5	&	$	1.81	$	&	**34567 &	-38.121	&	31.934	&	$	-0.19	\pm	0.48	$	&	$	1.51	\pm	0.80	$	&4.13&	$		$	\\
N.W. 6	&	$	1.24	$	&	****567 &	-38.838	&	30.803	&	$	-1.28	\pm	0.06	$	&	$	-0.21	\pm	0.11	$	&4.14&	$		$	\\
N.W. 7	&	$	0.86	$	&	**3*567 &	-41.625	&	24.954	&	$	-0.77	\pm	0.12	$	&	$	0.87	\pm	0.31	$	&2.22&	$		$	\\
N.W. 8	&	$	-0.02	$	&	****567 &	-45.198	&	19.198	&	$	-0.73	\pm	0.42	$	&	$	1.98	\pm	0.21	$	&2.46&	$		$	\\
\hline
S.E. 1	&	$	-6.69	$	&	123456* &	162.415	&	-232.163	&	$	1.13	\pm	0.14	$	&	$	-1.59	\pm	0.04	$	&4.83&	$		$	\\
S.E. 2	&	$	-6.20	$	&	1234567	&	156.283	&	-237.916	&	$	1.93	\pm	0.03	$	&	$	-1.16	\pm	0.03	$	&12.87&	$	0.479 \pm 0.022	$	\\
S.E. 3	&	$	-6.71	$	&	1**4567	&	147.805	&	-222.823	&	$	1.38	\pm	0.19	$	&	$	-0.83	\pm	0.11	$	&1.77&	$		$	\\
S.E. 4	&	$	-4.93	$	&	**34567	&	134.663	&	-228.13 	&	$	0.84	\pm	0.04	$	&	$	-1.19	\pm	0.08	$	&5.31&	$		$	\\
S.E. 5	&	$	-5.36	$	&	***4567 &	117.507	&	-230.753	&	$	0.36	\pm	0.04	$	&	$	-1.08	\pm	0.02	$	&2.72&	$		$	\\
S.E. 6	&	$	-6.91	$	&	****567	&	160.06	&	-235.217	&	$	1.97	\pm	0.70	$	&	$	-1.67	\pm	0.26	$	&2.52&	$		$	\\
S.E. 7	&	$	-4.08	$	&	**3****	&	121.616	&	-227.994	&	$	-		$	&	$	-		$	&1.41&	$		$	\\
S.E. 8	&	$	-7.05	$	&	*****6*	&	131.305	&	-229.803	&	$	-		$	&	$	-		$	&5.80&	$		$	\\
S.E. 9	&	$	-7.05	$	&	*****6*	&	159.065	&	-236.213	&	$	-		$	&	$	-		$	&0.88&	$		$	\\
S.E. 10	&	$	-5.77	$	&	**3****	&	54.856	&	-143.514	&	$	-		$	&	$	-		$	&2.15&	$		$	\\
\hline
F.S. 1	&	$	-5.78	$	&	1234567	&	11.783	&	-1241.213	&	$	1.42	\pm	0.03	$	&	$	-0.43	\pm	0.07	$	&13.98&	$	0.450 \pm 0.029	$	\\
F.S. 2	&	$	-4.41	$	&	***4567	&	12.405	&	-1227.991	&	$	1.64	\pm	0.08	$	&	$	0.05	\pm	0.21	$	&17.80&	$		$	\\
F.S. 3	&	$	-5.15	$	&	123456*	&	13.319	&	-1236.244	&	$	1.71	\pm	0.04	$	&	$	0.71	\pm	0.18	$	&5.74&	$		$	\\
F.S. 4	&	$	-4.43	$	&	**34567	&	13.227	&	-1230.149	&	$	1.35	\pm	0.08	$	&	$	1.00	\pm	0.11	$	&4.67&	$		$	\\
F.S. 5	&	$	-5.54	$	&	1234567	&	11.416	&	-1225.99	&	$	1.13	\pm	0.09	$	&	$	1.44	\pm	0.23	$	&11.02&	$		$	\\
\hline
B	&	$	0.77	$	&	1*3456*	&	639.914		&	148.019	&	$	1.89	\pm	0.06	$	&	$	-0.61	\pm	0.04	$	&3.06&	$		$	\\
F.S.W.&	$	-3.04	$	&	**3456*	&	-3457.993	&	-1553.682&	$	-2.56	\pm	-0.46	$	&	$	0.01	\pm	0.18	$	&3.29&	$		$	\\
A	&	$	-3.18	$	&	123456*	&	6.68	&	-40.77	&	$	-		$	&	$	-		$	&1.27&	$		$	\\
&&&&&&& Simultaneous fit & $0.467 \pm 0.010$ \\
\hline
\end{tabular}
\begin{tablenotes}
\item{\footnotesize{Column (2)}: Line of sight velocities, with respect to the local standard of rest, are quoted as those measured for the first detection.}
\item{\footnotesize{Column (3)}: Numbers indicate detection in the corresponding epoch, while asterisk represents non-detection.}
\item{\footnotesize{Column (8)}: Integrated fluxes are given as the average of measurements from all epochs.}
\end{tablenotes}
\end{center}
\end{table*}

In short, data reduction involves solving frequency dependant phase terms (group delay) using bright calibrators in the wider 15 IF data set, while time dependant phase terms are solved using the narrow but bright line emission of the maser. Both sets of solutions are then applied to the reference source, J0533+3451, inversely giving the angular separation of the maser from the phase center. Assuming that the reference source is fixed at its ICRF coordinates we then obtain the astrometric position of the maser from the derived separations. Since phase solutions determined using the maser are applied to itself we also produce high quality, astrometrically accurate maser maps. Typically, maser maps achieved rms noise values of $100 - 200$ mJy beam$^{-1}$ while the wider band continuum data achieved typical rms values of $\sim 1$ mJy.

Maser maps were produced using the AIPS task CLEAN in automatic mode (DOTV = 0), with a synthesised beam of dimensions $1.0\times1.2$ mas. Masers were identified from the maps using an automated \emph{SAD} routine in AIPS, employing a detection signal to noise cutoff of 7. Following common nomenclature, a `spot' refers to an individual maser emission peak, imaged in a single spectral channel - while a maser `feature' refers to a collection of spots considered to emanate from the same physical maser cloud. Maser spots were grouped into a feature if they shared a spectral feature and collocated within a radius of 1 mas. The associations of spots to features, and the flux weighted astrometric positions of maser features, were determined using a basic FORTRAN code which also removes the offset introduced by the difference in coordinates between the phase tracking center and the reference maser used in each epoch, in addition to shifting the maser maps to the absolute reference frame.

\begin{figure*}
\begin{center}
\includegraphics[scale=1.42]{I05274_FullScale.dat2.eps}
\caption{\small Internal motions of water masers in AFGL 5142. Vector magnitudes indicate proper motion while colours indicate LSR velocity. Annotations are the same as those in Table~\ref{TAB}, with S.E. i and o referring to the `inner' and `outer' South West groups, respectively. The map origin is at the phase tracing center of the maser data.
\label{FULL}}
\end{center}
\end{figure*}

\null

\section{Results}

\subsection{Maser distribution and l.o.s. velocities}
A total of 27 independent maser features were detected with VERA during the observing calendar. Their coordinates, line of sight (l.o.s) velocities, detection frequencies and other properties are recored in Table~\ref{TAB}.
With regards to temporal variability, some maser features persisted through all epochs, while others exhibited more sporadic behaviour. No periodic behaviour was noted. Figure~\ref{spectra} shows the temporal changes in the overall emission spectrum, and the variability of individual maser features can be read from the detection frequencies in column 3 of Table~\ref{TAB}.  The distributions of all maser features detected by VERA are shown in Figure~\ref{FULL}.

\begin{figure}
\begin{center}
\includegraphics[scale=0.88]{I05274.dat2.eps}
\caption{\small Same as Figure~\ref{FULL} for the MM1 core region.
\label{MM1}}
\end{center}
\end{figure}

Figure~\ref{MM1} shows a magnified view of the MM1 region. Masers trace three arc structures contained in a region of about $300 \times 300$ mas. There exists one maser arc to the North West of the system center, and two proximate and similarly oriented arcs lying to the South East, one closer and one further from the system center. We refer to these as the `N.W.', the `S.E. inner', and `S.E. outer' maser arcs, respectively - and are hence categorised accordingly in column 1 of Table~\ref{TAB}. The N.W. masers arcs are redshifted- and the S.E. masers blueshifted with respect to the l.o.s velocity of the region ($-1.1$ km s$^{-1}$ \citealt{Zhang07}).

Maser feature A (see Table~\ref{TAB}) is also associated with MM1. This maser does not associate with the aforementioned arcs, and rather, stands in the midway - near the system center, at map position $(\alpha, \delta) = (7,-41)$ mas. This maser exhibits some remarkable properties which will be discussed in a later publication.

Our observations did not detect all of the masers reported by \citet{Goddi06}. In their VLBA observations they find more water masers associated with MM1 which are located further to the N.W. We include their results in our investigation of the maser distribution in AFGL 5142 MM1 (see Section~\ref{epi}).

Water masers associated with MM2 are arranged in a single arc (see Figure~\ref{FSBOW}). Due to their location `far south' of the phase tracking center we refer to these masers as the F.S. group in Table~\ref{TAB}. These masers likely correspond to maser features `IIb' of \citet{Goddi06}. All F.S. masers have slightly blueshifted l.o.s. velocities.

We found one maser to the North East (`B' in Table~\ref{TAB}) of the phase tracking center - possibly associated with MM1 and probably corresponding to maser feature `IIr' of \citet{Goddi06}. One maser feature was found in the far South West (`F.S.W.'), associated with AFGL 5142 MM6.

\subsection{Fitting parallax and proper motion}

The astrometry of water masers observed with VERA was used to measure the annual parallax and resulting trigonometric distance of the AFGL 5142 region.
To assure the best possible precision for the distance estimate we performed parallax fitting only using maser features which were detected in all 7 epochs - spanning just over one year. Four maser features passed this criterion, two of which were in the N.W. bowshock, one in the S.E. bowshock and one in the F.S. bow shock. Parallax fitting was done simultaneously for the four maser features, assuming a common distance, to average and smooth out noise-like and structure induced errors. Additional error floors were added and reduced iteratively until a $\chi ^2$ value of unity was reached. Required error floors were 0.04 and 0.27 mas in R.A and Dec., respectively. \

We estimate the annual parallax of AFGL 5124 to be $\pi=0.467 \pm 0.010$ mas, corresponding to a trigonometric distance of $D=2.14^{+0.051}_{-0.049}$ kpc, firmly placing it in the Perseus Arm. Fits are shown in Figure~\ref{PIE}. To confirm the reliability of our simultaneous fitting the parallaxes (\emph{see} Table~\ref{TAB}: column 7) and proper motions of each of the four maser features were then fit independently and were found to be consistent with each other and with the simultaneous fit.
Note that our consistent distances to MM1 and MM2 confirm that they are in close physical proximity.

Using Feature N.W. 1 as reference we obtained the relative proper motions of all other maser features that were detected in at least three epochs, of which there were 21. These relative motions were then converted to absolute proper motions (\emph{see} Table~\ref{TAB}: columns 6,7) using the absolute proper motion of Feature N.W. 1 from the astrometric fitting procedure described earlier. Proper motion errors were calculated as the quadrature sum of the standard deviation in the relative proper motion, and the proper motion uncertainty of Feature N.W. 1.

\begin{figure}
\vspace{-0.666cm}
\begin{center}
\includegraphics[scale=0.85]{I05274.eps}
\caption{\small Annual parallax motions of four maser features in AFGL 5142 in (\emph{above}) R.A. and (\emph{below}) Dec. directions. Maser features are coloured arbitrarily and error floors are represented by vertical bars.
\label{PIE}}
\end{center}
\end{figure}

\subsection{Systemic motion and internal motions}

Absolute (observed) proper motions are a composition of the internal motions of masers in the frame of the driving source, and the apparent systemic motion of the source across the sky as governed by Galactic dynamics. In the case where a source's internal motions are expected to be symmetric the sum of all motion vectors should balance to zero - the systemic proper motion of the region is derived from the residual of this summation. In AFGL 5124 we used the proper motions of masers associated with the N.W. and S.E. maser bowshocks in MM1 (8 and 6 maser features, respectively) since the masers in this region are known to trace a symmetric bipolar ejection \citep{Goddi07}. The resulting systemic motion was ($\mu_{\alpha}\cos\delta$, $\mu_{\delta}$) = ($ +0.32\pm0.27$, $-0.22\pm0.47$) mas yr$^{-1}$, where errors were calculated as the quadrature sum of the average errors of all maser motions in each bowshock.

Subtracting the systemic proper motion from all absolute proper motions reveals the internal motions in the frame of the driving source, these are shown in Figures~\ref{FULL} and~\ref{MM1}.

Masers in MM1 trace arc shaped shocks which are expanding from the center of the system at an average velocity of 15 km s$^{-1}$. We interpret these as bow shocks tracing the protostellar jet reported by \citet{Goddi06,Goddi11} and discuss these in further detail later. The masers associated with MM2 also appear to trace an expanding arc. However the detection of only Easterly propagating components i.e. without a symmetric counterpart in the West, precludes detailed analysis of its kinematics.

Maser feature A could not be assigned a proper motion since its motion in the sky-plane was strongly non-linear.

\begin{figure}
\begin{center}
\includegraphics[scale=0.7]{Elli_Scenarios2.eps}
\caption{\small Ellipses fit to the VERA data, assuming a common origin and indipendant ejections for N.W., S.E.i and S.E.o masers. VLBA data from 2004 \citep{Goddi06} are also shown.
\label{SINGLE}}
\end{center}
\end{figure}

\vspace{-0.1cm}

\subsection{Single vs Episodic ejection}
\label{epi}

To investigate the physical nature of ejections from MM1 we performed least squares fitting of elliptical jet shells to the maser distributions to evaluate their morphology and to determine the center position of the system. Since our maser data reveal two independent bowshocks in the S.E. lobe we must consider that ejections may occur in multiple episodes rather than a single ejection event, thus ellipses were fit under the condition that the three ejections, each traced by a bowshock, originate from a common center.
We find that the data are best fit to concentric ellipses with an origin at $(\alpha, \delta) = (33,-54)$ mas, shown as a star in Figure~\ref{SINGLE}. The center position coincides with the methanol disk reported by \citet{Goddi07} and is near to maser feature A which exhibits non-linear motions (\emph{see} Section~\ref{ALLVLBI}). The physical distances of the bowshocks from the system center were 240, 437 and 475 AU for the N.W., S.E. inner and S.E. outer masers, respectively, assuming that these bowshocks are close the plane of the sky. The shocked surfaces extend about $\sim 50$ AU in a direction perpendicular to the outflow direction, while their separations from the system center range from $240-475$ AU, indicating a high degree of collimation.

The difference in the ejection radii of the N.W. and S.E. bowshocks groups in our model can be interpreted either as separate asynchronus ejections expanding from the center at constant velocity, or as a quasi-single ejection event with non-equal propagation velocities - possibly due to a density gradient in the vicinity of the YSO. The former interpretation appears favoured considering the results of \citet{Goddi06} who find masers further in the N.W. direction than those observed with VERA (Figure~\ref{SINGLE}). These masers align well with the ellipses fit to the VERA data (eventhough the VLBA data were not used for fitting) and may thus be counterparts to the S.E. bow shocks.

Assuming expansion at constant velocity the ages of ejection events - i.e. the time it would have taken for ejecta to reach the positions of the observed bowshocks - can be estimated from their distance from the star divided by their outward motion. We measured average internal motions of 1.4 mas yr$^{-1}$ for each of the lobes, in opposite directions. For the ejections observed in AFGL 5142 MM1 we estimate ejection ages of 80, 145 and 159 years for the N.W., S.E. inner and S.E. outer bowshocks, respectively. In this scenario the S.E. bow shocks would be roughly twice the age of the N.W. bow shock, with 14 years interval between the ejections of the inner- and outer-S.E. ejections.

\section{Discussion}

\subsection{Physical properties of the jet}
\subsubsection{Outflow driving mechanism}
\label{ostriker}

Masers in AFGL 5142 MM1 trace arcs of shocked gas at the leading edge of a protostellar jet.
Using the proper motions of masers measured with VLBI we can investigate the physical properties of the jet by comparison with published outflow models. VLBI maser data are well suited for comparison with the models of \citet{Lee01} and \citet{Ostriker01} who describe the motions of gas in the shocked region between a protostellar outflow and the ambient gas. In their works they compare two models; a jet driven bowshock and a stellar wind.
Generally, the prevalence of each model can be determined by the amount of transverse motion observed near the head of the jet; in a jet driven outflow the transverse velocity is large as the jet material disperses at the bowshock while the transverse velocities in a stellar wind tend to zero at the tip since in a momentum driven stellar wind the gas propagates radially outward from the YSO, and thermal pressure is negligible.

\begin{figure*}
\begin{center}
\includegraphics[scale=1.42]{I05274_VERA_Goddi.eps}
\caption{\small Combined view of 22 GHz water masers (filled circles) observed with VERA in 2010 (\emph{this work}), 22 GHz water masers (asterisk) observed with the VLBA in 2004 \citep{Goddi06} and 6.7 GHz methanol masers (triangles) observed with the EVN in 2004 \citep{Goddi07}. The inset shows the trajectory of maser feature A, moving in a clockwise fashion. The trajectory of feature A and proper motions of other masers are all converted to the YSO frame. The black asterisk symbol indicates the approximate origin of the episodic ejections, estimated from least-squares fitting of ellipses to the VERA maser data.
\label{COMBINED}}
\end{center}
\end{figure*}

Using a jet position angle of $-35 ^{\circ}$ we measure transverse proper motions in a velocity range of $ 0 - 13$ km s$^{-1}$, with a mean value of $4$ km s$^{-1}$.
The large range of transverse motions indicate that the outflow in AFGL 5142 MM1 is dominated by jet-driven bowshock kinematics. Since masers in AFGL 5142 MM1 exclusively trace the leading surface of the outflow it is difficult to quantify any contribution from a stellar wind - whose influence is more easily seen along the body of the outflow.
Other similar works comparing VLBI determined maser proper motions with the models of \citet{Lee01} and \citet{Ostriker01} also find MYSO jets to be dominated by jet-driven bowshocks with some contribution from a lower velocity wind \citep{Sanna12,Burns16b}.

\subsubsection{Jet momentum rate}
\label{momrate}

In addition to investigating the driving mechanism, we can calculate the momentum rate (force) of the jet via the expression used by \citet{Goddi11}:

\vspace{-0.35cm}
\begin{eqnarray}
\dot{P}_{\rm outflow} = 1.5 \times 10^{-3}~ V^2_{10}~ R^2_{100}~ (\Omega / 4 \pi)~ n_8~ [{\rm M}_{\odot}~ {\rm km s}^{-1}~ {\rm yr}^{-1}].
\label{FORCE}
\end{eqnarray}

In which V$_{10}$ is the maser velocity in units of 10 km s$^{-1}$, R$_{100}$ is the radial distance of the masers from the driving source in units of 100 AU, $\Omega$ is the opening solid angle of the jet and $n_8$ is the volume density of ambient gas in units of $10^8$ cm$^{-3}$.
The following parameters used in the calculation were derived from our observations of the N.W. maser bowshock; $v=14$ km s$^{-1}$, $R = 240$ AU, and S.E. bowshock; $v=14$ km s$^{-1}$, $R = 475$ AU. Furthermore we take $\Omega = 1.1$ sr from the jet continuum observations of \citet{Goddi11}.
Regarding the ambient density, $n_8$, we consider values in the range $10^7-10^9$ cm$^{-3}$ which are simulated initial ambient densities suitable for producing masers in shocks \citep{Kauf96}. We arrive at a range of momentum rates between $\dot{P} = 10^{-4} ~{\rm to}~ 10^{-2}$ M$_{\odot}$ km s$^{-1}$ yr$^{-1}$ for both bowshocks. These ranges are consistent with \citet{Goddi11} and the extended molecular outflows in AFGL 5142 MM1 (see Table~\ref{OFP}).

If we instead derive the jet opening angle directly from the maser data we arrive at $\Omega = 0.07$ sr and momentum rate ranges of $\dot{P} = 10^{-6} ~{\rm to}~ 10^{-4}$ M$_{\odot}$ km s$^{-1}$ yr$^{-1}$. The fundamental property that maser emission arises in regions of special physical conditions, and that they are observed preferentially along surfaces tangental to the observer, can complicate maser-derived estimates of the jet opening angle - except for cases where the jet head and body are more fully traced. As such, we argue that the jet opening angle derived from the continuum observations of \citet{Goddi11} is preferred.

\subsubsection{Primary jet velocity}

When the ram pressure of the jet is sufficient to balance that of the swept up material the relationship between the velocities of the primary jet, and that of the swept up gas, is given as $v_{j} \approx v_{ws} (1+ \sqrt{n / n_{j}}) $, adapted from \citet{Chernin94}. Where $v_{j}$ and $v_{ws}$ are the velocities of the primary jet and working surface (traced by the bowshock velocity), and $n$ and $n_j$ are the ambient and jet number densities, respectively. For the ambient density we use $n = 10^8$ cm$^{-1}$, adopted from shock models \citet{Hollenbach13}. Regarding the jet density, $n_j$, \citet{Garatti15} find primary jet electron densities of $n_e \sim 10^3$ to $10^4$ cm$^{-3}$ in their survey of MYSOs. Employing an ionisation fraction of $x_e = n_e/n_H = 0.1$ \citep{Hartigan94} the MYSO primary jet densities fall in the range of $n = 10^4$ to $10^5$ cm$^{-3}$. Finally, we estimate primary jet velocities in the range of $v_j \sim 100$ to $1000$ km s$^{-1}$. The jet velocity derived here for AFGL 5142 MM1 is consistent with other values for MYSOs from the literature \citep{Marti95,Curiel06,Sanna16,Guzman16,Kamenetzky16}.

\subsubsection{bowshock velocity gradient}

As a closing remark on the observationally determined jet physical properties in AFGL 5142 MM1, we note a shallow yet well structured velocity gradient across the N.W. bowshock masers (see Figure~\ref{VELGRAD}). This could be interpreted as a flattened and expanding ejection or as  rotation of the entrained material about the jet axis - a phenomenon widely sought in the context of massive star formation as it would suggest that outflows are driven by magneto-centrifugal jets (\emph{see also} \citealt{Burns15a}). Our data show that AFGL 5142 MM1 has potential as a target for future rotating jet investigations.

\subsubsection{Combined VERA, VLBA and EVN view of AFGL 5142 MM1}
\label{ALLVLBI}
AFGL 5142 has been observed by three different VLBI arrays; the VLBA \citep{Goddi06}, the EVN \citep{Goddi07}, and VERA (this work). These are presented together in Figure~\ref{COMBINED}.

The combined view of VLBI maser observations of AFGL 5142 depicts a prototypical disk-jet system; a structure ubiquitous in low mass star formation and becoming increasingly encountered in MYSOs with well studied cases including Orion Source I \citep{Hirota14}, IRAS 20126+4104 \citep{Moscadelli11,Chen16}. Water masers trace a narrow episodic jet which emanates at an angle near perpendicular to the protostellar disk traced by 6.7 GHz methanol masers - a transition exclusive to MYSOs.

Maser feature A, detected in this work only, exhibits a curved trajectory and lies close to the likely position of the star as was inferred independently by the center of kinematics of the 6.7 GHz methanol masers \citep{Goddi07}, the center of kinematics of the VLBA water masers \citep{Goddi06} and the concentric ellipses fit to VERA water masers - presented in Section~\ref{epi}. Non-linear motion implies the application of force - whether that be gravitationally or magnetically dominated we expect that the maser must be in close proximity to the central object.
The trajectory of feature A is shown in an inset in Figure~\ref{COMBINED} where motion progresses from the upper right of the inset in a clockwise direction. The trajectory has been shifted to the frame of the star.
We pursue a full decomposition of the trajectory of this maser in a future publication targeting the central $<100$ AU region.

Overall, the features highlighted by the combination of VLBA, EVN and VERA observations reveal a system comprising a circumstellar disk from within which a narrow and collimated jet is launched. These features resemble those representative of low mass YSOs.
As we showed in Section~\ref{ostriker}, the outflow launching mechanism in AFGL 5142 MM1 also resembles those common to low-mass YSOs, while the outflow momentum rate of AFGL 5142 MM1 (Section~\ref{momrate}) is larger than those typically seen in low-mass stars \citep{Beuther02}. Episodic ejection, which in turn implies that accretion also occurs episodically, is another feature shared with low-mass disk-aided star formation. These points, taken together with the morphology and 3D kinematics of structures near the MYSO compiled from VLBI observations, reveals that AFGL 5142 MM1 depicts a `scaled up' picture of low-mass star formation.

\subsection{Molecular outflows in the AFGL 5142 region}
\label{ofs}

Previous attempts at locating the progenitors of the multiple outflows seen in AFGL 5142 struggled on account of there being more outflows (at least 4) than major millimeter cores (only two), and because most outflows intersect one or both of MM1 or MM2.
We propose that several of the outflows in this region could be explained by episodic ejections from a slowly precessing outflow system. In a recent paper \citet{Liu16} consider that a precessing jet may be the best way to produce the extremely wide-angle bipolar outflow which forms the subject of their investigation. We share this view and extend on justifications described as follows.

Outflows A and C from \citet{Zhang07}, the `compact CO outflow' of \citet{Hunter95} and the SiO and HCO$^+$ outflows of \citet{Hunter99} share a trend whereby the redshifted lobe of the outflow extends to the North, while the blueshifted lobes extend to the south. We evaluate this trend in Figure~\ref{PAVST} which shows the P.A. of the aforementioned outflows as a function of angular extent, which we take as analogous to the outflow age (more extended outflows are older). Given that the age of the more extended outflows are a few $10^4$ years old \citep{Hunter95,Hunter99}, the precession has an approximate period of one revolution per $10^4 - 10^5$ years.
For this analysis we only use the Northern, redshifted outflows since the southernly directed blueshifted outflows often appear to truncate near MM2, which is known to coincide with the densest part of the molecular envelope (see Figure 4 of \citealt{Zhang07}). Outflow details used in the analysis and references to the original works are given in Table~\ref{OFP}.

\begin{figure}
\begin{center}
\includegraphics[scale=0.68]{PA_vs_time.eps}
\caption{\small The angular extent of outflows from Table~\ref{OFP} as a function of position angle. A sinusoid represents expected values for a slowly precessing outflow, fit to the data.
\label{PAVST}}
\end{center}
\end{figure}

\begin{figure*}
\begin{center}
\includegraphics[scale=1.99]{I05274_FullScale.datSHOES.eps}
\caption{\small Overlay showing water maser activity in the context of the millimeter core and molecular outflows. Small arrow vectors are the same as in Figure~\ref{FULL}. Greyscale shows the 1.3 mm continuum emission in the range of 0 - 0.07 mJy with contours at -4, 4, 8, 12, 16, 20 and 24 times the rms noise of 2.8 mJy beam $^{-1}$, from \citet{Palau13}, with the millimeter cores labeled. The red contour indicates the $6 \sigma$ detection adopted by \citet{Palau13}. Larger arrows with letters annotated indicate the directions of molecular outflows described in previous works \citep{Zhang07,Palau13}.
\label{SHOES}}
\end{center}
\end{figure*}

Figure~\ref{PAVST} demonstrates the plausibility that the outflows listed above emanated from a single, precessing jet system - which we can confidently attribute to MM1 thanks to works at the highest angular resolution (\citealt{Goddi06,Goddi11} and this work).
To further test this hypothesis we consider the momentum rates of each outflow, which we calculate or re-evaluate from the literature using our trigonometric distance of 2.14 kpc and record in column 4 of Table~\ref{OFP}. All outflows under consideration have comparable momentum rates, supporting the hypothesis that they were driven by the same progenitor.

Regarding the remaining outflows, outflow D was likely driven by MM2 since its East-West orientation matches the cental axis of the bowshock in the F.S. masers (Figures~\ref{SHOES},\ref{FSBOW}), and also the proper motions of the F.S. and F.S.E masers (Figure~\ref{SHOES}). Note that both millimeter cores of MM1 and MM2 are elongated in the direction perpendicular to the assigned outflows, a known attribute of disk-jet systems \citep{Reid07}.

The velocity orientation and position angle of outflow B is opposite to those of outflows A and C, ruling out an association with MM1. Outflow B may be driven by MM3 which could be a low mass YSO \citep{Zhang07}. MM3 sits between the blue- and redshifted lobes of outflow B where \citet{Zhang07} detect weak 8.4 GHz continuum emission at the location of the millimeter core. The authors thus argue that MM3 may be of stellar nature, while the remaining millimeter cores, MM6, MM7, MM8 and MM9 can be explained as dust condensations swept up by molecular outflows, D, A/C, A and C, respectively.

\begin{table}
\vspace{-0.2cm}
\begin{center}
\scriptsize
\caption{\small Outflow parameters from this work and the literature. \label{OFP}}
\begin{tabular}{ccccc}
\hline
Tracer & $Length$ & P.A. & Momentum rate & Notes\\
       & $^{\prime \prime}$ & [$^{\circ}$] & [M$_{\odot}$ km s$^{-1}$ yr$^{-1}$] &\\
\hline
CO (2-1)     & 40 & 34 & $1.11 \times 10^{-4}$ & 1 \\
HCO$^+$ (1-0)& -  & -  & $4.52 \times 10^{-3}$ & 2 \\
SiO (2-1)    & 30 & 10 &                     - & 2 \\
CO (2-1)     & 30 & 5  & $2.85 \times 10^{-3}$ & `Outflow A', 3 \\
CO (2-1)     & 10  &-40 & $2.38 \times 10^{-3}$ & `Outflow C', 3, 4\\
H$_2$O $(6_{16}-5_{23})$&0.1&-35& $10^{-4} ~{\rm to}~ 10^{-2}$ & this work\\
\hline
\end{tabular}
\begin{tablenotes}
\item{Column (5)}: 1 - \citet{Hunter95}; 2 - \citet{Hunter99}; 3 - \citet{Zhang07}; 4 - \citet{Palau11}
\end{tablenotes}
\end{center}
\end{table}

\subsection{Episodic ejection in AFGL 5142}

Bowshocks trace the leading edges of new ejection events. As such, multiple bowshocks imply a history of multiple ejections. In Section~\ref{epi} we present evidence of episodic ejection operating on timescales of $10^1$ years, based on maser data. In Section~\ref{ofs} we argue that several of the molecular outflows in AFGL 5142 appear to stem from AFGL 5142 MM1.
Regarding the dynamical timescales of these molecular outflows, the older, more extended outflows from MM1 are a few $10^4$ years old \citep{Hunter95,Hunter99}, while there are outflows associated with MM1 on 4 different scales (Table~\ref{OFP}) - corresponding to ejection timescales of every few $10^4$ years.

Since outflow activity is expected to correlate with accretion activity \citep{Garatti15} the $10^4$ year timescale of ejections could indicate alternating periods of accretion bursts and quiescence in AFGL 5142 MM1, consistent with timescales predicted for episodic accretion in low- \citep{Stamatellos11} and high-mass stars \citep{Meyer17}. The short episodic ejection seen on a 10 year timescales (Section~\ref{epi}) would likely originate in the inner disk.

The concept of episodic accretion is only recently being considered in the framework of massive star formation, though its ability to deal with the "luminosity problem" and suppress the accretion inhibiting radiation field show great potential in explaining how massive stars accumulate their mass. AFGL 5142 presents a promising observational target for studies of episodic ejection and accretion in MYSOs.

\section{Conclusions}

Using VERA, a dual-beam VLBI array dedicated to astrometry, we measured an annual parallax for the AFGL 5142 massive star forming region of $\pi=0.467 \pm 0.010$ mas, corresponding to a trigonometric distance of D=$2.14^{+0.051}_{-0.049}$ kpc.

Using the latest VLBI H$_2$O maser data we analyse a recent ejection from the AFLG 5142 MM1 region which produced prototypical bi-polar bowshocks that expand from a common center in a circumstellar disk. We use the 3D kinematics of the maser bowshocks to investigate the physical properties of the protostellar jet.

One of the masers near the likely position of the star exhibits clear non-linear motion. The immediate region around the star requires further deep observations to provide vital context for this very unusual finding.

Combining our data with those of other published VLBI observations and through proper motion analysis we find that AFGL 5142 MM1 is forming as a disk-jet system with evidence of episodic ejection, giving it the appearance of `scaled up' low-mass star formation.

Finally, we reinterpret the outflow history of AFGL 5142 MM1 as formation via an episodic, precessing outflow. We were able to suggest progenitor allocations for all outflows in AFLG 5142 known from the current literature.

\section*{Acknowledgments}

R.B. would like to acknowledge the Ministry of Education, Culture, Sports, Science and Technology (MEXT) Japan for support as part of the Monbukagakusho scholarship.

We kindly thank the referee for providing helpful feedback and guidance during the review of this work.

\input{./mn-bib.tex}
\bibliographystyle{mn2e}