text:
\section{Introduction}

\subsection{The classical case}
The theory of de Branges-Rovnyak spaces of analytic functions (see for instance \cite{alpay2017beurling,dbhsaf1,dbr1,dbr2})
allows us to prove Beurling-Lax type theorem in a one complex variable framework when leaving the setting of the Hardy space.
An illustrative example is given in the following theorem.
To give the statement, we must first recall that an element of $\mathbb{C}^{n \times n}$ which is
both selfadjoint and unitary is called a signature matrix.
Furthermore, for $\alpha\in\mathbb C$ and $f$ a matrix-valued function analytic in a neighborhood of $\alpha$,
let $R_{\alpha}$ denote the resolvent operator at $\alpha$:
\begin{equation}
\label{resolventOp}
R_\alpha f(z) =
\begin{cases}
\,\,
\dfrac{f(z)-f(\alpha)}{z-\alpha},
\quad
z\not=\alpha\\
\,\,
f^\prime(\alpha),
\quad \quad \quad \quad
\,
z=\alpha.
\end{cases}
\end{equation}
It follows that the resolvent identity \cite[Theorem I]{dbhsaf1}
\[
R_{\alpha} - R_{\beta}= (\alpha - \beta)  R_{\alpha}R_{\beta}
\]
holds for any function analytic in a connected neighborhood of $\alpha$ and $\beta$.

\begin{theorem}
\label{1_1}
Let $\Omega$ be an open subset of the complex plane, symmetric with respect to the real line, and let $J$ be a signature matrix.
Let $\mathcal X$ be a reproducing kernel Hilbert space of functions analytic in $\Omega$.  Then the reproducing
kernel of $\mathcal{X}$  is of the form
\begin{equation}
\label{ktheta}
K_T(z,w)
=
\frac{J - T(z) J T(w)^*}{-i(z-\overline{w}) },
\end{equation}
where $T$ is a $\mathbb C^{n\times n}$-valued function analytic in $\Omega$,
if and only if the following two conditions are fulfilled:
\begin{enumerate}
\item $R_\alpha \mathcal X\subset\mathcal X$ for all $\alpha\in\Omega$.
\item The structure identity
\begin{equation}
\label{eq:structure1}
\innerProductReg{R_\alpha f}{g} -
\innerProductReg{f}{R_\beta g} -
(\alpha - \overline{\beta}) \innerProductReg{R_\beta f}{R_\beta g}
=
i g(\beta)^* J f(\alpha)
\end{equation}
holds for all $\alpha,\beta \in \Omega$ and $f,g \in \mathcal X$.
\end{enumerate}
\label{thm:struct}
\end{theorem}
See \cite[Theorems III and IV]{dbhsaf1}, where $\Omega\cap\mathbb R$ is assumed non-empty; this restriction was later removed in \cite{rov-66} and in
\cite{ball-contrac} (for the corresponding theorem in the case of the disk, see \cite{ad-jfa}).
\smallskip

We remark that $T$ is $J$-contractive, that is, $T$ satisfies $T(z)J T(z)^*\le J$ for all $z$ in the intersection of $\Omega$ and the upper half-plane $\mathbb C_+$.
The positivity of the kernel \eqref{ktheta} in an open subset of the upper half-plane implies that $T$
has a meromorphic extension to $\mathbb C_+$ (and in fact, by reflection, to $\mathbb C\setminus \mathbb R$),
for which $K_T(z,w)$ is still positive definite.
\smallskip

If we further assume that $\infty\in\Omega$ and $\ker R_\alpha = \{0\}$ for an arbitrary $\alpha \in \Omega$
(and henceforth, for all $\alpha\in\Omega$), then the space  does not contain nonzero constants,
and $T$ is of the form
\[
T(z) = I - iC(zI - A) ^ {-1} C^* J
\]
where
\begin{equation*}
(A \, f)(z)=zf(z)-\lim_{w\rightarrow\infty}wf(w)\quad\text{and }\quad Cf=\lim_{w\rightarrow\infty}wf(w).
\end{equation*}
Furthermore, the resolvent operator satisfies $R_{\alpha} = (A- \alpha I )^{-1}$.
\smallskip

In his fundamental paper \cite{MR0027954} Beurling introduced a characterization of
invariant subspaces under the shift operator in the Hardy space $H^2(\mathbb{D})$.
These subspaces are characterized as the ones of the form $T H^2(\mathbb{D})$, where $T$ is an inner function.
An inner function is, by definition, an analytic function on the unit disk such that $|T(z)| \leq 1$ for $|z|<1$ and
with non-tangential values (which exist almost everywhere since $T$ is bounded) having modulus one.
Important generalizations, the vector-valued case and the infinite dimensional case, were presented later by Lax \cite{MR0140931} and Halmos \cite{MR0152896}, respectively.
\smallskip

We note that \eqref{eq:structure1} is automatically satisfied in the Hardy space $\mathbf H^2(\mathbb C_+)$.
Therefore, applying Theorem \ref{1_1} (restricted to subspaces of $\mathbf H^2(\mathbb C_+)$)
on the orthogonal complement implies the Beurling-Lax theorem, under the hypothesis of a symmetric domain of analyticity.

\subsection{Beurling theorem on Riemann surfaces}
\label{sucSectionIntoCRS}

In the setting of finite bordered Riemann surfaces,
several generalizations of Beurling's theorem were presented.
Sarason, in \cite{sarason2}, studied the invariant subspaces of certain multiplication operators on ${\bf L}^2$ in the case of the annulus.
Invariant subspaces in Hardy spaces over a finite bordered Riemann surface were studied by
Voichick \cite{MR0160920,voichick1966invariant}, Voichick-Zalcman \cite{MR0183883}, Hasumi \cite{MR0190790} and Forelli \cite{MR0193534} in the scalar-valued case,
while Abrahamse and Douglas \cite{AbDo} considered the vector-valued case
(see also Abrahamse \cite{Abr1} and Ball \cite{MR526287} for the closely related interpolation problems).
A dedicated investigation of invariant subspaces in Hardy space over a multiply connected
planer domain, including corresponding operator theory and Riesz bases,
was carried out by
Fedorov and Pavlov \cite{MR1025160} and Fedorov \cite{MR1072300,MR1070485}.
In addition to lifting the problem to the universal covering,
they use extensively the representation of the domain as a
finite branched covering of the unit disk.
\smallskip

Neville \cite{MR0301206,MR0586558} and Hasumi \cite{MR0364647,MR0407283} obtained Beurling Lax type theorems
for certain classes of infinitely connected open Riemann surfaces.
The crucial condition is that any flat complex line bundle on the open Riemann suface
admits bounded holomorphic sections; such open Riemann surfaces
were charaterized completely by Widom \cite{Widom} and referred as Parreau-Widom type, see also \cite{MR0050023}.
Finally, we also mention the work of Kupin and Yuditskii \cite{MR1473264},
as they have studied interpolation problems on
infinitely connected open Riemann surfaces of Parreau-Widom type.

\subsection{Compact real Riemann surfaces}
In the present paper we prove a counterpart of Theorem \ref{thm:struct} in the setting of compact real Riemann surface
(i.e. compact Riemann surface equipped with an anti-holomorphic involution $\tau$);
for a theorem in the finite dimensional case see \cite[Theorem 5.1]{av3}.
Taking the compact real Riemann surface to be the double of a finite bordered Riemann surface $S$
we obtain, as a corollary, a Beurling Lax type characterization of invariant subspaces of a Hardy space on $S$ under
a pair of certain multiplication operators.

The counterpart of the kernel $\frac{1}{-i(z-\overline{w})}$
in the compact real Riemann surface case is given by
\begin{equation}
\label{theta111}
K_{\zeta}(z,w)
=
\frac{\vartheta[\zeta](\tauBa{w}-z)}{i\vartheta[\zeta](0)E(z,\tauBa{w})}
\end{equation}
(see \cite[Section 2.4]{av3}), while the counterpart of the kernel
$K_S(z,w)$
is given in \eqref{eqKernPos} below.
Here $\vartheta[\zeta]$ is the theta function with characteristic $\zeta$ and $E(\cdot,\cdot)$ is the prime form
(for more details see Section \ref{sub21} below or \cite{fay1}).
Note that the kernel \eqref{theta111} is not always positive definite on a subset of a given real Riemann surface.
The case where it is positive definite ($\zeta \in T_0$ and the compact real Riemann surface is of dividing type, see below) corresponds to the Hardy space case.
\smallskip

Our approach to prove de Branges type theorem (and hence Beurling Lax type theorem) is substantially different
compared to the approaches of the papers mentioned in Section \ref{sucSectionIntoCRS}.
Essentially, all of these papers are using the method of lifting the problem to the universal
covering and using the classical Beurling theorem to solve the problem.
We, on the other hand, use commutative operators vessels and operator model theory \cite{MR2043236,KLMV,MR1634421} for pairs of commuting nonselfadjoint operators
in order to build explicitly the (inner) function $T$.
The function $T$ (actually a multiplicative multi-valued function, i.e. a mapping of line bundles)
is then the transfer function of a commutative two-operator vessel.
The inner space of this vessel is the given reproducing kernel Hilbert space
(the orthogonal complement in $H^2$ to the invariant subspace in the case of Beurling Lax theorem)
and the main operators are the compressed multiplication operators by a pair of real
meromorphic functions on the Riemann surface.
\smallskip

Furthermore, all the papers mentioned in Section \ref{sucSectionIntoCRS} consider the setting of multiplicative functions.
They use a non-canonical approach in the sense it is required to choose a measure on the boundary of the Riemann surface,
resulting in more complicated calculations. Our approach (as in \cite{av2,av3}) is to consider the half-order multiplicative differentials setting.
This approach is canonical in the sense that no selection of a measure is required.

\subsection{Structure of the paper}
The paper consists of five sections besides the introduction.
In Section \ref{secPre}, which consists of preliminaries, we review some basic definitions and results related to compact real Riemann surfaces.
In addition, we survey the theory of vessels associated to pairs of commutative nonselfadjoint operators,
and the corresponding operator models.
\smallskip

The main theorem, namely the de Branges structure theorem, the counterpart of Theorem \ref{1_1} in the compact real Riemann surfaces setting, is presented in Section \ref{secMainThm}.
In Section \ref{secH2}, we specialize our results to the setting of the Hardy space and obtain versions of Beurling's theorem on finite bordered Riemann surfaces.
Section \ref{secMultOp} is dedicated to study the compressed multiplication operators  associated to a real Riemann surface and a pair of meromorphic functions.
To ease the presentation, the proof of de Branges structure theorem is given later in Section \ref{secProofStruct}.
\smallskip

{\it The authors wish to thank the referee for carefully reading the manuscript
and for raising valuable suggestions to improve the original draft.}

\section{Preliminaries}
\setcounter{equation}{0}
\label{secPre}

In the preliminaries section we review some necessary background and notions related to the content of this paper.
In the first part we survey the compact real Riemann surfaces and their Jacobians.
The second part is dedicated to the theory of vessels associated to pairs of commutative operators.
The model space associated to a contractive mapping between certain line bundles defined on a compact real Riemann surface is presented in the last part.

\subsection{Compact real Riemann surfaces}
\label{sub21}

It is a well-known fact (see \cite{livsic2,KLMV} and Section \ref{sub22} in the upcoming pages) that real algebraic curves and compact real Riemann surfaces play an important role in the theory of operators vessels.
A survey of the main needed tools (including the prime form and the Jacobian) can be found in \cite[Section 2]{av3},
the descriptions of the Jacobian variety of a real curve and the real torii is in \cite{vinnikov5}.
For general background, we refer the reader to \cite{fay1,GrHa,gunning2,mumford1,mumford2}.
\smallskip

A compact Riemann surface $X$ of genus $g$ is called {\it real}
if there exists an anti-holomorphic involution $\tau : u \rightarrow \tauBa{u}$ from $X$ into itself.
Let $X_{\mathbb R}$ be the set, assumed nonempty, of real points in $X$ (that is the set of points $x \in X$ such that $x = \tauBa{x}$).
The set $X_{\mathbb R}$ contains $k \geq 1$ disjoint connected components denoted by $X_1,...,X_{k}$.
Two cases should be distinguished. The first, called the dividing case, is when $X \setminus X_{\mathbb R}$ is not connected; then it is a union of two connected components, $X_{-}$ and $X_{+}$.
The second, the non-dividing case, is when $X \setminus X_{\mathbb R}$ is connected.
\smallskip

Repeating the constructions as have been presented in \cite{ahlfors,av3,vinnikov5}, let $A_1,...,A_g$, $B_1,...,B_g$ be the canonical basis of the homology group $H_1(X,\mathbb{Z})$ such that the complex conjugation is given by
$T = \left( \begin{smallmatrix} I& H \\ 0 &-I \end{smallmatrix} \right)$.
Let $\omega_1,...,\omega_g$ be the normalized basis of the space of holomorphic differential on $X$, where $g$ is the genus of $X$.
The period matrix $\Gamma$ is defined by $\Gamma_{ij} = \int_{A_i}{\omega_j}$, and $\Gamma$ is symmetric and satisfies ${\rm Im}\, \Gamma > 0$.
It is convenient to use the following decomposition $ \Gamma = \frac{1}{2}H + i Y^{-1}$.
\smallskip

One associates to $X$ the Jacobian variety $J(X) = \mathbb{C}^g/\Lambda$,
where $\Lambda$ is a lattice defined by $\Lambda = \mathbb Z^g+\Gamma\mathbb Z^g$.
The Abel-Jacobi mapping from $X$ to $J(X)$ defined by
\[
\mu : p \longrightarrow
\left(
    \begin{array}{c}
      \int_{p_0}^{p}{\omega_1} \\
     \vdots \\
      \int_{p_0}^{p}{\omega_g} \\
          \end{array}
  \right),
\]
is well-defined for an arbitrary base point $p_0 \in X$.
We note that $\Lambda$ is invariant under $\tau$ and hence the complex conjugation is defined on $J(X)$.
Furthermore, we always choose $p_0 \in X_\mathbb R$ and it follows that $\tau$ and the complex conjugation
are invariant under the Abel-Jacobi map, i.e. $\mu(\tauBa{p}) = \overline{\mu(p)}$.
The corresponding theta function is given by
\[
\vartheta(\lambda) =
\sum_{n \in \mathbb{Z}^g}{\exp\bracketsA{i \pi n ^ t \Gamma n + 2 i \pi n ^ t \lambda}},
\]
and is a quasi-periodic function with respect to the lattice $\Lambda$,
that is,
\[
\vartheta (\lambda + m) = \vartheta (\lambda)
\quad {\rm and } \quad
\vartheta (\lambda + \Gamma n) = \exp \bracketsA{ -i \pi n ^ t \Gamma n - 2 \pi i n \lambda} \vartheta (\lambda),
\]
where $n,m \in \mathbb Z ^g$.
Therefore, the theta functions define divisors in $J(X)$.
The theta function with characteristic $a$ and $b$ in $\mathbb R^g$ is defined by
\[
\vartheta\left[ \!
\begin{array}{c}        a \\       b     \end{array} \!
\right] (\lambda) =
\sum_{n \in \mathbb{Z}^g}
{\exp{ \bracketsA{i \pi (n + a) ^ t \Gamma (n + a) + 2 i \pi (n + a) ^ t (\lambda+b)}}}.
\]
In this paper (as in \cite{av2,av3}), we consider the framework of multiplicative half-order differentials.
In order to construct and define the half-order differentials, we choose an atlas $\bracketsB{V_j , z_j}{j\in \mathcal J}$ on $X$,
for which every nonempty intersection is assumed to be simply connected.
Then there exists a family of analytic square-roots (see \cite{hawley-schiffer}), $\bracketsB{\sqrt{dz_j/dz_i}}{i,j\in \mathcal J}$, such that the following cocycle condition
\begin{equation}
\label{cocycleCond}
\sqrt{\frac{dz_i}{dz_j}}
=
\sqrt{\frac{dz_i}{dz_k}}
\sqrt{\frac{dz_k}{dz_j}}
\end{equation}
holds on $V_j \cap V_i \cap V_k $ (whenever the intersection is not empty).
Among the line bundles defined by \eqref{cocycleCond}, the line bundle corresponds to $-\kappa \in J(X)$ (where $\kappa$ is the Riemann constant)
plays an important role and is denoted by $\Delta$.
The sections of $\Delta$ are referred as half-ordered differentials.
\smallskip

For $\zeta \in J(X)$ we define the following multipliers over the cycles $A_j$ and $B_j$ by
\[
\chi (A_j) = \exp \bracketsA{-2 \pi i a_j},
\qquad
\chi (B_j) = \exp \bracketsA{-2 \pi i b_j}
\]
where $j=1,...,g$ and $\zeta = b+ \Gamma a$ where $a,b \in \mathbb R^g$.
The corresponding unitary flat line bundle is denoted by $L_\zeta$.
\smallskip

A multiplicative half-order differential corresponding to $\zeta$
is a section of $L_\zeta \otimes \Delta$.
Thus, a multiplicative half-order differential
is a family of functions $\bracketsB{f_j}{j \in \mathcal J}$,
such that
\[
f_i(u) = \sqrt{\frac{dz_j}{dz_i}}f_j(u), \qquad u \in \widetilde{V}_i \cap \widetilde{V}_j.
\]
and
\[
f_j(u_2) = f_j(u_1) \exp \bracketsA{2 \pi i (b^tm - a^t n)},
\]
where $u_1,u_2$ are elements on the universal covering such that $\widetilde{\mu}(u_2)-\widetilde{\mu}(u_1) = n + \Gamma m$
($\widetilde{\mu}: \widetilde{X} \rightarrow \mathbb C ^g$ is the lifting of $\mu$ to the universal coverings and $n,m \in \mathbb Z$).
We always assume that $\vartheta(\zeta) \neq 0$ which is equivalent to the fact that
$L_\zeta \otimes \Delta$ does not contain global holomorphic sections
i.e. there is no non-trivial multiplicative half-order differential which is globally holomorphic.
\smallskip

The Cauchy kernel plays important role in this framework as it is the analogue in the compact real Riemann surface to
the kernel $\frac{1}{-i(z-\overline{w})}$ in the upper half--plane.
In the line bundle case, unlike the vector bundle case, the Cauchy kernel can be described explicitly by (see \cite{av3})
\begin{equation}
\label{Cauchy}
K_{\zeta}(u,v)
=
\frac
{\vartheta[\zeta](\overline{v}-u)}
{i\vartheta[\zeta](0)E(u,\overline{v})}.
\end{equation}
Here we identify the points $u,v \in X$ with the images $\mu(v),\mu(u)$ in $J(X)$.
The prime form $E(u,v)$ is a multiplicative differential of order $-\frac{1}{2}$ in each of the variables $u$ and $v$.
It is defined by
\[
E(u,v) = \frac{\vartheta[\delta](v-u)}{\sqrt{\xi_\delta(u)}\sqrt{\xi_\delta(v)}}.
\]
Here $\delta$ is a non-singular odd half-period and
$\xi_\delta$ is defined by $(\sqrt{\xi^2_\delta(u)})^2 = \sum_{j=1}^g{\frac{\partial \vartheta [ \delta ]}{\partial z_j}  (0)w_j(u)    } $,
for more details see \cite{fay1}.

Its main property is that $E(u,v) = 0$ if and only if $u = v$,
and thus can be considered as the analogue for the compact Riemann surface case of the difference between two numbers in $\mathbb C$.
In terms of local coordinates, the prime form satsfies
\[
E(u,v) = (t(v)-t(u)) + o((t(v)-t(u))^2).
\]
For more details on the prime form, we refer to \cite[section 2.3]{av3} and \cite{fay1,mumford2}. \smallskip

The Cauchy kernel is Hermitian (see \cite[Proposition 2.8]{av3} and \cite{vinnikov5})
as long as $\zeta$ belongs to disjoint union of the $g$-dimensional real torii, which is given by:
\begin{align}
    \nonumber
T_{\nu}
=
\{
\zeta \in J(X) : \zeta
= &
\frac{1}{4} \text{diag}(H) +
\frac{\nu_1}{2}e_{g-k+2}
+
\ldots
+
\frac{\nu_{k-1}}{2}e_{g}
+
\\
&
i a_1 \mathfrak{Im} ~ \Gamma_1 + \cdots + i a_g \mathfrak{Im} ~ \Gamma_g
\},
\label{eqRealTorii}
\end{align}
where $\Gamma_1,...,\Gamma_g$ are the columns of the period matrix $\Gamma$,
$a_1,...,a_{g-k+1} \in \mathbb R / 2 \mathbb Z$,
$a_{g-k+2},...,a_{g} \in \mathbb R / \mathbb Z$
and $\nu \in \{ 0,1\}^{k-1}$.
Furthermore, the kernel \eqref{Cauchy} is positive on $X_+$ and negative on $X_-$ if $X$ is dividing and $\zeta \in T_0$ (\cite[Theorem 2.1]{av3}),
and hence automatically satisfies $\vartheta(\zeta) \neq 0$.
\smallskip

The Cauchy kernels satisfy an important identity, which is referred as the \textit{collection formula} (see \cite{av2} and \cite{av3}) and is used repeatedly in the sequel.
First, it is convenient to define, using the notations from \cite{av2}, the following matrices.
\begin{definition}
Let $y$ be a meromorphic function on $X$ of degree $n$ with simple poles $(p^{(j)})_{j=1}^{n}$
and residues $(-c_j)_{j=1}^{n}$ given in terms of some fixed local coordinates at the poles.
Then for $\lambda_1,\lambda_2 \in \mathbb{C}$, we set
\begin{align}
\label{formulaK1}
\mathbb{K}(\lambda_1,\lambda_2)
& =
(\lambda_1 - \lambda_2)
\left(
\frac{1}{\sqrt{dy}(u^{(i)})}
K_\zeta(u^{(i)},v^{(j)})
\frac{1}{\sqrt{dy}(v^{(j)})}
\right)_{i,j=1}^n
\\
\label{formulaK2BLA}
\mathbb{K}(\lambda_1,\infty)
& =
-
\left(
\frac{1}{\sqrt{dy}(u^{(i)})} K_\zeta(u^{(i)},p^{(j)}) \frac{\sqrt{c_j}}{\sqrt{dt_j}(p^{(j)})}
\right)_{i,j=1}^n
\\
\label{formulaK3}
\mathbb{K}(\infty,\lambda_1)
& =
\left(
\frac{\sqrt{c_i}}{\sqrt{dt_i}(p^{(i)})}
K_\zeta(p^{(i)},u^{(j)})
\frac{1}{\sqrt{dy}(u^{(j)})}
\right)_{i,j=1}^{n}
\end{align}
\end{definition}
where $(u^{(i)})_{j=1}^{n}$ and $(v^{(j)})_{j=1}^{n}$ are the $n$ distinct pre-images in $X$ of
$\lambda_1$ and $\lambda_2$, respectively.
We set $\mathbb{K}(\lambda_1, \lambda_2)=I$ for $\lambda_1 = \lambda_2$ by continuity.
\smallskip

The matrices \refMult{formulaK1}{formulaK3} satisfy the following relations, also known as the
collection formulas (\cite[Section 4]{av2}):
\begin{align}
\label{CollFormLimitVer}
\mathbb{K}(\lambda_1,\infty) \mathbb{K}(\infty,\lambda_2) & =  \mathbb{K}(\lambda_1,\lambda_2)
\\
\label{CollFormLimitVer2BLA}
\mathbb{K}(\lambda_1,\lambda_3) \mathbb{K}(\lambda_3,\lambda_2) & =  \mathbb{K}(\lambda_1,\lambda_2)
\\
\label{CollFormLimitVer3BLA}
\mathbb{K}(\lambda,\lambda) & = I
.
\end{align}
An additional version of the collection formula can be found in \cite[Lemma 4.1]{av3}.

\subsection{Commutative Vessels over Hilbert space}
\label{sub22}
It is a well-known fact that the best way to study a nonselfadjoint bounded operator is to view it as an element of an underlying \textit{colligation} rather than studying
directly the operator itself. There is deep connection
between invariant subspaces of such operator and factorization of the colligation characteristic function (see \cite{MR20:7221}).
As soon as we consider several commuting nonselfadjoint operators, the colligation does not carry enough structure of the interaction between the operators.
It seems, see \cite{livsic1,KLMV,vinnikov4,MR1634421}, that the appropriate framework to study a pair of commuting nonselfadjoint operators is via the notion of
commutative two--operator vessel. In this paper we consider the case of a pair of commuting operators,
although a generalization to the case of $n$-tuple of commuting operators does exist, see \cite{KLMV,MR3609238}.
\smallskip

A commutative two-operator vessel is a collection
\begin{equation*}
\mathcal{V} =
\bracketsA{A_1 \, , \, A_2 \, ; \, H \, , \, \Phi \, ,  \, E \, ; \,  \sigma_1 \, , \, \sigma_2 \, , \, \gamma \, , \, \widetilde{\gamma}}
\end{equation*}
where $H$ ("the inner space")  and $ E$ ("the external space") are Hilbert spaces, and ${\rm dim}\,  E < \infty$.
The operators $A_1$ and $A_2$ are bounded in $H$, commute ($A_1A_2=A_2A_1$) and with finite
non-Hermitian rank, i.e.
\begin{equation}
\label{collCond}
\frac{1}{i} \left(A_k - A_k^* \right) = \Phi^* \sigma_k \Phi,\quad k=1,2,
\end{equation}

where
$\sigma_1$ and $\sigma_2$ are selfadjoint operators in $E$.
Finally,
$\gamma$ and $\widetilde{\gamma}$ are bounded selfadjoint operators on $E$
which satisfy the {\it vessel conditions}:
\begin{align}
\label{VesselCond1}
\sigma_1 \Phi A_2^* - \sigma_2 \Phi A_1^* & = \gamma \Phi,
\\
\label{VesselCond2BLA}
\sigma_1 \Phi A_2 - \sigma_2 \Phi A_1 & = \widetilde{\gamma} \Phi,
\\
\label{VesselCond3}
i(\sigma_1 \Phi \Phi^* \sigma_2 - \sigma_2 \Phi \Phi^* \sigma_1)
& =
\widetilde{\gamma} - \gamma.
\end{align}
Equations \refMult{VesselCond1}{VesselCond3} also called the input, output and linkage vessel conditions, respectively.
The \textit{complete characteristic function} (CCF) of a vessel is defined by \cite[Section 3.4]{KLMV},
\begin{equation*}
W(\xi_1,\xi_2,z) = I -i \Phi(\xi_1 A_1 + \xi_2 A_2 -zI)^{-1} \Phi ^* (\xi_1 \sigma _1 + \xi_2 \sigma_2)
\end{equation*}
and satisfies the following metric properties:
\begin{align*}
W(\xi_1,\xi_2,z)
\left( \xi \sigma \right)
W(\xi_1,\xi_2,z)^*
 & =
\left( \xi \sigma \right)
\text{ where {\rm Im}($z$)}=0
\\
W(\xi_1,\xi_2,z)
\left( \xi \sigma \right)
W(\xi_1,\xi_2,z)^*
 & \geq
\left( \xi \sigma \right)
\text{ where {\rm Im}($z$)}>0,
\end{align*}
where $\left( \xi \sigma \right)$ stands for $\xi_1\sigma_1 +\xi_2 \sigma_2$.
The complete characteristic function, for fixed $\xi_1$ and $\xi_2$,
is analytic for all $z$ lie outside the spectrum of $\xi_1 A_1 + \xi_2 A_2$.
\smallskip

Given a commutative vessel, we define a polynomial (that we assume not identically zero) of two complex variables,
called the \textit{(affine) discriminant polynomial}, by
\[
p(\lambda_1,\lambda_2) = \det (\lambda_1 \sigma_2 - \lambda_2 \sigma_1 + \gamma).
\]
The associated real (affine) plane curve $C_0$ is called the \textit{(affine) discriminant curve} associated to $\mathcal V$.
Writing the discriminant polynomial in homogeneous form
\[
P(\nu_0,\nu_1,\nu_2) = \det (\nu_1 \sigma_2 - \nu_2 \sigma_1 + \nu_0 \gamma),
\]
leads to a corresponding real algebraic curve $C$ in $\mathbb P ^2$.
We always assume that $\det{(\xi_1 \sigma_1 + \xi_2 \sigma_2)}$ is not identically zero,
so $C$ is the projective closure of $C_0$.
\smallskip

In order to continue, we set the principal subspace $\widehat{ H} \subseteq H$ of a vessel to be
\[
\widehat{ H} =
\bigvee _{m_1,m_2 = 0} ^{\infty}{A_1^{m_1}A_2^{m_2}\Phi^*(E)}=
\bigvee _{m_1,m_2 = 0} ^{\infty}{A_1^{*m_1}A_2^{*m_2}\Phi^*(E)},
\]
and a vessel is said to be {\it irreducible} when $\widehat{H}=H$.
\smallskip

Two vessels
\[
\bracketsA{A_1 ,A_2 \; H ,\Phi ,E;  \sigma_1 ,\sigma_2 , \gamma , \widetilde{\gamma}}
\quad {\rm and} \quad
\bracketsA{A_1^\prime ,A_2^\prime \; H^\prime ,\Phi^\prime ,E;  \sigma_1 ,\sigma_2 , \gamma , \widetilde{\gamma}}
\]
are said to be {\it unitary equivalent} if they share common external data
(that is, both have a common external space $E$ and matrices $\sigma_1,\sigma_2, \gamma $ and $ \widetilde{\gamma}$)
and there exists a unitary operator $U$ from
${ H}$ to ${ H} ^\prime$ such that
\begin{equation}
\label{eqVesselEq}
A_1^\prime   U =  U A_1  ,
\quad
A_2^\prime  U =  U A_2
\quad
{\rm and }
\quad
\Phi^\prime U = \Phi .
\end{equation}

The following theorems were presented by Liv\v{s}ic (see for instance \cite{livsic1,KLMV}).

\begin{theorem}[{\cite[Theorem 2]{livsic3}}]
The polynomial $ p(A_1,A_2)$ vanishes on the principal subspace of the colligation
$\widehat{H}$.
\end{theorem}
As a consequence, one may notice, that the joint spectrum of $A_1$ and $A_2$ restricted to $\widehat{H}$ lies on the curve $C_0$.
The discriminant polynomial can be described in terms of the selfadjoint matrices $\sigma_1,\sigma_2$ and $\widetilde{\gamma}$.
\begin{theorem}[{\cite[Corollary 1]{livsic2}}]
The following equality holds:
\[
\det (\lambda_1 \sigma_2 - \lambda_2 \sigma_1 + \gamma)
=
\det (\lambda_1 \sigma_2 - \lambda_2 \sigma_1 + \widetilde{\gamma}).
\]
\end{theorem}
Hence
$\nu_1 \sigma_2 - \nu_2 \sigma_1 + \nu_0\gamma$
and
$\nu_1 \sigma_2 - \nu_2 \sigma_1 + \nu_0\widetilde{\gamma}$
are (the input and output, respectively) determinantal representations of the discriminant curve $C$.
Then for each $\nu = (\nu_0,\nu_1,\nu_2)\in C$ we define non-trivial subspaces of $E$ by:
\begin{align*}
\mathcal{E}(\nu) & = \ker (\nu_1 \sigma_2 - \nu_2 \sigma_1 +\nu_0\gamma)\\
\mathcal{\widetilde{E}}(\nu) & = \ker (\nu_1 \sigma_2 - \nu_2 \sigma_1 +\nu_0\widetilde{\gamma}).
\end{align*}
The complete characteristic function is a function of two independent variables
and does not admit satisfactory factorization properties as in the single-operator colligation case.
However, when restricting the CCF to $\mathcal E (\lambda)$ it fits to the geometrical descriptions of the system associated to the vessel.
\smallskip

The \textit{joint characteristic function} (JCF) \cite[Section 10.3]{KLMV},
\begin{equation*}
S(\lambda)
= \restr{W(\xi_1,\xi_2,\xi_1 \lambda_1 + \xi_2 \lambda_2)}{\mathcal E (\lambda)},
\end{equation*}
determines a map between the subspaces $\mathcal{E}(\lambda)$ and $\widetilde{\mathcal{E}}(\lambda)$.
Furthermore, it is independent of the choice of $\xi_1$ and $\xi_2$ when
$\xi_1 \lambda_1 + \xi_2 \lambda_2$  does not belong to the spectrum of $\xi_1 A_1+ \xi_2 A_2$.
As a system theory interpretation, the joint characteristic function determines the input-output relation
of the corresponding two-dimensional system by $y_0 = S(\lambda) u_0$, where $y_0 \in \widetilde{\mathcal{E}}(\lambda)$
and $u_0 \in \mathcal{E}(\lambda)$ are the amplitudes of the double periodic wave functions
with frequency $(\lambda_1,\lambda_2) \in C_0$ (see \cite{MR97m:30051,MR1634421}).
\smallskip

In order to continue, we assume that the discriminant polynomial $P(\nu)$, defining the discriminant curve $C$,
is an irreducible homogeneous polynomial and we denote by $X$ the desingularizing Riemann surface of $C$.
Therefore, in this paper, we consider the line bundle setting.
It can be shown that the $\dim \widetilde{\mathcal E}(\nu) \leq s$ where $s$ is the multiplicity of $\nu$ in $C$, see \cite[Prop. 10.5.1]{KLMV}.
Therefore $\dim \widetilde{\mathcal E}(\nu)=\dim{\mathcal E}(\nu)=1$ for all non-singular $\nu \in c$.
Therefore, in this paper, as opposed to a more general vector bundle setting, we consider the line bundle setting.
Hence, we may use the theory of theta functions and Jacobian varieties
to express concretely the Cauchy kernels, see Section \ref{sub21} above
(for more details, see \cite{MR97m:30051,MR1704479}).
We assume furthermore that the input and the output determinantal representations are maximal and fully saturated
as we define in the next several paragraphs allowing ourself to consider line bundles on the
desingularizing Riemann surface rather than line bundles on the singular algebraic curve,
see \cite{MR2964618} in setting of determinantal representation
and \cite{MR918564,MR0024985} in general.
\smallskip

We say that the output determinantal representation $\nu_1 \sigma_2 - \nu_2 \sigma_1 + \nu_0 \widetilde{\gamma}$ is maximal
if any point $\nu \in C$ satisfies $\dim \widetilde {\mathcal E} (\nu) = s$.
A stronger condition is given below but we first recall that (see \cite{MR2962792}) for
an $m \times m$ matrix with entries containing homogeneous polynomial of degree $m-1$ we define $V = {\mathrm adj}[U]$,
i.e. the matrix whose $(k,j)$ entry is $(-1)^{k+j}$ times the determinant of the
matrix obtained from $U$ by removing the $j$-th column and the $k$-th row.
Then we say that the output determinantal representation is {\it fully saturated}
if all the entries of the adjoint matrix $V$ vanish on the adjoint divisor:
$(V_{i,j}) \geq D_{\mathrm sing}$ for all $i, j = 1, . . . , m$.
Here the divisor of singularities is given by $D_{\mathrm sing} = (\omega)$, where $\omega$ is a fixed holomorphic differential defined by
$\omega = \frac{d \lambda_1}{\sfrac{\partial p}{\partial \lambda_2}} =- \frac{d \lambda_2}{\sfrac{\partial p}{\partial \lambda_1}}$
(i.e. $D_{\mathrm sing}$ has the poles of $\omega$ and zero at infinity of order $m-3$ at most).
As we will see below, under certain conditions, maximality of the output determinantal representation
implies the maximality of the input determinantal representation.
\smallskip

As follows from the theory of determinantal representations (see \cite{MR97m:30051,MR2043236,MR1634421}),
there exist isomorphisms (up to some twists) between the kernel bundles $\mathcal E$ and $\widetilde{\mathcal E}$
and a certain vector bundles of the form
$L_{\zeta} \otimes \Delta$ and $L_{\widetilde{\zeta}} \otimes \Delta$.
This isomorphisms can be given explicitly in terms of normalized sections of $\mathcal E$ and $\widetilde{\mathcal E}$,
denoted by ${\bf u}^{\times}$ and $\widetilde{\bf u}^{\times} $, respectively, as follows (see \cite[Equation 2-28]{MR1634421}):
\begin{equation}
\label{eqIsoExplicit}
f(p) \rightarrow \frac{1}{\omega(p)} f(p) \widetilde{\bf u}^\times(p).
\end{equation}
A normalized section ${\bf u}^{\times}$ is a meromorphic section of $\mathcal E \otimes \Delta$
with poles at the points of $C$ at infinity
and for each $p\in X$, ${\bf u}^{\times}(p)$ form a basis to $(\mathcal E \otimes \Delta)(p)$.
In the scalar case a normalized section can be explicitly described by the Cauchy kernels
(and the its derivatives in the case where the poles are not simple).
\smallskip

Under these isomorphisms, the joint characteristic function is translated into a mapping between flat vector bundles.
It is holomorphic in a neighborhood of the points of $C$ at infinity.
This mapping is called the \textit{Normalized joint characteristic function} (NJCF)
and is denoted in the sequel by $T(p)$.
An explicit relation between the two characteristic functions is given by (see \cite[Section 6]{MR1704479}):
\begin{equation}
\label{st}
S(p) {\bf u}^{\times}(p) = \widetilde{\bf u}^{\times}(p) T(p).
\end{equation}

We denote the left kernel bundles by $\mathcal E_l$ and $\widetilde{\mathcal E}_l$
and the associated matrices of normalized sections are denoted by
${\bf u}_{l} ^ {\times}$ and $\widetilde{\bf u}_{l} ^{\times} $, respectively.
The vectors ${\bf u}_{l} ^ {\times}(p)$ and $\widetilde{\bf u}_{l} ^{\times}(p) $ are line-vectors
containing the Cauchy kernels evaluated at $p$ and the $n$ poles of $y_1$ and $y_2$.
We refer the reader to \cite{MR1704479} for more details. There exists an isomorphism between the dual bundles of
$\mathcal E$ and $\widetilde{\mathcal E}$ and $\mathcal E_l$ and $\widetilde{\mathcal E}_l$
up to an appropriate twist
(more explicitly, $ \mathcal E_l \otimes \mathcal O (n-2)(-D_{\rm sing}) = (L_{\zeta} \otimes \Delta ) ^\vee  \otimes K$
where $D_{\rm sing}$ is the divisor of singularities of $C$ and $K$ is the canonical line bundle of $X$).
\smallskip

The following fundamental realization question arises, whether for a given line bundle mapping $T$ there exists a commutative vessel,
such that $T$ is its NJCF.
A solution to this question is given, in the line bundle setting, by the following realization theorem.
\begin{theorem}[{\cite[Theorem 2.3]{MR1634421}}, {\cite[Theorem 11.2.4]{KLMV}} and {\cite[Lemma 3.3]{MR1704479}}]
\label{realizationTh}
A multiplicative function $T(p)$ on $X$ with multipliers corresponding to $\zeta,\widetilde{\zeta} \in T_\nu$ is the normalized joint characteristic function of a
commutative two-operator vessel with discriminant polynomial $p(\lambda_1,\lambda_2)$ and maximal input and output determinantal representations corresponding to
$\zeta,\widetilde{\zeta} \in J(X)$
if and only if
$T(p)$ is holomorphic at the point of $C$ at infinity, meromorphic on $X \setminus X_{\mathbb R}$ and satisfies
$ T(p)\overline{T(\tauBa{p})}=1$ and the kernel
\begin{equation}
\label{eqKernPos}
\frac{\vartheta[\widetilde{\zeta}](\tauBa{q}-p)}{i \vartheta[\widetilde{\zeta}](0)E(p,\tauBa{q})}
-
T(p)
\frac{\vartheta[\zeta](\tauBa{q}-p)}{i \vartheta[\zeta](0)E(p,\tauBa{q})}
\overline{T(q)}
,
\end{equation}
is positive definite in $X \setminus X_{\mathbb R}$.
\end{theorem}

\begin{proposition}
\label{propCont}
Let $X$ be a compact Riemann surface of dividing type and let $\zeta,\widetilde{\zeta} \in T_0$.
Then $T$ is $(\zeta,\widetilde{\zeta})$-contractive (that is, the kernel \eqref{eqKernPos} is positive)
if and only if $T$ is contractive (i.e $\abs{T(p)}^2 \leq 1$).
\end{proposition}
\begin{pf}
Assume that the kernel \eqref{eqKernPos} is positive and let us choose $p=q\in X_+$.
Then, using the positivity of the Cauchy kernel on $X_+$, \eqref{eqKernPos} becomes
\[
\frac{\vartheta[\zeta](0)}{\vartheta[\widetilde{\zeta}](0)}
\frac{\vartheta[\widetilde{\zeta}](\tauBa{p}-p)}{\vartheta[\zeta](\tauBa{p}-p)}
\geq
\abs{T(p)}^2
\]
and then when one considers the limit as $p \rightarrow q_0 \in X_\mathbb R$.
Assuming $q_0 \in X_k$, $\tauBa{q_0}-q_0$ is equal to the $k$-th column of $\Gamma$.
Since $\zeta,\widetilde{\zeta} \in T_0$ and by the properties of the theta function, the limit of the left hand side is equal to one.
Hence, for $\varepsilon>0$, an element $p$ in $X_+$ sufficiently close to the boundary $X_\mathbb R$ satisfy $\abs{T(p)}^2 \leq 1 +\varepsilon$.
It remains to apply the maximum principle to conclude that $\abs{T(p)}^2 \leq 1$ for all $p \in X_+$.
\smallskip

Conversely, let us assume that $\abs{T(p)}^2 \leq 1$. Then $T$ is holomorphic on $X_+$ and the multiplication by $T$ is a contraction.
As a consequence, the result follows immediately since the Cauchy kernel
\eqref{Cauchy} is the reproducing kernel of the Hardy space $H^2(L_\zeta \otimes \Delta , X_+)$.
\end{pf}

The interplay between the joint characteristic function to the normalized joint characteristic function
can be described as follows.

\begin{Rk}
Let $T$ be the normalized joint characteristic function of a vessel $\mathcal V$, a mapping from $L_{\zeta}$ to $L_{\widetilde{\zeta}}$.
We associate the canonical determinantal representation and the matrix of normalized sections $\widetilde{\bf u}^{\times}(z)$ to the output bundle
and to the input bundle a (non-canonical) determinantal representation and ${\bf u}^{^\prime \times}(z)$.
Then, adjusting the input determinantal representation to a canonical representation
by $\alpha , \beta \in \mathbb C^{n \times n}$ according to
$\alpha (\lambda_1 \sigma_2 - \lambda_2 \sigma_1 + \lambda^\prime) \beta =  (\lambda_1 \sigma_2 - \lambda_2 \sigma_1 + \lambda)$
and furthermore consider the relation in \eqref{st} evaluated at infinity, lead to
\[
T(p^{(j)})
=
\widetilde{\bf u}^{\times}(p^{(j)})T(p^{(j)})
=
S(p^{(j)}) {\bf u}^{\prime \times}(p^{(j)})
=
S(p^{(j)}) \alpha {\bf u}^{ \times} (p^{(j)})
=
\alpha
.
\]
Hence, the adjustment is done by the multiplication of the values of $T$ at the points at infinity.
\end{Rk}

\subsection{The model space}
\label{secModelSpace}
We continue the discussion about commutative vessels,
by presenting the model space associated to a commutative two--operator vessel, see also \cite{av1,av3}.
For the analogue in the single operator case, the de Branges-Rovnyak operator model, see \cite{bc,dbr1,dbr2} and \cite{nf}.
\smallskip

Let $\zeta,\widetilde{\zeta} \in J(X)$ satisfying $\vartheta(\zeta) \neq 0$ and $\vartheta(\widetilde{\zeta}) \neq 0$
such that $\zeta$ and $\widetilde{\zeta}$ belong to the same real torii $T_\nu$
(or, equivalently, if $\zeta - \widetilde{\zeta} = \overline{\kappa} - \kappa$, where $\kappa$ is the Riemann constant)
and let $T(p)$ be a multiplicative function on $X$ with multiplier corresponding to $\zeta$ and $\widetilde{\zeta}$.

(or, equivalently, if $\zeta - \widetilde{\zeta} = \overline{\kappa} - \kappa$, where $\kappa$ is the Riemann constant)

Then $T(p)$ is called $(\zeta,\widetilde{\zeta})$-contractive if $T(p) \overline{T(\tauBa{p})} = 1$ and
the kernel \eqref{eqKernPos} is positive.
Then, the kernel \eqref{eqKernPos} defines a reproducing kernel Hilbert space denoted by $\mathcal H (T)$.
Its elements are sections of a line bundle of $L_{\widetilde{\zeta}} \otimes \Delta$ that are holomorphic on the domain of analyticity of $T$.
\smallskip

We continue with the definition of the operator model in the case of simple poles \cite[Equation 3-3]{MR1634421}
(see \cite{av3} and in Section \ref{secMultOp} below for the non-simple poles definition),
defined on the sections of the line bundle $L_{\widetilde{\zeta}}\otimes \Delta$ analytic in neighborhoods of the poles of $y$.
It is given by
\begin{equation}
M^{y}f(u)
\label{m_y}
=
y(u)f(u) + \sum_{m=1}^{n}{c_m f(p^{(m)})
\frac
{\vartheta[\widetilde{\zeta}](p^{(m)}-u)}
{\vartheta[\widetilde{\zeta}] (0)E(p^{(m)},u)}
}.
\end{equation}
where $y$ is a meromorphic function on $X$ of degree $n$ with distinct simple poles $p^{(1)},...,p^{(n)}$
and the corresponding residues (up to a sign) $c_1,...,c_n$.
Furthermore, for any pair of meromorphic functions $y_1,y_2 \in \mathcal M (X)$, the operators $M^{y_1}$ and $M^{y_2}$ commute
(see \cite{av3} and Section \ref{secMultOp} below for the theorem in the case of non-simple poles).
\smallskip

The counterpart of the resolvent operator \eqref{resolventOp}, denoted by $R_{\alpha}^{y}$, satisfies (in the neighborhood of infinity)
$R_{\alpha}^{y} = \left(M^{y} - \alpha I \right)^{-1}$ and is defined, for an $\alpha$ with $n$ distinct pre-images,
by (see \cite[Equation 3-4]{MR1634421})
\begin{equation}
\label{Ralpha}
R_{\alpha}^{y}f(u) =
\frac{f(u)}{y(u) - \alpha} -
\sum_{j=1}^{n}{\frac{f(u^{(j)})}{dy(u^{(j)})}\frac{\vartheta[\widetilde{\zeta}](u^{(j)}-u)}{\vartheta[\widetilde{\zeta}](0)E(u^{(j)},u)}},
\end{equation}
where $y(u^{(j)})=\alpha$ for $j=1,...,n$.
$R_{\alpha}^{y}$ is indeed the resolvent of $M^y$ since $R_{\alpha}^{y} = M^{\frac{1}{y-\alpha}}$ and Theorem \ref{rkModelOpAlg} below
(in the case that $y$ has only simple poles a direct calculation can be done).
The operators $R_{\alpha}^{y}$ and $M^y$ are independent of the choice of the local coordinates of the poles of $y$,
this point is discussed in Section \ref{secMultOp}.
Furthermore, the resolvent operators satisfy the resolvent identity (see \cite[Theorem 4.2]{av3})
\begin{equation}
\label{resId}
R^y_\alpha -R^y_\beta
=
(\alpha -\beta) R^y_\alpha R^y_\beta.
\end{equation}
The model (commutative two-operator) vessel, corresponding to contractive mapping $T$ and a pair of meromorphic functions generating $\mathcal M (X)$,
is the collection (see \cite[Theorem 3.1]{MR1634421}):
\begin{equation}
\label{eqModelVess}
\mathcal V_T =
\bracketsA{M^{y_1} \, , \, M^{y_2} \, ; \, \mathcal H(T) \, , \, \Phi \, ,  \, \mathbb C ^n \, ; \,  \sigma_1 \, , \, \sigma_2 \, , \, \gamma \, , \, \widetilde{\gamma}},
\end{equation}
where $\Phi$ is the evaluation operator at the poles of $y_1$ and $y_2$
(for the case of simple poles, for the case of high order poles, one should also take the derivatives, see Section \ref{secMultOp})
and $\sigma_1,\sigma_2,\gamma,\widetilde{\gamma}$ belong to $\mathbb C^{n \times n}$.
The mapping between the inner space $H$ of a vessel to its model space
does not depend on the direction $(\xi_1,\xi_2)$ and is given by,
see \cite[Equation 3-5]{MR1634421},
\begin{align}
\label{ModelMap0}
h \rightarrow  &
\frac{\xi_1 dy_1(z)+\xi_2 dy_2(z)}{\omega(z)}
P(\xi_1,\xi_2,z)
\Phi
\\ &
(\xi_1A_1+\xi_2A_2 - \xi_1y_1(z) -\xi_2y_2(z))^{-1}h
\nonumber
\end{align}
and defines a section of $\widetilde{\mathcal {E}} (p)$.
In \eqref{ModelMap0}, $h \in H$, $z \in X$, $\xi_1$ and $\xi_2$ are free parameters
and $\omega(z)$ is a meromorphic differential of order $(m-3)$ at the points of $C$ at infinity
and poles on the divisor of singularities.
The projection $\widetilde{P}$ of $E$ into the fiber $\widetilde{\mathcal {E}} (p)$,
is given by (see \cite[Equations 2.29]{MR97m:30051}):
\begin{equation}
\label{proj1}
\widetilde{P}(\xi_1,\xi_2,z)
=
\widetilde{\bf u}^{\times}(z) \widetilde{\bf u}^{\times}_l(z)
\frac{\xi_1 \sigma_1 + \xi_2 \sigma_2}{\xi_1 d y_1(z) + \xi_2 d y_2(z)}
,
\end{equation}
where $\widetilde{\bf u}^{\times} (z)$ and $ \widetilde{\bf u}^{\times} _l (z)$
are the normalized sections of $\widetilde{\mathcal E}(z)$
in column and row representations, respectively.
The mapping in \eqref{ModelMap0} may define a section of $L_{\widetilde{\zeta}}\otimes \Delta$ by pulling back according to the isomorphism \eqref{eqIsoExplicit}
and hence \eqref{ModelMap0} becomes:
\begin{align}
\label{ModelMap1}
h \rightarrow \varphi_h(z)
\overset{\text{def} } {=} &
{\bf \widetilde{u}}_l^\times(z)
(\xi_1 \sigma_1 + \xi_2 \sigma_2)
\Phi
\\ &
(\xi_1A_1+\xi_2A_2 - \xi_1y_1(z) -\xi_2y_2(z))^{-1}h.
\nonumber
\end{align}
\begin{Tm}
\label{vesselAreEq}
Let $\mathcal V$ be an irreducible commutative two-operator vessel and let
$y_1$ and $y_2$ be a pair of meromorphic functions generating $\mathcal M (X)$.
Then $\mathcal V$ is unitary equivalent to the model vessel $\mathcal V_T$ \eqref{eqModelVess}.
\end{Tm}

Since the proof of the statement in \cite{MR1634421} is rather sketchy,
we present the comprehensive proof of Theorem \ref{vesselAreEq} in Section \ref{secMultOp} below.
\smallskip

We note that two vessels sharing the same "external data" and the same characteristic functions
are unitarily equivalent to the same model vessel $\mathcal V_T$ \eqref{eqModelVess}.
Thus, we may conclude the following result.
\begin{Cy}
\label{coro12345}
Any two irreducible vessels with the same $E, \sigma_1, \sigma_2 , \gamma$ and $\widetilde{\gamma}$
are unitarily equivalent if and only if the pair of the associated normalized joint transfer functions are equal.
\end{Cy}
Corollary \ref{coro12345} has been proved in \cite{KLMV} in a substantially different method using the restoration formula.

\section{Statement of the structure theorem}
\label{secMainThm}

Before stating the main result of this section, the counterpart of Theorem \ref{1_1},
we fix some notations and conventions which are used in the upcoming sections.

\begin{Notation}
Let $y_1$ and $y_2$ be real (i.e. satisfying $y(p) = \overline{y(\tauBa{p})}$)
meromorphic functions of degree $n_1$ and $n_2$, respectively,
generating $\mathcal M(x)$ (the field of meromorphic functions on $X$).
Then
\begin{enumerate}[label=(A\arabic*)]
\item
For $\alpha$ and $\beta$ in $\hat{\mathbb{C}}$ and $k=1,2$,
we denote the $n_1$ and $n_2$ distinct pre-images of $\alpha$ and $\beta$ in $X$ of $y_k(\cdot)$ by
$\bracketsC{ \omega_k^{(j)} }{j= 1}{n_k}$ and $\bracketsC{ \nu_k^{(j)} }{j=1}{n_k} $.
\item
We denote by $\bracketsC{ p_k^{(j)} }{j=1}{n_k}$ the $n_k$ poles of $y_k(\cdot)$ where $k=1,2$.
\item
We denote by $C_0$ the embedding of $X$ by $(y_1(x),y_2(x))$ into a curve in $\mathbb C^2$ and by $C$ the projective closure of $C_0$ in $\mathbb P^2$.
\end{enumerate}
\end{Notation}
An important property of the model operator is given in the following lemma.
\begin{lem}
\label{lemMyBounded}
Let $X$ be a compact real Riemann surface and let $y$ be a meromorphic function on $X$.
Let ${\mathcal X}$ be a reproducing kernel Hilbert space of sections of $L_{\widetilde{\zeta}} \otimes \Delta$ (where $\widetilde{\zeta} \in T_{\nu}$)
analytic in an open and connected set $\Omega$ which invariant under $M^y$.
Then the model operator $M^y $ is bounded.
\end{lem}
\begin{pf}
    Since in a reproducing kernel Hilbert space, strong (or even weak)
    convergence implies pointwise convergence,
    we see that the operator $M^y$ is closed.
 Hence, since $M^y$ is closed and everywhere defined,
then, by the closed graph theorem, $M^y$ is bounded.
\end{pf}

We note that the following lemma, presented in \cite{av3} in the case where $\mathcal X$ is finite dimensional.
This result determines under which conditions $R^y_\alpha$ is the resolvent of the model operator $M^y$.
Since the model operator is bounded (Lemma \ref{lemMyBounded}) it follows that the resolvent exists at a neighborhood of infinity
and hence we may state the following result.

\begin{lem}[{\cite[Section 4]{av3}}]
\label{mainThLem}
Let $X$ be a compact real Riemann surface with
$y_1,y_2 \in \mathcal M(X)$ generating $\mathcal M (X)$.
Let $\Omega \subseteq X $ be an open set containing the poles of $y_1$ and $y_2$.
Let ${\mathcal X}$ be a reproducing kernel Hilbert space of sections of $L_{\widetilde{\zeta}} \otimes \Delta$ analytic in $\Omega$.
Then  $R^y_\alpha f$ and $M^y f$, $f \in {\mathcal X}$, are well-defined analytic sections on $\Omega$.
Furthermore,
if $\mathcal X$ is invariant under $M^y$ then $\mathcal X$
is invariant under $R^y_\alpha$ where $\alpha$ in a neighborhood of infinity, and under this condition
\begin{equation*}
R^y_\alpha = (M^y-\alpha I)^{-1}.
\end{equation*}
In particular, the kernel of the resolvent operator is trivial, i.e. $\ker R_\alpha^{y} = \{0\}$.
\end{lem}

The main result is presented in Theorem \ref{MainTh} below.
To prove Theorem \ref{MainTh}, we first state a simpler result in Theorem \ref{preTh}.
It contains the counterpart of the "if" part of Theorem \ref{1_1} under the assumption that the meromorphic functions
$y_1(\cdot)$ and $y_2(\cdot)$ have only simple poles. This assumption is dropped later in Theorem \ref{MainTh}.
In these statements, we recall that $T_\nu$ is given in \eqref{eqRealTorii}.

\begin{theorem}
\label{preTh}
Let $X$ be a compact real Riemann surface and let ${\mathcal X}$
be a reproducing kernel Hilbert space of sections of
$L_{\widetilde{\zeta}} \otimes \Delta$ ( where $\widetilde{\zeta} \in T_{\nu}$) analytic in an open and connected set $\Omega$.
We pick two meromorphic functions, $y_1$ and $y_2$, with simple poles generating $\mathcal M(X)$,
such that $\Omega$ contains the points above the singular points of $C$ and the poles of $y_1$ and $y_2$.
Furthermore, assume that for some $\alpha,\beta \in \mathbb C$ in the neighborhood of infinity
such that $\alpha \neq \beta$ and such that their $n$ pre-images lie within $\Omega$,
the following conditions hold:
\begin{enumerate}[label=(\roman*)]
\item
\label{assump1}
${\mathcal X}$ is invariant under $M^{y_1}$ and $M^{y_2}$.
\item
\label{assump3}
For every choice of $f,g \in {\mathcal X}$ it holds that
\begin{align}
\label{StructIdent4}
\innerProductReg{R_\alpha^{y_k}f}{g} & -
\innerProductReg{f}{R_\beta^{y_k}g} -
(\alpha - \overline{\beta})\innerProductReg{R_\alpha^{y_k}f}{R_\beta^{y_k}g}
=
\nonumber
\\ &
-i(\alpha - \overline{\beta}) \sum_{l,t=1}^{n}
{
\frac{f(\nu^{(l)}) \overline{g(\omega^{(t)})}}
{dy_k(\nu^{(l)})\overline{dy_k(\omega^{(t)})}}
\frac
{\vartheta[\widetilde{\zeta}](\nu^{(l)}-\tauBa{\omega^{(t)}})}
{\vartheta[\widetilde{\zeta}] (0)E(\nu^{(l)},\tauBa{\omega^{(t)}})}
}.
\end{align}
\end{enumerate}
Then the reproducing kernel of $\mathcal X$ is of the form
\begin{equation}
\label{RKstruct}
K_{\mathcal X}(p,q) =
{K_{\widetilde{\zeta}}(p,q) - T(p) K_{\zeta}(p,q) T(q)^*}
\end{equation}
for some $\zeta \in T_{\nu}$ and where $T(\cdot)$ is a $(\zeta, \widetilde{\zeta})$-contractive line bundles mapping.
\end{theorem}

\begin{Rk}
\label{LocCoor}
The structure identity does not depends on the choice of the local coordinates.
Note that the following holds:
\begin{align*}
f_i(p)
&=
\innerProductReg{f}{K_\zeta^i(\cdot,p)}
 =
\innerProductReg{f}{ g_{ij} (p) K_\zeta^j(\cdot,p)}
\\ & =
\overline{g_{ij}(p)}\innerProductReg{f}{ K_\zeta^j(\cdot,p)}
=
\overline{g_{ij}(p)} f_j(p),
\end{align*}
where $K_\zeta^j(q,p)$ is the Cauchy kernel in local coordinates $t_j$ and $g_{ij}$ is the corresponding transition function.
This is coherent with the fact that the collection $\left(f_i\right)_{i \in \mathcal{I}}$ defines a section of the line bundle $L_{\zeta} \otimes \Delta$.
\end{Rk}

\begin{Rk}
\label{rkBlaBla}
One may further consider the structure identity \eqref{StructIdent4} when $\alpha = \overline{\beta}$.
Then, the right hand side of \eqref{StructIdent4} is interpreted as the limit.
Unless specifically stated otherwise, this case is not considered in this upcoming pages.
\end{Rk}

The comprehensive counterpart (removing the simple poles assumption and adding the converse statement)
of Theorem \ref{1_1}, in the setting of compact real Riemann surfaces, is given below.
\begin{theorem}
\label{MainTh}
Let $X$ be a compact real Riemann surface and
let ${\mathcal X}$ be a reproducing kernel Hilbert space of sections
of $L_{\widetilde{\zeta}} \otimes \Delta$
(where $\widetilde{\zeta} \in T_{\nu}$).
\begin{enumerate}[label=(\alph*)]
\item
\label{stateA}
Assume that the elements of $\mathcal X$ analytic in an open and connected set $\Omega$.
Let $y_1$ and $y_2$ be two meromorphic functions generating $\mathcal M(X)$,
such that $\Omega$ contains the points above the singular points of $C$ and contains the
pre-images of some $(\beta_1,\beta_2) \in (\mathbb R \cup {\infty})^2$.
Furthermore, we assume that Assumptions \ref{assump1} and \ref{assump3} in Theorem \ref{preTh} are in force
for all points in some neighborhood of $(\beta_1,\beta_2)$.
Then the reproducing kernel of $\mathcal X$ is of the form
\begin{equation}
\label{rkhsKernel}
K_{\mathcal X}(p,q) =
{K_{\widetilde{\zeta}}(p,q) - T(p) K_{\zeta}(p,q) T(q)^*}
\end{equation}
for some $\zeta \in T_{\nu}$ and where $T(\cdot)$ is a $(\zeta, \widetilde{\zeta})$-contractive line bundles mapping.
\item
\label{stateB}
Conversely,
let $T$ be a $(\zeta,\widetilde{\zeta})$-contractive mapping and
$\mathcal X$ has a reproducing kernel of the form \eqref{rkhsKernel}.
Then all the elements of $\mathcal X$ are analytic on $\Omega_T$ (the region of regularity of $T$).
Furthermore, for any $y(\cdot)$, a real meromorphic function on $X$ such that all its poles are contained in $\Omega_T$,
$\mathcal X$ is $R_\alpha^{y}$--invariant and the structure identity \eqref{StructIdent4} holds.
\end{enumerate}
\end{theorem}
In the finite dimensional case, it has been proved \cite[Section 3]{av3} that $T(\cdot)$
has the form of a finite Blaschke product on a compact Riemann surface, that is,
a finite product of the Blaschke factors
\[
b_a(u) = \frac{E(u,a)}{E(u,\overline{a})}\exp \left(-2\pi (a -\overline{a})^t Y u\right).
\]
The proofs of Theorems \ref{preTh} and \ref{MainTh} are presented in Section \ref{secProofStruct}.
We now give the outline of the "only if" part of the proof of Theorem \ref{preTh}.
We start with the observation that by Lemma \ref{mainThLem}
the kernels of the operators $R_\alpha^{y_1}$ and $R_\beta^{y_2}$ are trivial and we have
\begin{equation}
\label{StructIdent5}
R_\alpha^{y_1} = (M^{y_1} - \alpha I) ^{-1}, \qquad R_\beta^{y_2} = (M^{y_2} - \beta I) ^{-1}.
\end{equation}
Then one proceeds as follows:
we start by presenting and constructing the two-operator vessel corresponding to the operators $M^{y_1}$ and $M^{y_2}$
and then we follow the next steps:
\begin{enumerate}[label=({\bf  Step {\arabic*}})]
\item
Prove that the colligation conditions of $M^{y_1}$ and $M^{y_2}$ are equivalent to the structure identities of $M^{y_1}$ and $M^{y_2}$, respectively.
\item
Show that the output vessel condition holds.
\item
Construct the matrix $\gamma$ such that the input vessel condition holds.
\item
Prove that, in our setting, the mapping \eqref{ModelMap1}, between the inner space of the vessel to the model space, is the identity mapping.
\item
Present the reproducing kernel in terms of the joint characteristic function.
\item
Show that the input and output determinantal representations are maximal.
\item
Conclude, by the reproducing kernel Hilbert space properties, that the reproducing kernel has the desired structure \eqref{RKstruct}.
\end{enumerate}
\section{Subspaces of \texorpdfstring{$H^2_{\widetilde{\zeta}}$}{ a } and a version of Beurling's Theorem}
\label{secH2}

In this section, we present two versions of Beurling's theorem on finite bordered compact Riemann surfaces.
For a short survey and related previous results, see Section \ref{sucSectionIntoCRS}.
\smallskip

Let $S$ be a open Riemann surface so that $S \cup \partial S$ is a finite bordered Riemann surface of genus $g_S$
whose boundary consists $k$ connected components, denoted by $X_1, ..., X_k$.
The double of $S$ is a compact Riemann surface $X$ with a natural antiholomorphic involution $\tau$, turning $X$ into a real compact  Riemann surface of genus $g = 2g_s	+k-1$.
The boundary $\partial S$ coincides with the set of fixed points of $\tau$ on $X$ (denoted by $X_{\mathbb R}$).
Furthermore, $X$ is a compact Riemann surface of dividing type since $X \setminus X_{\mathbb R}$ contains two connected components $X_+ = S$ and $X_-$.
\smallskip

As discussed in \cite[Section 2]{av2}, in a more general context,
a flat unitary line bundle on $S\cup\partial S$ can be uniquely extended to a
flat unitary line bundle on $X$  such that certain symmetry properties are fulfilled,
More precisely, to a line bundle $L_{\zeta}$ where $\zeta \in T_0$.
\smallskip

Let now $L_{\zeta}$ where $\zeta \in T_0$ be a flat unitary line bundle on $X$
and let $\Delta$ be a square root of the canonical line bundle as in Section \ref{sub21}.
Then, the corresponding Hardy space consists essentially of sections of the line bundle $L_{\zeta} \otimes \Delta$ analytic on $S$ and satisfying
(in the sense of non tangentional boundary values)
\[
\sum_{j=0}^{k-1}{\int_{X_j}{f(p)^*f(p)}} < \infty,
\]
it becomes a Hilbert space equipped with the inner product
\[
\innerProductReg{f}{g} =  \sum_{j=0}^{k-1}{\int_{X_j}{g(p)^*f(p)}}.
\]
See \cite{av2} for the precise definitions.
To simplify notations, we set $H^2_{\widetilde{\zeta}} = H^2 (X_+, L_{\widetilde{\zeta}}\otimes \Delta)$.
\smallskip

Before turning to the main theorem, let us recall the definition of dividing--type functions on a compact real Riemann surfaces.
\begin{definition}
\label{dividFunc}
A real meromorphic function $y$ on a compact real Riemann surface $X$ is of {\it dividing type}
if $u \in X_{\mathbb R}$ if and only if $y(u) \in {\mathbb R}$.
\end{definition}
A known result regarding functions of dividing type is presented below.
However, we note that only the first part is used in the sequel.
\begin{lem}[{\cite{MR0021108,ahlfors}, \cite{MR1072300,MR1070485} and \cite[Proposition 5.2]{shamovich2014livsic}}]
\label{elisLemma}
Let $y$ be a dividing function on $X$.
Then $y$ has only real simple poles and simple
zeros and its residues at the poles, with respect to a real local coordinate with positive orientation, are negative.
Conversely, if $X$ is dividing and $y$ is a real meromorphic function on $X$ with simple real poles and negative residues with respect to positive real
local coordinate, then $y$ is dividing.
\end{lem}
We now turn to state the first version of Beurling's theorem on finite bordered Riemann surfaces.
\begin{theorem}[Finite bordered Riemann surface version to Beurling's Theorem]
\label{BLth}
Let $S$ be a finite bordered Riemann surface and let $X$ be its double.
Let $H^2_{\widetilde{\zeta}}$ be the Hardy space on $X_+$ corresponding to $\widetilde{\zeta} \in T_0$
and let $y_1$ and $y_2$ be real meromorphic functions in $X$ of dividing type generating $\mathcal M (X)$.
Furthermore, assume that for $\mathcal H \subseteq H^2_{\widetilde{\zeta}}$ the following conditions hold:
\begin{enumerate}
\item $\mathcal H$ is a closed subspace of $H^2_{\widetilde{\zeta}}$ and is invariant under the multiplication operators
$\mathcal{M}_{\frac{1}{y_1(u)-\overline{\alpha}}}$ and $\mathcal{M}_{\frac{1}{y_2(u)-\overline{\alpha}}}$ for every $\alpha \in \mathbb C_+$.
\item The elements of $\mathcal H^{\perp}$ (the orthogonal complement of $\mathcal H$)
have analytic extensions with bounded point evaluations in a connected neighborhood of the poles of $y_1$ and $y_2$
and of the pre-images of the singular points of $C$ (the projective curve with coordinates $y_1$ and $y_2$).
\end{enumerate}
Then $\mathcal H$ is of the form
\[
\mathcal H = T H^2_{{\zeta}},
\]
where $T$ is a $(\zeta, \widetilde{\zeta})$-inner function for some $\zeta \in T_0$.
\end{theorem}

Notice that by Proposition \ref{propCont}, $\zeta, \widetilde{\zeta} \in T_0$ implies
$T$ is contractive and thus admits boundary values.
We say the $T$ is {\it inner} if $T$ equal to one almost everywhere on $\partial S$.
\smallskip

We note that the operator $\mathcal M_f$ denotes the conventional multiplication operator by a function $f$
while $M^y$, as before, denotes the model operator.
Furthermore, for $\alpha \in C_+$ and $y$ of dividing type, the operator $\mathcal{M}_{\frac{1}{y(u)-\overline{\alpha}}}$ is well-defined on $H^2_{\widetilde{\zeta}}$.
This follows since $\overline{\alpha} \in X_-$, the function $\frac{1}{y(u)-\overline{\alpha}}$ is meromorphic in $X$ and analytic in $S \cup \partial S$.
Hence, $\mathcal{M}_{\frac{1}{y(u)-\overline{\alpha}}}$ is bounded and sends $H^2_{\widetilde{\zeta}}$ to $H^2_{\widetilde{\zeta}}$.
\smallskip

Before heading to the proof, we present several preliminary results.
Our first goal is to show that the structure identity \eqref{StructIdent4} holds for all elements in $H^2_{\widetilde{\zeta}}$.
We start with the following result and we recall the proof for the sake of completeness.
\begin{lem}[{\cite[Theorem 4.3]{av3}}]
\label{R_alpha_evs}
Let $\alpha \in \mathbb C$ have $n$ distinct pre-images with respect to a real meromorphic function $y$
and let $\zeta \in J(X)$ such that $ \theta(\zeta) \neq 0$.
Then the Cauchy kernels are the eigenvectors of the resolvent operator
$R^y_\alpha$ with eigenvalue $\frac{1}{\overline{y(w)}-\alpha}$.
\end{lem}
\begin{pf}
We apply the resolvent operator $R^y_\alpha$ \eqref{Ralpha} on $K_{\widetilde{\zeta}}(z,w)$.
Then, a direct computation, using the collection formula \cite[Lemma 4.1]{av3}, yields the following
\begin{align*}
\big( R_\alpha^y \, & K_{\widetilde{\zeta}}(\cdot,w) \big) (z)
=
\frac{K_{\widetilde{\zeta}}(z,w)}{y(z)-\alpha}
-
\sum_{j=1}^{n}
{
\frac{K_{\widetilde{\zeta}}(u^{(j)},w)}{dy(u^{(j)})}
\frac
{\vartheta[\widetilde{\zeta}](u^{(j)}-z)}
{ \vartheta[\widetilde{\zeta}] (0)E(u^{(j)},z)}
}
\\
=
&
\frac{K_{\widetilde{\zeta}}(z,w)}{y(z)-\alpha}
-
\frac{-1}{i}
\sum_{j=1}^{n}
{
\frac{1}{dy(u^{(j)})}
\frac
{\vartheta[\widetilde{\zeta}](\overline{w} - u^{(j)})}
{ \vartheta[\widetilde{\zeta}] (0)E(\overline{w},u^{(j)})}
\frac
{\vartheta[\widetilde{\zeta}](u^{(j)}-z)}
{ \vartheta[\widetilde{\zeta}] (0)E(u^{(j)},z)}
}
\\
=
&
\frac{K_{\widetilde{\zeta}}(z,w)}{y(z)-\alpha}
-
\frac
{- \vartheta[\widetilde{\zeta}](\overline{w} - z)}
{ i \vartheta[\widetilde{\zeta}] (0)E(\overline{w},z)}
\left(
\frac{1}{y(z)-\alpha}
-
\frac{1}{\overline{y(w)}-\alpha}
\right)
\\
=
&
\frac{K_{\widetilde{\zeta}}(z,w)}{y(z)-\alpha}
-
K_{\widetilde{\zeta}}(z,w)
\left(
\frac{1}{y(z)-\alpha}
-
\frac{1}{\overline{y(w)}-\alpha}
\right)
\\
=
&
\frac{1}{\overline{y(w)}-\alpha}K_{\widetilde{\zeta}}(z,w)
.
\end{align*}
\end{pf}
Using the preceding lemma, we may conclude and prove that the structure identity holds on a dense subset of $H^2_{\widetilde{\zeta}}$.
\begin{lem}
\label{stIdH2}
Let $\alpha,\beta \in \mathbb C _+$. Then the structure identity \eqref{StructIdent4} holds on the linear span of Cauchy kernels inside $H^2_{\widetilde{\zeta}}$.
\end{lem}
\begin{pf}
Since $y$ is of dividing type, it maps $X_+$ to $\mathbb C_+$ and $X_-$ to $\mathbb C_-$. Thus, for $\alpha,\beta \in \mathbb C _+$, their pre-images are in $X_+$
and $R_\alpha^y$ and $R_\beta^y$ are well-defined operators in $H^2_{\widetilde{\zeta}}$.
\smallskip

It is enough to verify that \eqref{StructIdent4} holds on a pair of kernels
$f(z) = K_{\widetilde{\zeta}}(z,\omega)$ and $g(z) = K_{\widetilde{\zeta}}(z,\nu)$.
Then, starting with the left hand side of \eqref{StructIdent4} and using
Lemma \ref{R_alpha_evs}, we compute separately the three components.
    The first two components, using Lemma \ref{R_alpha_evs}, are given by:
\begin{align}
\label{firstcomp}
\innerProductReg{R_\alpha^y K_{\widetilde{\zeta}}(z,\omega)}{K_{\widetilde{\zeta}}(z,\nu)}
& =
\frac{1}{\overline{y(\omega)}-\alpha}
\innerProductReg{K_{\widetilde{\zeta}}(z,\omega)}{K_{\widetilde{\zeta}}(z,\nu)}
\\ & =
\nonumber
\frac{1}{\overline{y(\omega)}-\alpha}K_{\widetilde{\zeta}}(\nu,\omega),
\end{align}
and similarly,
\begin{equation}
\label{secondcomp}
\innerProductReg{K_{\widetilde{\zeta}}(z,\omega)}{R_\beta^y K_{\widetilde{\zeta}}(z,\nu)}
 =
\frac{1}{y(\nu)-\overline{\beta}}K_{\widetilde{\zeta}}(\nu,\omega).
\end{equation}
    The third element, applying Lemma \ref{R_alpha_evs} once again, is
\begin{align}
\label{thirdcomp}
\left( \alpha - \overline{\beta} \right)
\innerProductReg{R_\alpha^y K_{\widetilde{\zeta}}(z,\omega) }{R_\beta^y K_{\widetilde{\zeta}}(z,\nu)}
& =
\frac{\alpha - \overline{\beta} }{\overline{y(\omega)}-\alpha}
\innerProductReg{K_{\widetilde{\zeta}}(z,\omega)}{R_\beta^y K_{\widetilde{\zeta}}(z,\nu)}
\\
\nonumber
& =
\frac{\alpha - \overline{\beta}}
{(\overline{y(\omega)}-\alpha)(y(\nu)-\overline{\beta})}
K_{\widetilde{\zeta}}(\nu,\omega).
\end{align}
Summing all three components together, \eqref{firstcomp}, \eqref{secondcomp} and \eqref{thirdcomp}, one may conclude:
\begin{align}
\label{StructIdent123}
\innerProductReg{R_\alpha^y f}{g} - &
\innerProductReg{f}{R_\beta^y g} -
(\alpha - \overline{\beta})
\innerProductReg{R_\alpha^y f}{R_\beta^y g}
=
\nonumber
\\ = &
K_{\widetilde{\zeta}}(\nu,\omega)
\left(
\frac{1}{\overline{y(\omega)}-\alpha} -
\frac{1}{y(\nu)-\overline{\beta}} -
\frac
{\alpha - \overline{\beta}}
{(\overline{y(\omega)}-\alpha)(y(\nu)-\overline{\beta})}
\right)
\\
\nonumber
= &
K_{\widetilde{\zeta}}(\nu,\omega)
\frac{y(\nu) - \overline{y(\omega)}}
{(\overline{y(\omega)}-\alpha)(y(\nu)-\overline{\beta})}.
\end{align}
On the other hand, the right hand side (we use the notation $RHS$) of \eqref{StructIdent4}
can be simplified by using \cite[Lemma 4.1]{av3} twice.
The first part of the calculation is:
\begin{align}
\nonumber
RHS
=
&
-i(\alpha - \overline{\beta})
\sum_{l,t=1}^{n}
{
\frac
{K_{\widetilde{\zeta}}(u^{(l)},\omega) \overline{K_{\widetilde{\zeta}}(v^{(t)},\nu)}}
{dy(u^{(l)})\overline{dy(v^{(t)})}}
\frac
{\vartheta[\widetilde{\zeta}](u^{(l)}-\overline{v^{(t)}})}
{\vartheta[\widetilde{\zeta}] (0)E(u^{(l)},v^{(t)})}
}
\\
\nonumber
=
&
\sum_{t=1}^{n}
\frac{\overline{K_{\widetilde{\zeta}}(v^{(t)},\nu)}}{\overline{dy(v^{(t)})}}
\sum_{l=1}^{n}
\frac{\alpha - \overline{\beta}}{dy(u^{(l)})}
\frac
{\vartheta[\widetilde{\zeta}](\tauBa{\omega} - u^{(l)})}
{ \vartheta[\widetilde{\zeta}] (0)E(\tauBa{\omega},u^{(l)})}
\frac
{\vartheta[\widetilde{\zeta}](u^{(l)}-\overline{v^{(t)}})}
{\vartheta[\widetilde{\zeta}] (0)E(u^{(l)},\overline{v^{(t)}})}
\\
\nonumber
=
&
-i(\alpha - \overline{\beta})
\left(
\frac{1}{\overline{\beta}-\alpha}
-
\frac{1}{\overline{y(\omega)}-\alpha}
\right)
\sum_{t=1}^{n}
{
\overline{ \left( \frac{K_{\widetilde{\zeta}}(v^{(t)},\nu)}{dy(v^{(t)})} \right) }
\overline{K_{\widetilde{\zeta}}(\omega,\overline{v^{(t)}})}
}
\\
\label{Req2}
=
&
-i
\left(
1
+
\frac{\alpha - \overline{\beta}}{\overline{y(\omega)}-\alpha}
\right)
\overline{
\sum_{t=1}^{n}
{
\frac{1}{dy(v^{(t)})}
\frac
{\vartheta[\widetilde{\zeta}](\overline{\nu}-v^{(t)})}
{ \vartheta[\widetilde{\zeta}] (0)E(\overline{\nu},v^{(t)})}
\frac
{\vartheta[\widetilde{\zeta}]( v^{(t)} - \omega )}
{ \vartheta[\widetilde{\zeta}] (0)E(v^{(t)},\omega)}
}
}
.
\end{align}
Using the collection formula once again, we have the following equality
\begin{align}
\nonumber
&
\overline{
\sum_{t=1}^{n}
\frac{1}{dy(v^{(t)})}
\frac
{\vartheta[\widetilde{\zeta}](\overline{\nu}-v^{(t)})}
{ \vartheta[\widetilde{\zeta}] (0)E(\overline{\nu},v^{(t)})}
\frac
{\vartheta[\widetilde{\zeta}]( v^{(t)} - \omega )}
{ \vartheta[\widetilde{\zeta}] (0)E(v^{(t)},\omega)}
}
\\
\nonumber
=
&
\overline{
(-i)
\frac
{\vartheta[\widetilde{\zeta}](\overline{\nu}-\omega)}
{i\vartheta[\widetilde{\zeta}] (0)E(\omega,\overline{\nu})}
}
\left(
\frac{1}{\overline{y(\omega)}-\overline{\beta}}
-
\frac{1}{y(\nu)-\overline{\beta}}
\right)
\\
\label{Req1}
=
&
i
\left(
\frac{1}{\overline{y(\omega)}-\overline{\beta}}
-
\frac{1}{y(\nu)-\overline{\beta}}
\right)
K_{\widetilde{\zeta}}(\nu,\omega)
.
\end{align}
Substituting \eqref{Req1} in \eqref{Req2} leads to
\begin{align*}
RHS = &
-i
\left(
1
+
\frac{\alpha - \overline{\beta}}{\overline{y(\omega)}-\alpha}
\right)
i
\left(
\frac{1}{\overline{y(\omega)}-\overline{\beta}}
-
\frac{1}{y(\nu)-\overline{\beta}}
\right)
K_{\widetilde{\zeta}}(\nu,\omega)
\\
=
&
\frac{y(\nu)-\overline{y(\omega)}}
{(y(\nu)-\overline{\beta})(\overline{y(\omega)}-\alpha)}
K_{\widetilde{\zeta}}(\nu,\omega),
\end{align*}
as in \eqref{StructIdent123}.
\end{pf}

Moreover, using Lemma \ref{stIdH2}, we show below that the operator $R_{\alpha}^{y}$ is a bounded operator in $H^2_{\widetilde{\zeta}}$.

\begin{lem}
\label{R_alpha_evs1}
Let $\alpha \in \mathbb C _ +$ and let $y$ be a real meromorphic function of dividing type.
Then the resolvent operator $R_{\alpha} ^{y}$ is bounded in $H^2_{\widetilde{\zeta}}$.
\end{lem}
\begin{pf}
Using Lemma \ref{stIdH2}, the structure identity \eqref{StructIdent4}
holds for any linear combination of Cauchy kernels.
We use the structure identity in order to prove the boundedness of $R_{\alpha} ^{y}$.
In \eqref{StructIdent4}, we choose $f=g$ to be a linear combination of Cauchy kernels and we set $\beta = \alpha \in \mathbb{C}_{+}$.
Then we have:
\begin{align}
\label{StructIdent16}
2 & \left(
\mathfrak{Im} ~ \innerProductReg{R_\alpha^{y}f}{f} -
\imagg{\alpha} \norm{ R_\alpha^{y}f } ^2
\right)
=
\\ &
\nonumber
- 2 i \imagg{\alpha} \sum_{l,t=1}^{n}
{
\frac{f(\nu^{(l)})
\overline{f(\nu^{(t)})}}
{dy(\nu^{(l)})\overline{dy(\nu^{(t)})}}
\frac
{\vartheta[\widetilde{\zeta}](\nu^{(l)}-\overline{\nu^{(t)}})}
{\vartheta[\widetilde{\zeta}] (0)E(\nu^{(l)},\overline{\nu^{(t)}})}
}
.
\end{align}
The right hand side of Equation \eqref{StructIdent16}, in view of Remark \ref{LocCoor},
does not depend on the local coordinates choice.
As a result, we deduce from \eqref{StructIdent16} the following inequality:
\begin{equation}
\label{StructIdent977979}
\norm{ R_\alpha^{y}f } ^2
\leq
C_{\alpha} \ \norm{ R_\alpha^{y}f } \
\norm{f}
+
D_{\alpha}\norm{f}^2
,
\end{equation}
for some constants $C_{\alpha}$ and $D_{\alpha}$ (depend only on $\alpha$).
The inequality \eqref{StructIdent977979} is true for every $f$ in a densely defined subset of $H^2_{\widetilde{\zeta}}$.
Hence, dividing by $\norm{f}^2$ and taking the supremum over $f$ implies that the operator $R^{y}_{\alpha}$ is bounded.
\end{pf}

Combining the last two results, Lemma \ref{stIdH2} and Lemma \ref{R_alpha_evs1},
we conclude that since $R_\alpha^y$ is bounded and since the structure identity holds
on a dense subset, the structure identity holds in $H^2_{\widetilde{\zeta}}$.

\begin{corollary}
\label{stIdH3}
Let $\alpha,\beta \in \mathbb C _+$ and let $y(z)$ be a real meromorphic function of dividing type on a compact real Riemann surface $X$.
Then the structure identity \eqref{StructIdent4} holds in $H^2_{\widetilde{\zeta}}$.
\end{corollary}

The link between $R^y_\alpha$-invariant subspaces and subspaces which are
invariant under multiplication operators is illustrated in the following lemma.

\begin{lem}
\label{MM}
    Let $\mathcal H$ be a subspace of $H^2_{\widetilde{\zeta}}$ and $\alpha \in \mathbb C _+$.
    Then $\mathcal H$ is $R^y_\alpha$-invariant
    if and only if the orthogonal complement $\mathcal H ^{\perp}$
    is invariant under the operator of multiplication by $\frac{1}{y(\omega)-\overline{\alpha}}$.
\end{lem}
\begin{pf}
Since $R^y_\alpha$ is bounded, it is sufficient to prove that
\[
\innerProductReg{R^y_\alpha f }{ g } =
\innerProductReg{ f}{ \mathcal{M}_{\frac{1}{y(\cdot)-\overline{\alpha}}}g},
\]
where $\mathcal{M}$ is the multiplication operator and $f$ and $g$ are in $H^2_{\widetilde{\zeta}}$.
Let $f$ and $g$ be kernel elements,  $f=K_{\widetilde{\zeta}}(\cdot,\nu)$ and $g=K_{\widetilde{\zeta}}(\cdot,\omega)$.
Using the boundedness property of $R^y_\alpha$ in Lemma \ref{R_alpha_evs1}, one may obtain the following
\begin{align}
\label{eval1}
\innerProductReg{ K_{\widetilde{\zeta}}(\cdot,\nu)  }{ R^y_\alpha K_{\widetilde{\zeta}}(\cdot,\omega)}
& =
\innerProductReg{ \left(R^y_\alpha \right)^* K_{\widetilde{\zeta}}(\cdot,\nu) }{ K_{\widetilde{\zeta}}(\cdot,\omega) }
\\
& =
\nonumber
\left(\left( R^y_\alpha \right)^* K_{\widetilde{\zeta}}(\cdot,\nu)\right) (\omega).
\end{align}
On the other hand, using Lemma \ref{R_alpha_evs}, we have
\begin{align}
\label{eval2}
\nonumber
\innerProductReg{ K_{\widetilde{\zeta}}(\cdot,\nu) }{ R^y_\alpha K_{\widetilde{\zeta}}(\cdot,\omega)  }
& =
\innerProductReg{ K_{\widetilde{\zeta}}(\cdot,\nu) }{ \frac{1}{\overline{y(\omega)}-\alpha}  K_{\widetilde{\zeta}}(\cdot,\omega)  } \\
& =
\frac{1}{y(\omega)-\overline{\alpha}}  K_{\widetilde{\zeta}}(\omega,\nu).
\end{align}
Hence, combining \eqref{eval1} and \eqref{eval2} the desired result yields.
\end{pf}

To apply Theorem \ref{preTh},
one needs to show that the structure identity holds also on $\partial S$ in a neighborhood of the poles of $y_1$ and $y_2$.
This point is proved in Lemma \ref{MM1} below, but first, we present another preliminary result.
\begin{lem}
\label{MM13}
Let $X$ be a compact real Riemann surface and let $y$ be a dividing-type meromorphic function on $X$.
Furthermore, let $\mathcal H^\perp$ be a closed subspace of $H^2_{\widetilde{\zeta}}$
which is invariant under $R^y_\alpha$ where $\alpha \in \mathbb C_+$ and such that all elements of $H^2_{\widetilde{\zeta}}$
have analytic extensions
in a neighborhoods of the poles of $y$.
Then,
for any sequence $(\alpha_j)_{j \in \mathbb N} \subseteq \mathbb C_+$ converging to $\alpha_0 \in \mathbb R$,
for any $f \in \mathcal H^\perp$ and for all $p \in X_+$
the limit $\lim_{j \rightarrow \infty}  R_{\alpha_j}^y \, f(p)$
is well-defined and does not depend on the sequence $(\alpha_j)_{j \in \mathbb N}$.
\end{lem}
\begin{pf}
Let $f$ be an arbitrary element in $\mathcal H^\perp$
and let $(\alpha _j)_{j \in \mathbb N}$ be a sequence of elements in $\mathbb{C}_+$ converging to $\alpha_0 \in \mathbb R$.
We examine the limit $\lim_{j \rightarrow \infty}  R_{\alpha_j}^y f$ at a point $p \in X_{+}$
(a priori, the limit may depend on the sequence $(\alpha _j)_{j \in \mathbb N}$).
By definition, we have
\begin{equation}
\label{Ryalphafp}
(R_{\alpha_j}^y f)(p) =
\frac{f(p)}{y(p) - \alpha_j} -
\sum_{k=1}^{n}{\frac{f(u_j^{(k)})}{dy(u_j^{(k)})}\frac{\vartheta[\widetilde{\zeta}](u_j^{(k)}-p)}{\vartheta[\widetilde{\zeta}](0)E(u_j^{(k)},p)}},
\end{equation}
where $( u_j^{(k)} )_{k=1}^n$ are the points in $X_+$ such that $y(u_j^{(k)}) = \alpha_j$.
The first element on the right hand side of \eqref{Ryalphafp}, for a given $p$, converges to $\frac{f(p)}{y(p)-\alpha_0}$.
Moreover, by the continuity of $y$ and the analytic continuation property of $f$, it follows that
$f(u_j^{(k)}) \xrightarrow[ ]{j \rightarrow \infty}f(u_0^{(k)})$ where $u_0^{(k)}$ for $k=1,...,n$ are the points above $\alpha_0$.
Thus, for any $p \in X_+$, the point-wise limit exists in \eqref{Ryalphafp}
and does not depend on the choice of the sequence $(\alpha _j)_{j \in \mathbb N}$.
\end{pf}
\begin{lem}
\label{MM1}
Let $X$ be a compact real Riemann surface of dividing type,
let $y(\cdot)$ be a real meromorphic function of dividing type on $X$
and let $H^2_{\widetilde{\zeta}}$ be the Hardy space corresponding to $\widetilde{\zeta} \in T_0$.
Furthermore,
we assume that $\mathcal H$ is a closed subspace which is invariant under $\mathcal{M}_{\frac{1}{y(u)-\overline{\alpha}}}$
for $\alpha \in \mathbb C_+$
and the elements of $\mathcal H^\perp$ have analytic extensions with bounded point evaluations
in a neighborhood of the poles of $y(\cdot)$.
Then, for any $\alpha_0 \in \mathbb R$ in a neighborhood of infinity, the following hold:
\begin{enumerate}
\item
The operator $R_{\alpha_0}^y$ is $\mathcal H^\perp$-invariant.
\item
\label{lemAlpha0SI}
The structure identity may be extended to $\alpha_0 \neq \beta_0$
in $\mathbb R$ in a neighborhood of infinity
(see Remark \ref{rkBlaBla}, for the case where $\alpha_0=\beta_0\in \mathbb R$).
\end{enumerate}
\end{lem}
\begin{pf}
We fix $\alpha_0 \in \mathbb R$ and let $f$ be an element of $H^2_{\widetilde{\zeta}}$
with bounded evaluations in the neighborhood of the poles of $y$ and the fibre of $y$ above $\alpha_0$.
First, we show that $R_{\alpha_0} ^y \, f$ belongs $H^2_{\widetilde{\zeta}}$.
To see that this statement is indeed valid, we explicitly calculate the $H^2_{\widetilde{\zeta}}$ norm of $R_{\alpha_0} ^y \, f$.
The strategy is to divide the integration path $X_{\mathbb R}(\varepsilon)$ (a contour located at a distance $\varepsilon$ from $X_\mathbb R$) into two segments.
The first, is over the set of arcs around the pre-images of $\alpha_0$ (where $f$ has bounded evaluations).
We denote the first integration path by $X^{\alpha_0}_{\mathbb R}(\varepsilon)$.
The second, is $X^{\alpha_0}_{\mathbb R}(\varepsilon)^c \overset{\text{def} } {=} X_{\mathbb R}(\varepsilon) \backslash X^{\alpha_0}_{\mathbb R}(\varepsilon)$, where $\frac{1}{y(u)-\alpha}$ can be bounded.
Namely, we examine
\begin{equation}
\label{eqNormFdivide}
\normTreA{R^y_{\alpha_0} f}{H^2_{\widetilde{\zeta}}}{2}
=
\lim_{\varepsilon \rightarrow 0} \int_{X^{\alpha_0}_{\mathbb R}(\varepsilon)} \abs{R^y_{\alpha_0} f} ^2
+
\lim_{\varepsilon \rightarrow 0} \int_{X^{\alpha_0}_{\mathbb R}(\varepsilon)^c}  \abs{R^y_{\alpha_0} f} ^2
.
\end{equation}
Starting with the left summand, we move to zero-centered local coordinates $(t_j)_{j=1}^n$ and express
the resolvent operator at the neighborhood of $u^{j}_0$, a pre-image of $\alpha_0$ under $y$.
Using the analyticity of $f$ and $y$ at $u^{j}_0$, we consider
$f(t_j) = f(u^{j}_0) + f^{u^{j}_0}_1 t_j + f^{u^{j}_0}_2 t_j^2+O(\abs{t_j}^3)$
and
$y(t_j) =\alpha_0 + y^{u^{j}_0}_1 t_j +y^{u^{j}_0}_2 t_j^2 +O(\abs{t_j}^3)$
and using the prime form properties, we proceed as follows:
\begin{align}
\nonumber
\abs{R^y_{\alpha_0} f(t_j)} ^2
= &
\bigg|
\frac{f(u^{j}_0) + f^{u^{j}_0}_1 t_j + f^{u^{j}_0}_2 t_j^2+...}{(\alpha_0 + y^{u^{j}_0}_1 t_j +y^{u^{j}_0}_2 t_j^2 +...) - \alpha_0}
-
\frac{ f(u^{j}_0) }{y^{u^{j}_0}_1}\frac{1}{ t_j}
- \\ & \nonumber
\sum_{\substack{l=0 \\ l \neq j}}^n  \frac{f(u^{l}_0)}{dy(u^{l}_0)}
\frac{\vartheta[\zeta](u^{l}_0-u^{j}_0)}{\vartheta[\zeta](0)E(u^{l}_0,u^{j}_0)}
\bigg|^2
\\ \leq &
\abs{
\frac
{ y_1^{u^j} f^{u^{j}_0}_1 - y_2^{u^j} f(u^{j}_0) +O(\abs{t_j})  }
{ (y_1^{u^j})^2 + O(\abs{t_j})}
}^2
+
\lim_{u\rightarrow u^{j}_0}
\sum_{\substack{l=0 \\ l \neq j}}^n
\abs{
\frac{f(u^{l}_0)}{dy(u^{l}_0)}
\frac{\vartheta[\zeta](u^{l}_0-u)}{\vartheta[\zeta](0)E(u^{l}_0,u)}
}^2
.
\label{eqNormRf1}
\end{align}
Recall that $f$ has bounded evaluations at the pre-images of $y$ at $\alpha$
and hence the integration of the first term in \eqref{eqNormRf1} along
$X^{\alpha_0}_{\mathbb R}(\varepsilon)$, as $\varepsilon$ tends to zero, is finite.
Furthermore, since the Cauchy kernels in second term of \eqref{eqNormRf1} convergence to
$\frac{\vartheta[\zeta](u^{l}_0-u^{j}_0)}{\vartheta[\zeta](0)E(u^{l}_0,u^{j}_0)}$ can be (upper) bounded as long as $l \ne j$,
the integration of the second term in \eqref{eqNormRf1} is also finite.
\smallskip

$\frac{1}{y(u)-\alpha}$ is bounded from above by some constant $M$ on $X^{\alpha_0}_{\mathbb R}(\varepsilon)^c$
and the second summand of \eqref{eqNormFdivide} becomes
\begin{align}
\nonumber
\lim_{\varepsilon \rightarrow 0} \int_{X^{\alpha_0}_{\mathbb R}(\varepsilon)^c} \abs{R^y_{\alpha_0} h} ^2
= &
\lim_{\varepsilon \rightarrow 0}
\int_{X^{\alpha_0}_{\mathbb R}(\varepsilon)^c}
\abs{
\frac{f(u)}{y(u) -\alpha} -
\sum_{j=1}^n
h_{j}
\frac{\vartheta[\zeta](u^{j}_0-u)}{\vartheta[\zeta](0)E(u^{j}_0,u)}
}^2
\\ \leq &
M^2
\lim_{\varepsilon \rightarrow 0}
\int_{X^{\alpha_0}_{\mathbb R}(\varepsilon)^c}
\abs{
f(u)-
\sum_{j=1}^n
\frac{h_{j}}{M}
\frac{\vartheta[\zeta](u^{j}_0-u)}{\vartheta[\zeta](0)E(u^{j}_0,u)}
}^2,
\label{eqNormRf2}
\end{align}
where the coefficients $\left( h_{j} \right)_{j=1}^n$ depend on the pre-images of $\alpha_0$.
The statement follows since $f$ belongs to $H^2_{\widetilde{\zeta}}$ and furthermore
$\frac{\vartheta[\zeta](u^{j}_0-u)}{\vartheta[\zeta](0)E(u^{j}_0,u)}$ is bounded in $X^{\alpha_0}_{\mathbb R}(\varepsilon)^c$.
Thus, $f \in \mathcal H ^\perp$ implies that $R^y_{\alpha_0} f$ belongs to $H^2_{\widetilde{\zeta}}$.
\smallskip

In order to continue and conclude that $\mathcal H ^\perp$ is invariant under $R^y_{\alpha_0}$,
we fix a sequence $\left( \alpha_j \right)_{j=1}^{\infty} \subseteq \mathbb C_+$ converging to $\alpha_0 \in \mathbb R$ in a neighborhood of infinity
and show that $R^y_{\alpha_j} f$ converges to $R^y_{\alpha_0} f$ in norm.
Also here, we divide the integration path, now on $X_\mathbb R$, into two segments:
\begin{align}
\lim_{i \rightarrow \infty}
\normTreA{R^y_{\alpha_0} f  - R^y_{\alpha_i} f }{H^2_{\widetilde{\zeta}}}{2}
= &
\lim_{i \rightarrow \infty}
\int_{X^{\alpha_0}_{\mathbb R}(0)}
\abs{ R^y_{\alpha_0} f(u) - R^y_{\alpha_i} f(u)}^2
+ \nonumber \\ &
\lim_{i \rightarrow \infty}
\int_{X^{\alpha_0}_{\mathbb R}(0)^c}
\abs{R^y_{\alpha_0} f(u) - R^y_{\alpha_i} f(u) }^2.
\label{eqNormRf3}
\end{align}
The second summand tends to zero due to a similar analysis as in \eqref{eqNormRf2}.
To show that the first element also vanishes, we note that in general,
for any section $f$ analytic on a neighborhood $U$ of the fibre of $y$ over $\alpha_0$ (assuming the fibre is unramified),
$(R^y_{\alpha_0} f)(p)$ is jointly analytic in
$(\alpha,p) \in \{|\alpha-\alpha_0|<\varepsilon\} \otimes U$ (for sufficiently small $\varepsilon$).
The analyticity in $p$ is a consequence of Lemma \ref{MM13}, while the analyticity at $\alpha_0$
follows by (a similar computation appeared in \eqref{eqNormRf1}):
\[
R^y_{\alpha_0} f(t_j)
=
\frac
{ y_1^{u^j} f^{u^{j}}_1 - y_2^{u^j} f(u^{j}) +O(\abs{t_j})  }
{ (y_1^{u^j})^2 + O(\abs{t_j})}
+
\sum_{\substack{l=0 \\ l \neq j}}^n
\frac{f(u^{l})}{dy(u^{l})}
\frac{\vartheta[\zeta](u^{l}-u^{j})}{\vartheta[\zeta](0)E(u^{l},u^{j})}
,
\]
the pole at $u^j$ vanishes and the
fiber above $\alpha_0$ depends analytically on $\alpha_0$.
It then remains to apply the Hartog's Theorem, see for instance \cite[Theorem 2.2.8]{hormander1973introduction},
to conclude the local joint-analyticity property of $(R^y_{\alpha_0} f)(p)$.
\smallskip

As a consequence, $R^y_{\alpha_i} f(p)$ converges to $R^y_{\alpha_0} f(p)$ uniformly in $p$ on a neighborhood
of each $p^{(j)}$ in the fibre over $\alpha$, in particular on a small arc of the boundary near $p^{(j)}$ and so
the second limit in \eqref{eqNormRf3} is also zero.
Hence $\|R^y_{\alpha_i} f - R^y_{\alpha_0} f\|^2_{H^2_{\widetilde{\zeta}}}$ converges to zero,
i.e., $R^y_{\alpha_i} f$ converges to $R^y_\alpha f$ in $H^2_{\widetilde{\zeta}}$.
Thus, for any $f \in H^\perp$, by assumption also $R^y_{\alpha_i} f \in H^\perp$ for all $i \in \mathbb N$
and hence $R^y_{\alpha_0} f \in H^\perp$.
\smallskip

To show that \eqref{lemAlpha0SI} holds,
we choose a sequence $(\beta _j)_{j \in \mathbb N} \subset \mathbb C _ +$
converging to $\beta_0 \in \mathbb C_+$ or belongs to  $\mathbb R$ in a neighborhood of infinity such that $\alpha_0 \neq \beta_0$.
According to Lemma \ref{stIdH3}, the structure identity \eqref{resId} holds for any pair of elements of the sequences
$(\alpha _j)_{j \in \mathbb N}$ and $(\beta _j)_{j \in \mathbb N}$.
Considering $j \rightarrow \infty$, we obtain the following identity:
\begin{align}
\nonumber
\lim_{j\rightarrow \infty}
\big(
&
\innerProductReg{ R^{y}_{\alpha_j} f}{ g} -
\innerProductReg{ f}{  R^{y}_{\beta_j} g} -
(\alpha_j - \overline{\beta_j})\innerProductReg{  R^{y}_{\alpha_j} f }{  R^{y}_{\beta_j} g}
\big)
=
\\
\label{1234}
&
\lim_{j\rightarrow \infty}
\left(
-i(\alpha_j - \overline{\beta_j}) \sum_{l,t=1}^{n}
{
\frac{f(\nu_j^{(l)}) \overline{g(\omega_j^{(t)})}}
{dy_k(\nu_j^{(l)})\overline{dy_k(\omega_j^{(t)})}}
\frac
{\vartheta[\widetilde{\zeta}](\nu_j^{(l)}-\tauBa{\omega_j^{(t)}})}
{\vartheta[\widetilde{\zeta}] (0)E(\nu_j^{(l)},\tauBa{\omega_j^{(t)}})}
}
\right)
.
\end{align}
The limit on the right hand side exists due to the analytic continuation to neighborhoods of the fibers over $\alpha_0$ and $\beta_0$.
Since $R^y_{\alpha_j}$ converges to $R^y_{\alpha_0}$ in the operator norm,
the limit on the left hand side coincides with the operators $R^y_{\alpha_0}$ and $R^y_{\beta_0}$.
\smallskip

To complete the proof, by continuity, the limit on the right hand side of \eqref{1234} exists also for $\alpha_0 = \beta_0 \in \mathbb R$,
and we take the expression on the right hand side as its definition, see Remark \ref{rkBlaBla}.
\end{pf}

As we have gathered all the required preliminary results, we may present the proof of the first version of the Beurling's theorem.
\\

\begin{pf}[of Theorem \ref{BLth}]
Let $y_1(\cdot)$ and $y_2(\cdot)$, as in the statement, be meromorphic functions of dividing type.
Thus, using Lemma \ref{elisLemma}, the poles of $y_1$ and $y_2$ are real and simple. \smallskip

Applying Lemma \ref{MM}, the assumption that $\mathcal H$ is invariant under the multiplication operators
$\mathcal{M}_{\frac{1}{y_1(u)-\overline{\alpha}}}$ and $\mathcal{M}_{\frac{1}{y_2(u)-\overline{\beta}}}$
is translated to: $\mathcal H^{\perp}$ is invariant under the operators $R_{\alpha}^{y_1}$ and $R_{\beta}^{y_2}$ where $\alpha,\beta \in \mathbb C_+$.
By Corollary \ref{stIdH3}, the structure identity automatically holds in $H^2_{\widetilde{\zeta}}$ for all $\alpha,\beta \in \mathbb C _+$.
\smallskip

By assumption, the elements of $\mathcal H^{\perp}$ have analytic extension to neighborhood of the poles of $y_1$ and $y_2$
with bounded point evaluations and the elements of $\mathcal H^{\perp}$ have analytic extensions at the pre-images of the singular points of $C$.
Then, using Lemma \ref{MM1}, $\mathcal H ^\perp$ is invariant under the bounded operator $R_{\alpha}^{y_k}$ where $\alpha \in \mathbb R$ at the
neighborhood of infinity.
Furthermore, the structure identity can be extended to $\alpha \in \mathbb R$ in a neighborhood of infinity.
\smallskip

Combining all the observations above, we can apply Theorem \ref{MainTh} on the orthogonal complement $\mathcal H ^{\perp}$.
Thus, $\mathcal H^{\perp}$ is a reproducing kernel Hilbert space with reproducing kernel of the form
\[
K_{\mathcal H^{\perp}}(p,q) = {K_{\widetilde{\zeta}}(p,q) - T(p) K_{\zeta}(p,q) T(q)^*},
\]
where $T$ is inner from $H^2_\zeta$ to $H^2_{\widetilde{\zeta}}$ for some $\zeta \in T_0$.
Therefore, we may write
\begin{equation*}
\mathcal H ^{\perp} = H^2_{\widetilde{\zeta}} \ominus T H^2_{\zeta},
\end{equation*}
and so
\begin{equation*}
\mathcal H =  T H^2_{\zeta}.
\end{equation*}
By the general theory of reproducing kernel spaces, see \cite{dbr2,sarason94},
it follows that $H^\perp$ is contractively included in $H^2$ and
its generalized orthogonal complement with respect to the range norm is $T H^2$, i.e. $(H^\perp)^{[\perp]} = T H^2$.
On the other hand, we know it is isometrically included in $H^2$ and hence the multiplication by $T$ is isometry and so $T$ is inner.
\end{pf}

The natural condition for a Beurling Lax theorem using de Branges structure theorem, see \cite{sarason94},
is to assume that the invariant subspace is a contractively included subspace of the Hardy space.
However, in our approach, we note that the structure identity does not automatically hold in the Hardy space
and Lemma \ref{MM} when replacing $H^{\perp}$ by the Brangesian generalized orthogonal complement does not hold.
Thus, we add these assumptions in order to state a Beurling's theorem when $T$ is contractive instead of inner.

\begin{theorem}[Finite bordered Riemann surface version to Beurling's Theorem II]
Let $S$ be a finite bordered Riemann surface where $X$ is its double and
let $y_1$ and $y_2$ be real meromorphic functions in $X$ of dividing type generating $\mathcal M (X)$.
Let $H$ be a contractively included subspace subspace of the Hardy space $H^2_{\widetilde{\zeta}}$ on $S=X_+$
such that its generalized orthogonal complement $\mathcal H^{[\perp]}$
is invariant under $R^{y_1}_\alpha$ and $R^{y_2}_\alpha$ for all $\imagg{\alpha}>0$
and the structure identity holds in $\mathcal H^{[\perp]}$.
Then $\mathcal H$ is of the form
\[
\mathcal H = T H^2_{\zeta},
\]
for some $(\zeta, \widetilde{\zeta})$ contractive mapping $T$.
\end{theorem}

The Beurling's theorems in \cite{MR1070485} and \cite{MR0183883}
assume that $\mathcal H$ is invariant under all multiplication operators belong to some algebra of functions
(in \cite{MR0183883}, the collection of functions analytic inside a multiply connected domain $R$ and continuous in $\overline {R}$).
Hence, we are motivated to formulate the following version of the Beurling's theorem.

\begin{Cy}[Finite bordered Riemann surface version to Beurling's Theorem III]
Let $S$ be a finite bordered Riemann surface and let $X$ be its double.
Let $\mathcal H$ be a closed subspace of the Hardy space $H^2_{\widetilde{\zeta}}$ on $S=X_+$
corresponding to $\widetilde{\zeta} \in T_0$ such that the following conditions hold:
\begin{enumerate}
\item
$\mathcal H$ is invariant under the multiplication operators
of all elements in the algebra of functions analytic in $S$ and continuous on $\partial S$.
\item
The elements of $\mathcal H^{\perp}$ have analytic extensions with bounded point evaluations on
an open set $U \subseteq \partial S$.
\item
There exists a pair of meromorphic function $y_1$ and $y_2$ of dividing type such that their poles and the
singular points are in $U$.
\end{enumerate}
Then $\mathcal H$ is of the form
\[
\mathcal H = T H^2_{\zeta},
\]
where $T$ is a $(\zeta, \widetilde{\zeta})$-inner for some $\zeta \in T_0$.
\end{Cy}
\begin{pf}
The functions
$f_1(u) = {\frac{1}{y_1(u)-\overline{\alpha}}}$
and
$f_2(u) = {\frac{1}{y_2(u)-\overline{\beta}}}$
(where $\alpha,\beta \in \mathbb C_+$)
are analytic in $S$ and continuous on $\partial S$.
Therefore, by assumption, $\mathcal H$ is invariant under the multiplications by $f_1$ and $f_2$.
It remains to apply Theorem \ref{BLth}.
\end{pf}

We note that every meromorphic dividing function has poles on every boundary components
and hence a necessary condition is that $U \cap X_j \neq \emptyset$ for $j=0,...k$.
Moreover, it is natural to assume that it is a sufficient condition for
the existence of a pair of meromorphic functions of dividing type,
see \cite{MR1072300} for related results in the case of planar domains.

\section{Compressed multiplication operators on compact real Riemann surfaces}
\label{secMultOp}

We begin by giving an alternative definition of the compressed multiplication operator $M^y$
which is applicable also when $y$ has arbitrary poles.
Let $S$ be a finite set of elements in $X$.
We denote by $\mathcal M_S(L_{\zeta} \otimes \Delta)$ the germs of meromorphic sections of $L_{\zeta} \otimes \Delta$
in a neighborhood of $S$ and with poles only in $S$.
Then we consider the decomposition
\[
\mathcal M_S(L_{\zeta} \otimes \Delta)
=
\mathcal M_{S,+}(L_{\zeta} \otimes \Delta)
\oplus
\mathcal M_{S,-}(L_{\zeta} \otimes \Delta),
\]
where $\mathcal M_{S,+}(L_{\zeta} \otimes \Delta)$ denotes the set of germs of holomorphic sections of $L_{\zeta} \otimes \Delta$ in a neighborhood of $S$
while $\mathcal M_{S,-}(L_{\zeta} \otimes \Delta)$ is the set of global meromorphic sections of $L_{\zeta} \otimes \Delta$ with poles in $S$.
By assumption, $L_{\zeta} \otimes \Delta$ has no global non-zero holomorphic section and therefore
$$\mathcal M_{S,+}(L_{\zeta} \otimes \Delta) \cap \mathcal M_{S,-}(L_{\zeta} \otimes \Delta)=\{0\}.$$
We fix a real meromorphic function $y$ such that all the poles of $y$ belong to $S$ and we define
\begin{align*}
\widetilde{M}^y:  \mathcal{M}_{S,+} & \rightarrow \mathcal{M}_{S,+} \\
f(p)   & \rightarrow  (y(p)f(p))_+.
\end{align*}
The definitions of $\widetilde{M}^y$ and $M^y$ coincide.
This is true since there are no non-zero global analytic sections and
the difference between $\widetilde{M}^y$ and $M^y$ is a global analytic section and hence equal to zero.
\smallskip

In other words, $M^y$ is given by
$$(M^y f )(u)= y(u)f(u)+c_f(u).$$
Here
$f$ is a holomorphic section of $L_{\widetilde{\zeta}} \otimes \Delta $ in a neighborhood of $S$
and $c_f$ is the unique global meromorphic section with
divisor of poles contained in the divisor of poles of $y$,
such that $y(u)f(u)+c_f(u)$ is analytic at the poles of $y$.
\smallskip

Explicitly, we can rewrite the formula of the model operator in \eqref{m_y}
for the case $y$ has also non-simple poles by
\begin{align}
\nonumber
M^{y}f(u)
= &
y(u)f(u) +
\\ &
i
\sum_{m=1}^{n}
\sum_{l=1}^{s^{m}}
\sum_{j=l}^{s^{m}}
a_{-j,m}
\frac{ f^{(j-l)}(p^{(m)})}{ (j-l)!}
\frac{(-1)^{l-1}}{(l-1)!}
K_{\widetilde{\zeta}}^{(l-1)}(p^{(m)},u)
,
\label{MyMatrix}
\end{align}
where the set of the poles of $y$ and their order are given by $\left( p^{(m)} \right)_{m=1}^{n}$ order by $\left( s^{m} \right)_{m=1}^{n}$
and where $a_{-j,m}$ is the $j$-th Laurent coefficient of $y$ at $p^{(m)}$ (see also \cite[Equation 4.21]{av3}).
\smallskip

In order to prove Theorem \ref{lemCollCodI} below, we first present some preliminary results.
The next important result appeared partially in \cite[Theorem 4.10]{av3}.
Here, we present a more general statement with a comprehensive and an alternative proof.

\begin{theorem}
\label{rkModelOpAlg}
Let $S$ be a finite set of points on $X$ and let $\mathcal A$ be the algebra of meromorphic functions on $X$ whose poles are contained in $S$.
The mapping $y \mapsto M^{y}$, where $y$ is a meromorphic function on $X$ whose poles are contained in $S$,
is an algebra homomorphism from $\mathcal {A}$ to $\mathcal M_{S,+}(L_{\zeta} \otimes \Delta)$.
\smallskip

In other words, for a two-variables polynomial $p(x_1,x_2)$ and pair of meromorphic functions
$y_1$ and $y_2$ with poles in $S$, the model operator satisfies
$$M^{p(y_1,y_2)} = p(M^{y_1},M^{y_2}).$$
In particular, $M^{y_1 \, y_2} = M^{y_1}M^{y_2}$ and $M^{y_1 + y_2} = M^{y_1} + M^{y_2}$.
\end{theorem}
\begin{pf}[of Theorem \ref{rkModelOpAlg}]
Let $f$ be an element in $M_{S,+}$, then, by definition, $M^{y_1 + y_2}f = M^{y_1}f+M^{ y_2}f$.
We set  $M^{y_1}M^{y_2} f = h_+ = (h_+ + h_-)_+ $ and we show that $h_+ = M^{y_1 \, y_2} f$.
We also set $y_2 \cdot f_+ = g_+ + g_-$ and then we have
\begin{align*}
h_+(p) = & y_1 (p) g_+(p) - h_-(p)
\\ = &
y_1(p)\big(y_2(p) f(p) - g_-(p)\big)- h_-(p),
\end{align*}
and hence
\[
y_1(p)y_2(p)f(p) = h_+(p) + \big(h_-(p)+y_1(p)g_-(p)\big).
\]
However, $h_-(p)+y_1(p)g_-(p)$ is a global meromorphic section with poles on $S$,
so $\big(y_1(p)y_2(p)f(p)\big)_+  = h_+(p)$ and $M^{y_1 y_2}f = h_+$ follows.
\end{pf}

The main result in this section is given in the next theorem.

\begin{theorem}
\label{lemCollCodI}
Let $y_1$ and $y_2$ be a pair of meromorphic functions on a compact real Riemann surface $X$
Let us assume that the structure identity, given in \eqref{StructIdent4},
holds for a pair of meromorphic functions $y_1$ and $y_2$
and a reproducing kernel Hilbert space $\mathcal X$.
Then, the structure identity holds for all elements in the
algebra of meromorphic functions generated by $y_1$ and $y_2$
\end{theorem}

We start with the following technical result is required later in the proof of Proposition \ref{myGeneral1}.

\begin{lem}
\label{lemLimitZero}
Let $a,b,d \in \mathbb N$ such that $0\leq a,b < d$ and let $c_0,...,c_d$
be a sequence of real numbers such that $c_0 \neq 0$.
Then, the equality
\begin{equation}
\label{eqLimitDelta}
\lim_{x,y \rightarrow 0}
\sum_{j=0}^{d-1-a}
\frac{
c_{d-(j+a+1)} }{ b! j !}
\frac{\partial^{j+b}}{\partial x ^j \partial y ^b}
\frac
{
\sum_{q=0}^{d-1} c_{q}
\sum_{t=1}^{d-q} x^{q+t-1} y^{d-t}
}
{\sum_{p=0}^{d}c_{p}x^{p}    \sum_{p=0}^{d}c_{p}y^{p}}
=
\delta_{b,a}
\end{equation}
holds, where $\delta$ stands for the Kronecker delta.
\end{lem}
\begin{pf}
Along the proof, we use the notations $f(x) \overset{\text{def} } {=} \sum_{p=0}^{d}c_{p}x^{p}$,
\[
g[n] \overset{\text{def} } {=} \frac{1}{n!} \restr{ \frac{d^n}{dx^n} \frac{1}{f(x)} }{x=0}
    \qquad
    {\rm and}
    \qquad
h[m] \overset{\text{def} } {=} \sum_{n=0}^{m} g[n] g[m-n].
\]
According to Leibniz product rule, the following identity holds
\begin{align}
\nonumber
\lim_{x,y \rightarrow 0}
\frac{\partial^{j+b}}{\partial  x ^j\partial y^b}
\frac
{
\sum_{q=0}^{d-1} c_{q}
\sum_{t=1}^{d-q} x^{q+t-1} y^{d-t}
}
{f(x) \, f(y)}
&
=
\sum_{k=0}^{j}  \sum_{l=0}^{b}
\frac{j! g[k]}{(j-k)!}
\frac{b! g[l]}{(b-l)!}
\sum_{q=0}^{d-1} c_{q}
\times \\
&
\lim_{x,y \rightarrow 0}
\frac{\partial ^{(j-k)+(b-l)}}{\partial  x ^ {j-k} \partial  y ^{b-l}}
\sum_{t=1}^{d-q} x^{q+t-1} y^{d-t}
\label{eqLeib}
.
\end{align}
The two-variable polynomial in the numerator on the LHS of \eqref{eqLeib}
has zero coefficients for all monomials of combined degree less than $d-1$.
Hence, all mixed derivatives of orders $j$ and $b$ such that $j+b<d-1$ are zero.
Since the outer sum on the LHS of \eqref{eqLimitDelta} is summed upto $d-1-a$,
\eqref{eqLeib} is identically zero whenever $b<a$.
\smallskip

We move to consider the case where $a \leq b$.
First, we note that
\begin{align}
\nonumber
\sum_{q=0}^{d-1} c_{q}
&
\lim_{x,y \rightarrow 0}
\frac{\partial ^{(j-k)+(b-l)}}{\partial x ^{j-k} \partial  y ^{b-l}}
\sum_{t=1}^{q} x^{q+t-1} y^{d-t}
= \\ &
\begin{cases}
(j-k)! (b-l)! c_{(j-k)+(b-l)-(d-1)}
&
d-1 \leq (j-k)+(b-l),
\\
\quad
0
&
\it{otherwise}.
\end{cases}
\label{eqDeriv}
\end{align}
Here we repeat the same argument (as used in the $b<a$ case),
we note that if $j<d-b-1$, the mixed derivatives of the two-variable polynomial in the numerator of \eqref{eqLimitDelta} are again zero.
Hence, whenever $d-1 \leq b+j$, the RHS of \eqref{eqLeib} becomes
\begin{align*}
RHS
= &
\lim_{x,y \rightarrow 0}
\sum_{k = 0}^{j}
\sum_{l = 0}^{b}
\frac{j! b! g[k] g[l]}{(j-k)!(b-l)!}
\frac{\partial ^{(j-k)+(b-l)}}{\partial  x ^ {j-k} \partial  y ^{b-l}}
\sum_{q=0}^{d-1} c_{q}
\sum_{t=1}^{d-q} x^{q+t-1} y^{d-t}
\\ = &
\sum_{s = 0}^{j+b-d+1}
\sum_{l = 0}^{s}
(s-l)! g[s-l]
l! g[l]
c_{j+b-d+1 - s}
.
\end{align*}
The last equality is the result of setting $s \overset{\text{def} } {=} k+l$.
Combining the last observation together with \eqref{eqLeib} and \eqref{eqDeriv},
the LHS of \eqref{eqLimitDelta} can be rewritten as follows (we just set $j'=d-1-a-j$ and $b^\prime = b-a$):
\begin{align}
LHS =
&
\nonumber
\sum_{j=d-1-b}^{d-1-a}
c_{d-j-a-1}
\sum_{s = 0}^{j+b-d+1}
c_{j+b-d+1 - s}
\sum_{l = 0}^{s}
g[s-l] g[l]
\\ = &
\nonumber
\sum_{j'=0}^{b-a}
\sum_{s = 0}^{j'}
c_{b-a-j'}
c_{j' - s}
h[s]
\\ = &
\sum_{s=0}^{b^\prime}
h[s]
\sum_{m=s}^{b^\prime}
c_{m-s}
c_{b^\prime-m}
\label{eqZero}
.
\end{align}
We note that $h[s]$ is just the $s$-Taylor coefficient at zero of $\frac{1}{f(x)}\frac{1}{f(x)}$ (recall that $c_0 \neq 0$).
On the other hand, $\sum_{m=s}^{b^\prime} c_{m-s} c_{b^\prime-m}$ is the $(b^\prime-s)$-Taylor coefficient at zero of $f(x) f(x)$.
Hence, \eqref{eqZero} is the $b^\prime$-Taylor coefficient at zero of the multiplication of $f(x) f(x)$ and $\frac{1}{f(x)}\frac{1}{f(x)}$.
The $b^\prime$-Taylor coefficient of the function identically equal to one is clearly equal to zero whenever $b^\prime>0$,
or equivalently, when $b>a$ and equal to $1$ if $b^{\prime}=0$ or $b=a$.
Hence the claim follows.
\end{pf}

\begin{Rk}
Let $y(z)$ be a real meromorphic function with $n$ poles of order $s^n$.
Then, the poles of $y$ are either real or appear as conjugate pairs
(recall that $y$ is not necessarily of dividing type).
Hence, we denote by $\left(p^{(j)}\right)_{j=1}^r$ the real poles,
and by $\left(p^{(j)}\right)_{j=r+1}^{r+2m}$ the $m$ pairs of non-real poles ($n = r + 2m$).
\end{Rk}

\begin{prop}
\label{myGeneral1}
Let $y$ be a real meromorphic function on a compact real Riemann surface $X$.
Let $\mathcal X$ be a reproducing kernel Hilbert space of section of
$L_{\widetilde{\zeta}}\otimes\Delta$ which is invariant under the operators
$R^y_\alpha$, $R^y_\beta$ and $M^y$ for fixed $\alpha,\beta \in \mathbb C$.
Then the structure identity \eqref{StructIdent4} is equivalent to the colligation condition \eqref{collCond}.
Here $\sigma$ is defined by
\[
\sigma =
 \begin{psmallmatrix}
  \mathbf{P}^1		& 				&  				&		&			&		$		$\\
  			&	\ddots 		&  				&  		&			&		$		$\\
  			&			 	& 	\mathbf{P}^r		    &  		&			&		$		$\\
  			&			 	& 				& 0		& \mathbf{P}^{r+2}	&		&		&\\
  			&		 		&  				& \mathbf{P}^{r+1}& 0		&		$		$\\
  			&		 		&  				&  		&			&\ddots	$		$\\
  			&		 		&  				&  		&			&		&0	& \mathbf{P}^{r+2m} \\
  			&		 		&  				&  		&			&		& \mathbf{P}^{r+2m-1}& 0\\
 \end{psmallmatrix},
\]
where $\mathbf{P}^k$ are upper-skew-triangular Hankel matrices of sizes $s^k$
where the $(\gamma,\delta)$-entry (where $\gamma,\delta = 0,...,s^k-1$ and $\gamma + \delta < s^k $)
is equal to $a_{k,-(\gamma + \delta + 1)}$.
\end{prop}
\begin{pf}
First, let us consider the following version of the collection formula
\begin{align}
\nonumber
\big( y(v) - y(w) \big)& K_{\widetilde{\zeta}}(v,w)
=
\\ &
\sum_{r=1}^{n}
\sum_{\delta,\gamma = 0} ^{s^r-1}
K_{\widetilde{\zeta}}^{(\gamma)}(v,p^{(r)})
\mathbf{A}_{(r,\gamma,\delta)}
K_{\widetilde{\zeta}}^{(\delta)}(p^{(r)},w)
.
\label{eqCollFor}
\end{align}
To see that \eqref{eqCollFor} holds, we fix $w=w_0 \in X$, such that $w_0 \notin \mathcal{P}(y)$ and $y(w_0)$ has $n$ distinct pre-images.
Then, both sides of \eqref{eqCollFor} are sections of $L_{\widetilde{\zeta}}\otimes\Delta$.
Furthermore, both sides of \eqref{eqCollFor} have poles at $(p^{(j)})_{j=1}^n$ of orders $(s^{j})_{j=1}^n$, respectively,
and the left hand side has an additional removable simple pole at $v=w_0$.
We set the entries of $\mathbf{A}$ to be the residues of $y$
at $(p^{(j)})_{j=1}^n$ such that both sides of \eqref{eqCollFor} share the same principal part.
More precisely, we move to a local parameter $t_r$ at $p^{(r)}$ centered at zero.
Then, the expansion of $y(\cdot)$ at $p^{(r)}$ is given by
\[
y(t_r(u)) = a_{r,-s^r}t_r(u)^{-s^r} +...+ a_{r,-1}t_r(u)^{-1}+...,
\]
while the expansion of the Cauchy kernel $K_{\widetilde{\zeta}}(p^{(r)},w_0)$ is
\[
K_{\widetilde{\zeta}}(p^{(r)},w_0)
=
\sum_{j=0}^{\infty}
K^{(j)}_{\widetilde{\zeta}}(p^{(r)},w_0) t_r(u)^j.
\]
Hence, the $l$-th residue (where $0 < l \leq s^r$) of the LHS of \eqref{eqCollFor} is
\begin{equation}
\label{eqLRes}
\sum_{j=l}^{s^r} a_{-j}K^{(s^r-j)}_{\widetilde{\zeta}}(p^{(r)},w_0)
=
\sum_{j=0}^{s^r-l} a_{-s^r+j}K^{(j)}_{\widetilde{\zeta}}(p^{(r)},w_0)
.
\end{equation}
Using the previous observation, we set the entries of the $r$-th block of $\mathbf{A}$ by
$$\mathbf{A}_{\gamma,\delta}=a_{r,-(\delta+\gamma+1)}$$
and then the residues at $p^{(r)}$ on both sides of \eqref{eqCollFor} coincide and equal to \eqref{eqLRes}.
Since, by assumption, there are no non-zero global holomorphic sections,
we may conclude that \eqref{eqCollFor} holds for any $v \in X$.
\smallskip

It is convenient to use the following identities (see Lemma \ref{mainThLem})
\begin{equation*}
F = (M^{y} - \alpha I) ^{-1} f
\qquad {\rm and} \qquad
G  = (M^{y} - \beta I) ^{-1} g.
\end{equation*}
Recall that using \eqref{StructIdent5}, $F$ and $G$ are the images of resolvent operators acting on $f$ and $g$, respectively.
Thus, the left hand side of the structure identity \eqref{StructIdent4} can then be rewritten as
\begin{align}
\label{StructIdent7}
\innerProductReg{(M^{y}-  \alpha I)^{-1} f}{g} & -
\innerProductReg{f}{(M^{y}-\beta I)^{-1}g} -
\nonumber
\\ &
(\alpha - \overline{\beta})
\innerProductReg{ (M^{y}-\alpha I)^{-1}f}{(M^{y}- \beta I)^{-1}g}
\nonumber
\\ = &
\innerProductReg{
(M^{y}-\beta I)^{-*}
\left( M^{y}-(M^{y})^* \right)
(M^{y}-\alpha I)^{-1}f}{g}
\nonumber
\\ = &
\innerProductReg{
\left( M^{y}-(M^{y})^* \right) F}{G}.
\end{align}

Substituting \eqref{eqCollFor} on the RHS and \eqref{StructIdent7} on the LHS of \eqref{StructIdent4} leads to the following expression
\begin{align}
\nonumber
\innerProductReg{
\left( M^{y}-(M^{y})^* \right) F}{G}
=
\sum_{r=1}^{n}
\sum_{\delta,\gamma = 0} ^{s^r-1}
&
\left(
\sum_{l=1}^{n}
\frac{f(v^{(l)})}{dy(v^{(l)})}
K_{\widetilde{\zeta}}^{(\gamma)}(v^{(l)},p^{(r)})
\right)
\mathbf{A}_{(r,\gamma,\delta)}
    \times \\ \times &
\left(
\sum_{t=1}^{n}
K_{\widetilde{\zeta}}^{(\delta)}(p^{(r)},\tauBa{w^{(t)}})
\frac{\overline{g(w^{(t)})}}
{\overline{dy(w^{(t)})}}
\right)
\label{eq95}
.
\end{align}
Noting that $y(v ^{(l)})=\alpha$, the evaluation of $f$ at $v ^{(l)}$ in terms of $F$
(an equivalent formula ties $G$ and $g$), can be written as follows:
\begin{align*}
f(v^{(l)}) = & (M^y - \alpha) F(v^{(l)})
\\ = &
\sum_{m=1}^{n}
\sum_{l=1}^{s^{m}}
\sum_{j=l}^{s^{m}}
a_{-j,m}
\frac{ F^{(j-l)}(p^{(m)})}{ (j-l)!}
\frac{(-1)^{l-1}}{(l-1)!}
K_{\widetilde{\zeta}}^{(l-1)}(p^{(m)},v^{(l)})
.
\end{align*}
We write one of the expressions on the right hand side of \eqref{eq95}, as follows:
\begin{align}
\nonumber
\sum_{l=1}^{n}
\frac{f(v^{(l)})}{dy(v^{(l)})}
K_{\widetilde{\zeta}}^{(\gamma)}(v^{(l)},p^{(s)})
=
&
\sum_{m=1}^{n}
\sum_{r=1}^{s^{m}}
\sum_{j=r}^{s^{m}}
a_{-j,m}
\frac{ F^{(j-r)}(p^{(m)})}{ (j-r)!}
\frac{(-1)^{r-1}}{(r-1)!}
\times \\ &
\sum_{l=1}^{n}
\frac{K^{(r-1)}_{\widetilde{\zeta}}(p^{(m)},v^{(l)})K_{\widetilde{\zeta}}^{(\gamma)}(v^{(l)},p^{(s)})}{dy(v^{(l)})}
\label{eq34}
.
\end{align}
Differentiating the (simple-poles) collection formula $j$ times with respect to $v$
and $k$ times with respect to $w$ leads to
\begin{align}
\nonumber
\sum_{l=1}^{n}
&
\frac{K^{(j)}_{\widetilde{\zeta}}(v,u^{(l)})K_{\widetilde{\zeta}}^{(k)}(u^{(l)},w)}{dy(u^{(l)})}
   =  \\ & =
\frac{\partial ^{j+k}}{\partial  v ^j \partial w ^k}
\left( \frac{1}{y(w) - \alpha} - \frac{1}{y(v) - \alpha} \right) K_{\widetilde{\zeta}}(v,w)
\label{eqRepad45}
.
\end{align}
Therefore, the inner sum in \eqref{eq34} vanishes whenever $p^{(s)} \neq p^{(m)}$ and
by changing the summation indices and their order, the RHS of \eqref{eq34} becomes
\begin{equation}
\label{eqRepad}
\sum_{j=0}^{s^{m}-1}
F^{(j)}(p^{(s)})
\sum_{r=0}^{s^{m}-1-j}
\frac{(-1)^{r} \, a_{-(j+r+1),s} }{j! r!}
\sum_{l=1}^{n}
\frac{K^{(r)}_{\widetilde{\zeta}}(p^{(s)},\nu^{(l)})K_{\widetilde{\zeta}}^{(\gamma)}(\nu^{(l)},p^{(s)})}{dy(\nu^{(l)})}
.
\end{equation}
Evaluating the RHS of \eqref{eqRepad45} at $w=v=p^{(s)}$, setting $d=s^{s}$
and moving to local parameter $t_s$ of $p^{(s)}$ centered at zero, leads to
\begin{align*}
\frac{\partial ^{j+m}}{\partial v ^j\partial u ^m}
\left( \frac{t_s(u)^{d}}{c_{-d}+...+\hat{c}_0 t_s(u)^{d}} - \frac{t_s(v)^{d}}{c_{-d}+...+\hat{c}_0 t_s(v)^{d}} \right)
\frac{1}{t_s(u)-t_s(v)},
\end{align*}
where $\hat{c}_0 \overset{\text{def} } {=} (c_0 - \alpha)$.
The last term, equal to
\begin{align}
\nonumber
\sum_{l=1}^{n}
&
\frac{K^{(j)}_{\widetilde{\zeta}}(p^{(s)},u^{(l)})K_{\widetilde{\zeta}}^{(m)}(u^{(l)},p^{(s)})}{dy(u^{(l)})}
= \\ & =
\frac{\partial ^{j+m}}{\partial v ^j \partial u^m}
\frac
{
\sum_{l=1}^{d} c_{-l}
\sum_{k=0}^{l-1} t_s(u)^{d-l+k} t_s(v)^{d-1-k} }
{(c_{-d}+...+\hat{c}_0 t_s(u)^{d})(c_{-d}+...+\hat{c}_0 t_s(v)^{d})}
.
\label{eqCollCondQ}
\end{align}
Substituting \eqref{eqCollCondQ} in the RHS of \eqref{eqRepad}, leads to
\begin{align}
\nonumber
\sum_{j=0}^{s^{m}-1}
F^{(j)}(p^{(s)})
\sum_{r=0}^{s^{m}-1-j}
&
\frac{(-1)^{r} \, a_{-(j+r+1),s} }{j! r!}
\\ &
\frac{\partial ^{r+\gamma}}{\partial  v ^r\partial u^ \gamma }
\frac
{
\sum_{l=1}^{d} c_{-l}
\sum_{k=0}^{l-1} t_s(u)^{d-l+k} t_s(v)^{d-1-k} }
{(c_{-d}+...+\hat{c}_0 t_s(u)^{d})(c_{-d}+...+\hat{c}_0 t_s(v)^{d})}
\label{eqRepadQ}
\end{align}
and, by applying Lemma \ref{lemLimitZero} on \eqref{eqRepadQ},
we conclude that \eqref{eqRepad} is just equal to $F^{(\gamma)}(p^{(s)})$.
Similarly, we have
\[
\sum_{t=1}^{n}
K_{\widetilde{\zeta}}^{(\delta)}(p^{(r)},w^{(t)})
\frac{\overline{g(w^{(t)})}}{\overline{dy(w^{(t)})}}
=
\overline{G^{(\delta)}(\tauBa{p^{(r)}})}
.
\]
Finally, we move back to evaluate \eqref{eq95} and we summarize:
\begin{align*}
\innerProductReg{
\left( M^{y}-(M^{y})^* \right) F}{G}
= &
\sum_{r =1}^{n}
\sum_{\gamma = 0}^{s^{r}}
\sum_{\delta = 0}^{s^{r}}
\overline{G^{(\delta)}(\tauBa{p^{(r)}})}
\mathbf{A}_{r, \gamma, \delta}
F^{(\gamma)}(p^{(r)})
\\ = &
\Phi^*_y (G)\, \Per{y} \,  \mathbf{A} \, \Phi_y (F)
\\ = &
\Phi^*_y (G) \, \sigma \, \Phi_y (F)
.
\end{align*}
Here $\Phi_y$
is the poles-evaluation operator associated to $y$ up to derivatives of the poles order minus one,
$\sigma$ is equal to $\mathbf{A}$ up to blocks permutation corresponding to conjugate poles and
we denote the permutation matrix corresponding to a real meromorphic function $y$ by $\Per{y}$.
\end{pf}

We note that the matrix $\sigma_k$ is selfadjoint. It is an immediate consequence of the
selfadjointness of the Hankel matrix $\mathbf{A}$ and the assumption that
$y$ is real (and hence the residues of conjugate poles conjugate).
\smallskip

In the upcoming proofs we extensively use the following notation.

\begin{Notation}
\label{notatMatrix1}
We denote by
$\mathcal P(y_k)$ the set of poles of $y_k$,
we set $n_k = \Card{\mathcal{P}(y_k)}$ for $k=1,2$ and
$n = \Card{\mathcal{P}(y_1) \cup \mathcal{P}(y_2)}$.
The set of poles is denoted by $\left(p^{(m)}\right)_{m=1}^{n}$
and we denote by $s_{k}^{m}$ the order of $y_k$ at a pole $p^{(m)}$
(we set $s_{k}^{m}=0$ when $p^{(m)} \notin \mathcal{P}(y_k)$), for $k=1,2$ and $m=1,...,n$.
Furthermore, we denote the cumulative order of a pole by $s^{m}=s_1^{m}+s_2^{m}$ and the maximum by $\hat{s}^{j} = \max(s_1^{j},s_2^{j})$.

The notations $\Phi_{1}$, $\Phi_{11}$, $\Phi_{12}$ and $\widetilde{\Phi}_{12}$
stand for evaluations at all the poles in
$\mathcal P (y_1)$, $\mathcal P (y_1) \setminus \mathcal P (y_2)$, $\mathcal P (y_1) \cap \mathcal P (y_2)$ and $\mathcal P (y_1) \cap \mathcal P (y_2)$,
respectively, of an $L_{\widetilde{\zeta}}\otimes \Delta$ section $f$ and its derivatives up to orders of
$s^{m}_1-1$, $s^{m}_1-1$, $\hat{s}^{j}-1$ and $s^{m}-1$, respectively.

The block matrix $\mathbf{K}_{\zeta}(P_{1},P_2)$ is defined by
(similarly, we may define the vector $\mathbf{K}_{\zeta}(P_{1},u)$):
\[
\left[ K_{\zeta}^{(k_1-1,k_2-1)}(p^{(n)}_1,p^{(m)}_2)\right]_{(n,k_1),(m,k_2)},
\quad
p^{(n)}_1 \in \mathcal P (y_1)  \, \, {\rm  and  } \,\,  p^{(m)}_2 \in \mathcal P (y_2)
,
\]
where $n=1,...,n_1$, $k_1=0,...,s_1^{n}$, $m=1,...,n_2$ and $k=0,...,s_2^{m}$.
In the case where $\mathbf{K}_{\zeta}(\cdot,\cdot)$ contains an entry
which corresponds to $p^{(n)}_1=p^{(m)}_2$, then this entry is set to be zero.
\end{Notation}

\begin{Rk}
Using Notation \ref{notatMatrix1}, the model operator \eqref{MyMatrix} can be written in the following form
\begin{align*}
M^{y_k}f(u)
=
y_k(u)f(u) + i \mathbf{K}_{\widetilde{\zeta}}(u,P_k) \, \mathbf{A}_{y_k} \, \Phi_k(f).
\end{align*}
Here $\mathbf{A}_{y_k}$ is equal to upper-skew-triangular Hankel blocks matrix
which is presented in the proof of Proposition \ref{myGeneral1}.
\end{Rk}

\begin{pf}[of Theorem \ref{lemCollCodI}]
Using Proposition \ref{myGeneral1}, it is sufficient to prove that for an arbitrary polynomial $g$,
the meromorphic function $z_1 = g(y_1,y_2)$ satisfies the colligation condition.
It is clearly enough to that if the colligation condition holds for $M^{y_1}$ and $M^{y_2}$
then it holds for $M^{y_1}+M^{y_2}$ and $M^{y_1}M^{y_2}$.
\smallskip

First, by summing the colligation conditions corresponding to $y_1$ and $y_2$, we simply have
\begin{equation*}
M^{y_1} + M^{y_2} - M^{{y_1}*} - M^{{y_2}*}
=
i \Phi_1 ^* \sigma_1 \Phi_1 + i \Phi_2 ^* \sigma_2 \Phi_2
=
i \Phi ^* \sigma \Phi,
\end{equation*}
where $\Phi = \Phi_{1} \oplus \Phi_{12} \oplus \Phi_{2}$, see Notation \ref{notatMatrix1},
and
$\sigma = \sigma_{11} \oplus \widehat{\sigma}_{12} \oplus \sigma_{22}$.
The matrix $\sigma_{11}$ contains the Laurent coefficients of $y_1$ at all poles in $\mathcal P (y_1) \setminus \mathcal P (y_2)$.
$\widehat{\sigma}_{12}$ is the summation of the Laurent coefficients of $y_1$ and $y_2$ along the joint poles.
As for the case where a joint pole does not belong to $\mathcal P ( y_1 + y_2)$,
the matrix $\widehat{\sigma}_{12}$ contains zeros at the corresponding entries.
Thus, the colligation condition of $y_1 + y_2$ follows.
\smallskip

Moving to examine $M^{y_1}M^{y_2}$, one can show, using the commutativity of $M^{y_1}$ and $M^{y_2}$  and the colligation conditions associated $M^{y_1}$ and $M^{y_2}$ (Proposition \ref{myGeneral1}), the following
\begin{align}
\nonumber
\frac{1}{i}
&
\big(
\innerProductReg{M^{y_1}M^{y_2}f}{g} - \innerProductReg{f}{M^{y_1}M^{y_2}g}
\big)
\\
\nonumber
&
=
\innerProductReg{\sigma_1 \Phi_1 M^{y_2} f }{\Phi_1 g} +
\frac{1}{i}\big(\innerProductReg{M^{y_2}f}{M^{y_1}g}- \innerProductReg{f}{M^{y_1}M^{y_2}g}\big)
\\
\nonumber
&
=
\innerProductReg{\sigma_1 \Phi_1 M^{y_2} f }{\Phi_1 g} +
\frac{1}{i}\big(\innerProductReg{M^{y_2}f}{M^{y_1}g}- \innerProductReg{f}{M^{y_2}M^{y_1}g}\big)
\\
\label{mymyColl}
&
=
\innerProductReg{\sigma_1 \Phi_1 M^{y_2} f }{\Phi_1 g} +
\innerProductReg{\sigma_2 \Phi_2 f }{\Phi_2 M^{y_1} g}
.
\end{align}
For the sake of simplicity, we first assume that $\mathcal P (y_1)\cap \mathcal P (y_2)=\emptyset$.
It follows that the expression in \eqref{mymyColl}, using the matrix representations as presented in \eqref{MyMatrix}, becomes
\begin{align}
\nonumber
\frac{1}{i}
\big(
\innerProductReg{M^{y_1}M^{y_2}f}{g} -
&
\innerProductReg{f}{M^{y_1}M^{y_2}g}
\big)
= \\ =  &
\nonumber
\Phi^*_{1}(g)
\sigma_1
\left(
\Phi_{1}(y_2 \cdot f)
+
i
\mathbf{K}_{\widetilde{\zeta}}(P_1,P_2)
\mathbf{A}_{y_2}
\Phi_{2}(f)
\right)
+
\\
&
\left(
\Phi^*_{2}(y_1 \cdot g)
-
i
\Phi^*_{2}(g)
\mathbf{A}_{y_1}^*
K^*_{\widetilde{\zeta}}(P_2,P_1)
\right)
\sigma_2
\Phi_{2}(f)
\label{eqCollMultE}
.
\end{align}
We note that the selfadjoint matrices $\sigma_k$ and $ \mathbf{A}_{y_k}$ are equal up to a blocks-permutation matrix corresponding to conjugate poles
(compare \eqref{MyMatrix} and Proposition \ref{myGeneral1}), i.e. $\sigma_k = \Per{y_k} \mathbf{A}_{y_k}$
and hence we note that
$\sigma_k \mathbf{K}_{\widetilde{\zeta}}(P_k,u) = \mathbf{A}_{y_k} \mathbf{K}_{\widetilde{\zeta}}(\tauBa{P_k},u)$ holds.
Moreover, using the Hermitian structure of the Cauchy kernels, we have
$K^*_{\widetilde{\zeta}}(P_2,\tauBa{P_1})=\mathbf{K}_{\widetilde{\zeta}}(P_1,\tauBa{P_2})$ (see Notation \ref{notatMatrix1}).
Then, one may show that two of the terms in \eqref{eqCollMultE} vanish
as follows from the next computation:
\begin{align}
\nonumber
\frac{1}{i}
\big(
&
\innerProductReg{M^{y_1}M^{y_2}f}{g}
-
\innerProductReg{f}{M^{y_1}M^{y_2}g}
\big)
\\ =  &
\nonumber
\Phi^*_{1}(g)
\sigma_1
\Phi_{1}(y_2 \cdot f)
+
\Phi^*_{2}( g \cdot y_1)
\sigma_2
\Phi_{2}(f)
\\
&
i
\Phi^*_{1}(g)
\left(
\mathbf{A}_{y_1}
\mathbf{K}_{\widetilde{\zeta}}(\tauBa{P_1},P_2)
\mathbf{A}_{y_2}
-
\mathbf{A}_{y_1}
\mathbf{K}_{\widetilde{\zeta}}(\tauBa{P_1},\tauBa{P_2})
\sigma_2
\right)
\Phi_{2}(f)
\nonumber
\\ =  &
\nonumber
\Phi^*_{1}(g)
\sigma_1
\Phi_{1}(y_2 \cdot f)
+
\Phi^*_{2}( g \cdot y_1)
\sigma_2
\Phi_{2}(f)
\\
&
i
\Phi^*_{1}(g)
\left(
\mathbf{A}_{y_1}
\mathbf{K}_{\widetilde{\zeta}}(\tauBa{P_1},P_2)
\mathbf{A}_{y_2}
-
\mathbf{A}_{y_1}
\mathbf{K}_{\widetilde{\zeta}}(\tauBa{P_1},P_2)
\mathbf{A}_{y_2}
\right)
\Phi_{2}(f)
\nonumber
\\ = &
\Phi^*_{1}(g)
\sigma_1
\Phi_{1}(y_2 \cdot f)
+
\Phi^*_{2}( y_1 \cdot g)
\sigma_2
\Phi_{2}(f)
\label{eqCollMultT}
.
\end{align}
Equivalently, using the decomposition $\Phi_1(y_2 \cdot f) = \Phi_1(y_2)\Phi_1(f)$, where $\Phi_1(y_2)$
is a (Hankel-block) matrix containing the values of $y_2$ and its derivatives at the poles of $y_1$,
\eqref{eqCollMultT} can be written as follows:
\[
\Phi^*_{1}(g)
\sigma_1
\Phi_1(y_2)
\Phi_{1}( f)
+
\Phi^*_{2}( g)
\Phi^*_2(y_1)
\sigma_2
\Phi_{2}(f)
=
\Phi^*(g)
\sigma_{y_1 \, y_2}
\Phi(f)
,
\]
where $\Phi = \Phi_1 \oplus \Phi_2$ and $\sigma_{y_1 \, y_2} = \sigma_1 \Phi_1(y_2) \oplus \Phi_2^*(y_1) \sigma_2$.
Hence, the colligation condition, associated to the model operator acting on the multiplication of
two meromorphic functions with disjoint and non-simple poles, follows.
\smallskip

We now turn to the case of a pair of meromorphic functions with common poles, and we start by considering the term $\Phi_1 M^{y_2} f  $.
Let us evaluate the $k$-th derivative of $M^{y_2}f(u)$ at a non-simple pole (of $y_2$), say $p^{(r)}$ of order $s^{r}$ ($k < s^{r}$).
Considering a local parameter $t$ at $p^{(r)}$ centered at zero, we have
\begin{align}
\nonumber
M^{y_2}f(t(u))
=
&
y(t(u))f(t(u))
+
i
\mathbf{K}_{\widetilde{\zeta}}(t(u),P_2)
\mathbf{A}_{y_2}
\Phi_2 f
\\ = &
\label{eq159}
T^{(0)}(p^{(r)})
+
\sum_{n=0}^{\infty}
t(u)^n
\sum_{j=0}^{n+s^r}
a_{n-j}
\frac{f^{(j)}(p^{(r)})}{j!}
,
\end{align}
where
\[
T^{(k)}(p^{(r)})= i
\sum_{\substack{m=1 \\ m \neq r }}^{n}
\sum_{l=1}^{s^{m}}
\sum_{j=l}^{s^{m}}
a_{-j,m} \frac{  f^{(j-l)}(p^{(m)}) }{ (j-l)!}
\frac{(-1)^{l-1}}{(l-1)! }
K^{(k,l-1)}_{\widetilde{\zeta}}(p^{(r)},p^{(m)})
\]
and thus the $k$-th derivative of \eqref{eq159} is
\begin{align}
\nonumber
\restr{\frac{d^k }{d\,u ^k}M^{y}f(u)}{u=p^{(r)}}
=
&
\frac{d^k }{d\,u ^k} y(t(u))f(t(u))
+
i
\frac{d^k }{d\,u ^k} \mathbf{K}_{\widetilde{\zeta}}(t(u),P_2)
\mathbf{A}_{y_2}
\Phi_2 f
\\ = &
\nonumber
k!
\sum_{j=0}^{s^r} a_{k-j,r}\frac{f^{(j)}(p^{(r)})}{j!}+
T^{(k)}(p^{(r)}).
\\ = &
\label{eq44}
k!
\sum_{j=0}^{s^r} a_{k-j,r}\frac{f^{(j)}(p^{(r)})}{j!}+
\widetilde{K}_{\widetilde{\zeta}}(p^{(r)},P_2)
\mathbf{A}_{y_2}
\Phi_2 f
.
\end{align}
As a consequence, we note that
$\Phi_1 M^{y_2} f$ (and similarly $\Phi_2 M^{y_1} g$), as expected,
contains evaluations of the derivatives of $f$ (and $g$)
up to order $s^{r}$(=$s_1^{r} + s_2^{r}$).
Accordingly, we set $\mathbf{Q}_2$ (and similarly $\mathbf{Q}_1$) to be a matrix of blocks of sizes
$s^{r}$ where $p^{(r)} \in \mathcal P (y_1) \cap \mathcal P (y_2)$ with entries extracted from \eqref{eq44}.
Then, using the same arguments as in the $\mathcal P (y_1)\cap \mathcal P (y_2)=\emptyset$ case, \eqref{eqCollMultE} is now given by
\begin{align*}
\frac{1}{i}
&
\big(
\innerProductReg{M^{y_1}M^{y_2}f}{g} -
\innerProductReg{f}{M^{y_1}M^{y_2}g}
\big)
\\ =  &
\Phi^*_{11}(g)
    \sigma_{11}
\Phi_{11}(y_2 \cdot f)
+
\Phi^*_{22}( g \cdot y)
    \sigma_{22}
\Phi_{22}(f)
+
\widetilde{\Phi}_{12}^* (\sigma_1 \mathbf{Q}_2 + \mathbf{Q}_1 \sigma_1^*) \widetilde{\Phi}_{12}
+
\\
&
i
\Phi^*_{1}(g)
\left(
\sigma_1
\mathbf{K}_{\widetilde{\zeta}}(P_1,P_2)
\mathbf{A}_{y_2}
-
\mathbf{A}_{y_1}
\mathbf{K}^*_{\widetilde{\zeta}}(P_2,P_1)
\sigma_2
\right)
\Phi_{2}(f)
\\ = &
\Phi^*_{11}(g)
    \sigma_{11}
\Phi_{11}(y_2 \cdot f)
+
\Phi^*_{22}( g \cdot y_1)
    \sigma_{22}
\Phi_{22}(f)
+
\widetilde{\Phi}_{12}^* (\sigma_1 \mathbf{Q}_2 + \mathbf{Q}_1 \sigma_1^*) \widetilde{\Phi}_{12}
.
\end{align*}
Equivalently, by using the notation $\Phi = \Phi_{11} \oplus \widetilde{\Phi}_{12} \oplus \Phi_{22}$,
\eqref{mymyColl} is just
\begin{align*}
M^{y_1y_2} & - M^{(y_1y_2)*}
=
i \Phi ^*
\left[
\sigma_{11} \Phi_{11}(y_2) \oplus
(\sigma_1 \mathbf{Q}_2 + \mathbf{Q}_1 \sigma_2^*)
\oplus \Phi^*_{22}(y_1) \sigma_{22}
\right]
\Phi
\end{align*}
and the colligation condition associated to $M^{y_1 y_2}$ follows by setting
$$\sigma_{y_1 y_2}
\overset{\text{def} } {=}
\sigma_{11} \Phi_{11}(y_2) \oplus
(\sigma_1 \mathbf{Q}_2 + \mathbf{Q}_1 \sigma_2)
\oplus \Phi^*_{22}(y_1) \sigma_{22}.$$
\end{pf}
To conclude this section, we present the proof of Theorem \ref{vesselAreEq}. \\
\begin{pf}[of Theorem \ref{vesselAreEq}]
To prove the validity of the right most equality in \eqref{eqVesselEq},
we first mention that the inverse of an upper-skew-triangular Hankel matrix $G$
(the result can be extended naturally to a skew-triangular Hankel blocks matrices)
of size $s$ with entries $G_{i,\ell} = g_{i+\ell+1}$ where $ i + \ell < s$,
is the following lower-skew-triangular Hankel matrix:
$$ \left(G^{-1}\right)_{i,\ell} =
\restr{ \frac{d^{i+\ell-s+1}}{dx^{i+\ell-s+1}} \frac{1}{(i+\ell-s)! g(x)} }{x = 0} \quad  0 \leq \ell, i <s-1, \quad \ell+ i \geq s -1,$$
where $g(x) = \sum_{i=1}^{s} g_{s-i} x^i$.
This is true since for $\alpha \leq \beta$ (otherwise the entry is zero) we have
\begin{align}
\nonumber
\left(G G^{-1}\right)_{\alpha,\beta}
= &
\sum_{i = s-\beta}^{s-\alpha}
G_{\alpha,i} G^{-1}_{i,\beta}
=
\sum_{i = 0}^{\beta - \alpha}
G_{\alpha,s-\beta+i-1} G^{-1}_{s-\beta+i-1,\beta}
\\ = &
\nonumber
\sum_{i = 0}^{\beta - \alpha}
g_{s-(\beta-\alpha -i)}
\restr{ \frac{d^i}{dx^{i}} \frac{1}{i! g(x)}}{x=0}
\\ = &
\sum_{i = 0}^{\beta - \alpha}
\restr{
\frac{d^{\beta-\alpha-i}}{dx^{\beta-\alpha-i}} \frac{g(x)}{(\beta-\alpha-i)!} \frac{d^i}{dx^{i}} \frac{1}{i! g(x)}}{x=0}
=
\delta_{\beta-\alpha}.
\label{eqHankelInv}
\end{align}

We define the unitary mapping $U$ by \eqref{ModelMap0}, therefore we have
\[
U \, \colon \, h \mapsto  \varphi_h(z) = {\bf u}_l^{\times} (z) (\xi \sigma)\Phi (\xi A - \xi y(z) )^{-1}h,
\]
where $\xi \sigma$, $\xi A$ and $\xi y(z)$ stand for $\xi_1 \sigma_1+\xi_2\sigma_2$, $\xi_1A_1+\xi_2A_2$ and $\xi_1y_1(z) +\xi_2y_2(z)$, respectively.
We fix an arbitrary element $h\in H$,
we consider a local coordinate $t_j(u)$ centered at $p^{(j)} \in \mathcal P (y_1 ) \cup \mathcal P (y_2 )$
and we fix $(\xi_1,\xi_2)$ such that $\xi y(z)$ are of order $\hat{s}^{j}$ (see Notation \ref{notatMatrix1})	.
Then, $(\xi A - \xi y(z) )^{-1} h $ has a zero of order $\hat{s}^{j}$ at $p^{(j)}$
and with the Laurent coefficients of $\frac{1}{\xi y(z)}$.
We recall that ${\bf u}_l^{\times} (z)$ has poles at $(p^{(j)})$ of order $\hat{s}^{j}$.
It follows that $\varphi_h(z)$ is analytic at the poles of $\xi y$ and $\Phi_{y}\varphi_h(z)$ contains the
Taylor coefficients of $\varphi_h(z)$ at the poles of $\xi y(z)$.
Then, using $(g_{k})_{k=1}^{\hat{s}^{j}}$ to denotes Laurent
coefficients of the singular part of $\xi  y$ at $p^{(j)}$, we have
\begin{equation}
\label{eq123asd}
\Phi_{y} \varphi(z)
=
U_l^\times \mathbf{A}_{y}^{-1} (\xi \sigma) \Phi h
.
\end{equation}
Here $U_l^\times$ denotes the $n\times n$ matrix (actually, it is a block-matrix of blocks of sizes $\hat{s}_j$)
that consists of the coefficients of the singular part of ${\bf u}_l^{\times} (z)$, that is,
\[
U_l^\times =  \col_{j=1,\ldots,n} \col_{k=1,\ldots,\hat{s}_j}({\bf u}_l^{\times})_{(-k, p^{(j)})}.
\]
Furthermore, by multiplying $\frac{1}{\xi y (z)}$ by the local coordinates of the singularities of ${\bf u}_l^{\times} (z)$
we get entries which are exactly
of the form $\restr{\frac{1}{{\xi y (z)}}}{z=0}$ and hence, by \eqref{eqHankelInv}, is equal to the inverse matrix of $\mathbf{A}_{y}$
(presented in \eqref{MyMatrix} above).
$\Phi_{y} \varphi(z)$ in \eqref{eq123asd} does not depend on the selection of $(\xi_1,\xi_2)$ since \eqref{ModelMap0}
does not depend on $(\xi_1,\xi_2)$.
\smallskip

We note that when one considers the canonical determinantal representation and where $y$ has simple poles,
then $U_l^\times$ and $(\xi \sigma) = \mathbf{A}_y$.
It follows that the operator $\Phi$ in \eqref{eq123asd} is the evaluation operator
at the poles of $y_1$ and $y_2$, i.e. $\Phi_y (\varphi_h) = \Phi h$.
\smallskip

To verify the validity of the two remaining identities in \eqref{eqVesselEq},
we illustrate and continue the calculation for the operator $A_1$ and the direction $\xi=(1,0)$:
\begin{align}
\nonumber
\varphi_{A_1 h}(z   )
= &
{\bf u}_l^{\times} (z)
\sigma_1
\Phi
( A_1 - y_1(z) )^{-1} A_1 \, h
\\ = &
{\bf u}_l^{\times} (z)
\sigma_1 \Phi
( A_1 - y_1(z) )^{-1}
\left( A_1 - y_1(z) +y_1(z) \right)  \, h
\nonumber
\\ = &
y_1(z)
\varphi_h(z)
+
{\bf u}_l^{\times} (z) \sigma_1 \Phi h.
\label{eq123asd1}
\end{align}
Substituting \eqref{eq123asd} where $\xi = (1,0)$ into \eqref{eq123asd1}, leads to
\begin{align*}
\varphi_{A_1 h}(z)
= &
y_1(z) \varphi_h(z)
+
i
\mathbf{K}_{\widetilde{\zeta}}(z,P_1)
\mathbf{A}_{y_1}
\Phi_{y_1}
\varphi_h
.
\end{align*}
This is, see \eqref{MyMatrix}, the definition of the model operator $M^{y_1} \varphi_h$ corresponding to $y_1$.
\smallskip

Finally, it remains to show that the Hilbert space $H_T = \{ \varphi_h \, : \, h \in H \}$
with the inner product inherited form $H$,
is a reproducing kernel Hilbert space with the reproducing kernel \eqref{eqKernPos}.
We use the properties of the characteristic functions of the vessel $\mathcal V$;
\cite[Equation N2]{MR97m:30051}, the identity \eqref{st} and
an identity which is presented in the next section,
see \eqref{eqsec5A} below, to show the following
\begin{align*}
K_T(p,q) = & K_{\widetilde{\zeta}}(p,q) - T(p) K_{\zeta}(p,q) T(q)^*
\\ = &
\frac{\widetilde{\bf u}^{\times} _l (p)(\xi\sigma){\widetilde{\bf u}}^{\times} _l (q)^* }
{-i  (\xi y(p) - \xi \overline{y(q)}) }
-
\frac
{ T(p) {\bf u}^{\times} _l (p) ( \xi \sigma ) {\bf u}^{\times} _l (q)^* T(q)^* }
{-i  (\xi y(p) - \xi \overline{y(q)}) }
\\ = &
\widetilde{\bf u}^{\times} _l (p)
\frac{(\xi\sigma) - S(\xi,p) ( \xi \sigma ) S(\xi,q)^* }
{-i  (\xi y(p) - \xi \overline{y(q)}) }
\widetilde{\bf u}^{\times} _l (q)^*
\\ = &
\widetilde{\bf u}^{\times} _l (p)
(\xi\sigma)
\Phi
(\xi A - \xi y(p))^{-1}
(\xi A - \xi y(q))^{-*}
\Phi^*
(\xi \sigma)^*
\widetilde{\bf u}^{\times} _l (q)^*
.
\end{align*}
Thus, we may conclude that
\begin{align*}
\innerProductTri{\varphi_h(p)}{ K_T(p,q)}{H_T}
= &
\innerProductTri{h}{
(\xi A - \xi y(q))^{-*}
\Phi^*
(\xi \sigma)^*
\widetilde{\bf u}^{\times} _l (q)^*
}{H}
\\ = &
\widetilde{\bf u}^{\times} _l (q)
(\xi \sigma)
\Phi
(\xi A - \xi y(q))^{-1}
h
\\ = &
\varphi_h(q).
\end{align*}
\end{pf}
\section{Proof of the structure theorem}
\label{secProofStruct}

We prove Theorem \ref{preTh} in a number of steps.
As mentioned in Section \ref{secMainThm},
our strategy is to embed the operators $M^{y_1}$ and $M^{y_2}$
in a commutative two-operator vessel of the form
\begin{equation}
\label{vessel123}
( \, M^{y_1} \, ,  \,  M^{y_2} \, ; \, \mathcal X \, , \,  \Phi \, , \, E \, ; \,  \sigma_1 \, , \,  \sigma_2 \, , \, \gamma  \, , \, \widetilde{\gamma} \, ),
\end{equation}
where $n = \Card{\mathcal{P} (y_1) \cup \mathcal{P}(y_2)}$, $E=\mathbb{C}^n$, and where
$\Phi$ is the evaluation operator from $\mathcal{X}$ to $E$ at $\infty$, namely,
\begin{equation*}
f \rightarrow \
\begin{pmatrix}
  \Phi_{p_1}f \\
  \vdots \\
    \Phi_{p_n}f \\
 \end{pmatrix}=
 \begin{pmatrix}
  f(p^{(1)}) \\
  \vdots \\
  f(p^{(n)}),
 \end{pmatrix}.
\end{equation*}
Recall that (see \cite{MR1704479}) since $y_1$ and $y_2$ are real meromorphic functions generating $\mathcal M(X)$,
there exists a bi-rational embedding of $X$ to a projective curve $C \subseteq \mathbb P ^2$ of degree $n$,
which is exactly the discriminant curve associated to the vessel \eqref{vessel123}.
\smallskip

Before proceeding to the proof, we make some remarks.
\begin{Rk}
\label{T0}
In our setting, $\widetilde{\zeta} \in J(X)$ does not necessarily belong to $T_0$.
Therefore, for each component $X_j$ where $0 \leq j \leq k-1$, one may attach a sign $(-1)^{\mu_j}$ with $\mu_j = 0$ or $1$
(when $\widetilde{\zeta} \in T_0$, $\mu_j=0$ for all $j$)
such that ${\rm sign} \, K_{\widetilde{\zeta}}(q,q) = (-1)^{\mu_j}$ where $q$ is
close to $X_j$ on the positive side with respect to the chosen orientation.
For more details we refer to \cite{av3}.
\end{Rk}
\begin{Rk}
The meromorphic functions $y_1(z)$ and $y_2(z)$ are not necessarily of dividing type
(see Definition \ref{dividFunc}; i.e. not necessarily maps $X^+$ to the upper half-plane and $X^-$ to the lower half-plane).
As a result, the matrices $\sigma_1$ and $\sigma_2$ in vessel \eqref{vessel123} are
$1 \times 1$ and $2 \times 2$ block--diagonal matrices (see \eqref{sigma123} below),
rather than having a strictly diagonal structure, see also \cite{vinnikov5}.
\end{Rk}
\begin{pf}[of Theorem \ref{preTh}]
The proof consists of seven steps.
We start by recalling the relation between the colligation condition and the structure identity.

\begin{Step}
\label{step1}
The colligation conditions \eqref{collCond} of the model vessel \eqref{vessel123} are equivalent to
the structure identities (condition \ref{assump3}, Theorem \ref{preTh}),
where $\sigma_k$ for $k=1,2$ are given in \eqref{sigma123} below.
\end{Step}

This result is just a special case of Proposition \ref{myGeneral1} and hence the following identity holds
\begin{equation}
\label{eqCollCondAA}
\frac{1}{i}
\innerProductTri{\left( M^{y_k}-(M^{y_k})^* \right)F}{G}{\mathcal X}
=
\innerProductTri{\sigma_k \Phi_k (F)}{\Phi_k (G)}{E_k}
\qquad
k=1,2
.
    \end{equation}
Here $E_k \simeq \mathbb C ^{\card{\mathcal{P}(y_k)}}$ and $\sigma_k$ is given by
\begin{equation*}
\sigma_k =
 \begin{psmallmatrix}
  c_k^1		& 				&  				&		&			&		$		$\\
  			&	\ddots 		&  				&  		&			&		$		$\\
  			&			 	& 	c_k^r		&  		&			&		$		$\\
  			&			 	& 				& 0		& c_k^{r+2}	&		&		&\\
  			&		 		&  				& c_k^{r+1}& 0		&		$		$\\
  			&		 		&  				&  		&			&\ddots	$		$\\
  			&		 		&  				&  		&			&		&0	& c_k^{r+2m} \\
  			&		 		&  				&  		&			&		& c_k^{r+2m-1}& 0\\
 \end{psmallmatrix},
\end{equation*}
where $c_k^j$ is the residue of $y_k$ at $p^{(j)}$ and $\Phi_k$ is the evaluation operator at the poles of $y_k$.
\smallskip

We note that one may define a single evaluation operator, $\Phi$ of size $n$,
at the poles of $y_1$ and $y_2$.
The matrices, abusing notations, $\sigma_1$ and $\sigma_2$ are extended to size $\mathbb C^{n \times n}$,
where now an entry $c_k^j$ is the residue of $y_k$ at $p^{(j)}$ which is set to be zero whenever $p^{(j)}$ is not a pole of $y_k$.
We note that the matrices $\sigma_1$ and $\sigma_2$ are selfadjoint, since $y_1$ and $y_2$ are real.
Hence, we may rewrite \eqref{eqCollCondAA} as
\[
\frac{1}{i}
\innerProductTri{\left( M^{y_k}-(M^{y_k})^* \right)F}{G}{\mathcal X}
=
\innerProductTri{\sigma_k \Phi (F)}{\Phi (G)}{E}
\qquad
k=1,2
\]
and the colligation conditions \eqref{collCond} associated to the model vessel \eqref{vessel123}
hold with the corresponding $\sigma_1$ and $\sigma_2$.
\smallskip

Furthermore, in the view of Remark \ref{T0},
$\widetilde{\zeta}$ does not necessarily belong $T_0$.
Then for each component $X_j$ we associate a sign $(-1)^{\mu_j}$ such that
${\rm sign }\, K_{\widetilde{\zeta}}(q,q) = (-1)^{\mu_j}$ where $q \in X_j$ and $\mu_j \in \{0,1\}$.
Hence by setting
\begin{equation}
\label{sigma123}
\sigma =
\begin{pmatrix}
\sigma_{\mathbb R} & 0_{r \times 2m} \\
0_{2m \times r} & \sigma_{\mathbb C}\\
\end{pmatrix},
\end{equation}
where
\begin{equation*}
\sigma_{\mathbb R} =
\begin{pmatrix}
  (-1)^{\mu_1}c^1			& 					\\
  			&	\ddots 		&  					\\
  			&			 	& 	(-1)^{\mu_r}c^r \\
\end{pmatrix},
\end{equation*}
and
\begin{equation*}
\sigma_{\mathbb C} =
\begin{psmallmatrix}
 0		& (-1)^{\mu_{r+1}}c^{r+2}	&		&		&\\
 (-1)^{\mu_{r+1}}c^{r+1}& 0		&		$		$\\
 				&  		&			&\ddots	$		$\\
  				&  		&			&		&0	& (-1)^{\mu_{r+m}}c^{r+2m} \\
  				&  		&			&		& (-1)^{\mu_{r+m}}c^{r+2m-1}& 0\\
\end{psmallmatrix},
\end{equation*}
the colligation condition \eqref{collCond}, holds also for $\widetilde{\zeta} \in T_{\nu}$ where $\nu \in \mathbb N$.

\begin{Step}
\label{step2}
Let $\widetilde{\gamma}$ be defined by
\begin{equation}
\label{tildeGamma}
\widetilde{\gamma}_{j,l} =
\begin{cases}
c_{2}^{l} h_{1}^{l} - c_{1}^{l} h_{2}^{l},
&
\tauBa{p^{(j)}} = p^{(l)}
\\
\bracketsA{\overline{c_{2}^{j}} c_{1}^{l} - \overline{c_{1}^{j}} c_{2}^{l} }
\frac
{\vartheta[\widetilde{\zeta}]( p^{(l)}- \tauBa{p^{(j)}})}
{\vartheta[\widetilde{\zeta}] (0)E(p^{(l)},\tauBa{p^{(j)}})},
&
{\text otherwise},
\end{cases}
\end{equation}
where $c_{k}^{j}$ and $h_{k}^{j}$, for a $p^{(j)} \in \mathcal{P}(y_k)$,
are given by the expansion of $y_k$ at $p^{(j)}$ using local coordinate $t^j$:
\[
y_k(u) = - \frac{c^{m}_{k}}{t^j(u)}+h^{m}_{k} + O(|t^j|),
\]
otherwise ($p^{(j)} \notin \mathcal{P}(y_k)$), we have $c_{k}^{j}=0$ and $h_{k}^{j} = y_k(p^{(j)})$.
Then, the output vessel condition
\begin{equation}
\label{eq4_3}
\sigma_2\Phi M^{y_1} - \sigma_1\Phi M^{y_2} = \widetilde{\gamma}\Phi
\end{equation}
holds for $\sigma_1$ and $\sigma_2$ as defined in \eqref{sigma123}.
\end{Step}

First, we note that $\mathbf{A}_k$ is a diagonal matrix
(block-diagonal in a non-simple setting)
and furthermore that \eqref{eq4_3} becomes
\begin{equation}
\label{eq4_3A}
\mathbf{A}_2 \Phi M^{y_1} - \mathbf{A}_1 \Phi M^{y_2} = \Per{y} \widetilde{\gamma}\Phi
.
\end{equation}
Below, for the sake of simplicity,
we begin by calculating the $j$-th entry of the left hand side of \eqref{eq4_3A}
for the case where $p^{(j)} \in \mathcal{P}(y_2) \setminus \mathcal{P}(y_1) $,
we use $h_1^j=y_1(p^{(j)})$ and $c_1^j=0$ to conclude:
\begin{align}
\nonumber
[LHS]_j
= &
c_2^j
\left(
y_1(p^{(j)})f(p^{(j)})
+
\sum_{ p^{(l)} \in \mathcal{P}(y_1)}{c_1^l f(p^{(l)})
\frac
{\vartheta[\widetilde{\zeta}](p^{(l)}-p^{(j)})}
{\vartheta[\widetilde{\zeta}] (0)E(p^{(l)},p^{(j)})}
}
\right)
-
\\ & \nonumber
c_1^j
\left(
y_1(p^{(j)})f(p^{(j)})
+
\sum_{ p^{(l)} \in \mathcal{P}(y_2)}{c_2^l f(p^{(l)})
\frac
{\vartheta[\widetilde{\zeta}](p^{(l)}-p^{(j)})}
{\vartheta[\widetilde{\zeta}] (0)E(p^{(l)},p^{(j)})}
}
\right)
\\ = &
c_2^j
h_1^j
f(p^{(j)})
+
c_2^j
\sum_{ p^{(l)} \in \mathcal{P}(y_1)}
c_1^l f(p^{(l)})
\frac
{\vartheta[\widetilde{\zeta}](p^{(l)}-p^{(j)})}
{\vartheta[\widetilde{\zeta}] (0)E(p^{(l)},p^{(j)})}
\label{step2_2B}
.
\end{align}
For the case where $p^{(j)} \in \mathcal{P}(y_2) \cap \mathcal{P}(y_1) $,
we move to local coordinate and using the prime form properties:
\begin{align}
[\mathbf{A}_2 \Phi M^{y_1}f]_j
\nonumber
= &
c_2^j
 f(p^{(j)})
\bracketsA{
\frac{-c^{j}_{1}}{t(u)}+h^{j}_{1} + o(|t|) +
\frac
{c_1^j}
{t(u)}
} +
\\
\nonumber
&
c_2^j
\sum_{ p^{(j)} \neq p^{(l)} \in \mathcal{P}(y_1)  }
c_1^l f(p^{(l)})
\frac
{\vartheta[\widetilde{\zeta}](p^{(l)}-p^{(j)})}
{\vartheta[\widetilde{\zeta}] (0)E(p^{(l)},p^{(j)})}
.
\end{align}
Therefore, the LHS of \eqref{eq4_3A}, in the case of a joint pole, is given by
\begin{align}
\nonumber
[LHS]_j
=
\left(
c_2^j h_1^j
-
c_1^j
h_2^j
\right) f(p^{(j)})
& +
c_2^j
\sum_{ \substack{p^{(l)} \in \mathcal{P}(y_1) \\  p^{(j)} \neq p^{(l)} } }
c_1^l f(p^{(l)})
\frac
{\vartheta[\widetilde{\zeta}](p^{(l)}-p^{(j)})}
{\vartheta[\widetilde{\zeta}] (0)E(p^{(l)},p^{(j)})}
\\ & -
c_1^j
\sum_{ \substack{p^{(l)} \in \mathcal{P}(y_2) \\  p^{(j)} \neq p^{(l)} } }
c_2^l f(p^{(l)})
\frac
{\vartheta[\widetilde{\zeta}](p^{(l)}-p^{(j)})}
{\vartheta[\widetilde{\zeta}] (0)E(p^{(l)},p^{(j)})}
\label{step2_1}
.
\end{align}
One may note,
using that $c_k^j = 0$ when $p^{(j)} \notin \mathcal{P}(y_k)$,
The terms in \eqref{step2_2B} and \eqref{step2_1} can be jointly written as:
\begin{align}
\nonumber
[LHS]_j
= &
f(p^{(j)})
\left(
c_2^j h_1^j
-
c_1^j h_2^j
\right)
+ \\ &
\sum_{ \substack{p^{(l)} \in \mathcal{P}(y_1) \cup \mathcal{P}(y_2) \\  p^{(j)} \neq p^{(l)} } }
f(p^{(l)})
\left(
c_2^j c_1^l - c_1^j c_2^l
\right)
\frac
{\vartheta[\widetilde{\zeta}](p^{(l)}-p^{(j)})}
{\vartheta[\widetilde{\zeta}] (0)E(p^{(l)},p^{(j)})}
\label{step2_VV}
.
\end{align}
By extracting the coefficients of $f(p^{(j)})$ in \eqref{step2_VV},
two cases should be distinguished.
The first case is the diagonal entries.
For the $j$-th entry it is the coefficient of $f(p^{(j)})$,
and are in both cases equal to
\begin{equation}
\label{step2_3}
c_{2}^{j} h_{1}^{j} - c_{1}^{j} h_{2}^{j}
.
\end{equation}

Otherwise, the coefficients of $f(p^{(j)})$ are just
\begin{equation}
\label{step2_4}
\bracketsA{c_{2}^{j} c_{1}^{l} - c_{1}^{j} c_{2}^{l} }
\frac
{\vartheta[\widetilde{\zeta}]( p^{(l)}- p^{(j)})}
{\vartheta[\widetilde{\zeta}] (0)E(p^{(l)},p^{(j)})}
.
\end{equation}
Combining \eqref{step2_3} and \eqref{step2_4} and furthermore recalling that \eqref{step2_3} and \eqref{step2_4}
are the entries of $\widetilde{\gamma}$ up to permutation of conjugate poles,
the definition of $\widetilde{\gamma}$ as given in \eqref{tildeGamma} follows.

\begin{Step}
We set $\gamma$ by
\begin{equation}
\label{input}
\gamma_{j,k} =
\Psi_{j,k}+
\begin{cases}
c_{2}^{k} h_{1}^{k} - c_{1}^{k} h_{2}^{k}
&
\tauBa{p^{(j)}} = p^{(k)}
\\
\bracketsA{\overline{c_{2}^{j}} c_{1}^{k} - \overline{c_{1}^{j}} c_{2}^{k} }
\frac
{\vartheta[\widetilde{\zeta}]( p^{(k)}- \tauBa{p^{(j)}})}
{\vartheta[\widetilde{\zeta}] (0)E(p^{(k)},\tauBa{p^{(j)}})}
,
&
{\text otherwise},
\end{cases}
\end{equation}
where $\Psi_{j,k}= (\overline{c_1^j}c_2^k- \overline{c_2^j}c_1^k) K _{\mathcal{X}}(p^{(k)},\tauBa{p^{(j)}})$.
Then the output vessel condition \eqref{VesselCond3} holds.
\end{Step}
Since $\mathcal{X}$ is a reproducing kernel space, it follows that
\begin{equation}
\label{kernelEval}
\Phi\Phi^* = \left[ K _ {\mathcal{X}}(p^{(i)},p^{(j)}) \right]_{i,j=1,...,n}.
\end{equation}
Then, Equation \eqref{input} is derived by substituting
\eqref{sigma123}, \eqref{tildeGamma} and \eqref{kernelEval} in the linkage condition \eqref{VesselCond3}.
\begin{Step}
\label{step4}
For any $h\in\mathcal{X}$,
for any $v \in \Omega$ such that $\pi(v)$ does not belong to $\Spec{M^{y_1},M^{y_2}}$,
the mapping
\begin{equation}
\label{ModelMap}
h
\rightarrow
\widetilde{\bf u}^{\times}_l(v)
(\xi_1 \sigma_1 + \xi_2 \sigma_2)
\Phi(\xi_1A_1+\xi_2A_2 - \xi_1y_1(v) -\xi_2y_2(v))^{-1}h,
\end{equation}
is the identity,
where $\pi$ is bi-rational embedding of $X$ into a curve $C=\overline{C_0}\subset \mathbb P^2$.
\end{Step}
We recall that the mapping \eqref{ModelMap} is independent of the choice of the direction $(\xi_1,\xi_2)$
(see \cite[Section 3]{MR1634421} and Section \ref{secModelSpace} above for more details).
Hence, for the sake of simplicity, we illustrate the calculation for the direction $\xi = (1 \ 0)$.
A direct computation leads to:
\begin{align}
\nonumber
h \longrightarrow
&
\widetilde{\bf u}^{\times}_l(v)
(1 \cdot \sigma_1 + 0 \cdot \sigma_2)
\Phi(M^{y_1} - y_1(v) )^{-1}h
\\
= &
\nonumber
\widetilde{\bf u}^{\times}_l(v) \sigma_1
\Phi R^{y_1}_{y_1(z)}(h)
\\
\nonumber
= &
\widetilde{\bf u}^{\times}_l(v) \sigma_1
\Phi
\left(
\frac{h(\cdot)}{y_1(\cdot) - y_1(v)} -
\sum_{j=1}^{n}
{
\frac{h(u^{(j)})}{dy_1(u^{(j)})}
\frac
{\vartheta[\widetilde{\zeta}](u^{(j)}-\cdot)}
{\vartheta[\widetilde{\zeta}] (0)E(u^{(j)},\cdot)}
}
\right)
\end{align}
and so
\begin{align}
\nonumber
h
\longrightarrow
&
\widetilde{\bf u}^{\times}_l(v) \sigma_1
\left[
\restr{
\frac{h(u)}{y_1(v) - y_1(v)}
-
\sum_{j=1}^{n}
{
\frac{h(u^{(j)})}{dy_1(u^{(j)})}
\frac
{\vartheta[\widetilde{\zeta}](u^{(j)}-u)}
{\vartheta[\widetilde{\zeta}] (0)E(u^{(j)},u)}
}
}{u=p^{(l)}}
\right]_{l=1}^n
\\
= &
-
\widetilde{\bf u}^{\times}_l(v) \sigma_1
\left[
\sum_{j=1}^{n}
{
\frac{h(u^{(j)})}{dy_1(u^{(j)})}
\frac
{\vartheta[\widetilde{\zeta}](u^{(j)}-p^{(l)})}
{\vartheta[\widetilde{\zeta}] (0)E(u^{(j)},p^{(l)})}
}
\right]_{l=1}^n
\label{eqMapToMS}
,
\end{align}
where $ \left(  u^{(j)}  \right)_{j=1} ^ {n}$ are the points in $\Omega$ such that
$y_1(u^{(j)}) = y_1(v)$.
In the line bundle case,
the matrix of normalized section at $u$,
is given by a vector of Cauchy kernels at $u$ and $\pi^{-1}({\infty})$,
namely, it is of form
(see also \cite{MR1704479})
\begin{equation}
\label{normSectionCK}
[\widetilde{\bf u}^{\times}_l(v)]_j
=
K_{\widetilde{\zeta}}(v,p^{(j)})
=
\frac
{\vartheta[\widetilde{\zeta}](\tauBa{p^{(j)}}-v)}
{i \vartheta[\widetilde{\zeta}] (0)E(v,\tauBa{p^{(j)}})}
.
\end{equation}
Now, using the $2 \times 2$ diagonal-block structure of $\sigma_1$ (as in \eqref{sigma123}),
substituting \eqref{normSectionCK} and changing the summation order, we have:
\begin{align}
\nonumber
h \longrightarrow
&
-
\left[
\frac
{\vartheta[\widetilde{\zeta}](\tauBa{p^{(m)}} - v)}
{i\vartheta[\widetilde{\zeta}] (0)E(v,\tauBa{p^{(m)}})}
\right]^T_{m=1,...,n}
\sigma_1
\times
\\ &
\nonumber
\left[
\sum_{j=1}^{n}
{
\frac{h(u^{(j)})}{dy_1(u^{(j)})}
\frac
{\vartheta[\widetilde{\zeta}](u^{(j)}-p^{(l)})}
{\vartheta[\widetilde{\zeta}] (0)E(u^{(j)},p^{(l)})}
}
\right]_{l=1,...,n}
\\ = &
-
\sum_{j=1}^{n}
{
\frac{h(u^{(j)})}{dy_1(u^{(j)})}
\sum_{l=1}^{n}
{
\frac
{\vartheta[\widetilde{\zeta}](\tauBa{p^{(l)}} - v)}
{i\vartheta[\widetilde{\zeta}] (0)E(v,\tauBa{p^{(l)}})}
\sigma_1^{\tauBa{p^{(l)}}}
}
\frac
{\vartheta[\widetilde{\zeta}](u^{(j)}-\tauBa{p^{(l)}})}
{\vartheta[\widetilde{\zeta}] (0)E(u^{(j)},\tauBa{p^{(l)}})}
}.
\label{eq12344}
\end{align}
We use a version of the collection formula as presented in \eqref{CollFormLimitVer} to simplify the second summation in \eqref{eq12344},
which is equal to $-dy_1(v)$ whenever $u^{(j)} = v$ and zero otherwise.
Hence, using $\sigma_1$ as given in \eqref{sigma123}, one may conclude that
\begin{align*}
h \longrightarrow
&
-
\sum_{j=1}^{n}
h(u^{(j)})
(- \delta_{u^{(j)},v} )
=
h(v),
\end{align*}
where $\delta_{u^{(j)},v}$ denotes the Kronecker delta, that is,
$\delta_{u^{(j)},v} = 1$ if and only if $v = u^{(j)}$.
\begin{Step}
\label{step6A}
The reproducing kernel of $\mathcal{X}$,
in terms of the joint characteristic function $S$ of the vessel \eqref{vessel123}, is given by
\begin{align}
\label{eqRkS}
K_{\mathcal X}(p,q)
= &
\frac{\widetilde{\bf u}^{\times} _l (p)(\xi\sigma)\widetilde{\bf u}^{\times} _l (q)^* }
{-i  (\xi y(p) - \xi \overline{y(q)}) }
-
\frac{\widetilde{\bf u}^{\times} _l (p) S(\xi,p) ( \xi \sigma ) S(\xi,q)^*\widetilde{\bf u}^{\times} _l (q)^*}
{-i (\xi  y(p) - \xi  \overline{y(q)}) }
\end{align}
where $\pi(p)$ and $\pi(q)$ lie outside $\Spec{M^{y_1},M^{y_2}}$.
\end{Step}
$\mathcal X$, by assumption, is a reproducing kernel Hilbert space and in particular
\begin{equation}
\label{qwe123A}
h(z)
=
\innerProductReg{h}{K_{\mathcal X}(\cdot,z)}_{\mathcal X}.
\end{equation}
On the other hand, the mapping to the model space is, by Step \ref{step4}, the identity.
Thus, we have
\begin{align}
\nonumber
\hat{h}(z)
= &
\widetilde{\bf u}^{\times} _l (p)
(\xi_1 \sigma_1 + \xi_2 \sigma_2)
\Phi
(\xi_1A_1+\xi_2A_2 - \xi_1 y_1(z) -\xi_2y_2(z))^{-1}h
\\ = &
\label{qwe123B}
h(z).
\end{align}
Combining equations \eqref{qwe123A} and \eqref{qwe123B}, one may conclude that the reproducing kernel
can be expressed explicitly in terms of the model space mapping by:
\begin{align}
\label{eqsec5B}
K_{\mathcal X}(p,q) =
&
\widetilde{\bf u}^{\times} _l (p) (\xi \sigma)
\Phi(\xi A - \xi y(p))^{-1}\times
\\
&
\nonumber
\hspace{21mm}
(\xi A - \xi y(q))^{-*}\Phi^*
(\xi \sigma)^*
\widetilde{\bf u}^{\times} _l (p)^*
.
\end{align}
A classical computation in the single-operator colligation setting
(a similar computation can be found in \cite[Chapter 10]{KLMV})
yields the following relation:
\begin{align}
\nonumber
(\xi\sigma)
\Phi
(\xi A - \xi y(p))^{-1}
&
(\xi A - \xi y(q))^{-*}
\Phi^*
(\xi \sigma)^*
\\ & =
\frac
{( \xi \sigma ) - S(\xi,p) ( \xi \sigma ) S(\xi, q)^*  }
{-i(\xi y(p) - \xi \overline{y(q)})},
\label{eqsec5A}
\end{align}
where $S(\xi,z)$ is given by
\begin{equation*}
S(\xi,z)
=
I - i ( \xi \sigma ) \Phi(\xi A - \xi y(z))^{-1}\Phi^* .
\end{equation*}
Finally, we substitute \eqref{eqsec5A} in \eqref{eqsec5B} to conclude and get the expression in \eqref{eqRkS}.
\begin{Step}
\label{stepMax}
The input and the output determinantal representations of the vessel \eqref{vessel123} are maximal.
\end{Step}

For the definition of maximality, we refer to Section \ref{sub22} above, see also \cite{MR97m:30051,MR1634421}.
We start by showing the maximality of the output determinantal representation. \smallskip

The canonical determinantal representation constructed out from
a line bundle and a pair of meromorphic functions that
gives a birational embedding of the compact Riemann surface in
${\mathbb P}^2$ using the Cauchy kernels,
is always maximal.
The proof is given in \cite[Theorem 5.1]{MR1704479} under the assumption that the singularities of $C$ are ordinary multiple points (i.e. nodes).
\smallskip

We move to prove maximality of the input determinantal representation,
but first we present several preliminary results.
\smallskip

We recall that there are several definitions for the notion of the joint spectrum of a pair of commuting operator.
Here we mention the definition of Harte \cite{harte1972spectral}, Taylor \cite{MR42:3603}, and the one given in \cite{MR1634421}.
These various definitions coincide when the operators have compact imaginary parts, and in particular in our case (see \cite{MR1320541}).

\begin{lem}
The joint spectrum of $M^{y_1}$ and $M^{y_2}$ lies in $C_0$.
\end{lem}
\begin{pf}
Let us assume that $(\lambda_1,\lambda_2)\in \Spec{M^{y_1},M^{y_2}}$.
Then, using the spectral mapping theorem twice, we have
\[
P(\lambda_1,\lambda_2) \in P(\Spec{M^{y_1},M^{y_2}}) = \Spec{P(M^{y_1},M^{y_2})},
\]
where $P$ is the discriminant polynomial of the vessel.
Furthermore, using Theorem \ref{rkModelOpAlg} and the Cayley-Hamilton Theorem \cite[Theorem 4.1.2]{KLMV},
we may conclude the following
\[
M^{P(y_1,y_2)} = P(M^{y_1},M^{y_2}) = 0.
\]
Thus, $P(\lambda_1,\lambda_2)=0$ and $(\lambda_1,\lambda_2)\in C_0$ follows.
\end{pf}
\begin{lem}
\label{gExists}
Let $X$ be a compact Riemann surface of genus $g$ which is the normalization of
a curve $C \subset \mathbb P ^2$ with the bi-rational embedding $\pi: X \rightarrow \mathbb P^2$.
Let $\Omega$ be an open subset containing the pre-images of the singular points of $C$ and let $p \in \Omega$.
Then there exists a polynomial $h \in \mathbb C [z_1,z_2]$ such that $h(\pi(p)) = 0$ and all the zeros of $h \circ \pi$ on $X$ are in $\Omega$.
\end{lem}
\begin{pf}
We denote by $\mu_k$ the Abel-Jacobi mapping from $X^{(k)}$ to $J(X)$,
sending an effective divisors of degree $k$ to the Jacobian
and we recall that the curve $C$ is of degree $m$.
Let us assume, without loss of generality, that $p_0$, the base point of the Abel-Jacobi map, belongs to $\Omega$ and,
furthermore, we choose $p_0$ which does not belong to the finite set of the Weierstrass points.
\smallskip

Let $k>0$ be the degree of a polynomial $h \in \mathbb C [z_1,z_2]$ and
let $D$ be the effective divisor of $h$.
We know that $D$ is linear equivalent to $k \cdot L$, where $L$ is a divisor of straight line, that is,
\begin{equation}
\label{gDivisor}
D \equiv \deg (h) \, \cdot \, L.
\end{equation}
Conversely, let $D$ be an effective divisor of degree $k \, m$ containing the divisor of singularities, denoted by $D_{\rm sing}$.
Then, for sufficiently large $k$, $D$ is the divisor of some homogeneous polynomial of degree $k$
(by the completeness of the linear system of the adjoint curves, see e.g. \cite{fulton}).
By assumption, $\Omega$ contains the pre-images of the singularities of $C$ and
the existence of homogeneous polynomial of degree $k$ follows.
\smallskip

Let $D$ be a divisor that contains $p$ and the singular points, that is $D \geq p +D_{\rm sing}$.
We set an effective divisor $D^\prime$, supported on $\Omega$, by
\[
D^\prime = D - p - D_{\rm sing}.
\]
Recall that by the Abel-Jacobi theorem, \eqref{gDivisor} holds if and only if
\begin{equation}
\label{gDivisor12}
\mu_{k\,m}(D) = \deg (h) \, \cdot \, \mu_{m} (L),
\end{equation}
and hence, using \eqref{gDivisor12}, we have
\begin{align*}
\mu_{Deg(D^\prime)}(D^\prime)
=
k \, \mu_m(L) - \mu_{(1+deg(D)-g)}(p + D_{\rm sing}).
\end{align*}

We turn to consider the divisor $g \cdot p_0$.
This divisor, since $p_0$ is the base point of the Abel-Jacobi mapping, is zero under the mapping $\mu_{g}(\cdot)$.
Furthermore, since $p_0$ is not a Weierstrass point, $\mu_g(\cdot)$ is invertible at $g \cdot p_0$.
Therefore, $\mu_g(\Omega^{(g)})$ contains an open ball around zero
(where $\Omega^{(\ell)}$ stands for the set of effective divisors of degree $\ell$ with support in $\Omega$).
\smallskip

For any $j>0$ we have
$ j \cdot \mu_g(\Omega^{(g)}) \subseteq  \mu_{j \cdot g}(\Omega^{(jg)}) $
and hence for a $j_0$ sufficiently large, the equality $\mu_{j_0 \cdot g}(\Omega^{(j_0 g)}) = J(X)$ holds.
As a consequence, $\mu_{l}(\Omega^{(l)}) = J(X)$ for $l \geq j_0 g$.
It remains to choose $k$ such that
$k m - 1 - \deg \left(D_{\rm sing}\right) \geq j_0 g$.
\end{pf}

The following result appeared in \cite[Section 2]{MR97m:30051}) without a proof.

\begin{lem}
\label{lemMaxIO}
Let us assume that the singular points of $C$ lie outside the spectrum of $\left( A_1 , A_2 \right)$.
Then the following statements hold:
\begin{enumerate}
\item If one of the determinantal representations is maximal then so does the other.
\item If one of the determinantal representations is fully saturated then so does the other.
\end{enumerate}
\end{lem}
\begin{pf}
Assume that the output determinantal representation is maximal and all $\lambda \in C$ lie outside the joint spectrum of $\left( A_1 , A_2 \right)$.
Then, $S(\lambda)$, for any $\lambda$ outside the spectrum of $\left( A_1 , A_2 \right)$,
is a linear mapping from $\mathcal E(\lambda)$ to $\widetilde{\mathcal E}(\lambda)$.
Similarly, the JCF of the adjoint vessel $\widetilde{S}(\lambda)$ (see \cite{MR2043236}),
for any $\lambda$ outside the spectrum of $\left( A_1 , A_2 \right)$,
is a linear mapping from $\widetilde{\mathcal E}_l(\lambda)$ to $\mathcal E_l(\lambda)$.
Hence, $\mathcal E(\lambda)$ and $\widetilde{\mathcal E}(\lambda)$ share the same dimension.
\smallskip

Let us now assume that the output determinantal representation is fully saturated (see Subsection \ref{sub22}).
Then, for a singular affine point $\lambda_0 \in C_0$, we choose $\xi_1$ and $\xi_2$ such that
$\xi_1 \lambda^0_1 + \xi_2 \lambda^0_2$ does not belong to spectrum of $\xi_1 A_1 + \xi_2 A_2$
and by continuity it is also true in a neighborhood of $(\lambda^0_1,\lambda^0_2)$.
Then, one may consider the well-known identity (see for instance \cite{KLMV})
\[
E_l (\lambda)
W(\xi_1,\xi_2, {\xi_1\lambda_1 + \xi_2 \lambda_2} )
=
\widetilde{W}(\xi_1,\xi_2, \xi_1\lambda_1 + \xi_2 \lambda_2)
E(\lambda)
\]
and hence we have:
\begin{equation}
\label{eqBLa535}
E_l (\lambda)
=
\widetilde{W}(\xi_1,\xi_2, \xi_1\lambda_1 + \xi_2 \lambda_2 )
E(\lambda)
W(\xi_1,\xi_2, \xi_1 \overline{\lambda_1} + \xi_2 \overline{\lambda_2} )^*.
\end{equation}
Furthermore, since the adjoint operator is multiplicative,
that is, it satisfies ${Adj \,} (AB) = {Adj \,} (A) {Adj \,} (B)$, we have:
\[
\widetilde{V}(\lambda^0)
=
{Adj \,} S_l(\lambda^0) \, V(\lambda^0) \, {Adj \,} S(\overline{\lambda^0})^*.
\]
Then, since all the functions in \eqref{eqBLa535} are analytic in a neighborhood of $\lambda^0$,
the entries of $V(\lambda^0)$ and $\widetilde{V}(\lambda^0)$ share the same order of zeros
and if one is fully saturated, so does the other.
\smallskip

Finally, to understand the behavior at infinity, we note that
either $(1,0)$ or $(0,1)$ lies outside the spectrum.
It remains to repeat the same arguments for the affine points above.
\end{pf}

Using the previous results, we may conclude the main argument of this step.

\begin{pf}[of Step \ref{stepMax}]
Using Lemma \ref{lemMaxIO}, it is enough to show that for any $p_0 \in \Omega$ the point $(y_1(p_0),y_2(p_0))$
lies outside the joint spectrum of $M^{y_1}$ and $M^{y_2}$.
Although that for any $\alpha$ satisfying $y^{-1}(\alpha) \subseteq \Omega$ the operator $(M^y - \alpha I)^{-1}$ is well defined,
this does not mean that $\mathcal X$ is invariant under $(M^y - \alpha I)^{-1}$.
Thus, our approach is different. We associate to each $\alpha$ such that $y^{-1}(\alpha) \subseteq \Omega$,
a pair of commuting operators and build a corresponding vessel and study the mapping to its model space.
The image of $\mathcal X$ under the model mapping is a reproducing kernel Hilbert space of sections
defined on the complement of the joint spectrum of the pair of operators.
Furthermore, the pair of operators have finite non-Hermitian rank, thus their spectrum consists of real points and possibly isolated non-real points.
Hence, an analytic extension to the neighborhood of a point can be always considered as an extension from infinity (which is of bounded point evaluations).
\smallskip

First, in \eqref{eqRkS}, let $\alpha$ be a point in the neighborhood of infinity and we denote the $n$ distinct
pre-images of $\alpha$ under $\pi$ by $q_1,...,q_n$.
As a consequence, $q_1,...,q_n$ belong to $\Omega$ and the columns of ${\bf u}^{\times}(q)$
are linearly independent and form an invertible matrix.
Using Lemma \ref{gExists}, there exists a polynomial $g$ satisfying $g(\pi(p_0))=0$ such that the entire fiber of $p_0$ belongs to $\Omega$.
We define a new meromorphic function $w_1(z) = g(y_1(z),y_2(z))$ and set $w_2(z) = h(y_1(z),y_2(z))$ for some two-variable polynomial $h$
such that $w_1$ and $w_2$ generate $\mathcal M (X)$.
It follows, by construction, that $(y_1(p_0),y_2(p_0))$ lies outside the joint spectrum of $M^{w_1}$ and $M^{w_2}$.
\smallskip

Then, by Theorem \ref{lemCollCodI}, the colligation conditions for any element in the algebra of meromorphic functions generated by $y_1$ and $y_2$,
in particular, it is true for $w_1$ and $w_2$.
Applying Proposition \ref{myGeneral1},
there exist $\sigma_{w_1}$ and $\sigma_{w_2}$ such that $M^{w_1}$ and $M^{w_2}$ satisfy the colligation conditions.
Hence, repeating Steps (\ref{step2}-\ref{step6A}) but now with $M^{w_1}$ and $M^{w_2}$ implies that the collection
\[
( \, M^{w_1} \, ,  \,  M^{w_2} \, ; \, \mathcal X \, , \,  \Phi \, , \, E \, ; \,  \sigma_{w_1} \, , \,  \sigma_{w_2} \, , \, \gamma  \, , \, \widetilde{\gamma} \, ),
\]
is again a commutative two-operator vessel.
\smallskip

Assuming $p_0\in\Omega$ is regular, the fiber containing $p$ has $n$ distinct elements (assuming all belong to $\Omega$).
The characteristic function of the vessel can be extended to a neighborhood of $y(p_0)$.
Let us assume that
$(y_1(p_0),y_2(p_0))$ belongs to $\Spec{M^{y_1},M^{y_2}}$.
Then, by the spectral mapping theorem for a pair of commuting operators \cite{MR1320541}, the following holds:
\begin{align*}
(g(y_1(p_0),y_2(p_0)),h(y_1(p_0),y_2(p_0)))  \in& \Spec{g(M^{y_1},M^{y_2}),h(M^{y_1},M^{y_2})} \\ =  & \Spec{M^{w_1},M^{w_2}}.
\end{align*}
However, since $w_1(p_0) =  g(y_1(p_0),y_2(p_0))=0$, $(y_1(p_0),y_2(p_0))$
lies outside the joint spectrum of $M^{w_1}$ and $M^{w_2}$, a contradiction.
Hence $(y_1(p_0),y_2(p_0))$ does not belong to the joint spectrum of $M^{y_1}$ and $M^{y_2}$.
\smallskip

We turn now to the case where $p_0 \in \Omega$ is singular and where $\Omega$ contains the entire fiber of $p_0$.
The colligation corresponding to $M^{w_1}$ has finite non-Hermitian rank and hence the characteristic function has only isolated singularities.
We then set $\alpha_0 = w_1(p_0)$ and we show that $\alpha_0$ does not belong to the spectrum of $M^{w_1}$.
Note that using previous result, all elements in a punctured neighborhood of $\alpha_0$ lie outside the spectrum.
We use the following well-known fact (see for instance \cite{brod}): assume the characteristic function of a colligation can be extended analytically to a (full) neighborhood of $\alpha_0$,
then $\alpha_0$ lies outside the spectrum of the colligation.
We consider the reproducing kernel formula as given in \eqref{eqRkS} where
the matrix $[{\bf u}^{\times}_l(p_0)]$, containing the (line) normalized sections of the per-images of $y(p_0)$
(similarly, we define $[{\bf u}^{\times}(p_0)]$).
This matrix is invertible and it is given by
\[
[{\bf u}^{\times}_l(p_0)]^{-1}
=
\frac{\sigma [{\bf u}^{\times}(p_0)] }{dy(p_0)}.
\]
We choose a local parameter $t$ at $p_0$ centered at zero and then $y(t)=t^r$ where $r$ is the order of the singularity at $p_0$.
We then consider the $r$ points $t, \varepsilon t,..., \varepsilon^{r-1} t$, where $\varepsilon$ is the $r$-th root of unity.
All elements in $[{\bf u}^{\times}(p_0)]$ are analytic, i.e. we may use $\phi(t)=\sum_{k=0}^{\infty} \phi_k  t^k$, and then we have
\begin{align*}
\sum_{j=0}^{r-1} \frac{\phi(\varepsilon^j t)}{y'(\varepsilon^j t)}
= &
\sum_{j=0}^{r-1} \frac{\sum_{k=0}^{\infty} \phi_k (\varepsilon^j t)^k}{r(\varepsilon^j t)^{r-1}}
=
\frac{1}{r t^{r-1}} \sum_{k=0}^{\infty} \phi_k t^k \sum_{j=0}^{r-1}  (\varepsilon^j)^{-(r-1)}  (\varepsilon^j)^k t^k
\\ = &
\frac{1}{r} \sum_{k=0}^{\infty} \phi_k t^{k-r+1} \sum_{j=0}^{r-1}  (\varepsilon^{k+1})^j
.
\end{align*}
The inner summation vanishes as long as $\varepsilon^{k+1} \neq 1$,
as it is the sum of all $r$ unit roots.
In particular, the inner summation vanishes whenever $k < r-1$
and thus the negative-index coefficients are zero.
It follows by \eqref{eqRkS} that the characteristic function can be extended to $\alpha_0$ and hence to a (full) neighborhood of $\alpha_0$.
Therefore, $\alpha_0$ does not belong to the joint spectrum of $M^{w_1}$ and $M^{w_2}$
and hence $\alpha_0$ lies outside the joint spectrum of $M^{y_1}$ and $M^{y_2}$.
\end{pf}
\begin{Step}
The reproducing kernel of $\mathcal{X}$ is equal to
\begin{align*}
K_{\mathcal{X}}(z,q)
= &
K_{\widetilde{\zeta}}(z,p) -
T(z)
K_{\zeta}(z,p)
T(q)^*,
\end{align*}
where $T$ is the normalized joint characteristic function
associated to the model vessel \eqref{vessel123}.
\end{Step}
Using Step \ref{stepMax}, the input and output determinantal representations are maximal and hence
we turn to examine the NJCF (the normalized joint characteristic function, which is related to the JCF by \eqref{st}).
Before proceeding, we note that the matrices of normalized sections associated to $\widetilde{\mathcal E} (p)$ and $\mathcal E (p)$,
i.e. $\widetilde{\bf u}^{\times} (p)$ and $ \widetilde{\bf u}^{\times} _l (p)$, satisfy the identities:
\begin{equation}
\label{eqC}
\widetilde{\bf u}_{l}^\times (\tauBa{p})  = \widetilde{\bf u}^{\times} (p)^*,
\qquad
{\rm and}
\qquad
{\bf u}_{l}^\times (\tauBa{p})  = {\bf u}^{\times} (p)^*.
\end{equation}
Using \eqref{eqC}, $T(p) T (\tauBa{p})^{*} = 1$ and $S(\xi,\tauBa{z})^{*} ( \xi \sigma ) S(\xi,z) =  ( \xi \sigma) $
(see \cite[Equation 2-33]{MR1634421} and \cite[Equation 10.84]{KLMV}, respectively),
\eqref{st} may be translated to
\begin{equation}
\label{st2}
T(\tauBa{p}) {\bf u}_{l}^\times(p) = \widetilde{\bf u}_{l}^\times (p)  S(p).
\end{equation}
Note that the output representation has canonical form ($\widetilde{\bf u}^{\times}(p)={\bf u}_{\widetilde{\zeta}}^{\times}(p)$),
while the input is not necessarily canonical.
However, it is equivalent to the canonical form (see \cite{vinnikov1} and \cite[Section 6]{MR1704479}),
that is, there exists $\rho \in GL(M^{mr},\mathbb C)$ satisfying
$\rho \left( \xi \sigma \right) \rho^* = \left( \xi \sigma \right)$ such that
\begin{equation}
\label{eqRho}
{\bf u}^{\times} _l (z) = {\bf u}^{\times} _{l,\zeta} (z) \rho^* .
\end{equation}

Substituting (see \cite[Equation N2]{MR97m:30051}),
\begin{equation}
\label{EqN2}
K_{\widetilde{\zeta}} (p,q)
=
\widetilde{\bf u}^{\times}_l(p)
\frac{\xi_1\sigma_1 + \xi_2 \sigma_2}
{\xi_1 \left(y_1(p)-y_1(q) \right)  + \xi_2 \left(y_2(p)-y_2(q) \right)}
\widetilde{\bf u}^{\times}(q),
\end{equation}
\eqref{st2} and \eqref{eqRho} in \eqref{eqRkS},
we proceed to the following final calculation:
\begin{align*}
K_{\mathcal X}(p,q)
= &
K_{\widetilde{\zeta}}(p,q)
-
\frac
{T(p) {\bf u}^{\times} _l (p) ( \xi \sigma ) {\bf u}^{\times} _l (q)^*  T(q)^* }
{-i  ( \xi y(p) - \xi \overline{y(q)}) }
\\ = &
K_{\widetilde{\zeta}}(p,q)
-
T(p)
\frac
{ {\bf u}^{\times} _{l} (p) (\rho \xi \sigma \rho^* ) {\bf u}^{\times} _{l} (q)^* }
{-i( \xi y(p) - \xi \overline{ y(q)})}
T(q)^*
\\
= &
K_{\widetilde{\zeta}}(p,q)
-
T(p)
\frac
{ {\bf u}^{\times} _{l,\zeta} (p) ( \xi \sigma ) {\bf u}^{\times} _{l,\zeta} (q)^* }
{-i( \xi y(p) - \xi \overline{ y(q)})}
T(q)^*
\\ = &
K_{\widetilde{\zeta}}(p,q)
-
T(p)
K_{\zeta}(p,q)
T(q)^*
.
\end{align*}
\end{pf}

We now turn to present the proof of Theorem \ref{MainTh}.\\

\begin{pf}[of Theorem \ref{MainTh}]
Let $y_1(\cdot)$ and $y_2(\cdot)$ be two meromorphic functions,
not necessarily with simple poles, generating $\mathcal M (X)$.
Hence, there exists a bi-rational embedding of $X$ into a curve $C$
defined by the compactification to $\mathbb P ^2$ of the curve $C_0$ and is given by
\begin{align*}
\nu \colon X & \rightarrow  C_0 \subseteq {\mathbb C}^2
\\
x & \to (y_1(x),y_2(x)).
\end{align*}
Let us consider the pair of meromorphic functions defined by
\[
\widetilde{y}_k (z) \overset{\text{def} } {=} \frac{1}{y_k(z) - \alpha_k}
\qquad
\alpha_k \in \mathbb R
\qquad {\rm for} \qquad k=1,2.
\]
Since $y_k(z)$ has finite number of poles and
$\nu$ has finite number of critical points,
we can choose $\alpha_1$ and $\alpha_2$ such that the the pre-images
of $\alpha_1$ and $\alpha_2$ are non critical points of $\nu^\prime$,
the functions $\widetilde{y}_1 (z) $ and $\widetilde{y}_2 (z)$ have only simple poles
and
the structure identity holds for $\alpha_k$ and $\widetilde{y}_k (z) $ for $k=1,2$.
\smallskip

The functions $\widetilde{y}_1 (z)$ and $\widetilde{y}_2 (z)$ define an additional
embedding of $X$ into $\mathbb P ^2$, but to a different curve, denoted by $C^{\prime}$,
defined by the compactification of $C^{\prime}_0$:
\begin{align*}
\nu^{\prime} \colon X & \rightarrow C^{\prime}_0 \subseteq {\mathbb C}^2
\\
x & \to ((y_1(x)-\alpha_1)^{-1},(y_2(x)-\alpha_2)^{-1})
.
\end{align*}
The new critical points of $\nu^{\prime}$ (except the critical points inherited from $\nu$)
can be only the points above the poles of
$(y_1(x)-\alpha_1)^{-1}$ and $(y_2(x)-\alpha_2)^{-1}$.
However, by assumption, the fibers above $\alpha_1$ and $\alpha_2$ are regular.
Thus
$(\nu^{\prime})^{-1} C'_{\text{\rm sing}} \subseteq (\nu)^{-1} C_{\text{\rm sing}} $ and the regularity assumption of $(y_1(z),y_2(z))$
implies regularity of $(\widetilde{y}_1(z),\widetilde{y}_2(z))$.
\smallskip

We now can apply Theorem \ref{preTh}. Then, we embed the operators $M^{\widetilde{y}_1}$ and $M^{\widetilde{y}_2}$
into a commutative vessel of the form
\begin{equation*}
\mathcal V^{\prime} =
( \, M^{(y_1-\alpha_1)^{-1}} \, , \, M^{(y_2-\alpha_2)^{-1}} \, ;  \, \mathcal X \, , \,  \Phi \, , \, E \, ; \, \sigma_1 \, , \, \sigma_2 \, , \, \gamma, \, \widetilde{\gamma} \, ).
\end{equation*}
The reproducing kernel is of the form
\begin{equation*}
K(p,q)
=
K_{\widetilde{\zeta}}(p,q) -
T^{\prime}(p)
K_{\zeta}(p,q)
T^{\prime}(q)^*
\end{equation*}
where $T^{\prime}$ is the normalized joint characteristic function of the vessel $\mathcal V^{\prime}$.
We note that since $M^{\frac{1}{y(u)-\alpha}} = R^y_\alpha$,
the colligation conditions of $\mathcal V^{\prime}$ are derived from the structure identities of $\mathcal V$.
This completes the proof of part \ref{stateA}.
\smallskip

We now turn to prove part \ref{stateB}.
Let $T$ be a $(\zeta,\widetilde{\zeta})$-contractive mapping and let $y_1$ be a real meromorphic function.
We choose $y_2$, another real meromorphic function, such that all the poles of $y_2$ are contained in $\Omega$ and $y_1$ and $y_2$ generate $\mathcal M (X)$.
\smallskip

Then, one may associate to $\widetilde{\zeta}$ and $\zeta$ the canonical determinantal representations.
We multiply the canonical determinantal representations associated to $\zeta$ by the values of $T$ at infinity.
Thus, we keep maximality of the determinantal representation and force $T$ to be the identity at infinity.
Combining the last observations, we can apply the realization theorem for the NJCF (i.e. Theorem \ref{realizationTh}).
Thus, there exists a two-operator commutative vessel
\[
\mathcal V = ( \, A_1 \, , \, A_2  \, ; \,  H  \, , \,  \Phi \, , \, \mathbb C ^n \, ; \,  \sigma_1 \, , \, \sigma_2 \, , \, \gamma \, , \,  \widetilde{\gamma} \, )
\]
where $T$ is its normalized joint characteristic function (and where $H$ is a Hilbert space, $A_1$ and $A_2$ are bounded operators on $H$).
The associated model vessel is given by
\[
\mathcal V_T
= ( \, M^{y_1} \, , \, M^{y_2}  \, ; \,  \mathcal H(T) \, , \, \Phi_T \, , \, \mathbb C ^n \, ; \,  \sigma_1 \, , \, \sigma_2 \, , \, \gamma \, , \,  \widetilde{\gamma} \, )
,
\]
where $\mathcal H(T)$ is the reproducing kernel Hilbert space with the reproducing kernel \eqref{rkhsKernel} and $\Phi_T$ is the evaluation operator at the poles of $y_1$ and $y_2$.
Then $\mathcal V_T$ is an irreducible commutative two-operator vessel which is unitary equivalent,
on its principal subspace, to $\mathcal V$ and $T$ is its NJCF, see Theorem \ref{vesselAreEq}.
Furthermore, $\mathcal X$ is invariant under $M^{y_1}$ and then, by Lemma \ref{mainThLem},
is also $R_\alpha^{y_1}$-invariant for $\alpha$ in the neighborhood of infinity.
\smallskip

To complete the proof, we mention that by Proposition \ref{myGeneral1},
the structure identity is equivalent to the colligation condition.
Hence the colligation condition for $M^{y_1}$ in $\mathcal V_T$ implies the structure identity for $y_1(\cdot)$.
\end{pf}
\def\cfgrv#1{\ifmmode\setbox7\hbox{$\accent"5E#1$}\else
  \setbox7\hbox{\accent"5E#1}\penalty 10000\relax\fi\raise 1\ht7
  \hbox{\lower1.05ex\hbox to 1\wd7{\hss\accent"12\hss}}\penalty 10000
  \hskip-1\wd7\penalty 10000\box7} \def$'${$'$} \def$'${$'$}
  \def$'${$'$} \def\lfhook#1{\setbox0=\hbox{#1}{\ooalign{\hidewidth
  \lower1.5ex\hbox{'}\hidewidth\crcr\unhbox0}}} \def$'${$'$}
  \def$'${$'$} \def$'${$'$} \def$'${$'$} \def$'${$'$}
  \def$'${$'$}

redpajama_set_name: RedPajamaArXiv