\item Linear systems for multiple orthogonal polynomials and their duals
\begin{align}\label{flow.mop}
\begin{aligned}
\frac{\partial \A_{a'}}{\partial t_{j,a}}&=(B_{j,a}-\delta_{a,a'}x^j)\A_{a'},&  \frac{\partial \A_{a'}}{\partial
\bar t_{j,b}}&=(\bar B_{j,b})\A_{a'},&
\frac{\partial \bar \A_{b'}}{\partial t_{j,a}}&=-B_{j,a}^\top\bar\A_{b'},&  \frac{\partial \bar\A_{b'}}{\partial \bar
\end{aligned}
\end{align}
\item Lax equations
\begin{align}\label{laxeq}
\frac{\partial L_{a'}}{\partial t_{j,a}}&=[B_{j,a},L_{a'}], &\frac{\partial L_{a'}}{\partial \bar
t_{j,b}}&=[\bar B_{j,b},L_{a'}], &
\frac{\partial \bar L_{b'}}{\partial t_{j,a}}&=[B_{j,a}, \bar L_{b'}], &\frac{\partial \bar L_{b'}}{\partial
\bar t_{j,b}}&=[\bar B_{j,b}, \bar L_{b'}].
\end{align}
\item Zakharov--Shabat equations
\begin{align}\label{zseq}
\frac{\partial B_{j,a}}{\partial t_{j',a'}}-
\frac{\partial B_{j',a'}}{\partial t_{j,a}} +[B_{j,a},B_{j',a'}]&=0,\\
\frac{\partial \bar B_{j,b}}{\partial \bar t_{j',b'}}-
\frac{\partial \bar B_{j',b'}}{\partial \bar t_{j,b}} +[\bar B_{j,b},\bar B_{j',b'}]&=0,\\
\frac{\partial B_{j,a}}{\partial \bar t_{j',b'}}-
\frac{\partial \bar B_{j',b'}}{\partial t_{j,a}} +[B_{j,a},\bar B_{j',b'}]&=0.
\end{align}
\end{enumerate}
\end{theorem}
\begin{proof}
To prove \eqref{auxlinsys} we proceed as follows. In the first place we compute
\begin{align*}
\frac{\partial W_0}{\partial t_{j,a}}&=\Lambda_{1,a}^{j}W_0,  &\frac{\partial \bar W_0}{\partial \bar
t_{j,b}}&=(\Lambda_{2,b}^{\top})^{j}\bar W_0,
\end{align*}
and in the second place we observe that
\begin{align}
\frac{\partial W}{\partial t_{j,a}} &=\Big(\frac{\partial S}{\partial t_{j,a}}S^{-1}+L_a^{j}\Big)W,&
\frac{\partial \bar W}{\partial t_{j,a}} &=\Big(\frac{\partial \bar S}{\partial t_{j,a}}S^{-1}\Big)\bar W, \label{desfac1}\\
\frac{\partial  W}{\partial \bar t_{j,b}}&=\Big(\frac{\partial  S}{\partial \bar t_{j,b}}  S^{-1}\Big) W,&
\frac{\partial \bar W}{\partial \bar t_{j,b}}&=\Big(\frac{\partial \bar S}{\partial \bar t_{j,b}} \bar S^{-1}+
\bar L_b^{j}\Big) \bar W.\label{desfac2}
\end{align}
Now, using the factorization problem  we get
\begin{align*}
\frac{\partial S}{\partial t_{j,a}}S^{-1}+L_a^{j}
=\frac{\partial \bar S}{\partial t_{j,a}} \bar S^{-1},
\\
\frac{\partial \bar S}{\partial \bar t_{j,b}} \bar S^{-1}+
\bar L_b^{j}
=\frac{\partial S}{\partial \bar t_{j,b}} S^{-1},
\end{align*}
which, taking the $+$ part (upper triangular) and the $-$ part (strictly lower triangular) imply
\begin{align}\label{desfac3}
\frac{\partial S}{\partial t_{j,a}}S^{-1}&=-(L_a^{j})_-, &\frac{\partial \bar S}{\partial t_{j,a}} \bar
S^{-1}&=(L_a^{j})_+, & \frac{\partial \bar S}{\partial \bar
t_{j,b}} \bar S^{-1}&=-(\bar L_b^{j})_+,
& \frac{\partial S}{\partial \bar t_{j,b}} S^{-1}&=(\bar L_b^{j})_-,
\end{align}
so using \eqref{desfac3} into \eqref{desfac1} and \eqref{desfac2} with the definitions \eqref{defzs} we obtain
\eqref{auxlinsys}.
The linear system \eqref{linbaker} is obtained by inserting  \eqref{defbaker} into \eqref{auxlinsys}.

To obtain the Lax equations \eqref{laxeq} %is convenient to notice that
%  L_a&=W(t)\Lambda_{1,a}W(t)^{-1},&
%   \bar L_b&=\bar W(t)\Lambda_{2,b}^\top\bar W(t)^{-1},
we take derivatives of \eqref{deflax}
\begin{align*}
\frac{\partial L_{a'}}{\partial t_{j,a}}&=\Big[\frac{\partial S}{\partial t_{j,a}}S^{-1},L_{a'}\Big]=[B_{j,a},L_{a'}],&
\frac{\partial \bar L_{b'}}{\partial t_{j,a}}&=\Big[\frac{\partial \bar S}{\partial t_{j,a}}\bar S^{-1},\bar L_{b'}\Big]=[B_{j,a},\bar L_{b'}], \\
\frac{\partial L_{a'}}{\partial \bar t_{j,b}}&=\Big[\frac{\partial S}{\partial \bar t_{j,b}}S^{-1},L_{a'}\Big]=[\bar B_{j,b},L_{a'}],&
\frac{\partial \bar L_{b'}}{\partial \bar t_{j,b}}&=\Big[\frac{\partial \bar S}{\partial \bar t_{j,b}}\bar S^{-1},\bar L_{b'}\Big]=[\bar B_{j,b},\bar
\end{align*}
Finally, \eqref{zseq} are obtained as compatibility conditions for \eqref{auxlinsys}.
\end{proof}
All these equations provide us with different descriptions of  a multi-component integrable hierarchy of the 2D Toda
lattice hierarchy type that rules the
flows of the multiple orthogonal polynomials with respect to  deformed weights. This integrable hierarchy is
the Toda type extension of the
multi-component KP hierarchy considered in [9].

%This hierarchy  is a modification of the one considered in [28, 29].
%The variations are twofold, in the first place and in contrast with the Gauss factorization used  here inthe quoted papers we had a block Gauss factorization, difference that should  lead to modified dependent variables, as it happens with the KP and modified KP hierarchies; in the second place,  the blocks here are rectangular while in  mentioned papers where of square type.
\subsection{Darboux--Miwa discrete  flows}
We complete the previously considered  continuous flows with discrete flows, which we introduce through an iterated
application of  Darboux
transformations [5].

\begin{definition}
Given  sequences of complex numbers
\begin{align}
\begin{aligned}
\lambda_a&:=\{\lambda_a(n)\}_{n\in\Z}\subset\C, & a&=1,\dots,p_1,&
\bar\lambda_b&:= \{\bar \lambda_b(n)\}_{n\in\Z}\subset\C, & b&=1,\dots,p_2,
\end{aligned}
\end{align}
(where $\bar\lambda$ is not intended to denote the complex conjugate of $\lambda$) and two  vectors,
$(s_1,\dots,s_{p_1})\in\Z^{p_1}$ and $(\bar s_1,\dots,\bar s_{p_2})\in\Z^{p_2}$, we construct the following
\begin{align}
\begin{aligned}
D_0&:=\sum_{a=1}^{p_1}D_{0,a}, & D_{0,a}&:=\begin{cases}
\prod_{n=1}^{s_a}(\Lambda_{1,a}-\lambda_a(n)\Pi_{1,a}), &s_a>0,\\
\Pi_{1,a},& s_a=0,\\
\prod_{n=1}^{|s_a|}(\Lambda_{1,a}-\lambda_a(-n)\Pi_{1,a})^{-1},& s_a<0,
\end{cases}\\
%\Big(\delta_{s_a,0}\Pi_{1,a}+(1-\delta_{s_a,0})\prod_{n=1}^{|s_a|}(\Lambda_{1,a}-\lambda_a(n\sgn s_a)\Pi_{1,a})^{\sgn s_a}\Big),\\
\bar   D_0^{-1}&:=   \sum_{b=1}^{p_2} \big(\bar D_{0}^{-1}\big)_b & \big(\bar D_{0}^{-1}\big)_{b}&
\prod_{n=1}^{\bar s_b}(\Lambda_{2,b}^\top-\bar \lambda_b(n)\Pi_{2,b}), &\bar s_b>0,\\
\Pi_{2,b},&\bar s_b=0,\\
\prod_{n=1}^{|\bar s_b|}(\Lambda_{2,b}^\top-\bar \lambda_b(-n)\Pi_{2,b})^{-1},&\bar s_b<0,
\end{cases}
%   \sum_{b=1}^{p_2}\prod_{n=1}^{|\bar s_b|}(\Lambda_{2,b}^\top-\bar\lambda_b(n\sgn s_b)\Pi_{2,b})^{\sgn b},
\end{aligned}
\end{align}
where $s:=\{s_a,\bar s_b\}_{\substack{a=1,\dots,p_1\\b=1,\dots,p_2}}$ denotes the set of discrete times,
%, and as usual $\sgn s=1$ if $s\geq 0$ and $\sgn s=-1$ if $s<0$.
in terms of which we define the  deformed moment matrix
\begin{align}
g(s)=D_0(s)g\bar D_0(s)^{-1}.
\end{align}
\end{definition}

\begin{pro}
The moment matrix $g(s)$ has the same form as the moment matrix $g$ but with new weights
\begin{align}\label{d.evolved}\begin{aligned}
w_{1,a}(s,x)&=\Ds_{a}(x,s_a)w_{1,a}(x),& \Ds_a&:=\begin{cases}
\prod_{n=1}^{s_a}(x-\lambda_a(n)), &s_a>0,\\
1,& s_a=0,\\
\prod_{n=1}^{|s_a|}(x-\lambda_a(-n))^{-1},& s_a<0,
\end{cases}\\
%  (x-\lambda_a(\sgn s_a))^{\sgn s_a}\cdots(x-\lambda_a(s_a))^{\sgn s_a}w_{1,a}(x),\\
w_{2,b}(s,x)&=\bar \Ds_b(x,\bar s_b)^{-1}w_{2,b}(x),& \bar \Ds_b^{-1}&:=\begin{cases}
\prod_{n=1}^{\bar s_b}(x-\bar \lambda_b(n)), &\bar s_b>0,\\
1,&\bar s_b=0,\\
\prod_{n=1}^{|\bar s_b|}(x-\bar \lambda_b(-n))^{-1},&\bar s_b<0,
\end{cases}\end{aligned}
%(x-\bar\lambda_b(\sgn \bar s_b))^{\sgn \bar s_b}\cdots(x-\bar\lambda_b(\bar s_b))^{\sgn \bar s_b}w_{2,b}(x).
\end{align}
\end{pro}

Thus, the proposed discrete evolution introduces new zeroes and poles in the weights at the points  defined by
sequences of
$λ$'s. For example,  in the $a$-th direction, the $s_a$ flow in the positive direction, $s_a→s_a+1$, introduces
a new zero at the point
$λ_a(s_a+1)$, while  if we move in the negative direction, $s_a→s_a-1$, it introduces a simple pole at
$λ_a(s_a-1)$. Let us stress that for the time being we have not ensured the reality and positiveness/negativeness  of the evolved
weights, this will be considered later on.

\subsubsection{Miwa transformations}
Here we show that the discrete flows just introduced can be reproduced with the aid of Miwa shifts in the continuous
\begin{definition}
We consider two types of Miwa transformations:
\begin{enumerate}
\item We introduce the following time
\begin{align}
t \rightarrow t\mp[z^{-1}]_a:=\Big\{t_{j,a'}  \mp\delta_{a',a}\frac{1}{jz^j},\;\bar t_{j,b'}\Big\}_{\substack{j=1,2,\dots,\\a'=1,\dots,p_1,\\b'=1,\dots,p_2}},
\end{align}
\item Dual   time shifts are
\begin{align}
t\rightarrow t\pm\overline{[z^{-1}]}_b:=\Big\{t_{j,a'} ,\; \bar t_{j,b'}
\pm\delta_{b',b}\frac{1}{jz^j}\Big\}_{\substack{j=1,2,\dots,\\a'=1,\dots,p_1,\\b'=1,\dots,p_2}}.
\end{align}
\end{enumerate}
\end{definition}
\begin{pro}
The Miwa transformations produce the following effect on the weights
\begin{align}
\label{vertex weight}
w_{1,a'}(x,t\mp [z^{-1}]_a,s)&=  \Big(1-\frac{x}{z}\Big)^{\pm\delta_{a,a'}}w_{1,a'}(x,t,s),&
w_{2,b'}(x,t\mp [z^{-1}]_a,s)&=  w_{2,b'}(x,t,s),\\
\label{vertex dualweight}
w_{2,b'}(x,t\pm\overline{[z^{-1}]}_b,s)&= \Big(1-\frac{x}{z}\Big)^{\pm\delta_{b,b'}}w_{2,b'}(x,t,s).
\end{align}
\end{pro}
\begin{proof}
When we consider what happens to the evolutionary factors under these shifts we find
\begin{align}
\exp\Big(\sum_{j}t_{j,a'}x^j\Big)\rightarrow
\exp\Big(\sum_{j}\big(t_{j,a'}\mp\delta_{a',a}\frac{x^j}{jz^j}\Big)\Big)=
\Big(1-\frac{x}{z}\Big)^{\mp\delta_{a',a}}\exp\Big(\sum_{j}t_{j,a'}x^j\Big),
\end{align}
and therefore the weights transform according to
\begin{align}\label{vertex weight evol}
w_{1,a'}(x,t,s)\rightarrow   \Big(1-\frac{x}{z}\Big)^{\pm\delta_{a,a'}}w_{1,a'}(x,t,s),
\end{align}
which is like the Darboux transformations considered previously.
For the dual Miwa shifts  we consider what happens to the evolutionary factors under these shifts
\begin{align}
\exp\Big(-\sum_{j,b'}\bar t_{j,b'}x^j\Big)\rightarrow   \exp\Big(-\sum_{j,b'}\Big(\bar
\Big(1-\frac{x}{z}\Big)^{\pm 1}\exp\Big(-\sum_{j,b'}\bar t_{j,b'}x^j\Big),
\end{align}
and the transformation for the weights is

\begin{align}\label{vertex dual weight evol}
w_{2,b'}(x,t,s)\rightarrow   \Big(1-\frac{x}{z}\Big)^{\pm\delta_{b,b'}}w_{2,a'}(x,t,s).
\end{align}
\end{proof}
Thus, a comparison of \eqref{d.evolved}, \eqref{vertex weight} and \eqref{vertex dualweight} leads to
\begin{pro}
Miwa transformations and discrete flows can be identified as follows
\begin{align}
\label{miwa-darboux}\begin{aligned}
c_a w_{1,a}(x,t,s_a)&=
\begin{cases}
w_{1,a}(x,t-\sum_{n=1}^{s_a}\big[\lambda_a(n)^{-1}\big]_a), & s_a>0,\\
w_{1,a}(x,t),& s_a=0,\\
w_{1,a}(x,t+\sum_{n=1}^{|s_a|}\big[\lambda_a(-n)^{-1}\big]_a), & s_a<0,
\end{cases}&
\prod_{n=1}^{s_a}(-\lambda_a(n))^{-1},&s_a>0\\
\prod_{n=1}^{|s_a|}(-\lambda_a(-n)),&s_a<0,
\end{cases}\\
\bar c_b w_{2,b}(x,t,\bar s_b)&=
\begin{cases}
w_{2,b}(x,t+\sum_{n=1}^{\bar s_b}\overline{\big[\bar\lambda_b(n)^{-1}\big]}_b,x), & \bar s_b>0,\\
w_{2,b}(x,t),& \bar s_b=0,\\
w_{2,b}(x,t-\sum_{n=1}^{|\bar s_b|}\overline{\big[\bar \lambda_b(-n)^{-1}\big]}_b), & \bar s_b<0,
\end{cases}&
\bar c_b:=\begin{cases}
\prod_{n=1}^{\bar s_b}(-\bar \lambda_b(n))^{-1},&\bar s_b>0\\
1, & \bar s_b=0\\
\prod_{n=1}^{|\bar s_b|}(-\bar\lambda_b(-n)),&\bar s_b<0.
\end{cases}
\end{aligned}
\end{align}
\end{pro}
As a conclusion, the discrete flows and Miwa shifts in the continuous flows are the very same thing, and therefore we could  work with continuous flows and Miwa transformations or with continuous/discrete flows. This discussion justifies the
Miwa part in the name we gave to these discrete flows.

\subsubsection{Bounded from below measures}
Of course, in order to preserve the link with multiple orthogonal polynomials, these discrete flows must  preserve the
reality, regularity and sign constance  of the
weights, which generically is not the case.  When the support of the weights is bounded from below, i.e. there are
finite real numbers  $K_a$ and $K_b$,
such that $supp (w_1,aμ̣)⊂[K_a,∞)$ and $supp (w_2,bμ̣)⊂[K̅_b,∞)$,
a possible solution is to place all the new
zeroes and poles in the real line but outside the corresponding support, $λ_a(n)<inf(supp
(w_1,aμ̣))$ and
$λ̅_b(n)<inf(supp (w_2,bμ̣))$. A different approach, which will be considered in Appendix \ref{III}, is
to arrange the zeroes in complex conjugate

To analyze the consequence of the discrete flows  on the integrable hierarchy we introduce two sets of shifts
\begin{definition}
\begin{enumerate}
\item Let us consider the sets of shift operators $\{T_a\}_{a=1}^{p_1}$ and $\{\bar T_b\}_{b=1}^{p_2}$, where
$T_a$ stands for the shift $s_a\mapsto s_a+1$ and $\bar T_b$ stands for $s_b\mapsto \bar s_b+1$. The rest of the variables $\{s_{a'},\bar s_{b'}\}$ will
remain constant.
\item  We introduce
\begin{align}
\begin{aligned}
q_a&:=  \I-\Pi_{1,a}(\I+\lambda_a(s_a+1))+\Lambda_{1,a},\\
\bar q_b&:=   \I  -\Pi_{2,b}(\I+\bar \lambda_b(\bar s_b+1))+\Lambda_{2,b}^\top.
\end{aligned}
\end{align}

\item We also define the operators
\begin{align}\label{delta}
\delta_a&:=Sq_aS^{-1}=\I-C_{a}(\I+\lambda_a(s_a+1))+L_a,&
\bar\delta_b&:=\bar S \bar q_b\bar S^{-1}=\I-\bar C_{b}(\I+\bar\lambda_b(\bar s_b+1))+\bar L_b,\\
C_a&:=S\Pi_{1,a}S^{-1},& \bar C_b&:=\bar S\Pi_{2,b}\bar S^{-1}.\notag
\end{align}
Here the matrices $\delta_a$ and $\bar\delta_b$  are called lattice resolvents.

\item  Finally the semi-infinite wave matrices
\begin{align}\label{discrete-wave}
W&:=SD_0,& \bar W&:=\bar S\bar D_0.
\end{align}
\end{enumerate}
\end{definition}
Observe that
\begin{align}
\begin{aligned}
(T_a D_0)D_0^{-1}&=q_a,&\bar D_0^{-1} (T_a\bar D_0)&=\I,\\
\bar D_0(\bar T_b\bar D_0^{-1})&=\bar q_b,& (\bar T_b  D_0) D_0^{-1}&=\I.&
\end{aligned}
\end{align}
When we   assume that the semi-infinite matrices $δ_a$ and $δ̅_b$ are $LU$ factorizable as in
\eqref{facto}, i.e. all their principal minors do not
vanish, we can write
\begin{align}\label{fac.delta}
\delta_a&=\delta_{a,-}^{-1}\delta_{a,+},& \bar \delta_b&=\bar \delta_{b,-}^{-1}\bar\delta_{b,+},
\end{align}
where $δ_a,-$ and $δ̅_b,-$ are  lower matrices as is  $S$  in \eqref{facto}, and $δ_a,+$ and
$δ̅_b,+$ are upper matrices as
$S̅$ in \eqref{facto}. We now show  that when the deformed moment matrix $g(s)$ is factorizable, and therefore the
multiple orthogonality  makes sense, the following  holds
\begin{pro}
If the deformed moment matrix $g(s)$ is factorizable for all values of $s$ then  so is $\delta_a$ and $\bar
\delta_b$ with
\begin{align}
\begin{aligned}
\delta_{a,+}&= (T_a\bar S)\bar S^{-1}, &   \delta_{a,-}&=(T_a S)S^{-1},\\
\bar\delta_{b,+}&=(\bar T_b\bar
S)\bar S^{-1}, &  \bar\delta_{b,+}&=(\bar T_b S)S^{-1}.
\end{aligned}
\end{align}
\end{pro}
\begin{proof}
When we apply the discrete shifts to the Gauss--Borel factorization problem $g(s)=S^{-1}(s)\bar S(s)$ we get
\begin{align*}
T_a(S^{-1})T_a(\bar S)&=T_a g(s)=(T_aD_0)D_0^{-1} g(s)=q_ag(s) &&\Rightarrow &((T_a S)S^{-1})^{-1}(T_a\bar S)\bar
S^{-1}&=\delta_a ,\\
\bar T_b(S^{-1})\bar T_b(\bar S)&=\bar T_b g(s)=g(s)\bar D_0(T_b\bar D_0^{-1})=g(s)\bar q_b &&\Rightarrow & ((\bar
T_b S)S^{-1})^{-1}((\bar T_b\bar
S)\bar S^{-1}&=\bar \delta_b,
\end{align*}
and the desired result follows.
\end{proof}
Therefore, we can consider the following
\begin{definition}
The semi-infinite matrices $\omega_a$ and $\bar\omega_b$ are given by
\begin{align}\label{omega}
\omega_a&:=\delta_{a,-} \delta_a=\delta_{a,+},&
\bar \omega_b&:=\bar \delta_{b,-}=\bar \delta_{b,+} \bar \delta_b^{-1},
\end{align}
\end{definition}
and show that
\begin{pro}
\begin{enumerate}
\item    The following auxiliary linear systems
\begin{align}\label{d.wave}
\begin{aligned}
T_aW &=\omega_a W,&T_a\bar W &=\omega_a \bar  W,\\
\bar  T_bW &=\bar \omega_b W,&\bar T_b\bar W &=\bar \omega_b \bar  W,
\end{aligned}
\end{align}
are satisfied.
\item The Lax matrices fulfill the following relations
\begin{align}\label{d.lax}
\begin{aligned}
T_aL_{a'}&=\omega_aL_{a'}\omega_a^{-1},  &  T_a\bar L_{b}&=\omega_a\bar L_{b}\omega_a^{-1},\\
\bar T_bL_{a}&=\bar\omega_bL_{a}\bar\omega_b^{-1},&  \bar T_b\bar L_{b'}&=\bar\omega_b\bar
L_{b'}\bar\omega_b^{-1}.        \end{aligned}
\end{align}
\item The following discrete Zakharov--Shabat compatibility conditions hold
\begin{align}\label{d.ZS2}
( T_a\omega_{a'})\omega_{a}&=( T_{a'}\omega_{a})\omega_{a'},&
( T_a\bar \omega_{b})\omega_{a}&=( \bar T_{b}\omega_{a})\bar \omega_{b},&
( \bar T_b\bar \omega_{b'})\bar\omega_{b}&=( \bar T_{b'}\bar\omega_{b})\bar \omega_{b'}.
\end{align}
\item When the discrete and continous flows are considered simultaneously, the
following equations
\begin{align}\label{d.ZS1}\begin{aligned}
\bar
T_{a}\bar B_{j,b}&=(\bar\partial_{j,b}\omega_{a})\omega_{a}^{-1}+\omega_{a}\bar
\bar  T_{b'}\bar B_{j,b}=(\bar\partial_{j,b}\bar\omega_{b'})\bar\omega_{b'}^{-1}+\bar\omega_{b'}\bar
\end{aligned}
\end{align}
are obtained.
\end{enumerate}
\end{pro}
\begin{proof}
We compute
\begin{align*}
T_a W&=(T_a S)(T_a D_0)=(T_a S) S^{-1} S q_a S^{-1} S D_0=\delta_{a,-} \delta_a W=\delta_{a,+} W,\\
T_a \bar W&=(T_a\bar S)\bar D_0=(T_a\bar S)\bar S^{-1} \bar S \bar D_0=\delta_{a,+}\bar W,\\
\bar  T_b  W&=(  \bar  T_b S) D_0=( \bar  T_b S) S^{-1}  S  D_0=\bar \delta_{b,-}  W,\\
\bar  T_b \bar W&=(\bar T_b \bar S)(T_b \bar D_0)=(T_b \bar S) \bar S^{-1} \bar S \bar q_b^{-1} \bar S^{-1}
\bar S \bar D_0=\bar \delta_{b,+} \bar
\delta_b^{-1} \bar W=\bar \delta_{b,-}\bar W,
\end{align*}
from where we deduce  \eqref{d.wave}, which in turn imply \eqref{d.lax} and \eqref{d.ZS2}.

\end{proof}
The simultaneous consideration of continuous and discrete flows leads to the replacement  $W_0→W_0D_0$ and $W̅_0→W̅_0D̅_0$, and the
corresponding modification of the weight's flows  is achieved by the multiplication of the continuous and discrete
evolutionary factors, in this context we
also have \eqref{d.ZS1}. These discrete flows could be understood as a sequence of Darboux transformations of $LU$ and $UL$ types in the
terminology of [5], which motivates the Darboux part in name we give to these discrete flows.
In fact,  we have that the lattice resolvents satisfy
\begin{align*}
\delta_a&=\delta_{a,-}^{-1}\delta_{a,+}&&\Rightarrow&
\delta_{a,+}\delta_{a,+}^{-1}=\delta_{a,+}\delta_{a,-}^{-1},\\
\bar\delta_b&=\bar\delta_{b,-}^{-1}\bar\delta_{b,+}&&\Rightarrow&
\bar T_b\bar\delta_b&=\bar\omega_b\bar\delta_b\bar\omega_b^{-1}=\bar\delta_{b,-}\bar\delta_{b,-}^{-1}
\bar \delta_{b,+}\bar \delta_{b,-}^{-1}=\bar\delta_{b,+}\bar\delta_{b,-}^{-1},
\end{align*}
which amounts to the typical permutation of the $LU$ factorization to the $UL$ factorization. When there is only one
component we have $δ=L+λ$ and
$δ̅=λ̅+L̅$ and the shift corresponds to the classical $LU$ or $UL$ Darboux transformations.

If  $A^(k)_a(x,s),A̅_b^(k)(x,s)$ are the multiple orthogonal polynomials and  dual polynomials in the
$x$ variable  corresponding to the  discrete evolution of the  weights \eqref{d.evolved} respectively we have
the discrete version of Proposition \ref{conti}
\begin{pro}
The wave and adjoint wave functions \eqref{discrete-wave} are
\begin{align}
\Psi_a^{(k)}(z,s)&=A^{(k)}_a(z,s)\Ds_a(z,s_a)&
( \bar \Psi_b^*)^{(k)}(z,s)&=\bar A_b^{(k)}(z,s)\bar \Ds_b(z,\bar s_b)^{-1},
\end{align}
and the  the linear forms
\begin{align}\label{linear.form.baker.2}
Q^{(k)}(x,s)&=\sum_{a=1}^{p_1}A_a^{(k)}(x,s)w_{1,a}(x,s), &   \bar Q^{(k)}(x,s)&=\sum_{b=1}^{p_2}(\bar
\end{align}
associated with the weights  $w_{1,a}(x,s), w_{2,b}(x,s)$, can be expressed as
\begin{align}\label{linear.form.baker.1}
Q^{(k)}(x,s)&:=\sum_{a=1}^{p_1}\Psi_a^{(k)}(x,s)w_{1,a}(x), &   \bar
\end{align}
in terms of which we have the equations
\begin{align}\label{cauchy.baker}
\bar  \Psi_b^{(k)}(z,s)&=\int \frac{Q^{(k)}(x,s)}{z-x}w_{2,b}(x)\d \mu(x),&
( \Psi_a^*)^{(k)}(z,s)&=\int \frac{\bar Q^{(k)}(x,s)}{z-x}w_{1,a}(x)\d \mu(x).
\end{align}
\end{pro}
Here the Cauchy transforms must be interpreted in exactly the same terms as in Proposition \ref{pro:cauchy  tr}.
Observe that \eqref{cauchy.baker} do not correspond to the functions of the second kind
\begin{align}%\label{cauchy.baker}
\bar  C_b^{(k)}(z,s)&:=\int \frac{Q^{(k)}(x,s)}{z-x}w_{2,b}(x,s)\d \mu(x),&
C_a^{(k)}(z,s)&:=\int \frac{\bar Q^{(k)}(x,s)}{z-x}w_{1,a}(x,s)\d \mu(x).
\end{align}
Notice also that from \eqref{delta}-\eqref{fac.delta}, relations that hold true for any $g$ and not only
for the moment matrix, we get
\begin{lemma}\label{lemon}
We have that
\begin{align}\label{form.omega}\begin{aligned}
\omega_a&=\omega_{a,0}\Lambda^{|\vec n_1|-n_{1,a}+1}+\omega_{a,1}\Lambda^{|\vec
\bar \omega_b&=\bar\omega_{b,0}(\Lambda^\top)^{|\vec
n_2|-n_{2,b}}+\dots+\bar\omega_{b,|\vec n_2|-n_{2,b}+1},\\
\omega_a^\top&=\rho_{a,0}(\Lambda^\top)^{|\vec n_1|-n_{1,a}+1}+\rho_{a,1}(\Lambda^\top)^{|\vec
\bar \omega_b^\top&=\bar\rho_{b,0}\Lambda^{|\vec n_2|-n_{2,b}+1}+\bar\rho_{b,1}\Lambda^{|\vec
\end{aligned}
\end{align}
for some diagonal semi-infinite matrices \begin{align}
\begin{aligned}
\omega_{a,j}&=\diag(\omega_{a,j}(0),\omega_{a,j}(1),\dots),\\
\bar\omega_{b,j}&=\diag(\bar\omega_{b,j}(0),\bar\omega_{b,j}(1),\dots),\\
\rho_{a,j}&=\diag(\rho_{a,j}(0),\rho_{a,j}(1),\dots),\\
\bar\rho_{b,j}&=\diag(\bar\rho_{b,j}(0),\bar\rho_{b,j}(1),\dots),
\end{aligned}
\end{align}
\begin{align}
\begin{aligned}
\rho_{a,j}(k)&:=\omega_{a,j}(k-|\vec n_1|+n_{1,a}-1+j),\\
\bar\rho_{b,j}(k)&:=\bar\omega_{b,j}(k+|\vec n_2|-n_{2,b}+1-j).
\end{aligned}
\end{align}
\end{lemma}
that with
\begin{definition}
We define
\begin{align}\label{gamma}\begin{aligned}
\gamma_{a,a'}(s,x)&:=  (1-\delta_{a,a'}(1+\lambda_a(s_a+1)-x)),\\
\gamma_{b,b'}(s,x)&:=  (1-\delta_{b,b'}(1+\bar\lambda_b(\bar s_b+1)-x)),
\end{aligned}
\end{align}
\end{definition}
leads to
\begin{pro}
The following equations
\begin{align}\label{multi.variante}
\begin{aligned}
(T_{a'} A_a^{(k)})\gamma_{a,a'}&=\omega_{a',0}(k)A_a^{(k+|\vec n_1|-n_{1,a'}+1)}+\dots+\omega_{a',|\vec
n_1|-n_{1,a'}+1}(k) A_a^{(k)},\\
\bar T_{b'}A_a^{(k)}&=\bar\omega_{b',0}(k)A_a^{(k-|\vec n_2|+n_{2,b'}-1)}+\dots+\bar\omega_{b',|\vec n_2|-n_{2,b'}+1}
(k) A_a^{(k)}.
\end{aligned}\\
\label{dual.multi.variante}
\begin{aligned}
\rho_{a',0}(k) (T_{a'} \bar A_b^{(k-|\vec n_1|+n_{1,a'}-1)})+\dots+ \rho_{a',|\vec n_1|-n_{1,a'}+1}(k) (T_{a'} \bar
A_b^{(k)})&= \bar A_b^{(k)},\\
\Big( \bar\rho_{b',0} (k)(\bar T_{b'} \bar A_b^{(k+|\vec n_2|-n_{2,b'}+1)})+
\dots+ \bar\rho_{b',|\vec n_2|-n_{2,b'}+1} (k) (\bar T_{b'} \bar A_b^{(k)})\Big)\bar\gamma_{b,b'}&= \bar A_b^{(k)},
\end{aligned}
\end{align}
are fulfilled.
\end{pro}
\begin{proof}
For \eqref{multi.variante} recall the discrete auxiliary systems for $W$, while for \eqref{dual.multi.variante} just
consider that \begin{align*}
\omega_a^\top T_a((\bar W^{-1})^\top)&= (\bar W^{-1})^\top,&
\bar  \omega_b^\top \bar T_b((\bar W^{-1})^\top)&= (\bar W^{-1})^\top.
\end{align*}
\end{proof}
Notice that relations \eqref{multi.variante} and \eqref{dual.multi.variante} are among multiple orthogonal
polynomials in the same ladder but with different weights, they link the polynomials for the weights $w_1,a,
w_2,b$ with those with $T_a'w_1,a,T_a' w_2,b$ or  $T̅_b'w_1,a,T̅_b' w_2,b$.

\subsection{Symmetries, recursion relations and string equations }

We now return to the discussion of the symmetry of the moment matrix that we started in \S \ref{sec:recursion} but
with  evolved weights and the use of Lax matrices. The first observation is the following
\begin{pro}
The $j$-th power of the evolved Jacobi type matrix introduced in \S \ref{sec:recursion} is related with Lax matrices
through what we call a string equation:  \begin{align}\label{bigradedj}
J^{j}&=\sum_{a=1}^{p_1}L_a^j=\sum_{b=1}^{p_2}\bar L_b^j,& j=1,2,\dots,
\end{align}
and the multiple orthogonal polynomials are  eigen-vectors:
\begin{align}\label{bigradedjp}
(J^{j})^\top\bar\A_{b'}&=x^j\bar \A_{b'},
\end{align}
for $a'=1,\dots,p_1$ and $b'=1,\dots,p_2$.
\end{pro}
\begin{proof}
Using \eqref{sym} it can be proven by induction on $j$ that for any $j\geq 1$ the following equation holds
\begin{align}\label{symj}
\Lambda_{1,a}^j g  \Pi_{2,b}=
\Pi_{1,a}g(\Lambda_{2,b}^\top)^j,
\end{align}
so that
\begin{align}\label{stringj}
L_{a}^j\bar C_{b}= C_{a}\bar L_{b}^j.
\end{align}
Summing over $a,b$ we deduce \eqref{bigradedj}. Moreover \eqref{bigradedjp} is obtained as follows
\begin{align}
(J^j)^\top\bar\A_{a'}&=(\bar S^{-1})^\top\sum_{b=1}^{p_2}\Lambda_{2,a}^j\bar S^\top  (\bar
S^{-1})^\top\chi_{2,b'}=x^j\bar \A_{b'}.
\end{align}
\end{proof}

We are ready to show that the symmetry  \eqref{sym} induces a  corresponding invariance on Lax matrices and multiple
orthogonal polynomials
\begin{pro}\label{proposition: symmetry}
The following relations hold for $j = 1,2,\dots$
\begin{align}
\big(\sum_{a=1}^{p_1}\frac{\partial}{\partial t_{j,a}}+\sum_{b=1}^{p_2}\frac{\partial}{\partial \bar t_{j,
\big(\sum_{a=1}^{p_1}\frac{\partial}{\partial t_{j,a}}+\sum_{b=1}^{p_2}\frac{\partial}{\partial \bar t_{j,
\big(\sum_{a=1}^{p_1}\frac{\partial}{\partial t_{j,a}}+\sum_{b=1}^{p_2}\frac{\partial}{\partial \bar t_{j,
\big(\sum_{a=1}^{p_1}\frac{\partial}{\partial t_{j,a}}+\sum_{b=1}^{p_2}\frac{\partial}{\partial \bar t_{j,
\end{align}
\end{pro}
\begin{proof}
See Appendix \ref{II}.
\end{proof}

\subsection{Bilinear equations and $\tau$-functions}
The proof of the bilinear identity needs three lemmas. For the first one, let  $W_n⃗_1,n⃗_2,W̅_n⃗_1,n⃗_2$ be  the wave matrices associated with the moment matrix $g_n⃗_1,n⃗_2$; so that, $W_n⃗_1,n⃗_2g_n⃗_1,n⃗_2=W̅_n⃗_1,n⃗_2$. Then, we have
\begin{lemma}
\label{lemma: matrix bilinear}
The wave matrices associated with different compositions and times satisfy
\begin{align}\label{bilinear-wave}
W_{\vec n_1,\vec n_2}(t,s)\pi_{\vec n_1',\vec n_1}^\top W_{\vec n_1',\vec n_2'}(t',s')^{-1}=\bar W_{\vec n_1,\vec
n_2}(t,s)\pi_{\vec n_2',\vec n_2}^\top\bar W_{\vec n_1',\vec n_2'}(t',s')^{-1},
\end{align}
\end{lemma}
\begin{proof}
We consider simultaneously the following equations
\begin{align*}
W_{\vec n_1,\vec n_2}(t,s)g&=\bar W_{\vec n_1,\vec n_2}(t,s),\\
W_{\vec n_1',\vec n_2'}(t',s')\pi_{\vec n_1',\vec n_1} g\pi_{\vec n_2',\vec n_2}^\top&=\bar W_{\vec n_1',\vec
\end{align*}
where $g=g_{\vec n_1,\vec n_2}$, and we get
\begin{align*}
W_{\vec n_1,\vec n_2}(t,s)^{-1}\bar W_{\vec n_1,\vec n_2}(t,s)= \pi_{\vec n_1',\vec n_1}^\top W_{\vec n_1',\vec
n_2'}(t',s')^{-1} \bar W_{\vec n_1',\vec n_2'}(t',s')\pi_{\vec n_2',\vec n_2}=g,
\end{align*}
and the result becomes evident.
\end{proof}
For the second one, let  $(·)_-1$ denote the   coefficient in $z^-1$ in the Laurent expansi\'{o}n around
$z=∞$ (place where the Cauchy transforms make sense).
\begin{lemma}\label{chi-id}
For the vectors $\chi_a$ the following formulae hold
\begin{align*}
\big(\sum_{a=1}^{p}\chi_a(\chi_a^*)^\top\big)_{-1}=\big(\sum_{a=1}^{p}\chi^*_a\chi_a^\top\big)_{-1}=\I,
\end{align*}
\end{lemma}
and therefore
\begin{lemma}
For any couple of semi-infinite matrices $U$ and $V$ we have
\begin{align}\label{residue1}
UV&=\Big(  \sum_{a=1}^{p_1}(U\chi_{1,a}) \big(V^\top\chi_{1,a}^*\big)^\top\Big)_{-1}\\\label{residue2}
&=\Big(  \sum_{b=1}^{p_2}(U\chi_{2,b}^*)\big(V^\top\chi_{2,b}\big)^\top\Big)_{-1},
\end{align}
\end{lemma}
\begin{proof}
It follows easily from Lemma \ref{chi-id}:
\begin{align*}
\Big(  \sum_{a=1}^{p_1}(U\chi_{1,a}) \big(V^\top\chi_{1,a}^*\big)^\top\Big)_{-1}&=
U    \Big(\sum_{a=1}^{p_1}\chi_{1,a}(\chi_{1,a}^*)^\top \Big)_{-1}V=UV,\\
\Big(  \sum_{b=1}^{p_2}(U\chi_{2,b}^*)\big(V^\top\chi_{2,b}\big)^\top\Big)_{-1}&=
U \Big(\sum_{b=1}^{p_2}\chi_{2,b}^*\chi_{2,b}^\top\Big)_{-1}V=UV.
\end{align*}
\end{proof}
We have the following
\begin{theorem}\label{theorem: bilinear}
\begin{enumerate}
\item The wave functions and their companions satisfy
\begin{align*}
\sum_{a=1}^{p_1}\oint_{\infty}\Psi^{(k)}_{\vec n_1,\vec n_2,a}(z,t,s)(\Psi^*_{\vec n_1',\vec
n_2',a})^{(l)}(z,t',s')\d z=
\sum_{b=1}^{p_2}\oint_{\infty}\bar\Psi^{(k)}_{\vec n_1,\vec n_2,b}(z,t,s)(\bar\Psi^*_{\vec n_1',\vec
n_2',b})^{(l)}(z,t',s')\d z.
\end{align*}
\item Multiple orthogonal polynomials, their duals and the corresponding second kind functions  are linked by
\begin{align}\label{bilinear2}
\sum_{a=1}^{p_1}\oint_{\infty}
A^{(k)}_{\vec n_1,\vec n_2,a}(z,t,s)\bar C_{\vec n_1',\vec n_2',a}^{(l)}(z,t',s')E_a(z)\d z=
\sum_{b=1}^{p_2}\oint_{\infty}C^{(k)}_{\vec n_1,\vec n_2,b}(z,t,s)\bar A_{\vec n_1',\vec
n_2',b}^{(l)}(z,t',s')\bar E_b(z)\d z,
\end{align}
\begin{align*}
\bar E_b&:=(\bar \Es_b\bar\Ds_b)(z,t,s)((\bar\Es_b\bar\Ds_b)(z,t',s'))^{-1}.
\end{align*}
\end{enumerate}
\end{theorem}
\begin{proof}
\begin{enumerate}
\item   If we set in  \eqref{residue1} $U=     W_{\vec n_1,\vec n_2}(t,s)$ and $V=\pi_{\vec n_1',\vec
n_1}^\top W_{\vec n_1',\vec n_2'}(t',s')^{-1}$
and in  \eqref{residue2} we put $U=\bar W_{\vec n_1,\vec n_2}(t,s)$ and $V= \pi_{\vec n_2',\vec n_2}^\top\bar
W_{\vec n_1',\vec n_2'}(t',s')^{-1}$
attending to
\eqref{bilinear-wave}, recalling that $\Psi_{\vec n_1,\vec n_2,a}=W_{\vec n_1,\vec n_2}\chi_{\vec n_1,a}$,
$\bar \Psi_{\vec n_1,\vec n_2,b}=\bar W_{\vec n_1,\vec n_2}\chi^*_{\vec n_2,b}$ and observing that
$\Psi_{\vec n_1',\vec n_2',a}^*=(W_{\vec n_1',\vec n_2'}^{-1})^\top\pi_{\vec n_1',\vec n_1}\chi^*_{\vec n_1,a}$
$\bar \Psi_{\vec n_1',\vec n_2',b}^*=(\bar W_{\vec n_1',\vec n_2'}^{-1})^\top\pi_{\vec n_2',\vec n_2}\chi_{\vec
n_2,b}$  we get  the stated  bilinear equation for the wave functions.\footnote{The reader familiarized with
Toda bilinear equations should
notice  that in the right hand term we are working at $z=\infty$ instead of, as customary, at $z=0$; the reason
is that for the definition of $\bar\A_b$ we
have used $\chi_2$ instead of $\chi_2^*$, in order to get polynomials in $z$, while normally one gets
polynomials in $z^{-1}$. See footnote 2. Moreover,
due to the redefinition of the wave functions there is no $\frac{\d z}{2\pi\text{i}z}$ factor}
\item
We can   write
\begin{multline*}
W_{\vec n_1,\vec n_2}(t,s)\pi_{\vec n_1',\vec n_1}^\top W_{\vec n_1',\vec n_2'}(t',s')^{-1}=(S_{\vec n_1,\vec
n_2}(t,s)W_{0,\vec n_1}(t,s)\pi_{\vec n_1',\vec n_1}^\top (W_{0,\vec n_1'}(t',s'))^{-1}\pi_{\vec n_1',\vec
n_1})\pi_{\vec n_1',\vec n_1}^\top S_{\vec n_1',\vec n_2'}(t',s')^{-1},
\end{multline*}
which strongly suggests to consider in  \eqref{residue1}
\begin{align*}
U&=S_{\vec n_1,\vec n_2}(t,s)W_{0,\vec n_1}(t,s)\pi_{\vec n_1',\vec n_1}^\top(W_{0,\vec
n_1'}(t',s'))^{-1}\pi_{\vec n_1',\vec n_1},& V&=\pi_{\vec n_1',\vec n_1}^\top S_{\vec n_1',\vec
\end{align*}
\begin{multline*}
\bar W_{\vec n_1,\vec n_2}(t,s)\pi_{\vec n_1',\vec n_1}^\top \bar W_{\vec n_1',\vec n_2'}(t',s')^{-1}=(\bar
S_{\vec n_1,\vec n_2}(t,s)\bar W_{0,\vec n_1}(t,s)\pi_{\vec n_1',\vec n_1}^\top (\bar W_{0,\vec
n_1'}(t',s'))^{-1}\pi_{\vec n_1',\vec n_1})\pi_{\vec n_1',\vec n_1}^\top \bar S_{\vec n_1',\vec
\end{multline*}
suggest to set in   \eqref{residue2}
\begin{align*}
U&=\bar S_{\vec n_1,\vec n_2}(t,s)\bar W_{0,\vec n_2}(t,s)\pi_{\vec n_2',\vec n_2}^\top(\bar W_{0,\vec
n_2'}(t',s'))^{-1}\pi_{\vec n_2',\vec n_2},& V&=\pi_{\vec n_2',\vec n_2}^\top \bar S_{\vec n_1',\vec
\end{align*}
The application of \eqref{residue1},\eqref{residue2} and \eqref{bilinear-wave} gives the alternative bilinear
relations \eqref{bilinear2}
where we have used  the evolved Cauchy transforms \eqref{evolved.ct} and  introduce the evolutionary factors
\begin{align*}
\bar E_b&:=(\bar \Es_b\bar\Ds_b)(z,t,s)((\bar\Es_b\bar\Ds_b)(z,t',s'))^{-1}.
\end{align*}
\end{enumerate}
\end{proof}
The factors involved in this definition were introduced in \eqref{evolved} and \eqref{d.evolved}, so that we assume
the discrete flows within the bounded from
below support scenario, while if we consider the two-step discrete flows the replacement of the $$-factors by the
$'$- factors \eqref{d.evolved2} is

It can be shown that for certain weights, for which we can use the Fubini and Cauchy theorems, and when one only
considers a finite number of continuous flows that the r.h.s and l.hs. in this bilinear relations are proportional to
$∫_ Q_n⃗_1,n⃗_2^(k)(x,t)Q̅_n⃗_1',n⃗_2'^(l)(x,t')μ̣(x)$. This is a direct
consequence of
\begin{pro}\label{proposition: bilinear fubini-cauchy}
We have the following identity
\begin{align}
\int_{\R} Q_{\vec n_1,\vec n_2}(x,t,s)\bar Q_{\vec n_1',\vec n_2'}^\top(x,t',s')\d \mu(x)&= W_{\vec n_1,\vec
n_2}(t,s)\pi_{\vec n_1',\vec n_1}^\top (W_{\vec n_1',\vec n_2'}(t',s'))^{-1}\\&=\bar W_{\vec n_1,\vec
n_2}(t,s)\pi_{\vec n_2',\vec n_2}^\top (\bar W_{\vec n_1',\vec n_2'}(t',s'))^{-1}.
\end{align}
\end{pro}
\begin{proof}
See Appendix \ref{II}.
\end{proof}
%=(W(t)W(t')^{-1})_{k,l}=(\bar W(t)\bar W(t')^{-1})_{k,l}$, or $2\pi\text{i}\int_{\R} Q(x,t)\bar Q(x,t')^\top \d\mu=W(t)W(t')^{-1}=\bar W(t)\bar W(t')^{-1}$ --which follows from the \eqref{g.evol} and \eqref{compact.g}.

%Recalling the expressions of the wave functions in terms of the multiple orthogonal polynomials, the duals and the corresponding Cauchy transforms we get that the previous bilinear reads
%\d z\int  \d\mu   \sum_{\substack{a=1,\dots,p_1\\
%\frac{A^{(k)}_a(z,t,s)\bar A^{(l)}_b(x,t',s')-A^{(k)}_a(x,t,s)\bar A^{(l)}_b(z,t',s')}{z-x}w_{1,a}(x,t,s)w_{2,b}(x,t',s')=0

%For example, if $\vec n_1=(2,2,2)$ and $\vec n_2=(1,1)$, we have the following table
%  \hline
%  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
%  $l$ & 0 & 1 & 2 & 3 & 4 & 5  & 6 &7&8&9&10\\ \hline
%  $\vec\nu_1$ & (0,0,0) & (1,0,0) & (2,0,0) & (2,1,0) & (2,2,0) & (2,2,1)& (2,2,2) & (3,2,2)&(4,2,2)& (4,3,2)&(4,4,2)\\
%  $\vec \nu_2$ & (0,0) & (1,0) & (1,1) & (2,1) & (2,2) & (3,2) & (3,3) &(4,3)&(4,4)& (5,4)&(5,5)\\
%  \hline
%Now $n\bar n=2$, $n\bar np_1=6$ and $n\bar np_2=4$, and for example  $x^2A_{(3,2,2),(3,3),a}$, $a=1,2,3$, i.e. $l=6$, is expressed as a linear combination of the 11 polynomials $A_{\vec\nu_1(l),\vec\nu_2(l-1),j}$, $l=1,..,11$.

Now,  we will perform a full characterization of the $τ$-functions associated with the multiple orthogonal polynomials
defined in this paper.
\begin{definition}
Let us define the following matrices
\begin{align}\label{def_g+a}
\left(\begin{BMAT}{ccccc}{cccc:c}
g_{1,0}&g_{1,1}&\cdots& g_{1,l-1}&g_{1,l}\\
\vdots&\vdots&&\vdots&\vdots \\
g_{l-1,0}&g_{l-1,1}&\cdots& g_{l-1,l-1}&g_{l-1,l}\\
g_{l_{+a},0}&g_{l_{+a},1}&\cdots& g_{l_{+a},l-1} &g_{l_{+a},l}
\end{BMAT}\right)&
\bar g^{[l+1]}_{+b}:&=
\left(\begin{BMAT}{cccc:c}{ccccc}
g_{0,0}&g_{0,1}&\cdots& g_{0,l-1} &g_{0,\bar l_{+b}}\\
g_{1,0}&g_{1,1}&\cdots& g_{1,l-1} &g_{1,\bar l_{+b}}\\
\vdots&\vdots&&\vdots&\vdots\\
g_{l-1,0}&g_{l-1,1}&\cdots& g_{l-1,l-1}&g_{l-1,\bar l_{+b}}\\
g_{l,0}&g_{l,1}&\cdots& g_{l,l-1}&g_{l,\bar l_{+b}}
\end{BMAT}\right).
\end{align}
\end{definition}
The matrix $g^[l+1]_+a$ is obtained from  $g^[l+1]$ replacing the last row (operation denoted by a dashed line)
\begin{align*}
\end{align*} and $g̅^[l+1]_+b$is obtained from  $g^[l+1]$ replacing the last column  by $(  g_0,l̅_+b,g_1,l̅_+b,…,g_l-1,l̅_+b,g_l,l̅_+b)^⊤$. It is clear that if $a_1(l)=a$ then $g^[l+1]_+a=g^[l+1]$ and if $a_2(l)=b$ then $g̅^[l+1]_+b=g^[l+1]$.
%(the first not yet appearing in the $b$-sector).

The minors of the these matrices \eqref{def_g+a}  will be denoted as $M^[l+1]_i,j=M̅^[l+1]_i,j$ for
$g^[l+1]$, $M^[l+1]_+a,i,j$ for $g^[l+1]_+a$ and $M̅^[l+1]_+b,i,j$ for $g̅^[l+1]_+b$.
Now we introduce the following determinants that are cofactors of the previously defined matrices
\begin{definition} \label{def_tau1}
The $\tau$-functions  are defined as follows
\begin{align} \label{def_taus}
\begin{aligned}
\tau_{+a,-a'}^{(l)}&:=(-1)^{l+l_{-a'}}M_{+a,l_{-a'},l}^{[l+1]}, &
\tau_{-b,-a}^{(l)}&:=(-1)^{\bar l_{-b}+l_{-a}}M_{l_{-a},\bar l_{-b}}^{[l+1]},
\end{aligned}\\
\label{def_dualtaus}
\begin{aligned}
\bar \tau_{+b,-b'}^{(l)}&:=(-1)^{l+\bar l_{-b'}}\bar M_{+b,l,\bar l_{-b'}}^{[l+1]},&
\bar \tau_{-a,-b}^{(l)}&:=(-1)^{l_{-a}+\bar l_{-b}}\bar M_{l_{-a},\bar l_{-b}}^{[l+1]}.
\end{aligned}
\end{align}
\begin{enumerate}
\item If $a_1(l)=a$ then we denote $\tau_{-a'}^{(l)}:=\tau_{+a,-a'}^{(l)}$ and $\bar \tau_{-b}^{(l)}:=\bar
\tau_{-a,-b}^{(l)}$.
\item If $a_2(l)=b$ then we denote $\tau_{-a}^{(l)}:=\tau_{-b,-a}^{(l)}$ and $\bar \tau_{-b'}^{(l)}:=\bar
\tau_{+b,-b'}^{(l)}$.
\item We also introduce $\tau^{(l)}=\bar \tau^{(l)} := \det g^{[l]}$ and
\begin{align*}
\tau^{(l+1)}_{+a}&:=\det   g^{[l+1]}_{+a},& \bar\tau^{(l+1)}_{+b}&:=\det  \bar  g^{[l+1]}_{+b}.
\end{align*}
\end{enumerate}
\end{definition}

If $a_1(l)=a$ then   $τ^(l+1)_+a=τ^(l+1)$, and if $a_2(l)=b$ then

Given a perfect combination $(μ,w⃗_1,w⃗_2)$ and the corresponding  set of multiple orthogonal polynomials
${A_[ν⃗_1;ν⃗_2],a}_a=1^p_1$, with degree vectors such that $|ν⃗_1|=|ν⃗_2|+1$,
there  exists a $(n⃗_1,n⃗_2)$ ladder  and an integer $l$ with $|ν⃗_1|=l+1$  and $|ν⃗_2|=l$  such
that the polynomials  ${A_a^(l)}_a=1^p_1$ coincide with ${A_[ν⃗_1;ν⃗_2],a}_a=1^p_1$.
The final result does not depend upon the particular $(n⃗_1,n⃗_2)$ ladder we choose to get up to the given
degrees in the ladder; however, the $τ$-functions do indeed depend on the ladder chosen through a global sign.
A simple sign-fixing rule is to choose the ladder $ n⃗_1=ν⃗_1$ and $n⃗_2=ν⃗_2+e⃗_p_2$.
We define
\begin{align*}
\tau_{[\vec\nu_1;\vec\nu_2]}&:=\tau^{(l)}_{\vec\nu_1,\vec\nu_2},& l&=|\vec\nu_1|-1=|\vec\nu_2|,
\end{align*}
and we deduce
\begin{pro}\label{tausss}
Given degree vectors $(\vec \nu_1,\vec \nu_2)$ such that $|\vec\nu_1|=|\vec \nu_2|+1$, a composition with  $ \vec
n_1=\vec \nu_1$ and $\vec n_2=\vec \nu_2+\vec e_{p_2}$ and $l=|\vec\nu_1|-1=|\vec\nu_2|$,
we have the following identities
\begin{gather*}% \label{def_taus}
\begin{aligned}
\tau_{+a,-a'}^{(l)}&=\varepsilon_{1,1}( a,a')\tau_{[\vec\nu_1-\vec e_{1,a'}+\vec e_{1,a};\vec\nu_2]},&
\bar \tau_{+b,-b'}^{(l)}&=\varepsilon_{2,2}(b,b')\tau_{[\vec\nu_1;\vec\nu_2-\vec e_{2,b'}+\vec e_{2,b}]},
\end{aligned}\\
\begin{aligned}
&\tau_{-b,-a}^{(l)}=\bar \tau_{-a,-b}^{(l)}&=\varepsilon_{2,1}(b, a)\tau_{[\vec\nu_1+\vec e_{1,p_1}-\vec
e_{1,a};\vec\nu_2+\vec e_{2,p_2}-\vec e_{2,b}]},
\end{aligned}
\end{gather*}
\begin{align*}
\varepsilon_{1,1}(a,a')&:=(-1)^{\sum_{j=1}^{a}\nu_{1,j}+\sum_{j=1}^{a'}\nu_{1,j}+\delta_{a,p_1}-1}, & a'& <a \\
\varepsilon_{1,1}(a,a')&:=(-1)^{\sum_{j=1}^{a}\nu_{1,j}+\sum_{j=1}^{a'}\nu_{1,j}+\delta_{a',p_1}},  & a'& > a \\
\varepsilon_{2,2}(b,b')&:=(-1)^{\sum_{j=1}^{b}\nu_{2,j}+\sum_{j=1}^{b'}\nu_{2,j}-1}, & b'& <b\\
\varepsilon_{2,2}(b,b')&:=(-1)^{\sum_{j=1}^{b}\nu_{2,j}+\sum_{j=1}^{b'}\nu_{2,j}}, & b'& >b\\
\varepsilon_{2,1}(b,a)&:=(-1)^{\sum_{j=1}^{b}\nu_{2,j}+\sum_{j=1}^{a}\nu_{1,j}+\delta_{b,p_2}}, \\
\varepsilon_{1,1}(a,a)&:=1=\varepsilon_{2,2}(b,b).
\end{align*}
In particular
\begin{align*}
\tau^{(l)}_{-a}&= \varepsilon_{1,1}(p_1,a)\tau_{[\vec\nu_1+\vec e_{1,p_1}-\vec e_{1,a};\vec\nu_2]},&
\bar   \tau^{(l)}_{-b}&=\varepsilon_{2,2}(p_2,b)\tau_{[\vec\nu_1;\vec\nu_2+\vec e_{2,p_2}-\vec e_{2,b}]},\\
\tau^{(l+1)}_{+a}&=\varepsilon_{1,1}(a,p_1)\tau_{[\vec \nu_1+\vec e_{1,a};\vec\nu_2+\vec e_{2,p_2}]},
& \bar\tau^{(l+1)}_{+b}&=\varepsilon_{2,2}(b,p_2)\tau_{[\vec \nu_1+\vec e_{1,p_1};\vec\nu_2+\vec e_{2,b}]}.
\end{align*}
\end{pro}
We now proceed to give the $τ$-function representation  of multiple orthogonal polynomials, their duals, second
kind functions and bilinear equations.
The $τ$-functions allow for compact expressions for the multiple orthogonal polynomials:
\begin{pro}\label{taumopis}
The mixed multiple orthogonal polynomials $A_a^{(l)}$, $A_{+a',a}^{(l)}$ and $A_{-b,a}^{(l)}$ have the following
$\tau$-function representation
\begin{align}\label{taumops}
\nu_2(l-1)],a}=z^{\nu_{1,a}(l)-1}\dfrac{\tau^{(l)}_{-a}(t-[z^{-1}]_a)}{\tau^{(l)}(t)}, & l\geq 1, \\
A_{+a',a}^{(l)}(z)&=A^{(\rII,a')}_{[\vec\nu_1(l-1)+\vec e_{1,a'};\vec
\nu_2(l-1)],a}=z^{\nu_{1,a}(l-1)+\delta_{a,a'}-1}\dfrac{\tau^{(l)}_{+a',-a}(t-[z^{-1}]_{a})}{\tau^{(l)}(t)}, & l\geq
\label{taumops-}
A_{-b,a}^{(l)}(z)&=A^{(\rI,b)}_{[\vec\nu_1(l);\vec \nu_2(l)-\vec e_{2,b}],a}=z^{
\nu_{1,a}(l)-1}\dfrac{\tau^{(l)}_{-b,-a}(t-[z^{-1}]_{a})}{\tau^{(l+1)}(t)}, & l\geq 1.
\end{align}
The dual  polynomials $\bar A_b^{(l)}$, $\bar A_{+b',b}^{(l)}$, and $\bar A_{-a,b}^{(l)}$ have the following
$\tau$-function representation
\begin{align}\label{taudualmops}
\bar A_b^{(l)}(z)&=\bar A^{(\rI,a_2(b))}_{[\vec\nu_2(l);\vec \nu_1(l-1)],b}=z^{\nu_{2,b}(l)-1}\dfrac{\bar
\tau_{-b}^{(l)}(t+\overline{[z^{-1}]}_b)}{ \tau^{(l+1)}(t)}, & l\geq 1,\\
\label{taudualmops+}
\bar A_{+b',b}^{(l)}(z)&=\bar A^{(\rII,b')}_{[\vec\nu_2(l-1)+\vec e_{2,b'};\vec \nu_1(l-1)],b}=z^{ \nu_{2,b} (l-1)+
\delta_{b,b'}-1} \dfrac{\bar \tau_{+b',-b}^{(l)}(t+\overline{[z^{-1}]}_{b})}{ \tau^{(l)}(t)}, & l\geq 1,\\
\label{taudualmops-}
\bar A_{-a,b}^{(l)}(z)&=\bar A^{(\rI,a)}_{[\vec\nu_2(l);\vec \nu_1(l)-\vec e_{1,a}],b}=z^{ \nu_{2,b} (l)-1}\dfrac{\bar
\tau^{(l)}_{-a,-b}(t+\overline{[z^{-1}]}_{b})}{ \tau^{(l+1)}(t)}, & l\geq 1.
\end{align}
\end{pro}
\begin{proof}
See Appendix \ref{II}
\end{proof}

Observe that in the simple ladder defined above $(ν⃗_1,ν⃗_2+e⃗_2,p_2)$ with
$l=|ν⃗_2|=|ν⃗_1|-1$ we have
\begin{align*}
\vec\nu_1(l)&=\vec \nu_1, &\vec\nu_2(l-1)&=\vec\nu_2,\\
\vec\nu_1(l-1)&=\vec \nu_1-\vec e_{1,p_1}, &\vec\nu_2(l)&=\vec\nu_2+\vec e_{2,p_2}.
\end{align*}
From  Proposition \ref{taumopis} we get
\begin{align}
\begin{aligned}
A^{(\rII,p_1)}_{[\vec\nu_1;\vec \nu_2],a}(z)&=\varepsilon_{1,1}(p_1,a)z^{\nu_{1,a}-1}\dfrac{\tau_{[\vec\nu_1+\vec
e_{1,p_1}-\vec e_{1,a}a;\vec \nu_2]}(t-[z^{-1}]_a)}{\tau_{[\vec\nu_1;\vec \nu_2]}(t)}, \\
A^{(\rII,a')}_{[\vec\nu_1-\vec e_{1,p_1}+\vec e_{1,a'};\vec \nu_2],a}(z)&=\varepsilon_{1,1}(
\dfrac{\tau_{[\vec\nu_1-\vec e_{1,a}+\vec e_{1,a'};\vec\nu_2]}(t-[z^{-1}]_{a})}{\tau_{[\vec\nu_1;\vec \nu_2]}(t)},\\
A^{(\rI,b)}_{[\vec\nu_1;\vec \nu_2+\vec e_{2,p_2}-\vec e_{2,b}],a}(z)&=\varepsilon_{2,1}(b, a)z^{
\nu_{1,a}-1}\dfrac{\tau_{[\vec\nu_1+\vec e_{1,p_1}-\vec e_{1,a};\vec\nu_2+\vec e_{2,p_2}-\vec
e_{2,b}]}(t-[z^{-1}]_{a})}{\tau_{[\vec\nu_1+\vec e_{1,p_1};\vec \nu_2+\vec e_{2,p_2}]}(t)},
\end{aligned}
\end{align}
\begin{align}
\begin{aligned}
\bar A^{(\rI,p_2)}_{[\vec\nu_2+\vec e_{2,p_2};\vec \nu_1-\vec
e_{2,p_2}-\vec e_{2,b}]}(t+\overline{[z^{-1}]}_b)}{\tau_{[\vec\nu_1+\vec e_{1,p_1};\vec \nu_2+\vec e_{2,p_2}]}(t)},\\
\bar A^{(\rII,b')}_{[\vec\nu_2+\vec e_{2,b'};\vec \nu_1-\vec e_{1,p_1}],b}&=\varepsilon_{2,2}(b,b')z^{ \nu_{2,b}+
\delta_{b',b}-1} \dfrac{\tau_{[\vec\nu_1;\vec\nu_2-\vec e_{2,b}+\vec e_{2,b'}]}(t+\overline{[z^{-1}]}_{b})}{
\tau_{[\vec\nu_1;\vec \nu_2]}(t)},\\
\bar A^{(\rI,a)}_{[\vec\nu_2+\vec e_{2,p_2};\vec \nu_1-\vec e_{1,a}],b}&=\varepsilon_{2,1}(b, a)z^{ \nu_{2,b}
+\delta_{b,p_2}-1}\dfrac{\tau_{[\vec\nu_1+\vec e_{1,p_1}-\vec e_{1,a};\vec\nu_2+\vec e_{2,p_2}-\vec
e_{2,b}]}(t+\overline{[z^{-1}]}_{b})}{\tau_{[\vec\nu_1+\vec e_{1,p_1};\vec \nu_2+\vec e_{2,p_2}]}(t)}.
\end{aligned}
\end{align}

We now present the $τ$-representation of the Cauchy transforms of the linear forms.
\begin{pro}\label{proposition: tau cauchy}
The Cauchy transforms have the following $\tau$-function representation
\begin{align}\label{tau3}
\bar C^{(l)}_a=z^{-\nu_{1,a}(l-1)-1}\frac{\tau_{+a}^{(l+1)}(t+[z^{-1}]_a)}{\tau^{(l+1)}(t)},\\
\label{tau4}
C^{(l)}_b=z^{-\nu_{2,b}(l-1)-1}\frac{\bar \tau_{+b}^{(l+1)}(t-\overline{[z^{-1}]}_b)}{\tau^{(l)}(t)}.
\end{align}
\end{pro}
\begin{proof}
See Appendix \ref{II}
\end{proof}

We have the representation
\begin{align*}
C^{(\rI,p_2)}_{[\vec\nu_2+\vec e_{2,p_2};\vec \nu_1-\vec
e_{1,p_1}],a}(z)&=\varepsilon_{1,1}(a,p_1)z^{-\nu_{1,a}-1+\delta_{a,p_1}}\frac{\tau_{[\vec \nu_1+\vec
e_{1,a};\vec\nu_2+\vec e_{2,p_2}]}(t+[z^{-1}]_a)}{\tau_{[\vec \nu_1+\vec e_{1,p_1};\vec\nu_2+\vec
C^{(\rII,p_1)}_{[\vec\nu_1;\vec \nu_2],b}(z)&=\varepsilon_{2,2}(b,p_2)z^{-\nu_{2,b}-1}\frac{\tau_{[\vec \nu_1+\vec
e_{1,p_1};\vec\nu_2+\vec e_{2,b}]}(t-\overline{[z^{-1}]}_b)}{\tau_{[\vec \nu_1;\vec\nu_2]}(t)}.
\end{align*}

%Taking into account the correspondence between Miwa shifts and Darboux--Miwa discrete flows \eqref{miwa-darboux} we get the following alternative
%expressions in terms of shifts on a discrete lattice where $\lambda_a(n)=\bar \lambda_b(n)=z$
%\label{las formulas}
%  A_a^{(l)}(z)&=z^{\nu_{1,a}(l)-1}\dfrac{\tau^{(l)}_{-a}(t-[z^{-1}]_a)}{\tau^{(l)}(t)}=\frac{T_a\tau_{-a}^{(l)}(t,s)}{\tau^{(l)}(t,s)},\\  \bar
%  A_b^{(l)}(z)&=z^{\nu_{2,b}(l)-1}\dfrac{\bar \tau_{-b}^{(l)}(t+\overline{[z^{-1}]}_b)}{ \tau^{(l+1)}(t)}=\frac{\bar T_b\bar
%  \tau_{-b}^{(l)}(t,s)}{\tau^{(l+1)}(t,s)},\\
%    \bar C_a^{(l)}(z)&=z^{-\nu_{1,a}(l-1)-1}\frac{\tau_{+a}^{(l+1)}(t+[z^{-1}]_a)}{\tau^{(l+1)}(t)}=
%    \frac{T_a^{-1}\tau_{+a}^{(l+1)}(t,s)}
%    {\tau^{(l)}(t,s)},\\
%      C_b^{(l)}(z)&=z^{-\nu_{2,b}(l-1)-1}\frac{\bar \tau_{+b}^{(l+1)}(t-\overline{[z^{-1}]}_b)}{\tau^{(l)}(t)}=\frac{\bar T_b^{-1}\bar
%      \tau_{-b}^{(l+1)}(t,s)}{\tau^{(l+1)}(t,s)}.
%\subsubsection{Bilinear equation in $\tau$-form and Fay identities}

Finally, we consider the $τ$-function representation of the bilinear equation
\begin{pro}
The $\tau$ functions fulfill the following bilinear relation
\begin{multline}\label{bilinear2?}
\sum_{a=1}^{p_1}\oint_{z=\infty}
z^{\nu_{1,a}(k)-\nu_{1,a}'(l-1)-2}\tau^{(k)}_{\vec n_1,\vec n_2,-a}(t-[z^{-1}]_a)
\tau_{\vec n_1',\vec n_2',+a}^{(l+1)}(t'+[z^{-1}]_a)E_a(z)\d z\\=
\sum_{b=1}^{p_2}\oint_{z=\infty}z^{\nu_{2,b}'(l)-\nu_{2,b}(k-1)-2}\bar \tau_{\vec n_1,\vec
\tau_{\vec n_1',\vec n_2',-b}^{(l)}(t'+\overline{[z^{-1}]}_b)\bar E_b(z)\d z.
\end{multline}
\end{pro}
\begin{proof}
Just consider \eqref{bilinear2} together with \eqref{taumops}, \eqref{taudualmops}, \eqref{tau3} and \eqref{tau4}.
\end{proof}
This bilinear relation can also be written as follows
\begin{multline}\label{bilinear2bis}
\sum_{a=1}^{p_1}\varepsilon_{11}(p_1,a)\varepsilon_{11}'(p_1,a)\oint_{z=\infty}
\tau_{[\vec \nu_1+\vec e_{1,p_1}-\vec e_{1,a};\vec \nu_2]}(t-[z^{-1}]_a)
\tau_{\vec \nu_1'+\vec e_{1,a};\vec \nu_2'+\vec e_{2,p_2}}^{(l+1)}(t'+[z^{-1}]_a)E_a(z)\d z\\=
\sum_{b=1}^{p_2}\varepsilon_{22}(p_2,b)\varepsilon_{22}'(p_2,b)\oint_{z=\infty}z^{\nu_{2,b}'+\delta_{b,p_2}-\nu_{2,b}-2}
\tau_{[\vec \nu_1+\vec e_{1,p_1};\vec \nu_2+\vec e_{2,b}]}(t-\overline{[z^{-1}]}_b)
\tau_{\vec \nu_1';\vec \nu_2'+\vec e_{2,p_2}-\vec e_{2,b}]}(t'+\overline{[z^{-1}]}_b)\bar E_b(z)\d z.
\end{multline}
That with the identification $m^*=ν⃗_1+e⃗_1,p_1$, $n^*=ν⃗_2$, $m=ν⃗_1'$ and $n=ν⃗_2'+e⃗_2,p_2$,  up to signs, is the bilinear relation (41) in [6].

%Following [4] we now derive Fay identities.
%  Given complex numbers
% \{ \lambda_{a,k},\mu_{a,l},\bar\lambda_{b,m},\bar\mu_{b,n}\}\subset\C
%with $a=1,\dots,p_1$,
%$b=1,\dots,p_2$, $k=1,\dots,q_1^a$, $l=1,\dots,r_1^a$, $m=1,\dots,q_2^b$, $n=1,\dots,r_2^b$,   the following Fay %identities are fulfill
%  \Fs=\bar\Fs
%\tau_{\vec n_1,\vec n_2,-a}^{(k)}(t-[\lambda_{a,k}^{-1}]_a)\tau_{\vec n_1',\vec %n_2',+a}^{(l+1)}(t'+[\lambda_{a,k}^{-1}]_a)  \lambda_{a,k}^{s_a+1}
%k}}^{q_1^a}(\lambda_{a,k}^{-1}-\lambda_{a,k'}^{-1})}\\&+\sum_{a=1}^{p_1}\frac{1}{s_a!}\frac{\partial^{s_a}}{\partial %x^{s_a}}\Big(
%\tau^{(k)}_{\vec n_1,\vec n_2,-a}(t-[x]_a)\tau^{(l+1)}_{\vec n_1',\vec %n_2',+a}(t'+[x]_a)\dfrac{\prod_{l=1}^{r_1^a}(x-\mu_{a,l}^{-1})}{\prod_{k=1}^{q_1 ^a}(x-\
%\bar \Fs&:=\sum_{b=1}^{p_2}\sum_{m=1}^{q_2^b}
%\bar\tau_{\vec n_1,\vec n_2,+b}^{(k+1)}(t-\overline{[\bar\lambda_{b,m}^{-1}]}_b)\bar\tau_{\vec n_1',\vec %n_2',-b}^{(l)}(t'+\overline{[\bar\lambda_{b,m}^{-1}]}_b)  \bar\lambda_{b,n}^{\bar
%m}}^{q_2^b}(\bar\lambda_{b,m}^{-1}-\bar\lambda_{b,m'}^{-1})}\\&+\sum_{b=1}^{p_2}\frac{1}{\bar %s_b!}\frac{\partial^{\bar s_b}}{\partial x^{\bar
%\bar \tau^{(k+1)}_{\vec n_1,\vec n_2,+b}(t-\overline{[x]}_b)\bar \tau^{(l)}_{\vec n_1',\vec %n_2',-b}(t'+\overline{[x]}_b)
%  s_a&:=\nu_{1,a}(k)-\nu_{1,a}(l-1)-1+r^a_1-q^a_1,&\bar s_b&:=\nu_{2,b}(l)-\nu_{2,b}(k-1)-1+r^b_2-q^b_2,
%with the understanding that $\partial^s/\partial x^s=0$ for $s<0$.
% The idea is the same as the one developed in [4]  for the derivation of Fay identities. We %put
%  t-t'=\sum_{\substack{a=1,\dots,p_1\\b=1,\dots,q_2}}\sum_{\substack{k=1,\dots,q_1^a\\l=1,\dots,r_1^a\\
%  %m=1,\dots,q_2^b\\n=1,\dots,r_2^b}}\Big([\lambda^{-1}_{a,k}]_a-[\mu_{a,l}^{-1}]_a+\overline{[\bar\lambda^{-1}_{b,m}]}_b
%  -\overline{[\bar\mu_{b,n}^{-1}]}_b\Big)
%and recall that in this case we have
%  %E_a&=\dfrac{\prod_{l=1}^{q_1^a}\Big(1-\dfrac{z}{\mu_{a,l}}\Big)}{\prod_{k=1}^{r_1^a}\Big(1-\dfrac{z}{\lambda_{a,k}}\Big)},&
%  \bar %E_b&=\dfrac{\prod_{n=1}^{q_2^b}\Big(1-\dfrac{z}{\bar\mu_{b,n}}\Big)}{\prod_{m=1}^{r_2^b}\Big(1-\dfrac{z}{\bar\lambda_{b,m}}
%  \Big)}.&
%Observe that these  $\partial/\partial x$-derivative terms come from the pole at

%We now discuss two types of different Fay identities
% \item Let us suppose that
%  l&=k-1,&t&=t^{(0)}+\overline{[\bar \lambda^{-1}]}_b-\overline{[\bar\mu^{-1}]}_b,& t'&=t^{(0)}-[z^{-1}]_a
%so that
%  t-t'&=[z^{-1}]_a+\overline{[\bar \lambda^{-1}]}_b-\overline{[\bar\mu^{-1}]}_b,&
%  s_a&=\nu_{1,a}(k)-\nu_{1,a}(k-2)-2,& \bar s_b&=-1.
%Observe that $s_a\in\{-2,-1,0\}$ and $s_a=0\Leftrightarrow a=a_1(k)=a_1(k-2)$ and therefore $n_{1,a}\geq 2$.

%We have different possibilities
%  \item Assume that  $s_a\neq 0$ so that
%  \begin{multline}\label{fay1}
%\tau_{-a}^{(k)}(t-[z^{-1}]_a+\overline{[\bar \lambda^{-1}]}_b-\overline{[\bar\mu^{-1}]}_b)\tau_{+a}^{(k)}(t)
%\bar\tau_{+b}^{(k+1)}(t-\overline{[\bar\mu^{-1}]}_b)\bar\tau_{-b}^{(k-1)}(t-[z^{-1}]_a+\overline{[\bar\lambda^{-1}]}_b)  %\end{multline}
%\item Let $s_a=0$, then
% \begin{multline}\label{fay2}
% \tau^{(k)}_{-a}(t+\overline{[\bar \lambda^{-1}]}_b-\overline{[\bar\mu^{-1}]}_b-[z^{-1}]_a)\tau^{(k)}_{+a}(t)
%-\tau^{(k)}_{-a}(t+\overline{[\bar \lambda^{-1}]}_b-\overline{[\bar\mu^{-1}]}_b)\tau^{(k)}_{+a}(t-[z^{-1}]_a)
%\bar\tau_{+b}^{(k+1)}(t-\overline{[\bar\mu^{-1}]}_b)\bar\tau_{-b}^{(k-1)}(t-[z^{-1}]_a+\overline{[\bar\lambda^{-1}]}_b)  %\end{multline}
%\item Let us suppose that
%  l&=k,&t&=t^{(0)}+[\lambda_1^{-1}]_a+\overline{[\bar\lambda_1^{-1}]}_b,& %t'&=t^{(0)}-[\lambda_2^{-1}]_{a}-\overline{[\bar\lambda_2^{-1}]}_{b},
%so that
%  %t-t'&=[\lambda_1^{-1}]_a+[\lambda_2^{-1}]_{a}+\overline{[\bar\lambda_1^{-1}]}_b+\overline{[\bar\lambda_2^{-1}]}_{b},\\
%  s_a&=\nu_{1,a}(k)-\nu_{1,a}(k-1)-3,& \bar s_b&=\nu_{2,b}(k)-\nu_{2,b}(k-1)-3.
%Observe that $s_a,\bar s_b\in\{-3,-2\}$.% and $s_a=0\Leftrightarrow a=a_1(k)=a_1(k-2)$ and therefore $n_{1,a}\geq 2$. %and
%and the Fay identity is
%\begin{multline} \label{fay3}
%\frac{ \lambda_1^{s_a+1}}{\lambda_1^{-1}-\lambda_2^{-1}}\\
%\frac{\lambda_2^{s_a+1}}{\lambda_2^{-1}-\lambda_1^{-1}} \\
%\frac{\bar\lambda_1^{\bar s_b+1}}{\bar\lambda_1^{-1}-\bar\lambda_2^{-1}}\\
%\frac{\bar\lambda_2^{\bar s_b+1}}{\bar\lambda_2^{-1}-\bar\lambda_1^{-1}}.

%If in \eqref{fay3} we assume that given a partition $(\vec n_1, \vec n_2)$ and an positive integer $l$ $s_a=\bar %s_b=-2$ for that $l$, then setting $\bar \lambda_1, \bar \lambda_2 \rightarrow \infty$ we obtain

%\begin{multline} \label{fay4}
%\frac{ \lambda_1^{-1}}{\lambda_1^{-1}-\lambda_2^{-1}}
%\frac{\lambda_2^{-1}}{\lambda_2^{-1}-\lambda_1^{-1}} \\

%With the same idea, using \eqref{fay2} we can set $\mu \rightarrow \infty $ and obtain

% \tau^{(l)}_{-a}(t+\overline{[\bar \lambda^{-1}]}_b-[z^{-1}]_a)\tau^{(l)}_{+a}(t)
%-\tau^{(l)}_{-a}(t+\overline{[\bar \lambda^{-1}]}_b)\tau^{(l)}_{+a}(t-[z^{-1}]_a)

%\bar\tau_{+b}^{(l+1)}(t)\bar\tau_{-b}^{(l-1)}(t-[z^{-1}]_a+\overline{[\bar\lambda^{-1}]}_b)  \end{multline}

%Expressions \eqref{fay4} and \eqref{fay5} can be generalized using different partitions for each $\tau$ function. For %instance, given $(\vec\nu_1,\vec\nu_2)$ and $(\vec\nu'_1,\vec\nu'_2)$ with $|\vec\nu_1|=|\vec\nu_2|+1$ and %$|\vec\nu'_1|=|\vec\nu'_2|+1$ and using the simple ladders defined above we have

%\begin{multline} \label{fay6}
%\varepsilon_{1,1}(p_1,a)\varepsilon_{1,1}'(a,p_1)\tau_{[\vec\nu_1+\vec e_{1,p_1}-\vec %e_{1,a};\vec\nu_2]}(t)\tau_{[\vec \nu_1'+\vec e_{1,a};\vec\nu_2'+\vec %e_{2,p_2}]}(t+[\lambda_1^{-1}]_a-[\lambda_2^{-1}]_a)
%\frac{ \lambda_1^{-1}}{\lambda_1^{-1}-\lambda_2^{-1}}+\\
%+\varepsilon_{1,1}(p_1,a)\varepsilon_{1,1}'(a,p_1)\tau_{[\vec\nu_1+\vec e_{1,p_1}-\vec %e_{1,a};\vec\nu_2]}(t+[\lambda_1^{-1}]_a-[\lambda_2^{-1}]_a)\tau_{[\vec \nu_1'+\vec e_{1,a};\vec\nu_2'+\vec %e_{2,p_2}]}(t)
%\frac{\lambda_2^{-1}}{\lambda_2^{-1}-\lambda_1^{-1}} \\
%\varepsilon_{2,2}(p_2,b)\varepsilon_{2,2}'(b,p_2)\tau_{[\vec \nu_1+\vec e_{1,p_1};\vec\nu_2+\vec %e_{2,b}]}(t+[\lambda_1^{-1}]_a)\tau_{[\vec\nu_1';\vec\nu_2'+\vec e_{2,p_2}-\vec e_{2,b}]}(t-[\lambda_2^{-1}]_a),

% \varepsilon_{1,1}(p_1,a)\varepsilon_{1,1}'(a,p_1)\tau_{[\vec\nu_1+\vec e_{1,p_1}-\vec %e_{1,a};\vec\nu_2]}(t+\overline{[\bar \lambda^{-1}]}_b-[z^{-1}]_a)\tau_{[\vec \nu_1'+\vec e_{1,a}-\vec %e_{1,p_1};\vec\nu_2']}(t)-\\
%-\varepsilon_{1,1}(p_1,a)\varepsilon_{1,1}'(a,p_1)\tau_{[\vec\nu_1+\vec e_{1,p_1}-\vec %e_{1,a};\vec\nu_2]}(t+\overline{[\bar \lambda^{-1}]}_b)\tau_{[\vec \nu_1'+\vec e_{1,a}-\vec %e_{1,p_1};\vec\nu_2']}(t-[z^{-1}]_a)
%\tau_{[\vec \nu_1+\vec e_{1,p_1};\vec\nu_2+\vec e_{2,b}]}(t+[\lambda_1^{-1}]_a)\tau_{[\vec\nu_1'-\vec %e_{1,p_1};\vec\nu_2'-\vec e_{2,b}]}(t-[z^{-1}]_a+\overline{[\bar\lambda^{-1}]}_b).  \end{multline}

%Expression \eqref{fay6} holds only for $a=p_1$ and $b=p_2$ and \eqref{fay7} holds for $a=p_1$.

%We want expressions for products of polynomials and Cauchy transforms. For example, if in \eqref{fay4} we set %$a=a_1(l)$ and $b=a_2(l)$ then we find

%A_a^{(l)}(\lambda)\bar C_a^{(l)} (\mu) &=\lambda^{\nu_{1,a}(l)-1} \mu^{-\nu_{1,a}(l-1)-1} %\dfrac{\tau^{(l+1)}(t+[\mu^{-1}]_a) \tau^{(l)}(t-[\lambda^{-1}]_a)}{\tau^{(l)}\tau^{(l+1)}} \\
%& = \frac{\lambda^{\nu_{1,a}(l)} \mu^{-\nu_{1,a}(l-1)}}{\lambda-\mu} \left( %\frac{\tau^{(l+1)}(t+[\mu^{-1}]_a-[\lambda^{-1}]_a)}{\mu \tau^{(l+1)}(t)}
%-\frac{\tau^{(l)}(t+[\mu^{-1}]_a-[\lambda^{-1}]_a)}{\lambda \tau^{(l)}(t)}\right)\\
% & \frac{\lambda^{\nu_{1,a}(l)} \mu^{-\nu_{1,a}(l-1)}}{\lambda-\mu} \left( %\frac{X_a(\mu^{-1}-\lambda^{-1})\tau^{(l+1)}(t)}{\mu \tau^{(l+1)}(t)}
%-\frac{X_a(\mu^{-1}-\lambda^{-1})\tau^{(l)}(t)}{\lambda \tau^{(l)}(t)}\right).

%In \eqref{fay5} we can set $a=a_1(l)$ and $b=a_2(l)$ and consequently we obtain

%A_a^{(l)}(\lambda)\bar A_b^{(l)} (\mu) &=\lambda^{\nu_{1,a}(l)-1} \mu^{\nu_{2,b}(l)-1} %\frac{\tau^{(l)}(t-[\lambda^{-1}]_a) \tau^{(l)}(t+\overline{[\mu^{-1}]}_b)}{\tau^{(l)}(t) \tau^{(l)}(t)} \\
%& = \lambda^{\nu_{1,a}(l)} \mu^{\nu_{2,b}(l)} %\left(\frac{\tau^{(l)}(t+\overline{[\mu^{-1}]}_b-[\lambda^{-1}]_a)}{\tau^{(l)}(t)}
%-\frac{\tau^{(l-1)}(t-[\lambda^{-1}]_a+\overline{[\mu^{-1}]}_b)}{\tau^{(l+1)}(t)}  \right) \\
%& = \lambda^{\nu_{1,a}(l)} \mu^{\nu_{2,b}(l)} \left(\frac{X_a(\lambda^{-1})\bar X _b(\mu^{-1})\tau^{(l)}(t)
%}{\tau^{(l)}(t)}-\frac{X_a(\lambda^{-1})\bar X _b(\mu^{-1})\tau^{(l-1)}(t)}{\tau^{(l+1)}(t)}  \right),
%where $X_a \bar X_b$ are vertex operators for Miwa shifts.

%We want to compute \eqref{las formulas}
%  A^{(l)}_{a_1(l)}(\lambda,t)\bar %A^{(l)}_{a_2(l)}(\mu,t)=\lambda^{k_1(l+1)-1}\mu^{k_2(l+1)-1}\dfrac{\tau^{(l)}(t-[\lambda^{-1}]_{a_1(l)})
%  \tau^{(l)}(t+\overline{[\mu^{-1}]}_{a_2(l)})}{\tau^{(l)}(t) \tau^{(l+1)}(t)}

\section*{Acknowledgements}
MM thanks economical support from the Spanish Ministerio de Ciencia e Innovaci\'{o}n, research
project FIS2008-00200 and UF thanks economical support from the SFRH / BPD / 62947 / 2009, Funda\c{c}\~{a}o para a Ci\^{e}ncia e a Tecnologia of Portugal,  PT2009-0031, Acciones Integradas Portugal, Ministerio de Ciencia e Innovaci\'{o}n de Espa\~{n}a, and
MTM2009-12740-C03-01. Ministerio de Ciencia e Innovaci\'{o}n de Espa\~{n}a. MM reckons different and clarifying discussions with Dr. Mattia Caffasso, Prof. Pierre van Moerbeke,  Prof. Luis Mart\'{\i}nez Alonso and Prof. David G\'{o}mez-Ullate. The authors of this paper are in debt with Prof. Guillermo L\'{o}pez Lagomasino who carefully read the manuscript and whose suggestions improved the paper indeed.

\begin{appendices}

\section{Proofs}\label{II}

\paragraph{Proof Proposition \ref{proposition: expressions for linear forms}}
The orthogonality relations can be recast into two alternative forms
\begin{align}\label{S.1}
\begin{pmatrix}
\end{pmatrix} \begin{pmatrix}
\vdots&\vdots&&\vdots\\
\end{pmatrix}&=- \begin{pmatrix}
\end{pmatrix},& l&\geq 1,\\
\label{S.2}
\begin{pmatrix}
\end{pmatrix} \begin{pmatrix}
\vdots&\vdots&&\vdots\\
\end{pmatrix}&=\underbrace{\begin{pmatrix}
0&0&\cdots&0& \bar S_{l,l}
\end{pmatrix}}_{\text{$l+1$ components}},& l&\geq 0.
\end{align}
From \eqref{linear forms S} we get
\begin{align}
Q^{(l)}&=\sum_{k=0}^l S_{l,k}\xi_{1}^{(k)}\notag
\\&=\xi_{1}^{(l)}-\begin{pmatrix}
\end{pmatrix}(g^{[l]})^{-1}\xi_1^{[l]}& \text{use \eqref{S.1}}\label{S.1.1}\\
&=\bar S_{l,l}\begin{pmatrix}
0 &0 &\cdots &0& 1
\end{pmatrix}
(g^{[l+1]})^{-1}\xi^{[l+1]}_1. &\text{use \eqref{S.2}}\label{S.1.2}
\end{align}
Cramer's method solves \eqref{S.1} as follows
\begin{align}\label{cofactor}
g^{[l]}}\sum_{j=0}^{l-1}g_{l,j}(-1)^{i+j+1}M^{(l)}_{i,j}=\frac{(-1)^{i+l}M^{(l+1)}_{i,l}}{\det g^{[l]}},
\end{align}
where $M_{i,j}^{(l)}$ is the $(i,j)$-minor of the truncated moment matrix $g^{[l]}$ defined in \eqref{gl}.
\begin{align*}
Q^{(l)}&=%\xi_1^{(l)}+\frac{1}{\det g^{(l)}}\sum_{i,j=0}^{l-1}g_{l,j}(-1)^{i+j+1}M^{(l)}_{i,j}
\frac{1}{\det g^{[l]}}\sum_{i=0}^l (-1)^{i+l} M_{i,l}^{(l+1)}  \xi_1^{(i)}\\&=
\frac{1}{\det g^{[l]}}\det
\left(\begin{BMAT}{cccc|c}{cccc|c}
\vdots &\vdots&            &\vdots&\vdots\\
\end{BMAT}\right),& l&\geq 1.
\end{align*}

The orthogonality relations for the dual system can be written also in two alternative forms
\begin{align}\label{dS.1}
\begin{pmatrix}
\vdots&\vdots&&\vdots\\
\end{pmatrix}\begin{pmatrix}
\bar S'_{0,l}\\\bar S'_{1,l}\\\vdots\\ \bar S_{l-1,l}
\end{pmatrix}&=-(\bar S_{l,l})^{-1}\begin{pmatrix}
\end{pmatrix},& l&\geq 1,\\ \begin{pmatrix}
\vdots&\vdots&&\vdots\\
\end{pmatrix}
\begin{pmatrix}
\bar S'_{0, l}\\
\bar S'_{1, l} \\
\vdots \\
\bar S'_{l, l}
\end{pmatrix}
\begin{pmatrix}
0 \\
0 \\
\vdots \\0\\
\end{pmatrix},& l&\geq 0.\label{dS.2}
\end{align}
As before,  \eqref{linear forms S} leads to the following expressions for the dual linear forms
\begin{align}
\bar  Q^{(l)}&=\sum_{k=0}^l \bar S'_{k,l}\xi_{2}^{(k)}
\notag\\&=(\bar S_{l,l})^{-1}\Big(\xi_{2}^{(l)}-(\xi_2^{[l]})^\top(g^{[l]})^{-1}\begin{pmatrix}
\end{pmatrix}\Big)& \text{use \eqref{dS.1}}\label{dS1.1}\\
0 \\0\\\vdots \\0\\ 1
\end{pmatrix}. &\text{use \eqref{dS.2}}\label{dS1.2}
\end{align}

From \eqref{dS.2} we obtain
\begin{align}\label{s'}
\bar S'_{j,l}=\big({g^{[l+1]}}^{-1}\big)_{j,l}=\frac{(-1)^{l+j}M^{(l+1)}_{l,j}}{\det(g^{[l+1]})}, \quad
\end{align}
and consequently
\begin{align*}
\bar Q ^{(l)}&=\sum_{j=0}^{l} \bar S'_{j,l} \xi_2^{(j)}=\frac{1}{\det
\\&=\frac{1}{\det g^{[l+1]}}\det
\left(\begin{BMAT}{cccc|c}{cccc|c}
g_{0,0}&g_{0,1}&\cdots&g_{0,l-1}& g_{0,l}\\
g_{1,0}&g_{1,1}&\cdots&g_{1,l-1}& g_{1,l} \\
\vdots &\vdots&        &\vdots&\vdots\\
g_{l-1,0}&g_{l-1,1}&\cdots&g_{l-1,l-1}&g _{l-1,l}\\
\xi_2^{(0)}& \xi_2^{(1)} & \cdots & \xi_2^{(l-1)} & \xi_2^{(l)}
\end{BMAT}\right),& l\geq 0.
\end{align*}

\paragraph{Proof Proposition \ref{proposition: expressions for ct}}
%If in \eqref{ct} we introduce \eqref{cofactor} we get
We have
\begin{align*}
C_b^{(l)}=\frac{1}{\det g^{[l]}}
\sum_{k=0}^l(-1)^{k+l}M^{(l+1)}_{k,l}\sum_{k_2=\nu_{2,b}(l-1)}^{\infty}z^{-k_2-1}\int
x^{k_1(k)}w_{1,a_1(k)}(x)w_{2,b}(x)x^{k_2}\d \mu(x),
\end{align*}
which according to  \eqref{Gamma}  recasts into
\begin{align}\notag
C_b^{(l)}&=\frac{1}{\det g^{[l]}}
\sum_{k=0}^l(-1)^{k+l}M^{(l+1)}_{k,l}\bar\Gamma_{k,b}^{(l)},\\
&=\frac{1}{\det g^{[l]}}\det
\left(\begin{BMAT}{cccc|c}{cccc|c}
\vdots &\vdots&            &\vdots&\vdots\\
\end{BMAT}\right),& l&\geq 1.&&\label{ctdet}
\end{align}
We also obtain
\begin{align*}
\bar C^{(l)}_a=\frac{1}{\det g^{[l+1]}}\sum_{k=0}^l
(-1)^{l+k}M^{(l+1)}_{l,k}\sum_{k_1=\nu_{1,a}(l-1)}^{\infty}z^{-k_1-1}\int  x^{k_1}w_{1,a}(x)w_{2,a_2(k)}x^{k_2(k)}\d
\mu(x),
\end{align*}
which can be written as \eqref{Gamma}
\begin{align}\notag
\bar C_a^{(l)}&=\frac{1}{\det g^{[l+1]}}
\sum_{k=0}^l(-1)^{k+l}M^{(l+1)}_{l,k}\Gamma_k^{(l)},\\
&=\frac{1}{\det g^{[l+1]}}\det
\left(\begin{BMAT}{cccc|c}{cccc|c}
\vdots &\vdots&            &\vdots&\vdots\\
\Gamma^{(l)}_{0,a}& \Gamma^{(l)}_{1,a}&\cdots& \Gamma^{(l)}_{l-1,a}& \Gamma^{(l)}_{l,a}
\end{BMAT}\right),& l&\geq 1.&&\label{ctdet-dual}
\end{align}

\paragraph{Proof Proposition \ref{proposition: CD type relations}}
From \eqref{mop-s} and \eqref{cauchy-S} we deduce that
\begin{align*}
(\bar\Cs_a(z))^\top\A_{a'}(z')&=(\chi_{1,a}^*(z))^\top\chi_{1,a'}(z')=\frac{\delta_{a,a'}}{z-z'},& |z'|<|z|,\\
(\Cs_b(z))^\top\bar\A_{a'}(z')&=(\chi_{2,b}^*(z))^\top\chi_{2,b'}(z')=\frac{\delta_{b,b'}}{z-z'},& |z'|<|z|,\\
(\bar\Cs_a(z))^\top\Cs_{b}(z')&=(\chi_{1,a}^*(z))^\top g\chi^*_{2,b}(z').
\end{align*}
The two first relations imply the corresponding equations in the Proposition. For the third we observe that from
\eqref{compact.g} we get
\begin{align*}
(\chi_{1,a}^*(z))^\top g\chi^*_{2,b'}(z')&=\int     (\chi_{1,a}^*(z))^\top
\xi_1(x)(\xi_2(x))^\top\chi^*_{2,b'}(z')\d\mu(x)\\
&=\int     (\chi_{1,a}^*(z))^\top
\chi_{1,a}(x)(\chi_{2,b}(x))^\top\chi^*_{2,b'}(z')w_{1,a}(x)w_{2,b}(x)\d\mu(x)\\
&=\int     \frac{1}{(z-x)(z'-x)}w_{1,a}(x)w_{2,b}(x)\d\mu(x)\\
&=-\frac{1}{z-z'}\int     \Big(\frac{1}{z-x}-\frac{1}{z'-x}\Big)w_{1,a}(x)w_{2,b}(x)\d\mu(x).
\end{align*}

\paragraph{Proof Proposition \ref{associated forms as mop}}
Using Definition \ref{associated linear forms} for the linear forms $Q^{(l)}_{+a}$ and multiplying by
$(\xi_2^{[l]}(x))^{\top}$ we have
\begin{align*}
Q_{+a}^{(l)}(x) (\xi_2^{[l]}(x))^{\top}=\xi_1^{(l_{+a})}(x)(\xi^{[l]}_2(x))^{\top}-
\begin{pmatrix}
g_{l_{+a},0}&g_{l_{+a},1}&\cdots& g_{l_{+a},l-1}
\end{pmatrix}(g^{[l]})^{-1}\xi_1^{[l]}(x)(\xi^{[l]}_2(x))^{\top},
\end{align*}
integrating both sides we get
\begin{align*}
\int  Q_{+a}^{(l)}(x) (\xi^{[l]}_2(x))^{\top} \d \mu(x)
& =\int  \xi_1^{(l_{+a})}(x)(\xi^{[l]}_2(x))^{\top} \d \mu(x)\\&\quad\quad\quad\quad-
\begin{pmatrix}
g_{l_{+a},0}&\cdots& g_{l_{+a},l-1}
\end{pmatrix}(g^{[l]})^{-1}\int  \xi_1^{[l]}(x)(\xi^{[l]}_2(x))^{\top} \d \mu(x)  \\
& =\int  \xi_1^{(l_{+a})}(x)(\xi^{[l]}_2(x))^{\top} \d \mu(x)-
\begin{pmatrix}
g_{l_{+a},0}& g_{l_{+a},1} &\cdots& g_{l_{+a},l-1}
\end{pmatrix}(g^{[l]})^{-1}g^{[l]} \\
& = \begin{pmatrix}
g_{l_{+a},0}&g_{l_{+a},1}&\cdots& g_{l_{+a},l-1}
\end{pmatrix}-
\begin{pmatrix}
g_{l_{+a},0}&g_{l_{+a},1}&\cdots& g_{l_{+a},l-1}
\end{pmatrix}\\
& =0,
\end{align*}
that written componentwise gives the following orthogonality relations
\begin{align*}
\int_{\R} Q^{(l)}_{+a}(x) w_{2,a_2(k)}(x) x^{k_2(k)} \d \mu(x),  \quad k=0,\dots,l-1,
\end{align*}
or equivalently
\begin{align*}
\int_{\R} Q^{(l)}_{+a}(x) w_{2,b}(x) x^k \d \mu(x)&=0,  &0&\leq k\leq \nu_{2,b}(l-1)-1,&b&=1,\dots,p_2.
\end{align*}
Notice that, $A^{(l)}_{+a,a}$ is monic and $\deg A^{(l)}_{+a,a}(x)=k_1(l_{+a})$ but $A^{(l)}_{+a,a'}$ with $a \neq a'$
satisfy $\deg A^{(l)}_{+a,a'} \leq k_1((l-1)_{-a'})$.
This means that the set of polynomials $A^{(l)}_{+a,a'}(x)$ have degrees determined by $\vec \nu_1(l-1)+\vec e_{1,a}$
and a normalization with respect to the $a$-th component of type  II; i.e,
$Q_{+a}^{(l)}=Q^{(\rII,a)}_{[\vec\nu_1(l-1)+\vec e_{1,a};\vec \nu_2(l-1)]}$.

In a similar way, the associated linear forms $Q^{(l)}_{-b}(x)$ solve a mixed multiple orthogonal problem that can be
obtained as follows. From Definition \ref{associated linear forms} and multiplying by $(\xi_2^{[l]}(x))^{\top}$ we
\begin{align*}
\end{align*}
integrating both sides
\begin{align*}
\int_{\R} Q_{-b}^{(l)}(x)(\xi_2^{[l+1]}(x))^{\top} \d \mu(x)=e_{\bar
l_{-b}}^\top(g^{[l+1]})^{-1}\int_{\R}\xi_1^{[l+1]}(x)(\xi_2^{[l+1]}(x))^{\top} \d \mu(x) = e_{\bar l_{-b}}^{\top},
\end{align*}
and written componentwise
\begin{align*}
\int_{\R} Q_{-b}^{(l)}(x) x^{k_2(k)}w_{2,a_2(k)}(x) \d \mu(x)= \delta_{k,\bar l_{-b}}, \quad k=0, \cdots, l,
\end{align*}
that is equivalent to
\begin{align*}
\int_{\R} Q^{(l)}_{-b}(x) w_{2,b}(x) x^k \d \mu(x)&=\delta_{k,\bar l_{-b}},  &0&\leq k\leq
\nu_{2,b}(l)-1,&b&=1,\dots,p_2.
\end{align*}
Hence, the set $A^{(l)}_{-b,a'}$ is a type I normalized to the $b$-th component solution for a mixed multiple
orthogonality problem; the degrees satisfy $\deg A^{(l)}_{-b,a'} \leq l_{-a'}$. Moreover, the fact that the last
orthogonality condition in the $b$-th component is missing gives the identification
$Q_{-b}^{(l)}=Q^{(\rI,b)}_{[\vec\nu_1(l);\vec \nu_2(l)-\vec e_{2,b}]}$.

Using Definition $\ref{associated linear forms}$ and multiplying by $\xi^{[l]}_1(x)$ we have
\begin{align*}
\bar Q_{+b}^{(l)}(x)\xi^{[l]}_1(x)&=\Big(\xi_1^{[l]}(x)\xi_2^{(\bar
\begin{pmatrix}
g_{0,\bar l_{+b}}\\g_{1,\bar l_{+b}}\\\vdots\\ g_{l-1,\bar l_{+b}}
\end{pmatrix}\Big),
\end{align*}
and integrating both sides
\begin{align*}
\int _\R\bar Q_{+b}^{(l)}(x)\xi_1(x)^{[l]} \d \mu(x) &= \Big(\int  \xi^{[l]}_1(x)\xi_2^{(\bar l_{+b})}(x)\d
\mu(x)-\int  \xi^{[l]}_1(x)(\xi_{2}^{[l]}(x))^\top \d \mu(x)(g^{[l]})^{-1}
\begin{pmatrix}
g_{0,\bar l_{+b}}\\g_{1,\bar l_{+b}}\\\vdots\\ g_{l-1,\bar l_{+b}}
\end{pmatrix}\Big) \\
& =\Big(\int  \xi^{[l]}_1(x)\xi_2^{(\bar l_{+b})}(x)\d \mu(x)-(g^{[l]})(g^{[l]})^{-1}
\begin{pmatrix}
g_{0,\bar l_{+b}}\\g_{1,\bar l_{+b}}\\\vdots\\ g_{l-1,\bar l_{+b}}
\end{pmatrix}\Big)\\
& =  \begin{pmatrix}
g_{0,\bar l_{+b}}\\g_{1,\bar l_{+b}}\\\vdots\\ g_{l-1,\bar l_{+b}}
\end{pmatrix}
\begin{pmatrix}
g_{0,\bar l_{+b}}\\g_{1,\bar l_{+b}}\\\vdots\\ g_{l-1,\bar l_{+b}}
\end{pmatrix}\\
\end{align*}
that componentwise leads to the following orthogonality relations
\begin{align*}
\int_{\R} \bar Q^{(l)}_{+b}(x) w_{1,a_1(k)}(x) x^{k_1(k)} \d \mu(x)=0,  \quad k=0,\dots,l-1,
\end{align*}
or alternatively
\begin{align*}
\int  \bar Q_{+b}^{(l)}(x)w_{1,a}(x)x^k\d \mu(x)&=0,&0&\leq k\leq \nu_{1,a}(l-1)-1,&a&=1,\dots,p_1.
\end{align*}
Notice that, $\bar A^{(l)}_{+b,b}$ is monic and $\deg \bar A^{(l)}_{+b,b}=k_2(\bar l_{+b})$ but $\bar A^{(l)}_{+b,b'}$
with $b \neq b'$ satisfy $\deg A^{(l)}_{+b,b'} \leq k_1(\overline{(l-1)}_{-b'})$.
This means that the polynomials $\bar A^{(l)}_{+b,b'}$ have degrees determined by $\vec \nu_2(l-1)+\vec e_{2,b}$ and a
normalization with respect to the $b$-th component of type  II; i.e,   $\bar Q_{+b}^{(l)}=\bar
Q^{(\rII,b)}_{[\vec\nu_2(l-1)+\vec e_{2,b};\vec \nu_1(l-1)]}$.

Finally, we obtain the orthogonality relations for the linear forms $\bar Q^{(l)}_{-a}(x)$. From the definition we
\begin{align*}
\xi_{1}^{[l+1]}(x)\bar Q^{(l)}_{-a}(x)=\xi_{1}^{[l+1]}(x)(\xi_{2}^{[l+1]}(x))^\top(g^{[l+1]})^{-1}e_{l_{-a}},
\end{align*}
and integrating both sides
\begin{align*}
\int_{\R} \xi_1^{[l+1]}(x) \bar Q_{-a}^{(l)}(x) \d \mu(x) =\Big(\int_{\R}\xi_{1}^{[l+1]}(x)(\xi_{2}^{[l+1]}(x))^\top
\d \mu(x) \Big) (g^{[l+1]})^{-1}e_{l_{-a}}=e_{l_{-a}},
\end{align*}
and componentwise that means
\begin{align*}
\int_{\R} \bar Q_{-a}^{(l)}(x) x^{k_1(k)}w_{1,a_1(k)} (x) \d \mu(x)= \delta_{k,l_{-a}}, \quad k=0, \cdots, l,
\end{align*}
or equivalently
\begin{align*}
\int_{\R} \bar Q_{-a}^{(l)}(x) x^{k}w_{1,a} (x) \d \mu(x)&=\delta_{k,l_{-a}}, &0&\leq k\leq
\nu_{1,a}(l)-1,&a&=1,\dots,p_1,
\end{align*}
so the set $\bar A^{(l)}_{-a,b'}$ is a type I normalized to the $a$-th component solution for a mixed multiple
orthogonality problem. The degrees satisfy $\deg A^{(l)}_{-a,b'} \leq l_{-b'}$; and therefore we conclude that $\bar
Q_{-b}^{(l)}=Q^{(\rI,b)}_{[\vec\nu_2(l);\vec \nu_1(l)-\vec e_{1,a}]}$.

\paragraph{Proof of Proposition \ref{proposition: symmetry}}
Taking $+$ and $-$ parts in \eqref{bigradedj} we obtain
\begin{align*}
\sum_{a=1}^{p_1}L_a^j=\sum_{a=1}^{p_1}(L_a^j)_++\sum_{a=1}^{p_1}(L_a^j)_-=\sum_{a=1}^{p_1}(L_a^j)_++\sum_{b=1}^{p_2}(\bar
L_b^j)_-=\sum_{a=1}^{p_1}B_{j,a}+\sum_{b=1}^{p_2} \bar B_{j,b}=\sum_{b=1}^{p_2}\bar L_b^j.
\end{align*}

Using Lax equations and observing that $L_aL_{a'}=L_{a'}L_a$ and $\bar L_b\bar L_{b'}=\bar L_{b'}\bar L_{b}$ we have
the following symmetries for the Lax operators
\begin{align*}
\big(\sum_{a=1}^{p_1}\frac{\partial}{\partial t_{j,a}}+\sum_{b=1}^{p_2}\frac{\partial}{\partial \bar
t_{j,b}}\big)L_{a'}=\Big[\sum_{a=1}^{p_1}B_{j,a}+\sum_{b=1}^{p_2} \bar
B_{j,b},L_{a'}\Big]=\sum_{a=1}^{p_1}[L_a^j,L_{a'}]=0, \\
\big(\sum_{a=1}^{p_1}\frac{\partial}{\partial t_{j,a}}+\sum_{b=1}^{p_2}\frac{\partial}{\partial \bar
L_{b'}=\Big[\sum_{a=1}^{p_1}B_{j,a}+\sum_{b=1}^{p_2} \bar B_{j,b},\bar L_{b'}\Big]=\sum_{b=1}^{p_2}[\bar L_b^j,\bar L_{b'}]=0.
\end{align*}
From \eqref{flow.mop} we conclude that the multiple orthogonal polynomials and their duals are also invariant
\begin{align*}
\big(\sum_{a=1}^{p_1}\frac{\partial}{\partial t_{j,a}}+\sum_{b=1}^{p_2}\frac{\partial}{\partial \bar
t_{j,b}}\big)\A_{a'}&=\Big(\sum_{a=1}^{p_1}B_{j,a}+\sum_{b=1}^{p_2} \bar
\big(\sum_{a=1}^{p_1}\frac{\partial}{\partial t_{j,a}}+\sum_{b=1}^{p_2}\frac{\partial}{\partial \bar
-\Big(\sum_{a=1}^{p_1}B_{j,a}+\sum_{b=1}^{p_2} \bar
\end{align*}

\paragraph{Proof of Proposition \ref{proposition: bilinear fubini-cauchy}}
We just follow the following chain of identities
\begin{align*}
W_{\vec n_1,\vec n_2}(t,s)\pi_{\vec n_1',\vec n_1}^\top  (W_{\vec n_1',\vec n_2'}(t',s'))^{-1}&=
W_{\vec n_1,\vec n_2}(t,s)\pi_{\vec n_1',\vec n_1}^\top g_{\vec n_1',\vec n_2'}(\bar W_{\vec n_1',\vec
n_2'}(t',s'))^{-1} \\
&    \hspace*{-4cm}=S_{\vec n_1,\vec n_2}(t,s)    W_{0,\vec n_1}(t,s)\pi_{\vec n_1',\vec n_1}^\top\Big(\int
\xi_{\vec n_1'}(x)\xi_{\vec n_2'}^\top(x)\d \mu(x)\Big)(W_{0,\vec n_2'}(t',s'))^{-1}
(\bar S_{\vec n_1',\vec n_2'}(t',s') )^{-1} \\
&    \hspace*{-4cm}=S_{\vec n_1,\vec n_2}(t,s)    W_{0,\vec n_1}(t,s)\Big(\int  \xi_{\vec n_1}(x)\xi_{\vec
n_2'}^\top(x)\d \mu(x)\Big)(W_{0,\vec n_2'}(t',s'))^{-1}
(\bar S_{\vec n_1',\vec n_2'}(t',s') )^{-1} \\
&    \hspace*{-4cm}=S_{\vec n_1,\vec n_2}(t,s)\Big(\int  \xi_{\vec n_1}(x,t,s)\xi_{\vec n_2'}^\top(x,t',s')\d
\mu(x)\Big)
(\bar S_{\vec n_1',\vec n_2'}(t',s') )^{-1} \\
&    \hspace*{-4cm}=\int  (S_{\vec n_1,\vec n_2}(t,s)(\xi_{\vec n_1}(x,t,s))(\bar S_{\vec n_1',\vec
n_2'}^\top(t',s')) ^{-1} \xi_{\vec n_2'}(x,t',s') )^\top \d \mu(x)\\
&\hspace*{-4cm}=\int_{\R} Q_{\vec n_1,\vec n_2}(x,t,s)\bar Q_{\vec n_1',\vec n_2'}^\top(x,t',s')\d \mu(x),
\end{align*}
where $\xi_{\vec n_1}(x,t,s)$ and $\xi_{\vec n_2'}(x,t',s')$ represent the vectors of weighted monomials but
with evolved weights.

\paragraph{Proof of Proposition \ref{taumopis}}

To find the $\tau$-function of the multiple orthogonal representation we first need two lemmas
\begin{lemma}\label{lemarows}
Let $R^{(j)}$ be the $j$-th row of $\tau^{(l)}(t)$ and $R^{(j)}_z$ the $j$-th row of $\tau^{(l)}(t-[z^{-1}]_a)$, then
\begin{equation}\label{eqrows}
\end{equation}
where $j'=j+1$ if $r_1(j)<n_{1,a}-1$, but $j'=j+(|\vec n_1|-n_{1,a})+1$ if $r_1(j)=n_{1,a}-1$. This is also valid for
$\tau^{(l)}_{-a}$ ,  $\tau^{(l)}_{+a,-a'}$ and for $\tau^{(l)}_{-b,-a'}$. \\

Let now be $C^{(j)}$ the j-th column of $\bar \tau^{(l)}$ and $C^{(j)}_z$ the $j$-th column of $\bar
\tau^{(l)}(t+\overline{[z^{-1}]}_b)$, then
\begin{equation}\label{eqcolumns}
\end{equation}
where  $j'=j+1$ if $r_2(l)<n_{2,b}-1$ but $j'=j+(|\vec n_2|-n_{2,b})+1$ if $r_2(j)=n_{2,b}-1$. This is also valid for
$\bar \tau^{(l)}_{-b}$, $\bar \tau^{(l)}_{+b,-b'}$ and for $\bar \tau^{(l)}_{-a,-b'}.$
\end{lemma}
\begin{proof}
It follows directly from \eqref{vertex weight} and \eqref{vertex
\end{proof}

Let us  recall the skew multi-linear character of  determinants and the consequent formulation in terms  of wedge
products of covectors. Observe that
\begin{lemma} \label{lemacovectors}
Given a set of covectors $\{r_1,\dots,r_n\}$ it can be  shown that
\begin{align}\label{covectors}
\bigwedge_{j=1}^n(zr_j-r_{j+1})=\sum_{j=1}^{n+1}(-1)^{n+1-j}z^{j-1} \;r_1\wedge r_2\wedge\dots\wedge \hat
r_j\wedge\dots \wedge r_{n+1},
\end{align}
where the notation $ \hat r_j$ means that we have erased the covector $r_j$ in the wedge product
$r_1\wedge\dots\wedge r_{n+1}$.
\end{lemma}
\begin{proof}
It can be done directly by induction.
\end{proof}

The proof of Proposition \ref{taumopis} relies on Lemma \ref{lemarows}, Lemma \ref{lemacovectors}, Corollary
\ref{cormops} and Proposition \ref{det-associated}. First let's focus on \eqref{taumops}; it is clear that
$z^{\nu_{1,a}(l)-1} \tau_{-a}^{(l)}(t-[z^{-1}]_a)$ expands in $z$ according to \eqref{eqrows} for $\tau ^{(l)}_{-a}$
and to \eqref{covectors}. Now $n=k_1(l_{-a})$ and the covectors  $r_j$ should be taken equal to those rows $R^{(j)}$
with $a_1(j)=a$. Observe that there are only $k_1( l_{-a})(=\nu_{1,a}(l)-1)$ rows which are non-trivially transformed.
In this form we get the identification of
\eqref{detmops} with  \eqref{taumops}, where the terms corresponding to the wedge with one covector deleted
corresponds to the minors
$M_{j,l}^{[l+1]}$. Now, looking to \eqref{taumops+} and \eqref{taumops-} we expand again in $z$ and use the same
technique based on \eqref{eqrows} for $\tau ^{(l)}_{+a,-a'}$ and $\tau ^{(l)}_{-b,-a}$ and \eqref{covectors}.
These allow to link \eqref{detmops+} to \eqref{taumops+} and \eqref{detmops-} to \eqref{taumops-}.

To prove \eqref{taudualmops} we proceed similarly. Looking at \eqref{eqcolumns} for $\bar \tau^{(l)}_{-b}$ observe
that there are only $k_2(\bar l_{-b})(=\nu_{2,b}(l)-1)$ columns which are non-trivially transformed. Now, recalling
\eqref{dualmops.2} and using \eqref{covectors} but with $r_j$ being the columns $C^{(j)}$, such that $a_2(j)=b$, and
$n=k_2(\bar l_{-b})$, we get the desired result. Finally for $\eqref{taudualmops+}$ and $\eqref{taudualmops-}$ we
expand around $z$ to see the equivalence between \eqref{detdualmops+} and \eqref{taudualmops+} and the equivalence
between \eqref{detdualmops-} and \eqref{taudualmops-}.

\paragraph{Proof of Proposition \ref{proposition: tau cauchy}}
We need the following two lemmas:
\begin{lemma}
Let $R^{(j)}$ be the $j$-th row of $g_{+a}^{[l+1]}$ and $R^{(j)}_z$ the $j$-th row of
$g_{+a}^{[l+1]}(t+[z^{-1}]_a)$,   we get
\begin{equation}
R^{(j)}_z=R^{(j)}+\delta_{a_1(j),a}\sum_{j'=1}^\infty z^{-k_1(j')}R^{(j+j')}\delta_{a_1(j+j'),a}.
\end{equation}
Let  $C^{(j)}$ be the $j$-th column of $\bar g_{+b}^{[l+1]}$ and
$C^{(j)}_z$ the $j$-th column of $\bar g_{+b}^{[l+1]}(t-\overline{[z^{-1}]}_b)$, then
\eqref{vertex dualweight} gives
\begin{equation}
C^{(j)}_z=C^{(j)}+\delta_{a_2(j),b}\sum_{j'=1}^\infty z^{-k_2(j')}C^{(j+j')}\delta_{a_2(j+j'),b}.
\end{equation}
\end{lemma}
\begin{proof}
For the first equality insert the expansion
\begin{align*}
\Big(1-\frac{x}{z}\Big)^{-1}=\sum_{k=0}^\infty\frac{x^k}{z^k}
\end{align*}
into \eqref{vertex dualweight}. The other equation is proven similarly.
\end{proof}

\begin{lemma}
The following  identity
\begin{align}\label{covectors.2}
\bigwedge_{j=1}^n
\Big(\sum_{i=0}^{\infty}r_{j+i}z^{-i}\Big)=r_1\wedge\dots\wedge r_{n-1}\wedge \Big(\sum_{i=0}^\infty
\end{align}
\end{lemma}
\begin{proof}
Use induction in $n$.
\end{proof}
Finally Proposition \ref{proposition: tau cauchy} is proven using  \eqref{ctdet-dual} and  \eqref{ctdet}.

\section{Discrete flows associated with binary Darboux transformations}\label{III}

When the supports of the measures are not bounded from below, (in which case  the new ``weights" \eqref{d.evolved} do
not have in general a definite sign and therefore
should not be considered as such), there is an alternative form of constructing discrete flows which preserve the
positiveness/negativeness of the measures.  The
construction is based in the previous one, but now the shift is the composition of two consecutive shifts associated
with the pair $\lambda_a(n)$ and $\lambda_a(n+1)$, being complex numbers conjugate to each other; i.e., we consider
\begin{definition}
We define a deformed moment matrix
\begin{align}
g(s)=D'_0(s)g(\bar D'_0(s))^{-1}.
\end{align}
\begin{align}
\begin{aligned}
D'_0&:=\sum_{a=1}^{p_1}D'_{0,a}, \\ D'_{0,a}&:=\begin{cases}
\prod_{n=1}^{s_a}\big(|\lambda_a(n)|^2\Pi_{1,a}-
2\operatorname{Re}(\lambda_a(n))\Lambda_{1,a}+\Lambda^{2}_{1,a}\big), &s_a>0,\\
\Pi_{1,a},& s_a=0,\\
\prod_{n=1}^{|s_a|}\big(|\lambda_a(-n)|^2\Pi_{1,a}-
2\operatorname{Re}(\lambda_a(-n))\Lambda_{1,a}+\Lambda^{2}_{1,a}\big)^{-1},& s_a<0,
\end{cases}
\end{aligned}
\end{align}
%\Big(\delta_{s_a,0}\Pi_{1,a}+(1-\delta_{s_a,0})\prod_{n=1}^{|s_a|}(\Lambda_{1,a}-\lambda_a(n\sgn s_a)\Pi_{1,a})^{\sgn s_a}\Big),\\
\begin{align}
\begin{aligned}
(\bar   D_0')^{-1}&:=   \sum_{b=1}^{p_2} \big((\bar   D_0(s)')^{-1}\big)_b, \\ \big((\bar   D_0(s)')^{-1}\big)_b
& :=\begin{cases}
\prod_{n=1}^{\bar s_b}(|\bar\lambda_b(n)|^2\Pi_{2,b}-2\operatorname{Re}(\bar\lambda_b(n))\Lambda_{2,b}^\top)
+(\Lambda_{2,b}^\top)^{2} )\big) &\bar s_b>0,\\
\Pi_{2,b},&\bar s_b=0,\\
(\prod_{n=1}^{\bar s_b}\big(|\bar\lambda_b(-n)|^2\Pi_{2,b}-2\operatorname{Re}(\bar\lambda_b(-n))\Lambda_{2,b}^\top)
+(\Lambda_{2,b}^\top)^{2} \big)^{-1},&\bar s_b<0.
\end{cases}
%   \sum_{b=1}^{p_2}\prod_{n=1}^{|\bar s_b|}(\Lambda_{2,b}^\top-\bar\lambda_b(n\sgn s_b)\Pi_{2,b})^{\sgn b},
\end{aligned}
\end{align}
\end{definition}
\begin{pro}
The previously defined deformed moment matrix corresponds to a moment matrix with the following  positive/negative evolved
\begin{align}\label{d.evolved2}\begin{aligned}
w_{1,a}(s,x)&=\Ds_{a}'(x,s_a)w_{1,a}(x),& \Ds'_a&:=\begin{cases}
\prod_{n=1}^{s_a}|x-\lambda_a(n)|^2, &s_a>0,\\
1,& s_a=0,\\
\prod_{n=1}^{|s_a|}|x-\lambda_a(-n)|^{-2},& s_a<0,
\end{cases}\\
%  (x-\lambda_a(\sgn s_a))^{\sgn s_a}\cdots(x-\lambda_a(s_a))^{\sgn s_a}w_{1,a}(x),\\
w_{2,b}(s,x)&=\bar \Ds'_b(x,\bar s_b)^{-1}w_{2,b}(x),& (\bar \Ds'_b)^{-1}&:=\begin{cases}
\prod_{n=1}^{\bar s_b}|x-\bar \lambda_b(n)|^2, &\bar s_b>0,\\
1,&\bar s_b=0,\\
\prod_{n=1}^{|\bar s_b|}|x-\bar \lambda_b(-n)|^{-2},&\bar s_b<0.
\end{cases}\end{aligned}
%(x-\bar\lambda_b(\sgn \bar s_b))^{\sgn \bar s_b}\cdots(x-\bar\lambda_b(\bar s_b))^{\sgn \bar s_b}w_{2,b}(x).
\end{align}
\end{pro}
Proceeding as in the previous case
\begin{definition}
We introduce
\begin{align}
\begin{aligned}
q'_a&:=%(T_a D'_0){D'_0}^{-1}=
\I-\Pi_{1,a}(\I-|\lambda_a(s_a+1)|^2)-2\operatorname{Re}(\lambda_a(s_a+1))
\Lambda_{1,a}+\Lambda_{1,a}^{2},\\
\bar q'_b&:=%    {\bar D}'_0(T_b\bar D'_0)^{-1}=
\I  -\Pi_{2,b}(\I-|\bar \lambda_b(\bar s_b+1)|^2)-2\operatorname{Re}(\bar\lambda_b(\bar
\end{aligned}
\end{align}
\begin{align}
\label{delta'}
\begin{aligned}
\delta_a'&:=\I-C_{a}(\I-|\lambda_a(s_a+1)|^2)-2\operatorname{Re}(\lambda_a(s_a+1))
\bar\delta_b'&:=\I-\bar C_{b}(\I-|\bar \lambda_b(\bar s_b+1)|^2)-2\operatorname{Re}(\bar\lambda_b(\bar s_b+1))\bar
L_b+\bar L_b ^2.
\end{aligned}
\end{align}
\end{definition}

The wave and adjoint wave functions  now have the form
\begin{align}
\Psi_a^{(k)}(z,s)&=A^{(k)}_a(z,s)\Ds'_a(z,s_a),&
( \bar \Psi_b^*)^{(k)}(z,s)&=\bar A_b^{(k)}(z,s)\bar\Ds'_b(z,\bar s_b)^{-1},
\end{align}
and the expressions \eqref{linear.form.baker.1}-\eqref{cauchy.baker} still hold.

If we introduce $\omega_a'$ and $\bar\omega_b'$ as in \eqref{omega} but replacing $\delta$ by $\delta'$,  the
equations \eqref{d.wave}-\eqref{d.ZS1} hold
true by replacement of $\omega$ by $\omega'$. Now, the form $\omega'$ differs from \eqref{form.omega} as now we have
\begin{align}\label{form.omega'}\begin{aligned}
\omega_a'&=\omega'_{a,0}\Lambda^{2(|\vec n_1|-n_{1,a}+1)}+\dots+\omega'_{a,2(|\vec
{ \bar \omega_b}'&=\bar\omega'_{b,0}(\Lambda^\top)^{2(|\vec n_2|-n_{2,b}+1)}+\dots+\bar\omega'_{b,2(|\vec
\end{aligned}
\end{align}

With the definition of
\begin{align}
\begin{aligned}
\gamma_{a,a'}'(s,x)&:=  (1-\delta_{a,a'}(1-|x-\lambda_a(s_a+1)|^2),\\
\gamma_{b,b'}'(s,x)&:=  (1-\delta_{b,b'}(1-|x-\bar\lambda_b(\bar s_b+1)|^2),
\end{aligned}
\end{align}
we have that
\begin{pro}
The present setting  \eqref{multi.variante} and \eqref{dual.multi.variante} are replaced by
\begin{align}%\label{multi.variante}
(T_{a'} A_a^{(k)})\gamma'_{a,a'}&=\omega_{a',0}'A_a^{(k+2(|\vec n_1|-n_{1,a}+1))}+\dots+
\omega_{a',2(|\vec n_1|-n_{1,a}+1)}(k) A_a^{(k)},\\
\bar T_{b'}A_a^{(k)}&=\bar\omega_{b,0}'(k)A_a^{(k-2(|\vec n_2|-n_{2,b'}+1))}+\dots+\bar\omega_{b,2(|\vec
n_2|-n_{2,b'}+1)}' A_a^{(k)},
\end{aligned}\right.\\
\rho_{a',0}'(T_{a'} \bar A_b^{(k-2(|\vec n_1|-n_{1,a'}+1))})+\dots+ \rho_{a',2(|\vec n_1|-n_{1,a'}+1)}'(k) (T_{a'} \bar
A_b^{(k)})&= \bar A_b^{(k)},\\
\Big( \bar\rho_{b',0}' (k)(\bar T_{b'} \bar A_b^{(k+2(|\vec n_2|-n_{2,b'}+1)})+
\dots+\bar\rho'_{b',2(|\vec n_2|-n_{2,b'}+1)} (\bar T_{b'} \bar A_b^{(k)})\Big)\bar\gamma_{b,b'}'&= \bar A_b^{(k)}.
\end{aligned}\right.
\end{align}
\end{pro}
\end{appendices}

\begin{thebibliography}{99}

[1]M. Adler and P. van Moerbeke, \emph{Group factorization, moment matrices and Toda lattices}, Int. Math.
Res. Notices \textbf{12} (1997) 556-572.

[2] M. Adler and P. van Moerbeke, \emph{Generalized orthogonal polynomials, discrete KP and
Riemann--Hilbert  problems}, Commun. Math. Phys. \textbf{207} (1999) 589-620.

[3] M. Adler and P. van Moerbeke, \emph{Vertex operator solutions to the discrete KP
hierarchy}, Commun. Math. Phys. \textbf{203} (1999) 185-210.

[4] M. Adler and P. van Moerbeke, \emph{The spectrum of coupled random matrices}, Ann.
Math. \textbf{149} (1999) 921-976.

[5] M. Adler and P. van Moerbeke, \emph{Darboux transforms on band matrices, weights and
associated polynomials}, Int. Math. Res. Not. \textbf{18} (2001) 935-984.

[6] M. Adler, P. van Moerbeke and P. Vanhaecke, \emph{Moment matrices and multi-component
KP, with applications to random matrix theory}, Commun. Math. Phys. \textbf{286} (2009) 1-38.

[7] C. \'{A}lvarez-Fern\'{a}ndez, U. Fidalgo and M. Ma\~{n}as,
\emph{The multicomponent 2D Toda hierarchy: generalized matrix orthogonal polynomials, multiple
orthogonal polynomials and Riemann--Hilbert problems}.
Inv. Prob. \textbf{26} (2010) 055009   (17 pp).

[8] R. Apery. \emph{Irrationalite de $\zeta(2)$ et $\zeta(3)$}, Ast\`{e}risque \textbf{61} (1979) 11-13.

[9] M. J. Bergvelt and A. P. E. ten Kroode, \emph{Partitions, vertex operators constructions and
multi-component KP equations}, Pacific J. Math. \textbf{171} (1995) 23-88.

[10] F. Beukers, \emph{Pad\'{e} approximation in number theory,} Lec. Not. Math. \textbf{888}, Springer
Verlag, Berlin, 1981, 90-99.

[11] P.M. Bleher and A.B.J. Kuijlaars, \emph{ Random matrices with external source and multiple orthogonal
polynomials}, Int. Math. Res. Not. \textbf{2004} (2004), 109-129.

[12] G. Carlet, \emph{The extended bigraded Toda hierarchy} J. Phys. A: Math.  Gen. \textbf{39} (2006)

[13] J. Coates, \emph{On the algebraic approximation of functions}, I, II, III. Indag. Math. \textbf{28}
(1966) 421-461.

[14] E. Daems and A. B. J. Kuijlaars,  \emph{Multiple orthogonal polynomials of mixed type and
non-intersecting Brownian motions}, J. of Approx. Theory \textbf{146} (2007) 91-114.

[15] E. Daems and A. B. J. Kuijlaars,  \emph{A Christoffel--Darboux formula for multiple
orthogonal polynomials}, J. of Approx. Theory \textbf{130} (2004) 188-200.

[16]  E. Date, M. Jimbo, M. Kashiwara,  and T. Miwa, \emph{Operator approach to the
Kadomtsev-Petviashvili equation. Transformation groups for soliton equations. III}  J. Phys. Soc. Jpn.
\textbf{50},  (1981) 3806-3812.

[17]   E. Date, M. Jimbo, M. Kashiwara,  and T. Miwa, \emph{Transformation groups for soliton equations.
Euclidean Lie algebras and reduction of the KP hierarchy},  Publ. Res. Inst. Math. Sci. \textbf{18} (1982)

[18]   E. Date, M. Jimbo, M. Kashiwara,  and T. Miwa, \emph{Transformation groups for soliton equations} in
\emph{Nonlinear Integrable Systems-Classical Theory and Quantum Theory}  M. Jimbo and T. Miwa (eds.) World
Scientific, Singapore, 1983, pp. 39-120.

[19] R. Felipe  F. Ongay,  \emph{Algebraic aspects of the discrete KP hierarchy}, Linear Algebra App.
\textbf{338} (2001) 1-17.

[20] U. Fidalgo Prieto and G. L\'{o}pez Lagomasino, \emph{Nikishin systems are perfect}, \texttt{arXiv:1001.0554}

[21] A. S. Fokas, A. R. Its and A. V. Kitaev, \emph{The isomonodromy approach to matrix models in 2D quatum
gravity}, Commun. Math. Phys. (1992) 395-430.

[22] Ch. Hermite, \emph{Sur la fonction exponentielle}, C. R. Acad. Sci. Paris 77 (1873), 18-24, 74-79,
226-233, 285-293; reprinted in his Oeuvres, Tome III, Gauthier-Villars, Paris, 1912, 150-181.

[23] H. Jager, \emph{A simultaneous generalization of the Pad\'{e} table,  I-VI}, Indag. Math. \textbf{26} (1964),

[24]    V. G. Kac and J. W. van de Leur, \emph{The n-component KP hierarchy and representation theory}, J.
Math. Phys. \textbf{44} (2003) 3245-3293.

[25] A. B. J. Kuijlaars,   \emph{ Multiple orthogonal polynomial ensembles}, \texttt{arXiv:0902.1058v1 [math.PR] } (2009).

[26]  I. G. Macdonald,  \emph{Symmetric Functions and Hall Polynomials}, Clarendon Press, Oxford, 1995.

[27] K. Mahler, \emph{Perfect systems}, Compos. Math. \textbf{19} (1968), 95-166.

[28] M. Ma\~{n}as, L. Mart\'{\i}nez Alonso and C. \'{A}lvarez-Fern\'{a}ndez, \emph{The multicomponent 2D
Toda hierarchy: discrete flows and string equations}, Inv. Prob. \textbf{25} (2009) 065007 (31 pp).

[29] M. Ma\~{n}as and L. Mart\'{\i}nez Alonso, \emph{The multicomponent 2D Toda hierarchy: dispersionless
limit}, Inv. Prob. \textbf{25}  (2009) 115020 (22 pp).

[30] M. Mulase, \emph{Complete integrability of the Kadomtsev--Petviashvili equation}, Adv. Math.
\textbf{54} (1984)

[31]
E. M. Nikishin,  \emph{On simultaneous Pad\'{e} approximants} Matem. Sb. {\bf 113} (1980), 499--519 (Russian); English
translation in Math. USSR Sb.  {\bf 41} (1982), 409-425.

[32] A. Yu. Orlov and  D. M. Scherbin, \emph{Fermionic representation for basic hypergeometric functions
related to Schur polynomials}, \texttt{arXiv:nlin/0001001v4 [nlin.SI]}.

[33]  M. Sato, \emph{Soliton equations as dynamical systems on infinite dimensional Grassmann manifolds},
Res. Inst. Math. Sci. Kokyuroku \textbf{439}, (1981) 30-46 .

[34]  J. A.  Shohat and J.D. Tamarkin, \emph{The problem of moments},   American Mathematical Society (1943).

[35] B. Simon, \emph{The Christoffel-Darboux Kernel},  Proceedings of Symposia in Pure Mathematics
\textbf{79}:``Perspectives in Partial Differential Equations, Harmonic Analysis and Applications: A Volume in
Honor of Vladimir G. Maz'ya's 70th Birthday'',  (2008) 295-336. \texttt{arXiv:0806.1528}

[36]  V.N. Sorokin, \emph{On simultaneous approximation of several linear forms}, Vestnik Mosk Univ., Ser
Matem \textbf{1} (1983) 44-47.

[37] K. Ueno and K. Takasaki, \emph{Toda lattice hierarchy}, in Group Representations and Systems
of Differential Equations, Adv. Stud. Pure Math. \textbf{4} (1984) 1-95.

[38] R. P. Stanley, \emph{Enumerative combinatorics}, Cambridge University Press, Cambridge (1998).

[39] W. Van Assche, J. S. Geromino and A. B. J. Kuijlaars, \emph{Riemann--Hilbert problems for multiple
orthogonal polynomials} in: Bustoz et al (eds.), Special Functions 2000: Current Perspectives and Future
Directions, Kluwer Academic Publishers, Dordrecht, 2001, pp 23-59.

\end{thebibliography}

\end{document}