text:
\begin{document}

\title{Prohibitions caused by nonlocality \\ for Alice-Bob Boussinesq-KdV type systems}
\author{S. Y. Lou$^{1,2}$}
\affiliation{
$^{1}$\footnotesize{Center for Nonlinear Science and Department of Physics, Ningbo University, Ningbo, 315211, China}\\
$^2$\footnotesize Shanghai Key Laboratory of Trustworthy Computing, East China Normal University, Shanghai 200062, China }

\begin{abstract}
  It is found that two different celebrate models, the Korteweg de-Vrise (KdV) equation and the Boussinesq equation, are linked to a same model equation but with different nonlocalities. The model equation is called the Alice-Bob KdV (ABKdV) equation which was derived from the usual KdV equation via the so-called consistent correlated bang (CCB) companied by the shifted parity (SP) and delayed time reversal (DTR). The same model can be called as the Alice-Bob Boussinesq (ABB) system if the nonlocality is changed as only one of SP and DTR. For the ABB systems, with help of the bilinear approach and recasting the multi-soliton solutions of the usual Boussinesq equation to an equivalent novel form, the multi-soliton solutions with even numbers and the head on interactions are obtained. However, the multi-soliton solutions with odd numbers and the multi-soliton solutions with even numbers but with pursuant interactions are prohibited. For the ABKdV equation, the multi-soliton solutions exhibit many more structures because an arbitrary odd function of $x+t$ can be introduced as background waves of the usual KdV equation.\\
\leftline{\pacs{02.30.Ik, }}
\end{abstract}

\maketitle

\section{Introduction}
Recently, the so-called nonlocal integrable systems (or the Alice-Bob systems) have attracted many attentions of scientists because of the first paper presented by
 Ablowitz and Musslimani \cite{AM} and the possible applications for two-place physics (Alice-Bob physics) \cite{AB,ABs}. The first two-place nonlocal model is related to the places $\{x,\ t\}$ and $\{-x,\ t\}$, respectively, for the
nonlinear Schr\"odinger (NLS) equation\cite{AM}
\begin{equation}
iA_{t}+A_{xx}\pm A^2B=0,
\ \ \  B=\hat{f}A=\hat{P}\hat{C}A=A^*(-x,t), \label{Eq1}
\end{equation}
where the operators $\hat{P}$ and $\hat{C}$ are the usual parity and charge conjugation.
It is clear that the equation system \eqref{Eq1} with its conjugate form is form invariant under the discrete transformation group with the generators $\{\hat{P},\ \hat{C},\ \hat{T}\}$, where $\hat{T}$ is the time reversal operator.
In literature, the nonlocal nonlinear Schr\"odinger equation \eqref{Eq1} is also called parity-time reversal (PT) symmetric. PT symmetry plays an important role in the quantum physics \cite{2} and many other areas of physics, such as the quantum chromodynamics \cite{3}, electric circuits \cite{4}, optics \cite{5,6}, Bose-Einstein condensates \cite{7}, atmospheric and oceanic dynamics \cite{AB,Tang} etc.

It is well known that there are various correlated and/or entangled events that may happen in different times and places. To describe some types of two-place physical problems, Alice-Bob (AB) systems \cite{AB,ABs} are proposed by using the shifted parity ($\hat{P}_{s},\ \hat{P}_{s}x=-x+x_0 $), delayed time reversal ($\hat{T}_{d},\ \hat{T}_{d}t=-t+t_0$) and charge conjugate ($\hat{C}$) symmetries. If one event (A, Alice event) is correlated/entangled to another (B, Bob event), we denote the correlated relation as $B=\hat{f}A$ for suitable $\hat{f}$ operators. Usually, the event $A=A(x,\ t)$ happened at $\{x,\ t\}$ and event $B=B(x',\ t')$ happened at $\{x',\ t'\}=\hat{f}\{x,\ t\}$. In fact, $\{x',\ t'\}$ is usually far away from $\{x,\ t\}$. Hence, the intrinsic two-place models or Alice-Bob systems are nonlocal. In addition to the nonlocal nonlinear Schr\"odinger equation \eqref{Eq1}, there are many other types of two-place nonlocal models, such as the nonlocal KdV systems \cite{9,AB,ABs}, the nonlocal modified KdV systems \cite{10,11,ABs}, the discrete nonlocal NLS systems \cite{12}, the coupled nonlocal NLS systems \cite{13}, the nonlocal Davey-Stewartson systems \cite{14,15,16}, general nonlocal NLS equation \cite{TL}  and the nonlocal peakon systems \cite{ABP} including AB Xia-Qiao-Zhou (ABXQZ), AB Cammasa-Holm (ABCH), AB Degasperis-Procesi (ABDP), AB Novikov (ABN), AB FORQ (ABFORQ) and so on.

In \cite{ABs}, we proposed a series of other types of integrable AB systems including the AB-KdV systems, AB-mKdV systems, AB-KP systems, AB-sine Gordon systems, AB-NLS systems, AB-Toda systems and AB-H1 systems. Furthermore, by using the $\hat{P}_{s},\ \hat{T}_{d}$ and $\hat{C}$ symmetries, their $\hat{P}_{s},\ \hat{T}_{d}$ and $\hat{C}$ invariant muti-soliton solutions are obtained in the elegant forms for all the AB systems listed in \cite{ABs}.
 In addition, we established a most general AB-KdV equation and presented its $\hat{P}_{s},\ \hat{T}_{d}$ and $\hat{C}$ invariant Painlev\'e II reduction and soliton-cnoidal periodic wave interaction solutions \cite{AB}.

 In section II of this paper, we write down an integrable AB real system with three different nonlocal properties. One non-locality is related to the ABKdV system while the other two are associated with the ABB systems. The Lax pairs for these two kinds of AB systems are explicitly given. The multi-soliton solutions for the ABB systems are studied in section III and the possible prohibitions on multi-soliton solutions are also discussed in this section. Section IV is devoted to investigating the multi-soliton solutions of the ABKdV system. The last section is a short summary and discussion.

\section{AB integrable systems come from same equation with different non-localities}

In Ref. \cite{AB}, a special ABKdV system (after some re-scaling transformations)
\begin{eqnarray}
&&A_t=2(A+B)A_x+A_{xxx}+G(A,B),\label{ABKdVG}\\
&& B=\hat{P}_s\hat{T}_dA=A(-x+x_0,\ -t+t_0),\label{Bxt}
\end{eqnarray}
is derived from the usual KdV equation by using the so-called consistent correlated bang (CCB) \cite{CCB} with an arbitrary $\hat{P}_s\hat{T}_d$ invariant functional $G(A,\ B)=\hat{P}_s\hat{T}_dG(A,\ B)$.

In this paper, we consider a special $\hat{P}_s\hat{T}_d$ invariant selection of $G(A,\ B)$,
\begin{eqnarray}
&&G(A,\ B)=\frac12\left(1-2A-2B\right)\left(A-B\right)_x+\frac12(B-A)_{xxx}.\label{G}
\end{eqnarray}
Under the selection \eqref{G}, \eqref{ABKdVG} becomes
\begin{eqnarray}
&&A_t=\frac12(A-B+(A+B)^2+A_{xx}+B_{xx})_x. \label{ABKdV}
\end{eqnarray}
The integrability of \eqref{ABKdV} with the nonlocal condition \eqref{Bxt} is trivial because it is derived from the usual KdV equation via CCB. In fact, the ABKdV system \eqref{ABKdV} with \eqref{Bxt} possesses the following Lax pair,
\begin{eqnarray}
&&\psi_{xx}=-\frac13(U-\Lambda)\psi,\  \label{Lx}\\
&&\psi_{t}=4\psi_{xxx}+\frac16\left(7U-V\right)_x\psi
+\frac16\left(10U+2V+3E\right)\psi_x,\label{Lt}\\
&&U\equiv\left(\begin{array}{cc}
A & B \\ B & A
\end{array}\right),\ \Lambda\equiv \left(\begin{array}{cc}
\lambda & \lambda \\ \lambda & \lambda
\end{array}\right),\ V\equiv \left(\begin{array}{cc}
B & A \\ A & B \end{array}\right),\ E\equiv \left(\begin{array}{cc}
1 & -1 \\ -1 & 1 \end{array}\right). \nonumber
\end{eqnarray}
Now, an important and interesting question is, are there any other types of non-localities such that \eqref{ABKdV} is still integrable? Fortunately, we can find that if we introduce the following two new real nonlocal conditions
\begin{eqnarray}
&&B=\hat{P}_sA=A(-x+x_0,\ t), \label{Bx}
\end{eqnarray}
and/or
\begin{eqnarray}
&&B=\hat{T}_dA=A(x,\ -t+t_0),\label{Bt}
\end{eqnarray}
\eqref{ABKdV} is really integrable because of the existence of the Lax pair,
\begin{eqnarray}
&&\psi_{xxx}=-\frac12(A+B)\psi_x-\frac1{4}\left[\frac{\sqrt{-3}}3(B-A)
+A_x+B_x-4\lambda\right]\psi,\label{Lax}
\end{eqnarray}
\begin{eqnarray}
&&\psi_t=\frac{\sqrt{-3}}3\left[3\psi_{xx}+(A+B)\psi\right],\ B=\hat{P}_s^{\alpha}\hat{T}_d^{1-\alpha}A, \ \alpha=0,1. \label{Lat}
\end{eqnarray}

It should be emphasized that though the nonlocal equation \eqref{ABKdV} is same for all three non-localities \eqref{Bxt},\ \eqref{Bx} and \eqref{Bt}, their integrable properties, Lax pairs, are quite different.

From the Lax pair \eqref{Lx} and \eqref{Lt}, we know that \eqref{ABKdV} with \eqref{Bxt} is a KdV type nonlocal system. However, \eqref{ABKdV} with \eqref{Bx} or \eqref{Bt} is a Boussinesq type nonlocal equation because its Lax pair possesses the form \eqref{Lax} and \eqref{Lat}.

\section{Mult-soliton solutions for the ABB system (\ref{ABKdV}) with (\ref{Bx}) and/or (\ref{Bt})}
The multi-linear approach, especially, Hirota's bilinear method, is a powerful method to looking for multi-soliton solutions for integrable nonlinear systems. By using the standard Hirota's bilinear operator,
\begin{equation*}
D_x^mD_t^nf\cdot f \equiv \left.(\partial_x-\partial_{x'})^m(\partial_t-\partial_{t'})^nf(x,t)f(x',t')\right|_{x'=x,t'=t},
\end{equation*}
it is straightforward to prove that the AB system \eqref{ABKdV} can be changed to the following eight-linear form, \begin{eqnarray}
&&2f^2g^2(g^2D_x^4f\cdot f+f^2D_x^4g\cdot g -2g^2D_t^2f\cdot f)+3(g^2D_xD_tf\cdot f-f^2g\cdot g)^2\nonumber\\
&&\quad -3(g^2D_x^2f\cdot f-f^2D_x^2g\cdot g)^2+2f^2g^2[g^2D_x^2f\cdot f-f^2D_x^2g\cdot g -g^2D_xD_t f\cdot f\nonumber\\
&&\quad +f^2D_xD_tg\cdot g+3 (D_xD_tf\cdot f)(D_x^2g\cdot g)-3 (D_xD_tg\cdot g)(D_x^2f\cdot f)=0, \label{MLABKdV}
\end{eqnarray}
after directly substituting the transformation
\begin{equation}
A=3 (\ln f)_{xx}+3 (\ln f)_{xt},\ B=3 (\ln g)_{xx}-3 (\ln g)_{xt}, \label{A}
\end{equation}
into the ABB equation \eqref{ABKdV} and integrating once for the variable $x$ with
\begin{equation}
g=f(-x+x_0,\ t) \quad \mbox{\rm or}\quad g=f(x,\ -t+t_0). \label{g}
\end{equation}
From the relations \eqref{Bx} or \eqref{Bt}, \eqref{A} and \eqref{g}, one can readily find
that
\begin{equation}
g=f f_1(t)f_2(x+t) \label{fg}
\end{equation}
where $f_1(t)$ and $f_2(x+t)$ are arbitrary functions of $t$ and $x+t$ respectively. By using the relation \eqref{fg}, \eqref{MLABKdV} is simplified to a very simple one
\begin{equation}
(D_{t}^2-D_x^4)f\cdot f=0 \label{BLBq}
\end{equation}
which is $f_1$ and $f_2$ independent. In fact, substituting \eqref{fg} into \eqref{A}, the result is also $f_1$ and $f_2$ independent. Thus, one can directly select $f$ as the shifted parity invariant or delayed time reversal invariant function, i.e., $f_1=f_2=1$ and
\begin{equation}
f(-x+x_0,\ t)=f(x,\ t)\quad \mbox{\rm or}\quad  f(x,\ -t+t_0)=f(x,\ t). \label{even}
\end{equation}
It is clear that the only remaining thing is to solve the bilinear equation \eqref{BLBq} with the invariant condition \eqref{even} to get multi-soliton solutions.

Eq. \eqref{BLBq} is just the bilinear form of the well known Boussinesq equation. The multi-soliton solutions of the usual bilinear Boussinsq system can be written as ($f=f_n$) \cite{ABs}
\begin{eqnarray}
f_n=\sum_{\{\nu\}}K_{\{\nu\}}\cosh\left(\sum_{i=1}^n\nu_i\xi_i\right), \ \xi_i=k_ix+2\delta_ik_i^2t+\xi_{i0},\ \delta_i^2=1,\label{nsol}
\end{eqnarray}
where the summation of $\{\nu\}=\{\nu_1,\ \nu_2,\ \ldots,\ \nu_n\}$ should be done for all
permutations of $\nu_i=1,\ -1,\ i=1,\ 2,\ \ldots,\ n$, and
\begin{equation}
K_{\{\nu\}}=\prod_{i<j}a_{ij},\quad a_{ij}^2=2k_i^2+2k_j^2-k_ik_j(\delta_i\delta_j+3\nu_i\nu_j).\label{K}
\end{equation}

Generally, the solution \eqref{nsol} is not shifted parity or delayed time reversal invariant. Before to give a general result, we consider the possible constraints on \eqref{nsol} because of the condition \eqref{even} for small $n$.

For $n=1$, we have
\begin{eqnarray}
f_1=\cosh\left(k_1x+2\delta_1k_1^2t+\xi_{10}\right), \  \delta_1^2=1.\label{sol1}
\end{eqnarray}
From \eqref{sol1}, the invariant condition \eqref{even} becomes
\begin{eqnarray}
\cosh\left[k_1\left(x_0-x\right)+2\delta_1k_1^2t+\xi_{10}\right]=
\cosh\left(k_1x+2\delta_1k_1^2t+\xi_{10}\right),\label{ev1x}
\end{eqnarray}
or
\begin{eqnarray}
\cosh\left[k_1x+2\delta_1k_1^2\left(t_0-t\right)+\xi_{10}\right]=
\cosh\left(k_1x+2\delta_1k_1^2t+\xi_{10}\right),\label{ev1t}
\end{eqnarray}
It is not difficult to see that there is no possible nontrivial solution of \eqref{ev1x} or \eqref{ev1t} for arbitrary $\{x,\ t\}$ by selecting parameters $k_1,\ \xi_{10},\ x_0$ and $t_0$ .

For $n=2$, the solution \eqref{nsol} becomes,
\begin{eqnarray}
f_2&=&\delta_+\sqrt{2k_1^2+2k_2^2-k_1k_2(\delta_1\delta_2+3)}
\cosh\left[(k_1+k_2)x+2(\delta_1k_1^2+\delta_2k_2^2)t+\xi_{10}+\xi_{20}\right]\nonumber\\
&&
+\delta_-\sqrt{2k_1^2+2k_2^2-k_1k_2(\delta_1\delta_2-3)}
\cosh\left[(k_1-k_2)x+2(\delta_1k_1^2-\delta_2k_2^2)t+\xi_{10}-\xi_{20}\right],\nonumber \\
&&
\delta_1^2=\delta_2^2=\delta_+^2=\delta_-^2=1,\label{sol2}
\end{eqnarray}
From the expression \eqref{sol2}, we know that if we take the parameter conditions,
\begin{eqnarray}
k_2=\pm k_1,\ \delta_2=\mp\delta_1,\ \xi_{10} = -\frac{k_1}2(2\delta_1k_1t_0+x_0),\ \xi_{20} = -\frac{k_2}2(x_0+2\delta_2k_2t_0),\label{kd2}
\end{eqnarray}
then two-soliton solution \eqref{sol2} is restricted to the simple form
\begin{eqnarray}
f_2&\sim & \delta_{\pm}
\cosh\left[2k_1\left(x-\frac{x_0}2\right)\right]+2\delta_{\mp}
\cosh\left[4k_1^2\left(t-\frac{t_0}2\right)\right] \label{sol21}
\end{eqnarray}
up to a neglected constant factor $\sqrt{2}k_1$.

It is not difficult to check that \eqref{sol21} is the only shifted parity or delayed time reversal invariant form of \eqref{sol2}.

The conditions $k_2=\pm k_1,\ \delta_2=\mp\delta_1$ for the two-soliton solution \eqref{sol2} implies that two solitons can exist only for head-on collisions with the same velocities and wave numbers. In other words, because of the nonlocal conditions \eqref{Bx} or \eqref{Bt}, the pursuant interactions between two solitons are prohibited. The head on collisions with different wave numbers (and then velocities) are also prohibited.

Fig. 1 exhibits the head on collision expressed by \eqref{A} with \eqref{sol21} for the field $A$ under the parameter selections
\begin{equation}
\delta_{\pm}=\delta_{\mp}=k_1=1.\label{p2}
\end{equation}

\input epsf
\begin{figure}
\centering\epsfxsize=7cm\epsfysize=5cm\epsfbox{BqFig01.eps}
\caption{Head on collision interaction between soliton and antisoliton for the field $A$ described by Eq. (\ref{A}) with \eqref{sol21} under the parameter
selections  (\ref{p2}).}\label{fig1}
\end{figure}

For $n=3$, the solution \eqref{nsol} possesses the form
\begin{eqnarray}
f_3&=&K_{\{\}}
\cosh\left(\xi_1+\xi_2+\xi_3\right)+K_{\{3\}}
\cosh\left(\xi_1+\xi_2-\xi_3\right)\nonumber\\
&&+K_{\{2\}}
\cosh\left(\xi_1-\xi_2+\xi_3\right)+K_{\{1\}}
\cosh\left(-\xi_1+\xi_2+\xi_3\right),\label{3sol}
\end{eqnarray}
where
\begin{eqnarray}
&&\xi_i= k_ix + 2\delta_ik_it +\xi_{i0}, \ \delta_i^2=1, \ i=1,\ 2,\ 3,\ a_{ij}^{\pm}=\sqrt{2k_i^2+2k_j^2-k_ik_j(\delta_i\delta_j\pm 3)}, \nonumber\\
&&K_{\{\}}=a_{12}^+a_{13}^+a_{23}^+, \ K_{\{3\}}=a_{12}^+a_{13}^-a_{23}^-, \ K_{\{2\}}=a_{12}^-a_{13}^+a_{23}^-, \ K_{\{1\}}=a_{12}^-a_{13}^-a_{23}^+. \label{KK}
\end{eqnarray}
In \eqref{3sol} and the following of this paper, the used notation $K_{\{i_1..i_j\}}$ means that for the permutation (the subscription of $K_{\{\nu\}}$), $\{\nu\}\equiv\{\nu_1,\ \nu_2,\ \ldots,\ \nu_n\}\equiv \{\nu_{i}=-1,\ i=i_{1},\ i_2,\ \ldots,\ i_j; \nu_{i}=1,\ i\neq i_{1},\ i_2,\ \ldots,\ i_j\} $. Thus, for $n=3$, $K_{\{\nu\}}$ in \eqref{3sol} possesses the forms $K_{\{\}}=K_{\{\nu_1=1,\nu_2=1,\nu_3=1\}},\ K_{\{1\}}=K_{\{\nu_1=-1,\ \nu_2=1,\ \nu_3=1\}}$ and so on.

After finishing some detailed analysis, one can find that there is no possible three soliton solution \eqref{3sol} because of the invariant condition \eqref{even}. Furthermore, all the soliton solutions with odd number $n$ will be prohibited by the shifted parity or the delayed time reversal nonlocality \eqref{even} or equivalently \eqref{Bx} and/or \eqref{Bt}.

For $n=4$, the solution \eqref{nsol} can be explicitly written as
\begin{eqnarray}
f_4&=&K_{\{\}}
\cosh\left(\xi_1+\xi_2+\xi_3+\xi_4\right)+K_{\{4\}}
\cosh\left(\xi_1+\xi_2+\xi_3-\xi_4\right)\nonumber\\
&&+K_{\{3\}}
\cosh\left(\xi_1+\xi_2-\xi_3+\xi_4\right)+K_{\{2\}}
\cosh\left(\xi_1-\xi_2+\xi_3+\xi_4\right)\nonumber\\
&&+K_{\{1\}}
\cosh\left(-\xi_1+\xi_2+\xi_3+\xi_4\right)+K_{\{34\}}
\cosh\left(\xi_1+\xi_2-\xi_3-\xi_4\right)\nonumber\\
&&+K_{\{24\}}
\cosh\left(\xi_1-\xi_2+\xi_3-\xi_4\right)+K_{\{23\}}
\cosh\left(\xi_1-\xi_2-\xi_3+\xi_4\right),\label{4sol}
\end{eqnarray}
where
\begin{eqnarray}
&&K_{\{\}}=a_{12}^+a_{13}^+a_{14}^+a_{23}^+a_{24}^+a_{34}^+, \ K_{\{4\}}=a_{12}^+a_{13}^+a_{14}^-a_{23}^+a_{24}^-a_{34}^-, \nonumber\\ &&K_{\{3\}}=a_{12}^+a_{13}^-a_{14}^+a_{23}^-a_{24}^+a_{34}^-, \ K_{\{2\}}=a_{12}^-a_{13}^+a_{14}^+a_{23}^-a_{24}^-a_{34}^+\nonumber\\
&&K_{\{1\}}=a_{12}^-a_{13}^-a_{14}^-a_{23}^+a_{24}^+a_{34}^+,\
K_{\{34\}}=a_{12}^+a_{13}^-a_{14}^-a_{23}^-a_{24}^-a_{34}^+\nonumber\\
&&K_{\{24\}}=a_{12}^-a_{13}^+a_{14}^-a_{23}^-a_{24}^+a_{34}^-,\
K_{\{23\}}=a_{12}^-a_{13}^-a_{14}^+a_{23}^+a_{24}^-a_{34}^-. \label{KK4}
\end{eqnarray}
After finishing some detailed analysis, for the four soliton solution \eqref{4sol}, there is only one independent selection
\begin{equation}
k_3=k_1\equiv\frac{k}2,\ k_4=k_2\equiv\frac{\kappa}2,\ \delta_3=-\delta_1=1,\ \delta_4=-\delta_2=1,\ \xi_{i0}=-\frac{k_i}2(x_0+2\delta_ik_i^2t_0). \label{kkdd}
\end{equation}
All other possible selections are equivalent to \eqref{kkdd}.

Under the restriction \eqref{kkdd}, the four soliton solution \eqref{KK4} becomes a $\hat{P}_s$ and $\hat{T}_d$ invariant form (up to a common factor $k\kappa$),
\begin{eqnarray}
&&f_4\sim 2(k^2-\kappa^2)\sqrt{(k^2+\kappa^2)^2-k^2\kappa^2}\left[\cosh\left(\eta_1\right)
+\cosh\left(\eta_2\right)+\cosh\left(\hat{T}_d\eta_1\right)
+\cosh\left(\hat{P}_s\eta_2\right)\right]\nonumber\\
&&\quad +\big(k^2+\kappa^2+k\kappa\big)(k+\kappa)^2
\cosh\left[\big(k-\kappa\big)\zeta\right] +\big(k^2+\kappa^2-k\kappa\big)(k-\kappa)^2
\cosh\left[\big(k+\kappa\big)\zeta\right] \nonumber\\
&&\quad +4\left[\big(k^2+\kappa^2\big)^2-k^2\kappa^2\right]
\cosh\left[\big(k^2-\kappa^2\big)\tau \right] +4\big(k^2-\kappa^2\big)^2
\cosh\left[\big(k^2+\kappa^2\big)\tau\right],\label{f4}\\
&&\eta_1=k\zeta
+\kappa^2\tau,\ \eta_2=\kappa \zeta
+k^2\tau,\ \zeta\equiv x-\frac{x_0}2,\ \tau\equiv t-\frac{t_0}2.
\end{eqnarray}

Fig. 2 is a plot of the head on collision interaction for the four soliton solution described by \eqref{A} with \eqref{f4} and the parameters are fixed as
\begin{equation}
x_0=t_0=0,\ k=2.5, \ \kappa=2. \label{p4}
\end{equation}

\input epsf
\begin{figure}
\centering\epsfxsize=7cm\epsfysize=5cm\epsfbox{BqFig02.eps}
\caption{Head on collision of two paired soliton-antisolitons exhibited in Eq. (\ref{A}) with \eqref{f4} and the parameter
selections  (\ref{p4}).}\label{fig2}
\end{figure}

It is clear that the solution \eqref{f4} is $\hat{P}_s$ and $\hat{T}_d$ invariant. The only independent selection \eqref{kkdd} for the four soliton solution implies that there are only soliton-antisoliton pairs with the same wave number but the opposite velocities.  While other types of soliton interactions such as the pursuant interaction and the head on collision with different wave numbers are prohibited by the shifted parity and or the delayed time reversal nonlocality conditions \eqref{even} or equivalently \eqref{Bx} and \eqref{Bt}. In fact, this conclusion is correct for all $n=2N$ soliton solution \eqref{nsol}.

For $n=2N$, the invariant condition of \eqref{nsol} becomes
\begin{equation}
k_{N+i}=k_i,\  \delta_{i}=-\delta_{N+i}=1,\
\xi_{i0}=-\frac{k_i}2(x_0+2\delta_ik_i^2t_0). \label{kN}
\end{equation}
Because of the paring condition \eqref{kN}, the even function property of cosh function in \eqref{nsol} and the summation is done for all possible permutations $\{\nu\}=\{\nu_1,\ \ldots,\ \nu_{2N}\}$ with $\nu_i=1,\ -1$, the solution $f_{2N}$ expressed by \eqref{nsol} with \eqref{kN} is always $\hat{O}$-invariant
\begin{equation}
\hat{O}\equiv \hat{P}_s^{\alpha}\hat{T}_d^{1-\alpha}, \ \alpha=0,\ 1. \label{O}
\end{equation}
 In other words, \eqref{nsol} with $n=2N$ and \eqref{kN} is even for both $\zeta\equiv \left(x-\frac{x_0}2\right)$ and $\tau\equiv \left(t-\frac{t_0}2\right)$.
For instance, for $N=3$ we have
\begin{eqnarray}
f_6&=&K_{\{23\}}\left[\cosh (y_{1+}) +\cosh \big(\hat{O}y_{1+}\big) \right]
+K_{\{35\}}\left[\cosh (y_{1-}) +\cosh \big(\hat{O}y_{1-}\big) \right]\nonumber\\
&&+ K_{\{13\}}\left[\cosh (y_{2+}) +\cosh \big(\hat{O}y_{2+}\big) \right]
+K_{\{34\}}\left[\cosh (y_{2-}) +\cosh \big(\hat{O}y_{2-}\big) \right]\nonumber\\
&&+K_{\{12\}}\left[\cosh (y_{3+}) +\cosh \big(\hat{O}y_{3+}\big) \right]
+K_{\{24\}}\left[\cosh (y_{3-}) +\cosh \big(\hat{O}y_{3-}\big) \right]\nonumber\\
&&+K_{\{3\}}\left[\cosh (z_{3+}) +\cosh \big(\hat{O}z_{3+}\big) \right]
+K_{\{146\}}\left[\cosh (z_{3-}) +\cosh \big(\hat{O}z_{3-}\big) \right]\nonumber\\
&&+K_{\{2\}}\left[\cosh (z_{2+}) +\cosh \big(\hat{O}z_{2+}\big) \right]
+K_{\{356\}}\left[\cosh (z_{2-}) +\cosh \big(\hat{O}z_{2-}\big) \right]\nonumber\\
&&+K_{\{1\}}\left[\cosh (z_{1+}) +\cosh \big(\hat{O}z_{1+}\big) \right]
+K_{\{346\}}\left[\cosh (z_{1-}) +\cosh \big(\hat{O}z_{1-}\big) \right]\nonumber\\
&&+K_{\{\}}\cosh\left[2(k_1+k_2+k_3)\zeta\right]
+K_{\{14\}}\cosh\left[2(k_1-k_2-k_3)\zeta\right]\nonumber\\
&&+K_{\{25\}}\cosh\left[2(k_2-k_1-k_3)\zeta\right]
+K_{\{36\}}\cosh\left[2(k_3-k_1-k_2)\zeta\right]\nonumber\\
&&+K_{\{456\}}\cosh\left[4(k_1^2+k_2^2+k_3^2)\tau\right]
+K_{\{156\}}\cosh\left[4(k_1^2-k_2^2-k_3^2)\tau\right]\nonumber\\
&&+K_{\{246\}}\cosh\left[4(k_2^2-k_1^2-k_3^2)\tau\right]
+K_{\{126\}}\cosh\left[4(k_3^2-k_1^2-k_2^2)\tau\right], \label{f6}
\end{eqnarray}
where
\begin{eqnarray}
y_{i\pm}&=&2k_i\zeta +4(k_{i+1}^2\pm k_{i+2}^2)\tau,\ i=1,2,3,\ \mod 3,\nonumber\\
z_{i\pm}&=&4k_i^2\tau+2(k_{i+1}\pm k_{i+2})\zeta, \ i=1,2,3,\ \mod 3. \label{yz}
\end{eqnarray}
Fig. 3 displays the three paired soliton-antisolitons head on collision interaction expressed by \eqref{A} with \eqref{f6} and parameter selections
\begin{equation}
x_0=t_0=0,\ k_1=-1.2, \ k_2=1.1,\ k_3=-1. \label{p6}
\end{equation}

\input epsf
\begin{figure}
\centering\epsfxsize=7cm\epsfysize=5cm\epsfbox{BqFig03.eps}
\caption{Plots of the head on collision interactions of six solitons expressed by Eq.  (\ref{A}) with the parameter
selections  (\ref{p6}).} \label{fig3}
\end{figure}

\section{Multi-soliton solutions for the ABKdV system (\ref{ABKdV}) with (\ref{Bxt})}
For the ABKdV system \eqref{ABKdV} with the nonlocal condition \eqref{Bxt}, we can make the bilinear transformation
\begin{equation}
A=3(\ln f)_{xx}+3(\ln g)_{xt},\ B=3(\ln f)_{xx}-3 (\ln g)_{xt}. \label{BLT}
\end{equation}
Substituting the transformation relation \eqref{BLT} into the ABKdV system \eqref{ABKdV} and its adjoint form (obtained by applying the operator $\hat{P}_s\hat{T}_d$ on \eqref{ABKdV}), we can find two separated bilinear equations
\begin{eqnarray}
&&(D_xD_t-D_x^4)f\cdot f=0,\label{BLf}\\
&&(D_xD_t-D_t^2)g\cdot g=0. \label{BLg}
\end{eqnarray}
It is interesting that the $f$ equation is just the bilinear form of the usual KdV equation while the general solution of \eqref{BLg} can be written as
\begin{eqnarray}
g=g_0(x)g_1(x+t)\sim g_1(x+t). \label{rg}
\end{eqnarray}
From the expression \eqref{BLT}, we know that $g_0(x)$ can be simply taken as $1$ without loss of generality.

Applying the nonlocal condition \eqref{Bxt}, we have
\begin{eqnarray}
\hat{P}_s\hat{T}_df=f(-x+x_0,-t+t_0)=
\frac{f\exp\left[f_1(t)x+f_0(t)\right]}{g_1(-x-t+x_0+t_0)g_1(x+t)} ,\label{rf}
\end{eqnarray}
where $f_0,\ f_1$ and $g$ are arbitrary functions of the indicated variables. Especially, $f$ may be $\hat{P}_s\hat{T}_d$-invariant if we take the following special selections
\begin{equation}
f_1(t)=f_0(t)=0,\ g_1(-x-t+x_0+t_0)=g_1(x+t)^{-1}. \label{fg1}
\end{equation}

In general, the solutions of \eqref{BLf} are not $\hat{P}_s\hat{T}_d$-invariant. Fortunately, because the expression \eqref{BLT} is invariant under the transformation
\begin{equation}
f\longrightarrow \exp(Kx+\Omega t+X)f \label{ff}
\end{equation}
with arbitrary constants $K,\ \Omega$ and $X$, the multi-soliton solutions of the ABKdV equation can be written as
\begin{equation}
A=3(\ln F_N)_{xx}+3(\ln G)_{xt}, \label{NS}
\end{equation}
where $G=G(x+t)$ is an arbitrary function with the condition
\begin{equation}
\hat{P}_s\hat{T}_dG=G^{-1}, \label{CG}
\end{equation}
and
\begin{equation}
F_N=\sum_{\{\nu\}}K_{\{\nu\}}\cosh\left(\sum_{i=1}^N\nu_i\xi_i\right),\ \xi_i=k_i\left(x-\frac{x_0}2\right)+4k_i^3\left(t-\frac{t_0}2\right) \label{FN}
\end{equation}
The summation of $\{\nu\}=\{\nu_1,\ \nu_2,\ \ldots,\ \nu_N\}$ in \eqref{FN} should be done for all permutations of $\nu_i=1,\ -1,\ i=1,\ 2,\ \ldots,\ N$, and
\begin{equation}
K_{\{\nu\}}=\prod_{i<j}(k_i-\nu_i\nu_j k_j).\label{K}
\end{equation}
It is clear that \eqref{FN} is really $\hat{P}_s\hat{T}_d$-invariant.

Different from the multi-soliton solutions of the ABB system discussed in the last section for the non-localities of \eqref{Bx} and \eqref{Bt}, there is no other prohibitions for the ABKdV system with the non-locality \eqref{Bxt} except for the locations of $N$ solitons.

Fig. 4 displays one soliton solution with a periodic background wave expressed by \eqref{NS} with \eqref{FN} for $N=k_1=1,x_0=t_0=0$ and $G(x+t)$ is taken as
\begin{equation}
G(x+t)=\exp(-0.2\sin(x+t)). \label{G1}
\end{equation}

\input epsf
\begin{figure}
\centering\epsfxsize=7cm\epsfysize=5cm\epsfbox{BqFig04.eps}
\caption{One soliton solution \eqref{NS} with \eqref{FN} for $N=k_1=1,x_0=t_0=0$ and the periodic background wave (\ref{G1}).} \label{fig4}
\end{figure}

Fig. 5 displays the interaction between the KdV soliton and the background soliton with few peaks expressed by \eqref{NS}, \eqref{FN} for $N=k_1=1,\ x_0=t_0=0$ and
\begin{equation}
G(x+t)=\exp\left[-\int\int \mbox{\rm sech}(x+t)\sin(2x+2t)\mbox{\rm dxdt}\right]. \label{G2}
\end{equation}

\input epsf
\begin{figure}
\centering\epsfxsize=7cm\epsfysize=5cm\epsfbox{BqFig05.eps}
\caption{Plot of the interaction between the usual KdV soliton and the background few cycle soliton expressed by Eq. (\ref{NS}), \eqref{FN} for $N=k_1=1,\ x_0=t_0=0$ and (\ref{G2}).} \label{fig3}
\end{figure}

The main difference between the usual KdV equation (\eqref{ABKdV} with $B=A$) and the ABKdV \eqref{ABKdV} with $B=\hat{P}_s\hat{T}_dA$) is that there is an arbitrary odd background wave
\begin{equation}
W(x+t)\equiv 3(\ln G(x+t))_{xt}=-\hat{P}_s\hat{T}_dW(x+t). \label{W}
\end{equation}

\section{summary and discussion}
It is shown that there are some really different phenomena in Alice-Bob systems (or namely nonlocal systems). From the result of this paper we can conclude that the following eight systems
\begin{eqnarray}
&&A_t=\frac12(A-B+(A+B)^2+A_{xx}+B_{xx})_x, \label{ABKdVA}\\
&&B=g_iA,\quad i=0,\ 1,\ 2,\ \ldots,\ 8, \label{ABKdVB}
\end{eqnarray}
where $g_i$ belongs to the eight order Abel group ${\cal G}$
\begin{eqnarray}
&&{\cal G}=\{g_0,\ g_1,\ \ldots,\ g_8\}\equiv \{I,\ \hat{P}_s,\ \hat{T}_d,\ \hat{C},\  \hat{P}_s\hat{T}_d,\ \hat{P}_s\hat{C},\  \hat{T}_d\hat{C},\ \hat{P}_s\hat{T}_d\hat{C}\},    \label{Group}
\end{eqnarray}
where $I=g_i^2$ for all $i$ is the identity of the group ${\cal G}$.  For $B=A$ and $B=\hat{C}A$, the models are just the real and complex local KdV equations. For $B=\hat{P}_s\hat{T}_d A$ and $B=\hat{P}_s\hat{T}_d\hat{C}A$, the models are related to the real and complex nonlocal KdV systems.
For $B=\hat{P}_sA$, $B=\hat{T}_dA$, $B=\hat{P}_s\hat{C}A$ and $B=\hat{T}_d\hat{C}A$, the models are related to the real and complex nonlocal Boussinesq systems.
For the (local and nonlocal, real and complex) KdV type equations, their integrability is guaranteed by the Lax pair \eqref{Lx} and \eqref{Lt}. For the (local and nonlocal, real and complex) Boussinesq type systems, their Lax pairs are the same as that given in \eqref{Lax} and \eqref{Lat}.

It is well known that in quantum physics, the existence of symmetries for a quantum system will lead to some prohibitions. It is interesting that the similar situation may also be found in classical physics. For the ABB systems, the model \eqref{ABKdVA} with $g_i \in \{\hat{P}_s,\ \hat{T}_d,\ \hat{P}_s\hat{C},\  \hat{T}_d\hat{C}\}$, because of the introduces  of the nonlocalities, some prohibitions are discovered. The number of solitons must be even. The odd number of solitons are prohibited. For the even number of solitons, all the solitons must be paired and every paired solitons must possess the same wave numbers, velocities but opposite moving directions. Only the head on interactions are allowed, which means the pursuant interactions are prohibited.

For the ABKdV cases, $g_i \in \{\hat{P}_s\hat{T}_d,\ \hat{P}_s\hat{T}_d\hat{C}\}$, different from the usual KdV case, for every known KdV solution, an arbitrary additional background wave described by an odd function of $x+t$ is allowed. For instance, the periodic wave and the few cycle soliton can be included as shown in Figs. 4 and 5.

Similar phenomena may be found for other kind of nonlinear nonlocal systems. For instance,
we may have a conjecture ($u\equiv A+B,\ v\equiv A-B,\ {\mbox i}=\sqrt{-1}$),
the following sixteen local and nonlocal models given by
\begin{eqnarray}
&&A_t+\frac12\left(2{\mbox i}uv-2u^3+u_{xx}-v\right)_x+2{\mbox i}uv_x=0, \label{mbq}\\
&&B=g_jA,\quad j=0,\ 1,\ 2,\ \ldots,\ 16,
\end{eqnarray}
are integrable,
where $g_j$ belongs to the sixteen order Abel group ${\cal G}$ with four generators, the shifted parity $\hat{P}_s$, the delayed time reversal $\hat{T}_d$, the charge conjugate $\hat{C}$ and the field reversal $\hat{F}$ defined by $\hat{F}A=-A$. The model \eqref{mbq} includes some interesting local and nonlocal systems related to the modified KdV and modified Boussinesq systems. The more about \eqref{mbq} and other types of possible nonlocal systems will be discussed elsewhere.

\section*{Acknowledgements}
The author is grateful  to thank Professors X. Y. Tang, D. J. Zhang, Z. N. Zhu, Q. P. Liu, X. B. Hu, Y. Q. Li and Y. Chen for their helpful discussions. The work was sponsored by the Global Change Research
Program of China (No.2015CB953904),  Shanghai Knowledge Service Platform for Trustworthy Internet of Things (No. ZF1213), the National Natural Science Foundations of China (No. 11435005) and K. C. Wong Magna Fund in Ningbo University.


text:
\begin{document}
\title{Computing recurrence coefficients of multiple orthogonal polynomials\footnote{GF is supported by MNiSzW Iuventus Plus grant Nr 0124/IP3/2011/71, MH and WVA are supported by FWO research project G.0934.13, KU Leuven research grant OT/12/073 and the Belgian Interuniversity Attraction Poles Programme P7/18.}}
\author{Galina Filipuk\footnote{Faculty of Mathematics, Informatics and Mechanics, University of Warsaw, Poland. \newline
filipuk@mimuw.edu.pl},
\ Maciej Haneczok\footnotemark[3],\ Walter Van Assche\footnote{Department of Mathematics, KU Leuven, Belgium. \newline
maciej.haneczok@wis.kuleuven.be; walter.vanassche@wis.kuleuven.be}}
\date{}
\maketitle

\begin{abstract}
Multiple orthogonal polynomials satisfy a number of recurrence relations, in particular there is a $(r+2)$-term recurrence relation connecting the type II multiple orthogonal polynomials near the diagonal (the so-called step-line recurrence relation) and there is a system of $r$ recurrence relations connecting the nearest neighbors (the so-called nearest neighbor recurrence relations). In this paper we deal with two problems. First we show how one can obtain the nearest neighbor recurrence coefficients (and in particular the recurrence coefficients of the orthogonal polynomials for each of the defining measures) from the step-line recurrence coefficients. Secondly we show how one can compute the step-line recurrence coefficients from the recurrence coefficients of the orthogonal polynomials of each of the measures defining the multiple orthogonality.
\end{abstract}

\section{Introduction}

Multiple orthogonal polynomials are polynomials of one variable that satisfy orthogonality conditions with respect to $r>1$ positive measures.
In this paper we will only consider positive measures on the real line. Let $\vec{n}=(n_1,\ldots,n_r) \in \mathbb{N}^r$ be a multi-index and
$|\vec{n}| = n_1 + \cdots + n_r$ its length and let $\mu_1,\ldots,\mu_r$ be $r$ positive measures on the real line. There are two types of multiple orthogonal polynomials \cite[Chapter 23]{Ismail}, \cite{ABVA}, \cite{ACVA}, \cite{ECVA}. Type I multiple orthogonal polynomials are such that $(A_{\vec{n},1},\ldots,A_{\vec{n},r})$ is a vector of $r$ polynomials, with $\deg A_{\vec{n},j} \leq n_j-1$, satisfying
\begin{equation}  \label{TypeIorth}
   \int x^k \sum_{j=1}^r A_{\vec{n},j}(x) \, d\mu_j(x) = 0, \qquad 0 \leq k \leq |\vec{n}|-2,
\end{equation}
with normalization
\[   \int x^{|\vec{n}|-1} \sum_{j=1}^r A_{\vec{n},j}(x) \, d\mu_j(x) = 1.  \]
Type II multiple orthogonal polynomials are monic polynomials $P_{\vec{n}}$ of degree $|\vec{n}|$ for which
\begin{equation} \label{TypeIIorth}
   \int x^k P_{\vec{n}}(x) \, d\mu_j(x) = 0, \qquad 0 \leq k \leq n_j-1
\end{equation}
holds for $1 \leq j \leq r$. The existence (and unicity) is not guaranteed, but if the type I and type II multiple polynomials exist with the above normalization, then they are unique and then the multi-index $\vec{n}$ is said to be a normal index. The measures $(\mu_1,\ldots,\mu_r)$ are a normal system
if all the multi-indices are normal.

Multiple orthogonal polynomials satisfy a number of recurrence relations. Let
$p_{rn}(x) = P_{(n,n,\ldots,n)}(x)$, $p_{rn+1}(x) = P_{(n+1,n,n,\ldots,n)}(x)$, and in general
$p_{rn+j}(x) = P_{(n+1,\ldots,n+1,n\ldots,n)}(x)$ for $0 \leq j \leq r$, where $(n+1,\ldots,n+1,n\ldots,n)$ has $j$ times the component $n+1$ and
$r-j$ times the component $n$, i.e.,
\[  (n+1,\ldots,n+1,n\ldots,n) =  (n,n,\ldots,n) + \sum_{i=1}^j \vec{e}_i, \]
where $\vec{e}_i$ are the standard unit vectors in $\mathbb{N}^r$. Then the \textit{step-line polynomials} $p_{m}(x)$ satisfy the following
$(r+2)$-term recurrence relation
\begin{equation}  \label{stepr-rec}
    xp_{m}(x) = p_{m+1}(x) + \sum_{j=0}^r \beta_{m,j} p_{m-j}(x),
\end{equation}
where $\beta_{m,j}$ are real recurrence coefficients and $p_0=1$ and $p_{-j}=0$ for $1 \leq j \leq r$. This recurrence relation corresponds
to the well known three term recurrence relation for orthogonal polynomials when $r=1$.
These step-line polynomials (the type II multiple orthogonal polynomials near the diagonal) are also known as $d$-orthogonal polynomials
(with $d=r$) and the orthogonality in \eqref{TypeIIorth} becomes
\[    \int p_m(x) x^k\, d\mu_j(x) = 0, \qquad 0 \leq k \leq \lfloor (m-j)/r \rfloor. \]
See, e.g., \cite{Maroni}, \cite{DouakMaroni}, \cite{BenCheikhDouak}. This recurrence relation only connects multiple orthogonal polynomials
near the diagonal $(n,n,\ldots,n)$.
All multiple orthogonal polynomials (of a normal system) are related by a system of $r$ recurrence relations relating $P_{\vec{n}}$ with
its nearest neighbors $P_{\vec{n}+\vec{e}_k}$ and $P_{\vec{n}-\vec{e}_j}$. The system of \textit{nearest neighbor recurrence relations}
is given by (see \cite{WVA})
\begin{equation}  \label{NNrec}
   x P_{\vec{n}}(x) = P_{\vec{n}+\vec{e}_k}(x) + b_{\vec{n},k} P_{\vec{n}}(x) + \sum_{j=1}^r a_{\vec{n},j} P_{\vec{n}-\vec{e}_j}(x), \qquad
    1 \leq k \leq r.
\end{equation}

In this paper we deal with two problems. For the first problem we assume that the recurrence coefficients $\beta_{m,j}$, $m \in \mathbb{N},\  0 \leq j \leq r$, in the step-line recurrence relation \eqref{stepr-rec} are given, and the goal is to find all the nearest neighbor recurrence coefficients
$a_{\vec{n},j}, b_{\vec{n},j}$ with $\vec{n} \in \mathbb{N}^r$ and $1 \leq j \leq r$. In particular this would give all the recurrence coefficients
of the (monic) orthogonal polynomials for each of the measures $\mu_j$:
\[   \int P_n(x;\mu_j)P_m(x;\mu_j)\, d\mu_j(x) = \gamma_n^{-2}(\mu_j) \delta_{m,n}, \]
which satisfy the three term recurrence relation
\[   xP_n(x;\mu_j) = P_{n+1}(x;\mu_j) + b_n(\mu_j) P_n(x;\mu_j) + a_n^2(\mu_j) P_{n-1}(x;\mu_j). \]
Indeed, one has
\[    b_n(\mu_j) = b_{n\vec{e}_j,j}, \quad   a_n^2(\mu_j) = a_{n\vec{e}_j,j}.  \]
In Section 2 we will show, for $r=2$, how one first can obtain the recurrence relation for the $r$ shifted step-line polynomials, i.e., the multiple orthogonal polynomials with a multi-index $k\vec{e}_j + \vec{n}$, with $k \in \mathbb{N}$ and $1\leq j \leq r$ fixed and $\vec{n}$ a multi-index on the step-line. Then, in a second step, we show how to obtain the nearest neighbor recurrence coefficients from the shifted step-line recurrence coefficients.
We explain the procedure for $r=2$ so that the reasoning is easy to follow and not obscured by the notation.

The second problem is the inverse of the first problem: we assume that the recurrence coefficients $a_{n+1}(\mu_j),b_n(\mu_j)$ of the orthogonal polynomials of each of the measures $\mu_j$, $1 \leq j \leq r$, are given, and we show how one can compute all the nearest neighbor recurrence
coefficients and the step-line recurrence coefficients from this input. In Section 3 we will explain the case $r=2$ in detail. We will also show
how to find the nearest neighbor recurrence coefficients from the marginal recurrence coefficients for general $r$.

In order to find the recurrence coefficients  $\beta_{m,j}$, $m \in \mathbb{N},\  0 \leq j \leq r$ in the step-line recurrence relation, one may use
either the Jacobi-Perron algorithm or the vector QD-algorithm. The Jacobi-Perron algorithm generates a vector continued fraction and was introduced by Jacobi in 1868 \cite{Jacobi} and studied in detail by Perron in 1907 \cite{Perron}. A modern version of the Jacobi-Perron algorithm and its relevance for simultaneous rational approximation of functions can be found in \cite{Parus}. Vector continued fractions and the Jacobi-Perron algorithm are quite popular in number theory to produce simultaneous Diophantine approximations to several real numbers, see the monographs of Schweiger \cite{Schweiger1} \cite{Schweiger2}. Another way to obtain the step-line recurrence coefficients from the moments of the measures $\mu_1,\ldots,\mu_r$ is to use a generalization of the QD-algorithm proposed by Van Iseghem \cite{Iseghem}.
The classical QD-algorithm of Rutishauser \cite{Rutis} can be used to find the recurrence coefficients $a_{n+1}(\mu_j),b_n(\mu_j)$ of the orthogonal polynomials of each of the measures $\mu_j$, but one can also use the (modified) Chebyshev algorithm as described in \cite[\S 2.1.7]{Gautschi}.
In this paper we assume that the step-line recurrence coefficients are given (for the first problem) or we assume that the recurrence coefficients
$a_{n+1}(\mu_j),b_n(\mu_j)$ are given for every measure $\mu_j$ with $1 \leq j \leq r$ (for the second problem).

\section{The recurrence coefficients along the step-line}

In this section we only consider multiple orthogonal polynomials with $r=2$. Hence the (monic) type II multiple orthogonal polynomials
$P_{n,m}$ depend on a multi-index $(n,m) \in \mathbb{N}^2$.
Let $p_{2n}(x) = P_{n,n}(x)$ and $p_{2n+1}(x) = P_{n+1,n}(x)$, then the recurrence relation along the step-line is
\begin{equation}  \label{step-line}
  xp_n(x) = p_{n+1}(x) + \beta_n p_n(x) + \gamma_n p_{n-1}(x) + \delta_n p_{n-2}(x).
\end{equation}
It is important to note that the step-line recurrence coefficients $(\beta_n)_{n \geq 0}, (\gamma_n)_{n\geq 1}$ and $(\delta_n)_{n \geq 2}$
do not determine the measures $\mu_1$ and $\mu_2$ in a unique way, even if we normalize the measures to be probability measures. The first measure
is determined uniquely as a probability measure, but for the second measure one can use any convex combination $\lambda \mu_1 + (1-\lambda) \mu_2$
because
\[   \lambda \int p_n(x) x^k\, d\mu_1(x) +  (1-\lambda) \int p_n(x) x^k d\mu_2(x) = 0, \qquad k \leq \lfloor \frac{n}{2} \rfloor -1, \]
since the first of these integrals vanishes for $k \leq \lfloor \frac{n-1}2 \rfloor $ and $\lfloor \frac{n}{2} \rfloor -1 \leq \lfloor \frac{n-1}2 \rfloor$, see \cite[Remark 2.2]{JC}.
This degree of freedom will
be reflected when we want to compute the recurrence coefficients of the orthogonal polynomials $p_n(x;\mu_1)$ and $p_n(x;\mu_2)$.

The nearest neighbor recurrence relations are
\begin{eqnarray}
   x P_{n,m}(x) &=& P_{n+1,m}(x) + c_{n,m} P_{n,m}(x) + a_{n,m} P_{n-1,m}(x) + b_{n,m} P_{n,m-1}(x), \label{NN1} \\
   x P_{n,m}(x) &=& P_{n,m+1}(x) + d_{n,m} P_{n,m}(x) + a_{n,m} P_{n-1,m}(x) + b_{n,m} P_{n,m-1}(x). \label{NN2}
\end{eqnarray}
As a consequence one has
\begin{equation}  \label{P-P}
    P_{n+1,m}(x) - P_{n,m+1}(x) = \kappa_{n,m} P_{n,m}(x),
\end{equation}
where $\kappa_{n,m} = d_{n,m}-c_{n,m}$.

\subsection{From step-line to shifted step-line}

We introduce for $j \geq 0$ the polynomials
\begin{equation}  \label{shiftj}
 p_{2n+j}^{(j,0)}(x) = P_{n+j,n}(x), \quad p_{2n+j+1}^{(j,0)}(x) = P_{n+j+1,n}(x),
\end{equation}
and for $k \geq 0$
\begin{equation}  \label{shiftk}
 p_{2n+k}^{(0,k)}(x) = P_{n,n+k}(x), \quad p_{2n+k+1}^{(0,k)}(x) = P_{n+1,n+k}(x).
\end{equation}
The polynomials $p_n^{(j,0)}$ are the multiple orthogonal polynomials on a shifted step-line with a shift $j$ in the direction of $\vec{e}_1 = (1,0)$. The polynomials $p_n^{(0,k)}$ are those on the shifted step-line with a shift $k$ in the direction of $\vec{e}_2 = (0,1)$.

\begin{figure}[ht]
\centering
\begin{tikzpicture}[domain=0:2]
\draw[->] (0,0) -- (4.5,0) node[below right] {$n$};
\draw[->] (0,0) -- (0,3.5) node[left] {$m$};
\draw[line width=0.2mm] (3,3) circle (0.7ex) node[above]{$P_{n,n}$};

\draw[line width=0.9mm, color=blue] (0,0) -- (0.5,0) -- (0.5,0.5) -- (1,0.5) -- (1,1) -- (1.5,1) -- (1.5,1.5)  -- (2,1.5) -- (2,2) -- (2.5,2) -- (2.5,2.5) -- (3,2.5) -- (3,3);

\draw[->, line width=0.4mm] (1.7,1) -- (2.3,1);
\draw[-, color=black, dashed] (1.5,0) -- (2,0) -- (2,0.5) -- (2.5,0.5) -- (2.5,1) -- (3,1) -- (3,1.5) -- (3.5,1.5) -- (3.5,2) -- (4,2) -- (4,2.5);
\draw[-, color=black, dashed] (2.5,0) -- (3,0) -- (3,0.5) -- (3.5,0.5) -- (3.5,1) -- (4,1) -- (4,1.5) -- (4.5,1.5) -- (4.5,2) -- (5,2) -- (5,2.5);

\draw[->, line width=0.4mm] (1,1.2) -- (1,1.8);
\draw[-, color=black, dashed] (0,1) -- (0,1.5) -- (0.5,1.5) -- (0.5,2) -- (1,2) -- (1,2.5) -- (1.5,2.5)  -- (1.5,3) -- (2,3) -- (2,3.5) --
(2.5,3.5) -- (2.5,4);

\draw[->] (3,3) -- (3.5,3) node[right] {$P_{n+1,n}$};
\end{tikzpicture}
\caption{Step-line and shifted step-lines}
\label{fig:step-lines}
\end{figure}

These shifted step-line polynomials again satisfy a four term recurrence relation, which we denote by
\begin{equation}  \label{step-line+j}
  xp_n^{(j,0)}(x) = p_{n+1}^{(j,0)}(x) + \beta_n^{(j,0)} p_n^{(j,0)}(x) + \gamma_n^{(j,0)} p_{n-1}^{(j,0)}(x) + \delta_n^{(j,0)} p_{n-2}^{(j,0)}(x),
\end{equation}
with initial conditions $\gamma _j^{(j,0)} = \delta _j^{(j,0)} = \delta _{j+1} ^{(j,0)}=0$ and in a similar way
\begin{equation}  \label{step-line+k}
  xp_n^{(0,k)}(x) = p_{n+1}^{(0,k)}(x) + \beta_n^{(0,k)} p_n^{(0,k)}(x) + \gamma_n^{(0,k)} p_{n-1}^{(0,k)}(x) + \delta_n^{(0,k)} p_{n-2}^{(0,k)}(x),
\end{equation}
with initial conditions $\gamma _k^{(0,k+1)} = \delta _k^{(0,k+1)} = \delta _{k+1} ^{(0,k+1)}=0$.
With this notation we have $\beta_n^{(0,0)} = \beta_n$, $\gamma_n^{(0,0)}=\gamma_n$, and $\delta_n^{(0,0)}=\delta_n$. We introduce two more sequences
$(c_n^{(j,0)})_{n\geq 0}$ and $(c_n^{(0,k)})_{n \geq 0}$ by
\begin{equation}  \label{c(j,0)}
      P_{n+j+1,n}(x) - P_{n+j,n+1}(x) = c_n^{(j,0)} P_{n+j,n}(x),
\end{equation}
and
\begin{equation}  \label{c(0,k)}
      P_{n+1,n+k}(x) - P_{n,n+k+1}(x) = c_n^{(0,k)} P_{n,n+k}(x),
\end{equation}
so that \eqref{P-P} gives
\begin{equation}  \label{kappa}
 c_n^{(j,0)} = \kappa_{n+j,n} = d_{n+j,n}-c_{n+j,n}, \quad c_n^{(0,k)} = \kappa_{n,n+k} = d_{n,n+k}-c_{n,n+k}.
\end{equation}

Our first result shows how one can obtain the recurrence coefficients $(\beta_n^{(j,0)},\gamma_n^{(j,0)},\delta_n^{(j,0)})$ from the recurrence coefficients along the step-line.

\begin{theorem}  \label{thm:(j,0)}
One has for all $j \geq 1$
\begin{eqnarray}  \label{2n+j}
   \beta_{2n+j}^{(j,0)} &=&  \beta_{2n+j}^{(j-1,0)} - c_n^{(j,0)}, \qquad n \geq 0,  \nonumber \\
   \gamma_{2n+j}^{(j,0)} &=& \gamma_{2n+j}^{(j-1,0)}, \qquad\qquad n \geq 1,  \\
   \delta_{2n+j}^{(j,0)} &=& \delta_{2n+j}^{(j-1,0)} - c_{n-1}^{(j,0)} \gamma_{2n+j}^{(j-1,0)}, \qquad n \geq 1. \nonumber
\end{eqnarray}
and
\begin{eqnarray}  \label{2n+j+1}
   \beta_{2n+j+1}^{(j,0)} &=&  \beta_{2n+j+1}^{(j-1,0)} + c_n^{(j,0)}, \qquad n \geq 0, \nonumber \\
   \gamma_{2n+j+1}^{(j,0)} &=& \gamma_{2n+j+1}^{(j-1,0)} + c_n^{(j,0)} (\beta_{2n+j}^{(j-1,0)} - \beta_{2n+j+1}^{(j,0)}), \qquad n \geq 0, \\
   \delta_{2n+j+1}^{(j,0)} &=& \delta_{2n+j+1}^{(j-1,0)} + c_n^{(j,0)} \gamma_{2n+j}^{(j-1,0)}, \qquad n \geq 0,  \nonumber
\end{eqnarray}
where $(c_n^{(j,0)})_{n \geq 0}$ is the solution of the Riccati type difference equation
\begin{equation}   \label{cn(j,0)-Riccati}
   c_n^{(j,0)} = \frac{c_{n-1}^{(j,0)}\delta_{2n+j+1}^{(j-1,0)}}{\delta_{2n+j}^{(j-1,0)} - c_{n-1}^{(j,0)} \gamma_{2n+j}^{(j-1,0)}},
\end{equation}
with initial condition
\begin{equation}   \label{c0(j,0)}
    c_0^{(j,0)} = - \frac{\delta_{j+1}^{(j-1,0)}}{\gamma_j^{(j-1,0)}} .
\end{equation}
\end{theorem}

\begin{proof}
Take the recurrence relation  \eqref{step-line+j} with $n$ replaced by $2n+j$
\[ xp_{2n+j}^{(j,0)}(x) = p_{2n+j+1}^{(j,0)}(x) + \beta_{2n+j}^{(j,0)} p_{2n+j}^{(j,0)}(x) + \gamma_{2n+j}^{(j,0)} p_{2n+j-1}^{(j,0)}(x)
+ \delta_{2n+j}^{(j,0)} p_{2n+j-2}^{(j,0)}(x), \]
use \eqref{shiftj} to find
\begin{equation}  \label{xp}
  xP_{n+j,n}(x) = P_{n+j+1,n}(x) + \beta_{2n+j}^{(j,0)} P_{n+j,n}(x) + \gamma_{2n+j}^{(j,0)} P_{n+j,n-1}(x)
 + \delta_{2n+j}^{(j,0)} P_{n+j-1,n-1}(x).
\end{equation}
Now use \eqref{c(j,0)} to replace $P_{n+j+1,n}(x)$ and (by changing $n$ to $n-1$) $P_{n+j,n-1}(x)$, to find
\begin{eqnarray}  \label{xp1}
 xP_{n+j,n}(x) &=& P_{n+j,n+1}(x) + \left(\beta_{2n+j}^{(j,0)} + c_n^{(j,0)}\right) P_{n+j,n}(x) +  \gamma_{2n+j}^{(j,0)} P_{n+j-1,n}(x) \nonumber \\
  & & +\ \left(\delta_{2n+j}^{(j,0)} + c_{n-1}^{(j,0)}\gamma_{2n+j}^{(j,0)}\right)  P_{n+j-1,n-1}(x).
\end{eqnarray}
On the other hand, we take the recurrence relation   \eqref{step-line+j} for $j-1$
\[   xp_n^{(j-1,0)}(x) = p_{n+1}^{(j-1,0)}(x) + \beta_n^{(j-1,0)} p_n^{(j-1,0)}(x) + \gamma_n^{(j-1,0)} p_{n-1}^{(j-1,0)}(x)
+ \delta_n^{(j-1,0)} p_{n-2}^{(j-1,0)}(x),  \]
replace $n$ by $2n+j$ and use  \eqref{shiftj} to find
\begin{eqnarray}  \label{xp2}
 xP_{n+j,n}(x) &=& P_{n+j,n+1}(x) + \beta_{2n+j}^{(j-1,0)} P_{n+j,n}(x) +  \gamma_{2n+j}^{(j-1,0)} P_{n+j-1,n}(x) \nonumber \\
   & & +\ \delta_{2n+j}^{(j-1,0)} P_{n+j-1,n-1}(x).
\end{eqnarray}
Comparing \eqref{xp1} and \eqref{xp2} gives
\begin{multline*}
   \left(\beta_{2n+j}^{(j,0)} + c_n^{(j,0)}\right) P_{n+j,n}(x) +  \gamma_{2n+j}^{(j,0)} P_{n+j-1,n}(x)
+ \left(\delta_{2n+j}^{(j,0)} + c_{n-1}^{(j,0)}\gamma_{2n+j}^{(j,0)}\right)  P_{n+j-1,n-1}(x) \\
  =  \beta_{2n+j}^{(j-1,0)} P_{n+j,n}(x) +  \gamma_{2n+j}^{(j-1,0)} P_{n+j-1,n}(x)
+ \delta_{2n+j}^{(j-1,0)} P_{n+j-1,n-1}(x).
\end{multline*}
The polynomials $P_{n+j,n}, P_{n+j-1,n}, P_{n+j-1,n-1}$ are linearly independent (since they have degrees $2n+j, 2n+j-1$ and $2n+j-2$), hence
one finds
\begin{eqnarray*}
  \beta_{2n+j}^{(j,0)} + c_n^{(j,0)} &=& \beta_{2n+j}^{(j-1,0)}, \\
  \gamma_{2n+j}^{(j,0)} &=& \gamma_{2n+j}^{(j-1,0)} , \\
  \delta_{2n+j}^{(j,0)} + c_{n-1}^{(j,0)}\gamma_{2n+j}^{(j,0)} & =& \delta_{2n+j}^{(j-1,0)},
\end{eqnarray*}
which gives the required relations \eqref{2n+j}.

In a similar way we start with the recurrence relation  \eqref{step-line+j} with $n$ replaced by $2n+j+1$
\[ xp_{2n+j+1}^{(j,0)}(x) = p_{2n+j+2}^{(j,0)}(x) + \beta_{2n+j+1}^{(j,0)} p_{2n+j+1}^{(j,0)}(x) + \gamma_{2n+j+1}^{(j,0)} p_{2n+j}^{(j,0)}(x)
+ \delta_{2n+j+1}^{(j,0)} p_{2n+j-1}^{(j,0)}(x), \]
and use \eqref{shiftj} to find
\[  xP_{n+j+1,n}(x) = P_{n+j+1,n+1}(x) + \beta_{2n+j+1}^{(j,0)} P_{n+j+1,n}(x) + \gamma_{2n+j+1}^{(j,0)} P_{n+j,n}(x)
 + \delta_{2n+j+1}^{(j,0)} P_{n+j,n-1}(x). \]
Use \eqref{c(j,0)} to replace $P_{n+j+1,n}(x)$ and $P_{n+j,n-1}(x)$, to find
\begin{eqnarray*}
  x P_{n+j,n+1}(x) + c_n^{(j,0)} x P_{n+j,n}(x) &=& P_{n+j+1,n+1}(x) + \beta_{2n+j+1}^{(j,0)} P_{n+j,n+1}(x) \\
   & & +\ \left(\gamma_{2n+j+1}^{(j,0)} + c_n^{(j,0)} \beta_{2n+j+1}^{(j,0)}\right) P_{n+j,n}(x) \\
   & & +\  \delta_{2n+j+1}^{(j,0)} P_{n+j-1,n}(x)
   + c_{n-1}^{(j,0)} \delta_{2n+j+1}^{(j,0)} P_{n+j-1,n-1}(x) .
\end{eqnarray*}
Use \eqref{xp2} to replace $xP_{n+j,n}(x)$ to find
\begin{eqnarray}  \label{xp3}
x P_{n+j,n+1}(x) &=&  P_{n+j+1,n+1}(x) + \left(\beta_{2n+j+1}^{(j,0)} - c_n^{(j,0)} \right)  P_{n+j,n+1}(x) \nonumber \\
   & & +\ \left(\gamma_{2n+j+1}^{(j,0)} + c_n^{(j,0)} \beta_{2n+j+1}^{(j,0)} - c_n^{(j,0)} \beta_{2n+j}^{(j-1,0)} \right) P_{n+j,n}(x) \nonumber \\
   & & +\  \left( \delta_{2n+j+1}^{(j,0)} - c_n^{(j,0)} \gamma_{2n+j}^{(j-1,0)} \right) P_{n+j-1,n}(x) \nonumber \\
   & & +\ \left( c_{n-1}^{(j,0)} \delta_{2n+j+1}^{(j,0)} - c_n^{(j,0)} \delta_{2n+j}^{(j-1,0)} \right) P_{n+j-1,n-1}(x) .
\end{eqnarray}
On the other hand, we take the recurrence relation \eqref{step-line+j} for $j-1$ with $n$ replaced by $2n+j+1$ and use
\eqref{shiftj} to find
\begin{eqnarray}  \label{xp4}
  xP_{n+j,n+1}(x) &=& P_{n+j+1,n+1}(x) + \beta_{2n+j+1}^{(j-1,0)} P_{n+j,n+1}(x) \nonumber \\
                  & & +\ \gamma_{2n+j+1}^{(j-1,0)} P_{n+j,n}(x) + \delta_{2n+j+1}^{(j-1,0)} P_{n+j-1,n}(x).
\end{eqnarray}
Comparing \eqref{xp3} and \eqref{xp4} then gives
\begin{multline*}
 \left(\beta_{2n+j+1}^{(j,0)} - c_n^{(j,0)} \right)  P_{n+j,n+1}(x) + \left(\gamma_{2n+j+1}^{(j,0)} + c_n^{(j,0)} \beta_{2n+j+1}^{(j,0)} - c_n^{(j,0)} \beta_{2n+j}^{(j-1,0)} \right) P_{n+j,n}(x) \\
  +  \left( \delta_{2n+j+1}^{(j,0)} - c_n^{(j,0)} \gamma_{2n+j}^{(j-1,0)} \right) P_{n+j-1,n}(x)
 + \left( c_{n-1}^{(j,0)} \delta_{2n+j+1}^{(j,0)} - c_n^{(j,0)} \delta_{2n+j}^{(j-1,0)} \right) P_{n+j-1,n-1}(x) \\
  = \beta_{2n+j+1}^{(j-1,0)} P_{n+j,n+1}(x) + \gamma_{2n+j+1}^{(j-1,0)} P_{n+j,n}(x) + \delta_{2n+j+1}^{(j-1,0)} P_{n+j-1,n}(x) .
\end{multline*}
The four polynomials $P_{n+j,n+1}, P_{n+j,n}, P_{n+j-1,n}, P_{n+j-1,n-1}$ are linearly independent, hence one finds
\begin{eqnarray}  \label{extra}
  \beta_{2n+j+1}^{(j,0)} - c_n^{(j,0)} & = & \beta_{2n+j+1}^{(j-1,0)}, \nonumber \\
  \gamma_{2n+j+1}^{(j,0)} + c_n^{(j,0)} \beta_{2n+j+1}^{(j,0)} - c_n^{(j,0)} \beta_{2n+j}^{(j-1,0)} & = & \gamma_{2n+j+1}^{(j-1,0)}, \nonumber \\
  \delta_{2n+j+1}^{(j,0)} - c_n^{(j,0)} \gamma_{2n+j}^{(j-1,0)} & =& \delta_{2n+j+1}^{(j-1,0)}, \nonumber \\
  c_{n-1}^{(j,0)} \delta_{2n+j+1}^{(j,0)} - c_n^{(j,0)} \delta_{2n+j}^{(j-1,0)} & = & 0.
\end{eqnarray}
The first three relations give \eqref{2n+j+1} and the last equation gives \eqref{cn(j,0)-Riccati} if we replace $\delta_{2n+j+1}^{(j,0)}$
by the third equation in \eqref{2n+j+1}. For $n=0$ the relation \eqref{xp3} becomes
\begin{eqnarray*}
x P_{j,1}(x) &=&  P_{j+1,1}(x) + \left(\beta_{j+1}^{(j,0)} - c_0^{(j,0)} \right)  P_{j,1}(x) \nonumber \\
   & & +\ \left(\gamma_{j+1}^{(j,0)} + c_0^{(j,0)} \beta_{j+1}^{(j,0)} - c_0^{(j,0)} \beta_{j}^{(j-1,0)} \right) P_{j,0}(x) \nonumber \\
   & &  -\  c_0^{(j,0)} \gamma_{j}^{(j-1,0)}  P_{j-1,0}(x) ,
\end{eqnarray*}
and if we compare the coefficient of $P_{j-1,0}(x)$ with the corresponding coefficient in \eqref{xp4} when $n=0$, then
\[   - c_0^{(j,0)} \gamma_{j}^{(j-1,0)} = \delta_{j+1}^{(j-1,0)}, \]
which gives \eqref{c0(j,0)}.
\end{proof}

The Riccati equation \eqref{cn(j,0)-Riccati} can be solved explicitly if all the step-line coefficients at shift $j-1$ are known.
The substitution $c_n^{(j,0)} = 1/d_n^{(j,0)}$ gives
\[    d_n^{(j,0)} = \frac{d_{n-1}^{(j,0)} \delta_{2n+j}^{(j-1,0)} - \gamma_{2n+j}^{(j-1,0)}}{\delta_{2n+j+1}^{(j-1,0)}}, \]
which is a first order linear recurrence relation. Its solution is
\begin{equation}  \label{d(j,0)}
d_n ^{(j,0)}  = \left( \sum_{i=1}^{n} \frac{\gamma _{2i+j} ^{(j-1,0)}} { \delta _{2i+j+1} ^{(j-1,0)} } \prod _{k=1}^{i} \frac {\delta _{2k+j+1}^{(j-1,0)}} {\delta _{2k+j}^{(j-1,0)}} + d_0^{(j,0)} \right) \prod _{k=1}^{n} \frac {\delta _{2k+j}^{(j-1,0)}}{\delta _{2k+j+1}^{(j-1,0)}} .
\end{equation}
However, this is only useful if one has explicit expressions for the step-line coefficients.

An algorithm for computing the recurrence coefficients for the shifted step-line is as follows. Assume that all the recurrence coefficients
$\beta_n^{(j-1,0)},\gamma_n^{(j-1,0)}, \delta_n^{(j-1,0)}$ are known, then one first computes the auxiliary sequence $(c_n^{(j,0)})_{n \geq 0}$
recursively by using \eqref{cn(j,0)-Riccati}. Once this is done, one uses the relations \eqref{2n+j} and \eqref{2n+j+1} to get the
recurrence coefficients for the shifted step-line with shift $j$. The following Maple procedure computes
$(\beta_n^{(j,0)})_{j \leq n \leq 2N-j+1}, (\gamma_n^{(j,0)})_{j+1 \leq n \leq 2N-j+1}$ and $(\delta_n^{(j,0)})_{j+2 \leq n \leq 2N-j+1}$ for $1 \leq j \leq J$.

\begin{verbatim}
shift_j:=proc(N,J)
  local n,j;
  for n from 0 to 2*N+J+1 do
    b[n](0):=beta0(n);
    g[n](0):=gamma0(n);
    d[n](0):=delta0(n);
  end do;
  for j from 1 to J do
    g[j](j):=0;
    d[j](j):=0;
    d[j+1](j):=0;
    c[0](j):=-d[j+1](j-1)/g[j](j-1);
    b[j](j):=b[j](j-1)-c[0](j);
    b[j+1](j):=b[j+1](j-1)+c[0](j);
    g[j+1](j):=g[j+1](j-1)+c[0](j)*(b[j](j-1)-b[j+1](j));
    for n from 1 to N-j do
      c[n](j):=c[n-1](j)*d[2*n+j+1](j-1)/(d[2*n+j](j-1)-c[n-1](j)*g[2*n+j](j-1));
      b[2*n+j](j):=b[2*n+j](j-1)-c[n](j);
      g[2*n+j](j):=g[2*n+j](j-1);
      d[2*n+j](j):=d[2*n+j](j-1)-c[n-1](j)*g[2*n+j](j-1);
      b[2*n+j+1](j):=b[2*n+j+1](j-1)+c[n](j);
      g[2*n+j+1](j):=g[2*n+j+1](j-1)+c[n](j)*(b[2*n+j](j-1)-b[2*n+j+1](j));
      d[2*n+j+1](j):=d[2*n+j+1](j-1)+c[n](j)*g[2*n+j](j-1);
    end do;
  end do;
end proc;
\end{verbatim}

\noindent\textbf{Remark:}
Observe that \eqref{cn(j,0)-Riccati} and the third equation in \eqref{2n+j} give
\[   c_n^{(j,0)} =c_{n-1}^{(j,0)} \frac{\delta_{2n+j+1}^{(j-1,0)}}{\delta_{2n+j}^{(j,0)}},  \]
and \eqref{extra} gives
\[   c_n^{(j,0)} =c_{n-1}^{(j,0)} \frac{\delta_{2n+j+1}^{(j,0)}}{\delta_{2n+j}^{(j-1,0)}}.  \]
Comparing both expressions gives
\[   \delta_{2n+j+1}^{(j,0)} = \frac{\delta_{2n+j+1}^{(j-1,0)}\delta_{2n+j}^{(j-1,0)}}{\delta_{2n+j}^{(j,0)}}.  \]
Hence one may therefore replace lines 16--24 in the Maple procedure by

\begin{verbatim}
    for n from 1 to N-j do
      d[2*n+j](j):=d[2*n+j](j-1)-c[n-1](j)*g[2*n+j](j-1);
      d[2*n+j+1](j):=d[2*n+j+1](j-1)*d[2*n+j](j-1)/d[2*n+j](j);
      c[n](j):=c[n-1](j)*d[2*n+j+1](j-1)/(d[2*n+j](j);
      b[2*n+j](j):=b[2*n+j](j-1)-c[n](j);
      g[2*n+j](j):=g[2*n+j](j-1);
      b[2*n+j+1](j):=b[2*n+j+1](j-1)+c[n](j);
      g[2*n+j+1](j):=g[2*n+j+1](j-1)+c[n](j)*(b[2*n+j](j-1)-b[2*n+j+1](j));
    end do;
\end{verbatim}

The formulas in Theorem \ref{thm:(j,0)} and the computations in the algorithm hold provided $\delta_{2n+j}^{(j,0)} \neq 0$.
If $\delta_{2n+j}^{(j,0)} = 0$ then we cannot compute $c_n^{(j,0)}$ and this quantity is needed in most of the other formulas
for the $j$-shifted step-line. The condition $\delta_{2n+j}^{(j,0)} \neq 0$ for all $j \geq 1$ and $n \geq 1$ is a sufficient condition
that implies that all the required recurrence coefficients can be computed and hence implies that the multi-indices $(n+k,n)$ are
normal.
\medskip

There is a similar result for the recurrence coefficients $(\beta_n^{(0,k)},\gamma_n^{(0,k)},\delta_n^{(0,k)})$.

\begin{theorem}  \label{thm:(0,k)}
One has for all $k \geq 0$
\begin{eqnarray*}
   \beta_{2n+k}^{(0,k+1)} &=&  \beta_{2n+k}^{(0,k)} + c_n^{(0,k)}, \qquad n \geq 0, \\
   \gamma_{2n+k}^{(0,k+1)} &=& \gamma_{2n+k}^{(0,k)}, \qquad\qquad n \geq 1,  \\
   \delta_{2n+k}^{(0,k+1)} &=& \delta_{2n+k}^{(0,k)} + c_{n-1}^{(0,k)} \gamma_{2n+k}^{(0,k)}, \qquad n \geq 1.
\end{eqnarray*}
and
\begin{eqnarray*}
   \beta_{2n+k+1}^{(0,k+1)} &=&  \beta_{2n+k+1}^{(0,k)} - c_n^{(0,k)}, \qquad n \geq 0, \\
   \gamma_{2n+k+1}^{(0,k+1)} &=& \gamma_{2n+k+1}^{(0,k)} - c_n^{(0,k)} (\beta_{2n+k}^{(0,k)} - \beta_{2n+k+1}^{(0,k+1)}), \qquad n \geq 0, \\
   \delta_{2n+k+1}^{(0,k+1)} &=& \delta_{2n+k+1}^{(0,k)} - c_n^{(0,k)} \gamma_{2n+k}^{(0,k)}, \qquad n \geq 1,
\end{eqnarray*}
where $(c_n^{(0,k)})_{n \geq 0}$ is the solution of the Riccati type difference equation
\begin{equation}   \label{cn(0,k)-Riccati}
   c_n^{(0,k)} = \frac{c_{n-1}^{(0,k)}\delta_{2n+k+1}^{(0,k)}}{\delta_{2n+k}^{(0,k)} + c_{n-1}^{(0,k)} \gamma_{2n+k}^{(0,k)}},
\end{equation}
with initial condition
\begin{equation}   \label{c0(0,k)}
    c_0^{(0,k)} =  \frac{\delta_{k+1}^{(0,k)}}{\gamma_k^{(0,k)}} , \qquad k \geq 1,
\end{equation}
and for $k=0$ the $c_0^{(0,0)}$ is a free parameter.
\end{theorem}
\begin{proof}
The proof is very similar to the proof of Theorem \ref{thm:(j,0)}, but one uses the recurrence relation \eqref{step-line+k} and
the relations \eqref{shiftk} and \eqref{c(0,k)}.
\end{proof}

An important difference is that one also needs $c_0^{(0,0)}$ which one can find by taking $n=0$ and $k=0$ in the relation
for $\beta_{2n+k}^{(0,k+1)}$, giving
\[  c_0^{(0,0)} = \beta_0^{(0,1)} - \beta_0^{(0,0)}. \]
Here $\beta_0^{(0,0)} = \beta_0$ is known as the first of the step-line recurrence coefficients, but $\beta_0^{(0,1)}$ can not
be obtained in terms of the step-line recurrence coefficients. Recall that the step-line recurrence coefficients do not determine the measures $\mu_1$ and $\mu_2$ but only determine $\mu_1$ and for the second measure any convex combination of $\mu_1$ and $\mu_2$ is possible. This degree of freedom
is reflected in $c_0^{(0,0)}$ being a free parameter.

Again, the Riccati equation \eqref{cn(0,k)-Riccati} can be solved explicitly. The substitution $c_n^{(0,k)} = 1/d_n^{(0,k)}$ gives
\[  d_n^{(0,k)} = \frac{d_{n-1}^{(0,k)} \delta_{2n+k}^{(0,k)} + \gamma_{2n+k}^{(0,k)}}{\delta_{2n+k+1}^{(0,k)}}, \]
and this first order linear recurrence has the following expression as solution
\begin{equation}  \label{d(0,k)}
d_n ^{(0,k)}  = \left( \sum_{i=1}^{n} \frac{\gamma _{2i+k} ^{(0,k)}} { \delta _{2i+k+1} ^{(0,k)} } \prod _{l=1}^{i} \frac {\delta _{2l+k+1}^{(0,k)}} {\delta _{2l+k}^{(0,k)}} + d_0^{(0,k)} \right) \prod _{l=1}^{n} \frac {\delta _{2l+k}^{(0,k)}}{\delta _{2l+k+1}^{(0,k)}} .
\end{equation}

The following Maple procedure computes $(\beta_n^{(0,k+1)})_{k \leq n \leq 2N-k+1}, (\gamma_n^{(0,k+1)})_{k+1 \leq n \leq 2N-k+1}$ and
$(\delta_n^{(0,k+1)})_{k+2 \leq n \leq 2N-k+1}$ for $0 \leq k \leq K$. It requires as extra input the first two moments of the measure $\mu_2$, i.e.,
$m(1,2)=m_1(\mu_2)$ and $m(0,2)=m_0(\mu_2)$.

\begin{verbatim}
shift_k:=proc(N,K)
  local n,k;
  for n from 0 to 2*N+K+1 do
    b[n](0) := beta0(n);
    g[n](0) := gamma0(n);
    d[n](0) := delta0(n);
  end do;
  c[0](0) := m(1,2)/m(0,2)-b[0](0);
  for k from 0 to K do
    g[k](k+1):=0;
    d[k](k+1):=0;
    d[k+1](k+1):=0;
    b[k](k+1):=b[k](k)+c[0](k);
    b[k+1](k+1):=b[k+1](k)-c[0](k);
    g[k+1](k+1):=g[k+1](k)-c[0](k)*(b[k](k)-b[k+1](k+1));
    for n from 1 to N-k do
      c[n](k):=c[n-1](k)*d[2*n+k+1](k)/(d[2*n+k](k)+c[n-1](k)*g[2*n+k](k));
      b[2*n+k](k+1):=b[2*n+k](k)+c[n](k);
      g[2*n+k](k+1):=g[2*n+k](k);
      d[2*n+k](k+1):=d[2*n+k](k)+c[n-1](k)*g[2*n+k](k);
      b[2*n+k+1](k+1):=b[2*n+k+1](k)-c[n](k);
      g[2*n+k+1](k+1):=g[2*n+k+1](k)-c[n](k)*(b[2*n+k](k)-b[2*n+k+1](k+1));
      d[2*n+k+1](k+1):=d[2*n+k+1](k)-c[n](k)*g[2*n+k](k);
    end do;
    c[0](k+1):=d[k+2](k+1)/g[k+1](k+1);
  end do;
end proc;
\end{verbatim}

\noindent\textbf{Remark:}
One has
\[   c_n^{(0,k)} = c_{n-1}^{(0,k)} \frac{\delta_{2n+k+1}^{(0,k)}}{\delta_{2n+k}^{(0,k+1)}}, \]
and in a way similar as before one has
\[   \delta_{2n+k+1}^{(0,k+1)} = \frac{\delta_{2n+k+1}^{(0,k)}\delta_{2n+k}^{(0,k)}}{\delta_{2n+k}^{(0,k+1)}}. \]
Hence one can replace the lines 16--24 in the above procedure by

\begin{verbatim}
    for n from 1 to N-k do
      d[2*n+k](k+1):=d[2*n+k](k)+c[n-1](k)*g[2*n+k](k);
      d[2*n+k+1](k+1):=d[2*n+k+1](k)*d[2*n+k](k)/d[2*n+k](k+1);
      c[n](k):=c[n-1](k)*d[2*n+k+1](k)/d[2*n+k](k+1);
      b[2*n+k](k+1):=b[2*n+k](k)+c[n](k);
      g[2*n+k](k+1):=g[2*n+k](k);
      b[2*n+k+1](k+1):=b[2*n+k+1](k)-c[n](k);
      g[2*n+k+1](k+1):=g[2*n+k+1](k)-c[n](k)*(b[2*n+k](k)-b[2*n+k+1](k+1));
    end do;
\end{verbatim}

A sufficient condition for normality of the multi-indices $(n,n+k)$ is now that $\delta_{2n+k}^{(0,k+1)} \neq 0$ for all $n \geq 1$ and $k \geq 0$.

\subsection{From shifted step-line to nearest neighbor recurrence coefficients}
Our next step is to find the nearest neighbor recurrence coefficients $a_{n,m},b_{n,m},c_{n,m},d_{n,m}$ in \eqref{NN1}--\eqref{NN2}
from the recurrence coefficients on the shifted step-lines.

\begin{theorem}  \label{thm:NN}
The coefficients of the nearest neighbor recurrence relations \eqref{NN1}--\eqref{NN2} are for $n\geq 0$ and $j \geq 1$ given by
\begin{eqnarray*}  \label{down}
  c_{n+j,n} &=& \beta_{2n+j}^{(j,0)}, \\
  d_{n+j,n} &=& c_{n+j,n} + c_n^{(j,0)}, \\
  a_{n+j,n} &=& - \frac{\delta_{2n+j+1}^{(j-1,0)}}{c_n^{(j,0)}}, \\
  b_{n+j,n} &=& \gamma_{2n+j}^{(j,0)} - a_{n+j,n} .
\end{eqnarray*}
and for $n\geq 1$ and $k \geq 0$
\begin{eqnarray*}
   c_{n,n+k} &=& \beta_{2n+k}^{(0,k)}, \\
   d_{n,n+k} &=& c_{n,n+k} + c_n^{(0,k)}, \\
   a_{n,n+k} &=& - \frac{\delta_{2n+k}^{(0,k)}}{c_{n-1}^{(0,k)}}, \\
   b_{n,n+k} &=& \gamma_{2n+k}^{(0,k)} - a_{n,n+k}.
\end{eqnarray*}
The initial coefficients are $a_{0,0}=b_{0,0}=0$, $c_{0,0} = \beta_0$ and $d_{0,0}$ is a free parameter.
\end{theorem}

Note that if $b_0(\mu_1)=m_1(\mu_1)/m_0(\mu_1)$ and $b_0(\mu_2)=m_1(\mu_2)/m_0(\mu_2)$ are the first recurrence coefficients
of the orthogonal polynomials for $\mu_1$ and $\mu_2$ respectively, then $c_{0,0} = b_0(\mu_1)$ and $d_{0,0}=b_0(\mu_2)$.
Hence the particular choice $d_{0,0}=m_1(\mu_2)/m_0(\mu_2)$ gives the nearest neighbor recurrence coefficients of the
multiple orthogonal polynomials for the measures $\mu_1$ and $\mu_2$.
Another choice of $d_{0,0}$ still gives $\mu_1$ as the first measure, but a linear combination of $\mu_1$ and $\mu_2$ for the second measure.

\begin{proof}
We start from \eqref{xp} and replace the polynomial $P_{n+j-1,n-1}$ by using \eqref{c(j,0)} to find
\begin{multline*}
   xP_{n+j,n}(x) = P_{n+j+1,n}(x) + \beta_{2n+j}^{(j,0)} P_{n+j,n}(x) + \gamma_{2n+j}^{(j,0)} P_{n+j,n-1}(x)  \\
 +\ \frac{\delta_{2n+j}^{(j,0)}}{c_{n-1}^{(j,0)}}  \left( P_{n+j,n-1}(x)-P_{n+j-1,n}(x) \right).
\end{multline*}
Now we compare this with the recurrence relation \eqref{NN1} with $n$ replaced by $n+j$ and $m=n$, then
\begin{eqnarray*}
       c_{n+j,n} &=& \beta_{2n+j}^{(j,0)}, \\
       b_{n+j,n} &=& \gamma_{2n+j}^{(j,0)} + \frac{\delta_{2n+j}^{(j,0)}}{c_{n-1}^{(j,0)}}, \\
       a_{n+j,n} &=& - \frac{\delta_{2n+j}^{(j,0)}}{c_{n-1}^{(j,0)}}.
\end{eqnarray*}
Use the last equation in (\ref{extra}) to obtain
\[ a_{n+j,n} = - \frac{\delta_{2n+j+1} ^{(j-1,0)} }{ c_n ^{(j,0)} }, \]
so that the formula is valid for all $n\geq 0$. For $d_{n+j,n}$ we use the relation \eqref{kappa} to find
\[  d_{n+j,n} = c_{n+j,n} + c_n^{(j,0)}. \]
 This gives the first part of the theorem.

For the second part of the theorem we use \eqref{step-line+k} with $n$ replaced by $2n+k$ and by \eqref{shiftk} this gives
\begin{equation}  \label{xp0}
  xP_{n,n+k}(x) = P_{n+1,n+k}(x) + \beta_{2n+k}^{(0,k)} P_{n,n+k}(x) + \gamma_{2n+k}^{(0,k)} P_{n,n+k-1}(x)
+ \delta_{2n+k}^{(0,k)} P_{n-1,n+k-1}(x).
\end{equation}
Replace the polynomial $P_{n-1,n+k-1}$ by using \eqref{c(0,k)} (but with $n$ replaced by $n-1$) to find
\begin{multline*}
  xP_{n,n+k}(x) = P_{n+1,n+k}(x) + \beta_{2n+k}^{(0,k)} P_{n,n+k}(x) + \gamma_{2n+k}^{(0,k)} P_{n,n+k-1}(x) \\
  +\   \frac{\delta_{2n+k}^{(0,k)}}{c_{n-1}^{(0,k)}} \left( P_{n,n+k-1}(x) - P_{n-1,n+k}(x) \right).
\end{multline*}
Comparing with \eqref{NN1} with $m=n+k$ then gives
\begin{eqnarray*}
        c_{n,n+k} &=& \beta_{2n+k}^{(0,k)}, \\
        a_{n,n+k} &=& - \frac{\delta_{2n+k}^{(0,k)}}{c_{n-1}^{(0,k)}}, \\
        b_{n,n+k} &=& \gamma_{2n+k}^{(0,k)} + \frac{\delta_{2n+k}^{(0,k)}}{c_{n-1}^{(0,k)}} .
\end{eqnarray*}
For $d_{n,n+k}$ we use \eqref{kappa} to find
\[   d_{n,n+k} = c_{n,n+k} + c_n^{(0,k)}. \]
This gives the second part of the theorem.
\end{proof}

In Maple one can compute the nearest neighbor recurrence coefficients using the following procedure:

\begin{verbatim}
nncoef:=proc(n,m)
   local j,k;
   if m<n then j:=n-m;
             shift_j(n,j);
             c(n,m):= b[2*m+j](j);
             d(n,m):= c(n,m)+c[m](j);
             a(n,m):= -d[2*m+j+1](j-1)/c[m](j);
             b(n,m):= g[2*m+j](j)-a(n,m);
          else k:=m-n;
             shift_k(m,k);
             c(n,m):= b[2*n+k](k);
             d(n,m):= c(n,m)+c[n](k);
             a(n,m):= -d[2*n+k](k)/c[n-1](k);
             b(n,m):= g[2*n+k](k)-a(n,m);
   end if;
   Vector([c(n,m), d(n,m), a(n,m), b(n,m)]);
end proc;
\end{verbatim}

\subsection{The recurrence coefficients of the marginal measures}
Now that we know the nearest neighbor recurrence coefficients, we can find the recurrence coefficients of the orthogonal polynomials
$P_n(x;\mu_1)$ for the measure $\mu_1$ and $P_n(x;\mu_2)$ for the measure $\mu_2$. The recurrence relations for these monic orthogonal
polynomials are
\[   xP_n(x;\mu_i) = P_{n+1}(x;\mu_i) + b_n(\mu_i) P_n (x;\mu_i) + a_n^2(\mu_i) P_{n-1}(x;\mu_i). \]
Observe that $P_n(x;\mu_1) = P_{n,0}(x)$ and $P_m(x;\mu_2) = P_{0,m}(x)$, so if we compare with the nearest neighbor recurrence relations
\eqref{NN1} and  \eqref{NN2} we find
\[    b_j(\mu_1) = c_{j,0}, \quad  a_j^2(\mu_1) = a_{j,0}, \]
and
\[   b_k(\mu_2) = d_{0,k},  \quad a_k^{2}(\mu_2) = b_{0,k}.  \]
If we use Theorem \ref{thm:NN} then this gives
\begin{equation}  \label{recmu1}
    b_j(\mu_1) = \beta_j^{(j,0)}, \quad  a_j^2(\mu_1) = \gamma_j^{(j-1,0)},
\end{equation}
where we used \eqref{c0(j,0)} to simplify the expression for $a_{j,0}$.
The recurrence coefficients for $\mu_2$ can be obtained more easily from \eqref{xp0} with $n=0$, which gives
\[  xP_{0,k}(x) = P_{1,k}(x) + \beta_k^{(0,k)} P_{0,k}(x) + \gamma_k^{(0,k)} P_{0,k-1}(x), \]
and if we replace $P_{1,k}$ by using \eqref{c(0,k)} with $n=0$, then
\[   xP_{0,k}(x) = P_{0,k+1}(x) + \left( c_0^{(0,k)} + \beta_k^{(0,k)} \right) P_{0,k}(x) + \gamma_k^{(0,k)} P_{0,k-1}(x), \]
so that
\begin{equation}  \label{recmu2}
    b_k(\mu_2) = \beta_k^{(0,k)} + c_0^{(0,k)}, \quad
      a_k^2(\mu_2) = \gamma_k^{(0,k)} .
\end{equation}

Observe that these results  allow us to find the recurrence coefficients of the orthogonal polynomials for the measures $\mu_1$ and $\mu_2$ if the
recurrence coefficients of the step-line multiple orthogonal polynomials are known. There are examples of multiple orthogonal polynomials for which the recurrence coefficients of the marginal orthogonal polynomials are not known. Such a situation occurs for instance in the case of multiple orthogonal polynomials associated to the modified Bessel functions of the first and second kind. In 1990 A.P. Prudnikov posed an open problem to find the orthogonal polynomials for the modified Bessel functions of the second kind $K_\nu(2\sqrt{x})$ on $[0,\infty)$ (see \cite[Problem 9 on pp. 239--241]{WVA-open}). It turned out that in this case it is more natural to consider multiple orthogonal polynomials for a pair of modified Bessel functions of the second kind $K_\nu(2\sqrt{x}),K_{\nu+1}(2\sqrt{x})$. This was shown by Van Assche and Yakubovich in \cite{WVA-Yaku} (see also \cite{BenCheikhDouak0}) and later for the modified Bessel functions of the first kind by Coussement and Van Assche in \cite{EC2} (see also \cite{Douak}).
Our algorithm allows us to find the recurrence coefficients of those polynomials, even though we are not able to find explicit expressions for them.
If we start from the step-line recurrence coefficients given in \cite[Thm.~4]{WVA-Yaku}
\begin{eqnarray*}
  \beta_n &=& (n+\alpha+1)(3n+\alpha+2\nu) - (\alpha+1)(\nu-1), \\
  \gamma_n &=& n(n+\alpha)(n+\alpha+\nu)(3n+2\alpha+\nu), \\
  \delta_n &=& n(n-1)(n+\alpha)(n+\alpha-1)(n+\alpha+\nu)(n+\alpha+\nu-1),
\end{eqnarray*}
for the multiple orthogonal polynomials with
\[  d\mu_1(x)=x^{\alpha+\nu/2} K_\nu(2\sqrt{x}), \quad d\mu_2(x) = x^{\alpha+(\nu+1)/2} K_{\nu+1}(2\sqrt{x}), \]
and put $\alpha=0$ and $\nu=0$, then Table \ref{K0} gives the results of our algorithm for the recurrence coefficients
$(a_{n},b_n)$ of the weight $K_0(2\sqrt{x})$.

\begin{table}[h]\footnotesize
\centering
\begin{tabular}{|r|c|c|}
\hline
$n$ & $a_n$ & $b_n$  \\
\hline
0 & --  & 1 \\
1 & 1.7320508075688772935 & 9.6666666666666666667 \\
2 & 8.5374989832437982487 & 28.186991869918699187 \\
3 & 20.265386777687130909 & 56.571895845674401834 \\
4 & 36.925214834648582674 & 94.823932737801348717 \\
5 & 58.518554562959399225 & 142.94410230778264607 \\
6 & 85.045955898223602580 & 200.93289913274452209 \\
7 & 116.50767686120789662 & 268.79060407933245800 \\
8 & 152.90385976282648737 & 346.51739199614374938 \\
9 & 194.23459164836084172 & 434.11337913848760712 \\
10 & 240.49992974325090503 & 531.57864673346522330 \\
\hline
\end{tabular}
\caption{Recurrence coefficients of orthogonal polynomials for $K_0(2\sqrt{x})$}
\label{K0}
\end{table}
Note that the values for the coefficients $a_n$ presented in Table \ref{K0} are the same as the ones that were computed in \cite{WVA-open}

using the moments.

\section{From marginal to  nearest neighbor}

In the previous section we started from the step-line recurrence coefficients and we showed how to find the nearest neighbor
recurrence coefficients and in particular the recurrence coefficients of the orthogonal polynomials for the marginal measures
$\mu_1$ and $\mu_2$. In this section we will investigate the inverse problem: suppose the recurrence coefficients
$(a_{n+1}^2(\mu_i),b_n(\mu_i))_{n \geq 0}$ are given for $i \in \{1,2\}$. How can one find all the nearest neighbor recurrence
coefficients $(a_{n,m}, b_{n,m}, c_{n,m}, d_{n,m})$ and the step-line recurrence coefficients $(\beta_n, \gamma_n,\delta_n)$?

\subsection{The nearest neighbor recurrence coefficients}

The nearest neighbor recurrence coefficients satisfy a system of non-linear partial difference equations, as was noted in \cite{WVA}.
We will briefly show how to find these partial difference equations. The nearest neighbor recurrence relations \eqref{NN1} and \eqref{NN2} can be written in a matrix form as
\[      Y_{n+1,m} = R_1(n,m) Y_{n,m}, \quad Y_{n,m+1} = R_2(n,m) Y_{n,m}, \]
where
\[   Y_{n,m} =   \begin{pmatrix} P_{n,m}(x) \\ P_{n-1,m}(x) \\ P_{n,m-1}(x) \end{pmatrix}  \]
and $R_1$ and $R_2$ are the two transfer matrices
\[    R_1(n,m) = \begin{pmatrix} x-c_{n,m} & -a_{n,m} & -b_{n,m} \\
                                 1 & 0 & 0 \\
                                 1 & 0 & d_{n,m-1}-c_{n,m-1}  \end{pmatrix},  \]
and
\[    R_2(n,m) = \begin{pmatrix} x-d_{n,m} & -a_{n,m} & -b_{n,m} \\
                                 1 & c_{n-1,m}-d_{n-1,m} & 0 \\
                                 1 & 0 & 0  \end{pmatrix}.  \]
Now there are two ways of finding $Y_{n+1,m+1}$ from $Y_{n,m}$ using these transfer matrices: one way is to first compute $Y_{n+1,m}$ and then to increase $m$ by one, which gives
\[    Y_{n+1,m+1} = R_2(n+1,m)R_1(n,m)Y_{n,m}. \]
Another way is to first compute $Y_{n,m+1}$ and then to increase $n$ by one, giving
\[    Y_{n+1,m+1} = R_1(n,m+1) R_2(n,m) Y_{n,m}.  \]
Comparing both expressions gives the matrix relation
\[     R_2(n+1,m)R_1(n,m) =  R_1(n,m+1) R_2(n,m).  \]
If one computes the entries of this matrix identity, then one finds the following partial difference relations (see also \cite[Thm.~3.1]{WVA}):
\begin{eqnarray}
   d_{n+1,m}-d_{n,m} &=& c_{n,m+1}-c_{n,m} , \label{cd} \\
   a_{n+1,m} + b_{n+1,m} - (a_{n,m+1} + b_{n,m+1}) &=& \det \begin{pmatrix} d_{n+1,m} & d_{n,m} \\ c_{n,m+1} & c_{n,m} \end{pmatrix}, \label{a+b} \\
   \frac{a_{n,m+1}}{a_{n,m}} &=& \frac{c_{n,m}-d_{n,m}}{c_{n-1,m}-d_{n-1,m}}, \label{acd} \\
   \frac{b_{n+1,m}}{b_{n,m}} &=& \frac{c_{n,m}-d_{n,m}}{c_{n,m-1}-d_{n,m-1}}.  \label{bcd}
\end{eqnarray}

We will show that these partial difference equations with boundary conditions
\begin{equation}  \label{BC1}
   a_{n,0} = a_n^2(\mu_1), \quad b_{n,0} = 0, \quad c_{n,0} = b_n(\mu_1), \qquad n \geq 0,
\end{equation}
and
\begin{equation}  \label{BC2}
   a_{0,m} = 0, \quad b_{0,m} = a_m^2(\mu_2), \quad d_{0,m} = b_m(\mu_2), \qquad m \geq 0,
\end{equation}
where $a_n^2(\mu_i),b_n(\mu_i)$ are the recurrence coefficients of the monic orthogonal polynomials for the measure $\mu_i$ $(i=1,2)$ (with
$a_0^2(\mu_1)=a_0^2(\mu_2)=0$),
can be solved recursively to find the nearest neighbor recurrence coefficients for the multiple orthogonal polynomials with the measures $(\mu_1,\mu_2)$.

\begin{theorem}  \label{thm:inverse}
Suppose $(a_n^2(\mu_i))_{n \geq 1}$ and $(b_n(\mu_i))_{n \geq 0}$ are the recurrence coefficients of the monic orthogonal polynomials
$P_n(x;\mu_i)$ for the measure $\mu_i$, i.e.,
\[    xP_n(x;\mu_i) = P_{n+1}(x;\mu_i) + b_n(\mu_i) P_n(x;\mu_i) + a_n^2(\mu_i) P_{n-1}(x;\mu_i), \qquad n \geq 0, \]
with $P_0(x;\mu_i)=1$ and $P_{-1}(x;\mu_i)=0$. Then the nearest neighbor recurrence coefficients for the type II multiple orthogonal polynomials
can be computed recursively by \eqref{cd}--\eqref{bcd}, using the boundary conditions \eqref{BC1}--\eqref{BC2}, provided $c_{n,m} \neq d_{n,m}$
for all $n,m \geq 0$.
\end{theorem}

\begin{proof}
We use induction on $k$, where $k=n+m$ is the length of the multi-index $(n,m)$ and show how to compute $a_{n,m},b_{n,m},c_{n,m},d_{n,m}$
when the recurrence coefficients are known for multi-indices of length less than $k$.
For $k=0$ we have that $a_{0,0}=b_{0,0}=0$ (these appear as coefficients of $P_{-1,0}$ and $P_{0,-1}$ and hence are not needed) and
$c_{0,0}=b_0(\mu_1)$, $d_{0,0}=b_0(\mu_2)$. Hence the case $k=0$ is settled. For $k=1$ we already have
$a_{1,0}=a_1^2(\mu_1),c_{1,0}=b_1(\mu_1), b_{0,1}=a_1^2(\mu_2),d_{0,1}=b_1(\mu_2)$, and we also defined $a_{0,1}=b_{1,0}=0$ (these
appear as coefficients for $P_{-1,1}$ and $P_{1,-1}$ and are not needed). This leaves only to determine $c_{0,1}$ and $d_{1,0}$.
If we use \eqref{cd} and \eqref{a+b} for $n=m=0$, then this gives the system of equations
\begin{eqnarray*}
      d_{1,0}-c_{0,1} &=& d_{0,0}-c_{0,0}, \\
      d_{1,0}c_{0,0} - c_{0,1}d_{0,0} &=& a_{1,0}+b_{1,0}-a_{0,1}-b_{0,1}.
\end{eqnarray*}
This is a linear system of two equations for $c_{0,1}$ and $d_{1,0}$. The determinant of the system is $d_{0,0}-c_{0,0}$ and hence this system has a unique solution whenever $b_0(\mu_1) \neq b_0(\mu_2)$.

\begin{figure}[ht]
\centering
\begin{tikzpicture}[domain=0:2]
\draw[->, color=black] (0,0) -- (3.5,0) node[below right] {$n$};
\draw[->, color=black] (0,0) -- (0,3.5) node[left] {$m$};

\draw[line width=0.5mm, color=blue, dashed] (0,0.5) -- (0.5,0);
\draw[line width=0.5mm, color=blue, dashed] (0,1) -- (1,0);
\draw[line width=0.5mm, color=blue, dashed] (0,1.5) -- (1.5,0);
\draw[->, line width=0.4mm] (1,1) -- (1.3,1.3);

\draw[line width=0.5mm, color=blue, dashed] (0,3) -- (3,0);

\draw[] (-0.3,-0.3)  node[right] {$0$};
\draw[] (0.2,-0.3)  node[right] {$1$};
\draw[] (0.7,-0.3)  node[right] {$2$};
\draw[] (1.2,-0.3)  node[right] {$3$};
\draw[] (1.7,-0.3)  node[right] {$\ldots $};
\end{tikzpicture}
\caption{Moving along the lines $m+n=k$}
\label{fig:moving}
\end{figure}

Suppose next that we know all the nearest neighbor recurrence coefficients for multi-indices $(n,m)$ of length $ \leq k$.
From \eqref{acd} we then find
\begin{equation*}
    a_{n,m+1} = a_{n,m} \frac{c_{n,m}-d_{n,m}}{c_{n-1,m}-d_{n-1,m}},
\end{equation*}
and from \eqref{bcd} we find
\begin{equation*}
    b_{n+1,m} = b_{n,m} \frac{c_{n,m}-d_{n,m}}{c_{n,m-1}-d_{n,m-1}}.
\end{equation*}
If we replace $n$ by $\ell$ and $m$ by $k-\ell$, then this gives
\begin{equation}  \label{anm1}
  a_{\ell,k-\ell+1} = a_{\ell,k-\ell} \frac{c_{\ell,k-\ell}-d_{\ell,k-\ell}}{c_{\ell-1,k-\ell}-d_{\ell-1,k-\ell}}, \qquad 1 \leq \ell \leq k.
\end{equation}
For $\ell=0$ and $\ell=k+1$ we use the boundary conditions
\[ a_{0,k+1} = 0, \quad a_{k+1,0} = a_{k+1}^2(\mu_1). \]
In a similar way
\begin{equation}   \label{bnm1}
  b_{\ell+1,k-\ell} = b_{\ell,k-\ell} \frac{c_{\ell,k-\ell}-d_{\ell,k-\ell}}{c_{\ell,k-\ell-1}-d_{\ell,k-\ell-1}}, \qquad 0 \leq \ell \leq k-1.
\end{equation}
For $\ell = -1$ and $\ell=k$ we use the boundary conditions
\[  b_{0,k+1} = a_{k+1}^2(\mu_2), \quad b_{k+1,0} = 0.  \]
The expressions on the right of \eqref{anm1}--\eqref{bnm1} contain coefficients of multi-indices of length $k$ and $k-1$ and hence they are known.
If we use \eqref{cd} and \eqref{a+b}, then
\begin{eqnarray*}
     d_{n+1,m} - c_{n,m+1} &=& d_{n,m} - c_{n,m} , \\
     d_{n+1,m}c_{n,m} - c_{n,m+1} d_{n,m} &=& a_{n+1,m}+b_{n+1,m}-a_{n,m+1}-b_{n,m+1},
\end{eqnarray*}
which is a linear system for $c_{n,m+1}$ and $d_{n+1,m}$ with determinant $d_{n,m}-c_{n,m}$. This system has a unique solution
whenever $c_{n,m}\neq d_{n,m}$. This solution is
\[   c_{n,m+1} = c_{n,m} + \frac{a_{n+1,m}+b_{n+1,m}-a_{n,m+1}-b_{n,m+1}}{c_{n,m}-d_{n,m}}, \]
and
\[   d_{n+1,m} = d_{n,m} + \frac{a_{n+1,m}+b_{n+1,m}-a_{n,m+1}-b_{n,m+1}}{c_{n,m}-d_{n,m}}.  \]
Replacing $n$ by $\ell$ and $m$ by $k-\ell$ then gives
\begin{equation} \label{cnm1}
    c_{\ell,k-\ell+1} = c_{\ell,k-\ell} + \frac{a_{\ell+1,k-\ell}+b_{\ell+1,k-\ell}-a_{\ell,k-\ell+1}-b_{\ell,k-\ell+1}}
   {c_{\ell,k-\ell}-d_{\ell,k-\ell}},  \qquad 0 \leq \ell \leq k,
\end{equation}
and for $\ell=k+1$ we use the boundary condition
\[     c_{k+1,0} = b_{k+1}(\mu_1). \]
Similarly we have
\begin{equation} \label{dnm1}
    d_{\ell+1,k-\ell} = d_{\ell,k-\ell} + \frac{a_{\ell+1,k-\ell}+b_{\ell+1,k-\ell}-a_{\ell,k-\ell+1}-b_{\ell,k-\ell+1}}
   {c_{\ell,k-\ell}-d_{\ell,k-\ell}},   \qquad 0 \leq \ell \leq k,
\end{equation}
and for $\ell=-1$ we use the boundary condition
\[      d_{0,k+1} = b_{k+1}(\mu_2).  \]
\end{proof}

We can implement this in Maple using the following procedure which computes $a_{n,m},$ $b_{n,m},c_{n,m},d_{n,m}$ for $n+m \leq N+M$. It requires
the input $a1(n) = a_n^2(\mu_1)$ and $b1(n)=b_n(\mu_1)$ for the first measure $\mu_1$ and $a2(n) = a_n^2(\mu_2)$ and $b2(n)=b_n(\mu_2)$
for the second measure $\mu_2$, for $0 \leq n \leq N+M$, where we set $a_0^2(\mu_1)=0=a_0^2(\mu_2)$. In particular it gives the
coefficients $c_{N,M},d_{N,M},a_{N,M},b_{N,M}$ which will be given in the output explicitly.

\begin{verbatim}
IP:=proc(N,M)
 local n,m,k;
 for n from 0 to N+M do
  c(n,0):=b1(n);
  a(n,0):=a1(n);
  b(n,0):=0;
 end do;
 for m from 0 to N+M do
  d(0,m):=b2(m);
  a(0,m):=0;
  b(0,m):=a2(m);
 end do;
 for n from 1 to N+M do
  for k from 1 to n-1 do
    a(k,n-k):=a(k,n-k-1)*(c(k,n-k-1)-d(k,n-k-1))/(c(k-1,n-k-1)-d(k-1,n-k-1));
    b(k,n-k):=b(k-1,n-k)*(c(k-1,n-k)-d(k-1,n-k))/(c(k-1,n-k-1)-d(k-1,n-k-1));
  end do;
  for k from 1 to n do
   c(n-k,k):=c(n-k,k-1)
       +(a(n-k+1,k-1)+b(n-k+1,k-1)-a(n-k,k)-b(n-k,k))/(c(n-k,k-1)-d(n-k,k-1));
  end do;
  for k from 1 to n do
   d(k,n-k):=c(k-1,n-k+1)-c(k-1,n-k)+d(k-1,n-k);
  end do;
 end do;
Vector([c(N,M),d(N,M),a(N,M),b(N,M)]);
end proc;
\end{verbatim}

\subsection{The step-line recurrence coefficients}
Next we will show how to compute the recurrence coefficients in the step-line recurrence \eqref{step-line} if one knows the nearest neighbor
recurrence coefficients.

\begin{theorem}  \label{thm:stepline}
Suppose that the nearest neighbor recurrence coefficients are given. Then the step-line recurrence coefficients in \eqref{step-line}
are given by
\begin{eqnarray*}
     \beta_{2n} &=& c_{n,n},  \qquad n \geq 0,\\
     \gamma_{2n} &=& a_{n,n} + b_{n,n}, \qquad n \geq 1,\\
     \delta_{2n} &=& a_{n,n} ( c_{n-1,n-1} - d_{n-1,n-1}), \qquad n \geq 1,
\end{eqnarray*}
and
\begin{eqnarray*}
     \beta_{2n+1} &=& d_{n+1,n}, \qquad n \geq 0, \\
     \gamma_{2n+1} &=& a_{n+1,n} + b_{n+1,n}, \qquad n \geq 0, \\
     \delta_{2n+1} &=& b_{n+1,n} (d_{n,n-1}-c_{n,n-1}), \qquad n \geq 1.
\end{eqnarray*}
\end{theorem}

\begin{proof}
From the nearest neighbor recurrence relation \eqref{NN1} we find
\[   xP_{n,n}(x) = P_{n+1,n}(x) + c_{n,n} P_{n,n}(x) + a_{n,n} P_{n-1,n}(x) + b_{n,n} P_{n,n-1}(x). \]
Use \eqref{P-P} to replace $P_{n-1,n}$ to find
\begin{multline*}   xP_{n,n}(x) = P_{n+1,n}(x) + c_{n,n} P_{n,n}(x)
    + (a_{n,n}+b_{n,n}) P_{n,n-1}(x)  \\
+ a_{n,n} (c_{n-1,n-1}-d_{n-1,n-1})P_{n-1,n-1}(x).
\end{multline*}
If we compare this with \eqref{step-line} with $n$ replaced by $2n$, then we find the relations for the even recurrence coefficients.
The proof for the odd recurrence coefficients is similar: use \eqref{NN2} to find
\[   xP_{n+1,n}(x) = P_{n+1,n+1}(x) + d_{n+1,n} P_{n+1,n}(x) + a_{n+1,n} P_{n,n}(x) + b_{n+1,n} P_{n+1,n-1}(x), \]
and replace $P_{n+1,n-1}$ using \eqref{P-P}, giving
\begin{multline*}
   xP_{n+1,n}(x) = P_{n+1,n+1}(x) + d_{n+1,n} P_{n+1,n}(x) + (a_{n+1,n} + b_{n+1,n})P_{n,n}(x) \\
+ b_{n+1,n} (d_{n,n-1}-c_{n,n-1}) P_{n,n-1}(x).
\end{multline*}
Comparing with the recurrence relation \eqref{step-line} with $n$ replaced by $2n+1$ gives the required result.
\end{proof}

\subsection{The nearest neighbor coefficients for general $r$}

For general $r$ there are more partial difference equations for the nearest neighbor recurrence coefficients. The nearest neighbor recurrence relations \eqref{NNrec} can be written as
\[    Y_{\vec{n} + \vec{e}_k} = R_{k}(\vec{n}) Y_{\vec{n}}, \qquad 1 \leq k \leq r, \]
where
\[   Y_{\vec{n}} = \begin{pmatrix}
                   P_{\vec{n}}(x) \\ P_{\vec{n}-\vec{e}_1}(x) \\ \vdots \\ P_{\vec{n}-\vec{e}_r}(x)
                  \end{pmatrix}  \]
and $R_{k}(\vec{n})$ are $(r+1)\times(r+1)$ transfer matrices given by
\[   R_{k}(\vec{n}) = \begin{pmatrix}
                      x - b_{\vec{n},k} & - a_{\vec{n},1} & \cdots & -a_{\vec{n},k} & \cdots & -a_{\vec{n},r} \\
                      1 & b_{\vec{n}-\vec e_1,1}-b_{\vec{n}-\vec e_1,k} & \cdots & 0 & \cdots & 0 \\
                      \vdots &   & \ddots &  & & \vdots \\
                      1 & 0 & \cdots & 0 & \cdots & 0 \\
                      \vdots &   & & & \ddots & \vdots\\
                      1 & 0 & \cdots & 0 & \cdots & b_{\vec{n}-\vec e_r,r}-b_{\vec{n}-\vec e_r,k}
                      \end{pmatrix}.  \]
Expressing that $Y_{\vec{n}+\vec{e}_i+\vec{e}_j}$ can be computed in two ways when $i \neq j$ is done by
\[    R_i(\vec{n}+\vec{e}_j)R_j(\vec{n}) = R_j(\vec{n}+\vec{e}_i)R_i(\vec{n}), \]
and this gives the following partial difference relations \cite[Thm.~3.2]{WVA}: for all $1 \leq i \neq j \leq r$  one has
\begin{eqnarray}
   b_{\vec{n}+\vec{e}_i,j}-b_{\vec{n},j} &=& b_{\vec{n}+\vec{e}_j,i}-b_{\vec{n},i},  \label{PD1} \\
   \sum_{k=1}^r a_{\vec{n}+\vec{e}_j,k} - \sum_{k=1}^r a_{\vec{n}+\vec{e}_i,k} &=& \det \begin{pmatrix}
                               b_{\vec{n}+\vec{e}_j,i} & b_{\vec{n},i} \\
                               b_{\vec{n}+\vec{e}_i,j} & b_{\vec{n},j} \end{pmatrix}, \label{PD2}  \\
   \frac{a_{\vec{n}+\vec{e}_j,i}}{a_{\vec{n},i}} &=&   \frac{b_{\vec{n},j}-b_{\vec{n},i}}{b_{\vec{n}-\vec{e}_i,j}-b_{\vec{n}-\vec{e}_i,i}}.
   \label{PD3}
\end{eqnarray}

\begin{theorem}
Suppose the recurrence coefficients $(a_n^2(\mu_i))_{n\geq 1}$ and $(b_n(\mu_i))_{n \geq 0}$ of the monic orthogonal polynomials for the measure
$\mu_i$ are known $(1 \leq i \leq r)$.
Then the nearest neighbor recurrence coefficients $a_{\vec{n},j}, b_{\vec{n},j}$ $(1 \leq j \leq r)$ can be computed from
\eqref{PD1}--\eqref{PD3} and the boundary conditions
\[   a_{n\vec{e}_j,j} = a_n^2(\mu_j), \quad b_{n\vec{e}_j,j} = b_n(\mu_j), \qquad n \geq 0,\ 1 \leq j \leq r, \]
where $a_0^2(\mu_j)=0$, and
\[  a_{n\vec{e}_i,j} = 0, \qquad n \geq 0,\ i \neq j, \]
provided that $b_{\vec{n},i} \neq b_{\vec{n},j}$ for all multi-indices $\vec{n} \in \mathbb{N}^r$ and $1 \leq i \neq j \leq r$.
\end{theorem}

\begin{proof}
We use induction on the length $N=|\vec{n}|$ of the multi-index $\vec{n}$.
For $|\vec{n}|=0$ we see that $a_{\vec{0},j}=0$ and $b_{\vec{0},j}=b_0(\mu_j)$ for $1 \leq j \leq r$. If $|\vec{n}|=1$ then
$\vec{n}=\vec{e}_i$ for some $i$ with $1 \leq i \leq r$. Therefore $a_{\vec{n},j}= a_{\vec{e}_i,j} = 0$ whenever $i \neq j$ and
$a_{\vec{e}_i,i}=a_1^2(\mu_i)$. Furthermore $b_{\vec{n},i} = b_{\vec{e}_i,i} = b_1(\mu_i)$. Using \eqref{PD1} we also find
\[   b_{\vec{e}_i,j} = b_{\vec{e}_j,i} + b_0(\mu_j) - b_0(\mu_i),  \]
and  \eqref{PD2} gives
\[   b_{\vec{e}_j,i}b_0(\mu_j) - b_{\vec{e}_i,j}b_0(\mu_i) =  \sum_{k=1}^r a_{\vec{e}_j,k} - \sum_{k=1}^r a_{\vec{e}_i,k}. \]
Solving this linear system gives for $i \neq j$
\[    b_{\vec{e}_j,i} = b_0(\mu_i) + \frac{\sum_{k=1}^r a_{\vec{e}_j,k} - \sum_{k=1}^r a_{\vec{e}_i,k}}{b_0(\mu_j)-b_0(\mu_i)}, \]
provided $b_0(\mu_i) \neq b_0(\mu_j)$. Hence all the nearest neighbor recurrence coefficients are known for $|\vec{n}|=1$.

Suppose that we know all the nearest neighbor recurrence coefficients with multi-indices of length $\leq N$. Let $\vec{m}$ be a multi-index
of length $N+1$. In order to compute $a_{\vec{m},i}$ we choose a $j \neq i$ such that $m_i \geq 1$ and write $\vec{m} = \vec{n} + \vec{e}_j$,
where $\vec{n}$ is a multi-index of length $N$.
Then from \eqref{PD3} we find that
\begin{equation}  \label{ar1}
  a_{\vec{m},i} =  a_{\vec{n}+\vec{e}_j,i} = a_{\vec{n},i} \frac{b_{\vec{n},j}-b_{\vec{n},i}}{b_{\vec{n}-\vec{e}_i,j}-b_{\vec{n}-\vec{e}_i,i}}.
\end{equation}
The coefficients on the right hand side have a multi-indices of length $N$ or $N-1$ and hence \eqref{ar1} allows us to compute
$a_{\vec{m},i}$ when $m_i=n_i \geq 1$. If $m_i=0$ then $a_{\vec{m},i}$ is the coefficient of the polynomial $P_{\vec{m}-\vec{e}_i}$ which is $0$, hence
we don't need this coefficient and we can set it equal to $0$. If $m_j=0$ for all $j \neq i$ then $\vec{m}=(N+1)\vec{e}_i$ and we have
$a_{(N+1)\vec{e}_i,i} = a_{N+1}^2(\mu_i)$.

If $i \neq j$ then the equations \eqref{PD1} and \eqref{PD2} are a linear system for the two unknowns $b_{\vec{n}+\vec{e}_i,j}$ and $b_{\vec{n}+\vec{e}_j,i}$.
The solution is
\begin{equation}  \label{br1}
    b_{\vec{n}+\vec{e}_j,i} = b_{\vec{n},i} + \frac{\sum_{k=1}^r a_{\vec{n}+\vec{e}_j,k}- \sum_{k=1}^r a_{\vec{n}+\vec{e}_i,k}}
                              {b_{\vec{n},j}-b_{\vec{n},i}},
\end{equation}
provided that $b_{\vec{n},j} \neq b_{\vec{n},i}$. Hence if there exists $j \neq i$ with $m_j \geq 1$, then $\vec{m}=\vec{n}+\vec{e}_j$, and we can compute $b_{\vec{m},i}$ from \eqref{br1}.
If $m_j=0$ for all $j \neq i$ then $\vec{m}=(N+1)\vec{e}_i$ and $b_{(N+1)\vec{e}_i,i} = b_{N+1}(\mu_i)$.
\end{proof}

config: arxiv


TITLE: Does $ a_{\lfloor n + \sqrt{n/2}\rfloor} \le a_n\left(1 + \frac{b}{n^{\frac{1}{2}+\varepsilon}}\right) $ imply that $a_n$ is bounded?
QUESTION [1 upvotes]: Let $(a_n)_{n\ge 1}$ be an increasing sequence of positive real numbers and suppose there is a constant $b > 0$ and $\varepsilon > 0$ such that
$$
a_{\lfloor n + \sqrt{n/2}\rfloor} \le a_n\left(1 + \frac{b}{n^{\frac{1}{2}+\varepsilon}}\right)
$$
for all $n$ sufficiently large. Does it follow that $a_n$ is bounded? i.e., is $\limsup a_n < \infty$?
If the left hand side of the inequality was $a_{n+1}$, then we could use induction to derive
$$
a_{n + 1} \le a_{n_0}\prod_{k=n_0}^{n}\left(1 + \frac{b}{k^{\frac{1}{2}+\varepsilon}}\right)
$$
and even then, it is not clear whether $a_n$ is bounded. Moreover, the form of the index on the left hand side of the inequality makes it difficult to apply this method. Is there a way to do this?
Any help or comments are welcome.

REPLY [1 votes]: Write $f(n) = \bigl\lfloor n + \sqrt{n/2} \bigr\rfloor$ and consider the sequence $(n_k)_{k\geq 0}$ defined by
$$ n_0 \geq 2 \qquad\text{and}\qquad n_{k+1} = f(n_k). $$
Since $f(n) \geq n+1$ for $n \geq 2$, we find that  $(n_k)_{k\geq 0}$ is strictly increasing. Moreover, for any sufficiently large $n$, we have
$$
\sqrt{f(n)}
\geq \biggl( n + \frac{\sqrt{n}}{2} \biggr)^{1/2}
= \sqrt{n} \biggl( 1 + \frac{1}{2\sqrt{n}} \biggr)^{1/2}
\geq \sqrt{n} \biggl( 1 + \frac{1}{5\sqrt{n}} \biggr)
= \sqrt{n} + \frac{1}{5}. \tag{*}
$$
So by choosing $n_0$ to be sufficiently large, it follows that
$$ \sqrt{n_k} \geq \sqrt{n_0} + \tfrac{1}{5}k \quad \text{for all} \quad k \geq 0 $$
Then
$$ a_{n_k}
\leq a_{n_0} \prod_{j=0}^{k-1} \biggl( 1 + \frac{b}{n_j^{0.5+\epsilon}} \biggr)
\leq a_{n_0} \prod_{j=0}^{k-1} \biggl( 1 + \frac{b}{(\sqrt{n_0} + \frac{1}{5}j)^{1+2\epsilon}} \biggr) $$
and this upper bound converges as $k\to\infty$. Therefore $(a_n)$ is bounded above.

Addendum. In fact we can give a precise range of $n$ for which the inequality $\text{(*)}$ holds.

To begin with, the first step of $\text{(*)}$ holds precise when $ \left\lfloor \sqrt{\frac{n}{2}} \right\rfloor \geq \frac{\sqrt{n}}{2} $, which turns out to be true if and only if $n \geq 18$.

Next, it follows that $\sqrt{1+x} \geq 1 + \frac{2}{5}x$ if and only if $0 \leq x \leq \frac{5}{4}$. So the third step of $\text{(*)}$ holds for any $n$.

Combining altogether, $\text{(*)}$ holds for any $n \geq 18$ and we may take $n_0 = 18$.


TITLE: Chain rule for implicit functions
QUESTION [0 upvotes]: Let $F_1(x_1,x_2,x_3)=f(x_1,f(x_1,x_2,x_3),x_3)$ and $F_2(x_1,x_2,x_3)=f(x_1,x_2, f(x_1,x_2,x_3))$. Find $\displaystyle \frac{\partial F_i}{\partial x_j}$ for all $i=1,2$ and $j=1,2,3$.

I know that here I have to use chain rule. But since they are not explicit expressions it is very confusing me to find them. Can somebody please explain me how to find them?

REPLY [1 votes]: To compute $\frac{\partial F_1}{\partial x_j}$ you have to derive each coordinate of $f$ and multiply by the derivative w.r.t. to $x_j$ of "what is inside" that coordinate, then  you add all. Let us call $y_1,y_2,y_3$
the variables of $f$ (i.e. $f(y_1,y_2,y_3))$, then
$$\frac{\partial F_1}{\partial x_1}(x_1,x_2,x_3)=\frac{\partial f}{\partial y_1}(x_1,f(x_1,x_2,x_3),x_3)\\+\frac{\partial f}{\partial y_2}(x_1,f(x_1,x_2,x_3),x_3)\cdot
\frac{\partial f}{\partial y_1}(x_1,x_2,x_3)
$$
$$\frac{\partial F_1}{\partial x_2}(x_1,x_2,x_3)=\frac{\partial f}{\partial y_2}(x_1,f(x_1,x_2,x_3),x_3)\cdot
\frac{\partial f}{\partial y_2}(x_1,x_2,x_3)
$$
$$\frac{\partial F_1}{\partial x_3}(x_1,x_2,x_3)=\frac{\partial f}{\partial y_2}(x_1,f(x_1,x_2,x_3),x_3)\cdot
\frac{\partial f}{\partial y_3}(x_1,x_2,x_3)+\frac{\partial f}{\partial y_3}(x_1,f(x_1,x_2,x_3),x_3)
$$
$$\frac{\partial F_2}{\partial x_1}(x_1,x_2,x_3)=\frac{\partial f}{\partial y_1}(x_1,x_2,f(x_1,x_2,x_3))+\\
\frac{\partial f}{\partial y_3}(x_1,x_2,f(x_1,x_2,x_3))\cdot \frac{\partial f}{\partial y_1}(x_1,x_2,x_3)
$$
$$\frac{\partial F_2}{\partial x_2}(x_1,x_2,x_3)=\frac{\partial f}{\partial y_2}(x_1,x_2,f(x_1,x_2,x_3))+\\
\frac{\partial f}{\partial y_3}(x_1,x_2,f(x_1,x_2,x_3))\cdot \frac{\partial f}{\partial y_2}(x_1,x_2,x_3)
$$
$$\frac{\partial F_2}{\partial x_3}(x_1,x_2,x_3)=\frac{\partial f}{\partial y_3}(x_1,x_2,f(x_1,x_2,x_3))\cdot \frac{\partial f}{\partial y_3}(x_1,x_2,x_3)
$$


TITLE: Complex number locus on Argand Plane
QUESTION [0 upvotes]: If the point $z$ in the complex plane describes a circle with radius $2$ with centre at the origin then the point $z+\frac{1}{z}$ describe...
Options
A Circle
A parabola
an ellipse
a hyperbola
What I did... Did some elementary manipulations on the given term and conditions but couldn't really equate my manipulations to anything.

REPLY [2 votes]: Putting $z=2e^{it}$ gives the $x$ and $y$ coordinates of the transformed
shape as $x=(5/2)\cos t$ and $y=(3/2)\sin t$. These parametrise the ellipse
$$\frac{x^2}{(5/2)^2}+\frac{y^2}{(3/2)^2}=1.$$


TITLE: Waveguide and ethernet
QUESTION [3 upvotes]: I am having a doubt about waveguides and the traditional Ethernet cable.
The newest Ethernet i.e. the 10Gbps Ethernet has an astounding speed of 10Gbps. Which consists of 4 lanes in Cu as physical media. Each wire pair carries 3.125 G baud as signal!
Now that's enormous even for Cu cable.
Does that somehow violate the waveguide principle, that the impedance of the carrier should be lesser than the impedance of the environment, and as the freq. of wave increases the impedance of carrier increases...
Some light in this might be helpful...
Regards
Vineet

REPLY [0 votes]: 1000BaseT and 10GBaseT both use amplitude modulation so that a cable with a maximum frequency of 500Mhz (cat6) can carry much higher amounts of data.


TITLE: Expected value of $x^t\Sigma x$ for multivariate normal distribution $N(0,\Sigma)$
QUESTION [3 upvotes]: For standard normal distribution, the expected value of $x^2$ is $1$. A natural question is that in the multivariate case, what is the expected value of $x^t\Sigma x$ for multivariate normal distribution $x \sim N(0,\Sigma)$? I have difficulty to carry out the integral, but would guess the result is related to the norm of $\Sigma$.

REPLY [6 votes]: Note that $x^T \Sigma x = tr(x^T \Sigma x) = tr( x x^T \Sigma)$ by invariance of trace to cyclic permutations and trace of a scalar is itself. Also, note that trace is linear (so we can bring expectations inside the trace). Since $x$ is zero mean, $E[x x^T ] = \Sigma$.
Then we have, by the aforementioned, $E[x^T \Sigma x] = E[ tr( x x^T \Sigma) ] = tr(E[x x^T \Sigma]) = tr(E[x x^T ] \Sigma) = tr(\Sigma^2)$.


TITLE: Distance between trunctated random walk and its normal form
QUESTION [4 upvotes]: I have $$X_i \sim N(0,1), \quad S_n=X_1+\cdots+X_n,$$
$$  \mathscr{S}_n (t, \omega) := \frac{1}{ \sqrt{n} } \sum_{i=1}^{n} \left[ S_{i-1} (\omega ) + n \left( t - \frac{i-1}{n} \right) X_{i}(\omega) \right] \textbf{1}_{ \left( \frac{i-1}{n}, \frac{i}{n} \right] } (t) $$ and let $$d(f,g) := \sup_{x \in [0,1]}| f(x) - g(x) | .$$
Then $$ Y_i^{(v)} := X_i \textbf{1}_{\{ |X_i| \leq v \} },
        X_i^{(v)} := \left( Y_i^{(v)} - E [ Y_1^{(v)} ] \right) $$ analogically we define
$$ S_n^{(v)} = \sum_{i=1}^n X_i \textbf{1}_{ \left\{ \left| X_i \right| \leq v \right\} } - E \left[ \sum_{i=1}^n X_i \textbf{1}_{ \left\{ \left| X_i \right| \leq v \right\} } \right] $$ and $\mathscr{S}_n^{(v)}$.
In the lemma 3.2 (lecture notes on Donsker's theorem Davar Khoshnevisan,p.6, https://www.math.utah.edu/~davar/ps-pdf-files/donsker.pdf) it's shown that
$$ \sup_{n \geq 1}  \| d \left( \mathscr{S}_n^{(v) }, \mathscr{S}_n \right) \|_2 \leq 2 \sqrt{E \left[ X_1^2; |X_1| > v \right] } \text{ for all } v>0.$$
But to prove the Donsker's theorem we also need that
$$ \lim_{v \to \infty} \sup_n P \left( d \left( \mathscr{S}_n^{(v) }, \mathscr{S}_n \right) > \lambda \right)=0$$ for all $\lambda>0$ which should follow from the inequality above and  Chebyshev’s inequality. I don't know how to show this.

REPLY [2 votes]: You have
$$
\begin{aligned}
&\sup_{n\ge1}\sqrt{E[d\left( \mathscr{S}_n^{(v) }, \mathscr{S}_n \right)^2]} \\
&=\sup_{n\ge1}\|d\left( \mathscr{S}_n^{(v) }, \mathscr{S}_n \right) \|_2 \leq 2 \sqrt{E \left[ X_1^2; |X_1| > v \right] } \text{ for all } v>0.
\end{aligned}
$$
So, by the Chebyshev/Markov inequality,
$$
\begin{aligned}
&\lim_{v \to \infty} \sup_n P \left( d \left( \mathscr{S}_n^{(v) }, \mathscr{S}_n \right) > \lambda \right) \\
&\le\lim_{v \to \infty} \sup_n
\frac{E[d\left( \mathscr{S}_n^{(v) }, \mathscr{S}_n \right)^2]}{\lambda^2} \\
&\le\lim_{v \to \infty}
\frac{4E \left[ X_1^2; |X_1| > v \right] }
{\lambda^2}
=0,
\end{aligned}
$$
because, by the dominated convergence theorem, $E \left[ X_1^2; |X_1| > v \right]\to0$ as $v\to\infty$.


TITLE: How to compare combination with power?
QUESTION [0 upvotes]: How can we compare combination with power when the both are very very large numbers?
For example how to compare $\binom{1000}{500}$ with $2^{500} , 3^{500} , 4^{500}$ ?

REPLY [2 votes]: First, note that if $\log(a) > \log(b)$, then $a > b$.
Stirlings approximation lets us write
$$\log(n!) = n\log(n)-n$$
So that since
$${1000 \choose 500} = \frac{1000!}{500! 500!}$$
we have
$$\log\left(\frac{1000!}{500! 500!}\right) = \log(1000!) - \log((500!)^2)\approx (1000\log(1000) - 1000) - 2(500\log(500) - 500) = 1000(\log(1000)-\log(500)) = 1000 \log(2) $$
we can compare this to, for example, $2^{500}$ by comparing their logarithms:
$$\log\left(\frac{1000!}{500! 500!}\right) \approx 1000 \log(2) > \log(2^{500}) = 500 \log(2)$$
$$1000 \log(2) > 500 \log(2)$$
so that ${1000 \choose 500} > 2^{500}$.