\section{Introduction}
\label{sec:intro}

Exactly solvable quantum mechanical systems of one degree of freedom 
(Schr\"odinger equation: $\mathcal{H}\phi_n(x)=\mathcal{E}_n\phi_n(x)$) can be
deformed keeping solvability by the Darboux transformations.
By such deformations, new types of orthogonal polynomials, exceptional or
multi-indexed polynomials, are obtained \cite{gkm08}--\cite{idQMcH}
for the Askey-scheme of hypergeometric orthogonal polynomials \cite{ismail,kls}.
They satisfy second order differential or difference equations and form
a complete set of orthogonal basis in an appropriate Hilbert space in spite of
missing degrees, by which the restrictions of Bochner's theorem \cite{ismail}
are avoided.
We have studied orthogonal polynomials based on quantum mechanical formulations:
ordinary quantum mechanics (oQM) and two kinds of discrete quantum mechanics
(dQM), dQM with pure imaginary shifts (idQM) and dQM with real shifts (rdQM)
\cite{os24}. The Schr\"odinger equation for oQM is a differential equation and
that for dQM is a difference equation. The coordinate $x$ for oQM and idQM is
continuous and that for rdQM is discrete.

Depending on the choice of seed solution, the Darboux transformation can be
divided into three types: isospectral, state-deleting and state-adding.
When the wavefunction of the virtual state $\tilde{\phi}_{\text{v}}(x)$,
eigenstate $\phi_n(x)$ and pseudo virtual state
$\tilde{\phi}^{\text{pv}}_{\text{v}}(x)$ is used as a seed solution,
the Darboux transformation is isospectral, state-deleting and state-adding,
respectively.
The wavefunctions $\tilde{\phi}_{\text{v}}(x)$ and
$\tilde{\phi}^{\text{pv}}_{\text{v}}(x)$ are obtained from $\phi_n(x)$ by
twisting parameters.
For oQM and idQM, $\tilde{\phi}_{\text{v}}(x)$ and
$\tilde{\phi}^{\text{pv}}_{\text{v}}(x)$ are solutions of the Schr\"odinger
equation and not square integrable \cite{os25,os27,idQMcH,os29,os30}.
On the other hand, for rdQM, they satisfy the Schr\"odinger equation except
for the boundary \cite{os26,os35,casoidrdqm}.
The Darboux transformations with $\tilde{\phi}_{\text{v}}(x)$ give the case-(1)
multi-indexed polynomials and those with $\phi_n(x)$ or
$\tilde{\phi}^{\text{pv}}_{\text{v}}(x)$ give the case-(2) multi-indexed
polynomials.
Here, the case-(1) is the case that the set of missing degrees of the
multi-indexed polynomials is $\{0,1,\ldots,\ell-1\}$, and the case-(2) is
otherwise.
Another type of seed solution is discussed for the oQM systems with a finite
number of eigenstates $\phi_n(x)$ ($n=0,1,\ldots,n_{\text{max}}$) \cite{os28}.
We call them overshoot eigenfunctions
$\tilde{\phi}^{\text{os}}_n(x)$,
which have the same form as the eigenstates $\phi_n(x)$ but with
$n>n_{\text{max}}$ and are not square integrable.
The overshoot eigenfunctions correspond to the virtual or pseudo virtual
states.

The Hamiltonian of a finite type rdQM system is a real symmetric matrix of order
$N+1$. The coordinate $x$ takes a value in $\{0,1,\ldots,N\}$ and the number
of the eigenstates (eigenvectors) $\phi_n(x)$ ($n=0,1,\ldots,N$) is $N+1$.
By using the pseudo virtual state $\tilde{\phi}^{\text{pv}}_{\text{v}}(x)$
obtained from $\phi_n(x)$ by twisting parameters, the state-adding Darboux
transformation and the Casoratian identities are studied in \cite{casoidrdqm}.
Recently another type of seed solution that gives the state-adding Darboux
transformation was found by Miki, Tsujimoto and Vinet \cite{mtv22}.
They studied single-indexed exceptional Krawtchouk polynomials and one of them
(their type (\romannumeral2)) corresponds to the state-adding Darboux
transformation.
Motivated by their work, Sasaki and the present author studied the $M$-step
state-adding Darboux transformations, whose seed solutions have the same form
as the eigenstate $\phi_n(x)$ but with $n>N$ \cite{os40}. This situation
corresponds to the overshoot eigenfunctions $\tilde{\phi}^{\text{os}}_n(x)$ in
the oQM systems with a finite number of eigenstates.
The overshoot eigenfunction has an infinite norm, but the seed solution here
has a ``zero norm.''
(The meaning of ``zero norm'' is the following.
The Schr\"odinger equation (with $x=0,1,\ldots,N$) is a matrix eigenvalue
problem, and it can be interpreted as a difference equation with a continuous
$x$. The function $\phi_n(x)$ with $n>N$ satisfies this difference equation
and vanishes at $x=0,1,\ldots,N$.
So, $\phi_n(x)$ with $n>N$ satisfies the Schr\"odinger equation and it is a
zero vector, namely zero norm.)
After the $M$-step state-adding Darboux transformations with the seed solutions
$\phi_n(x)$ ($n\in\mathcal{D}=\{d_1,d_2,\ldots,d_M\}$, $d_j>N$), the order of
the Hamiltonian becomes $N+M+1$, and the coordinate $x$ takes a value in
$\{-M,-M+1,\ldots,N\}$. The deformed Hamiltonian has $N+M+1$ eigenvectors.
It is easy to find $N+1$ eigenvectors which correspond to the original
eigenvectors. They are expressed as
$\text{W}_{\text{C}}[\phi_{d_1},\ldots,\phi_{d_M},\phi_n](x)\times(\cdots)$
with $n\in\{0,1,\ldots,N\}$.
However, it is difficult to find extra $M$ eigenvectors.
In \cite{os40}, the special case $\mathcal{D}=\{N+1,N+2,\ldots,N+M\}$ is
studied in detail.

In this paper we study the deformations of finite type rdQM systems by the
$M$-step state-adding Darboux transformations with the seed solutions
$\phi_n(x)$ ($n>N$) and obtain all eigenvectors.
The original systems are described by twelve orthogonal polynomials: ($q$-)Racah,
(dual, $q$-)Hahn, Krawtchouk and five types of $q$-Krawtchouk.
In \cite{mtv22}, which corresponds to $M=1$ case, one extra eigenvector is
derived by two methods: (1) solving the difference equation, (2) shifting $N$
to $N+\varepsilon$ and taking $\varepsilon\to 0$ limit.
For general $M$ case, the first method is difficult and we adopt the second
method. Extra $M$ eigenvectors are obtained from
$\text{W}_{\text{C}}[\phi_{d_1},\ldots,\phi_{d_M},\phi_n](x)\times(\cdots)$
in the ``$n\to d_i$ limit'', which is achieved by shifting $N$ to
$N+\varepsilon$ and taking $\varepsilon\to 0$ limit.
The eigenvectors are described by new multi-indexed orthogonal polynomials
$\check{Q}_{\mathcal{D}',n}(x)$.

This paper is organized as follows.
In section \ref{sec:DT} the finite type rdQM systems are recapitulated and
the multi-step Darboux transformations with seed solutions $\phi_n(x)$ are
discussed.
In section \ref{sec:KAmiop} the results obtained in \S\,\ref{sec:DT} are
applied to the case of seed solutions $\phi_n(x)$ with $n\leq N$, which
corresponds to the state-deleting Darboux transformations.
The eigenvectors are described by the Krein-Adler type multi-indexed
orthogonal polynomials $\check{P}^{\text{KA}}_{\mathcal{D},n}(x)$
($n\in\{0,1,\ldots,N\}\backslash\mathcal{D}$).
This case was studied in \cite{os22}, but the explicit forms of their
difference equations, orthogonal relations etc. are new results.
Section \ref{sec:newmiop} is the main part of the paper.
The results obtained in \S\,\ref{sec:DT} are
applied to the case of seed solutions $\phi_n(x)$ with $n>N$, which
corresponds to the state-adding Darboux transformations.
The eigenvectors are described by new multi-indexed orthogonal polynomials
$\check{Q}_{\mathcal{D}',n}(x)$ ($n\in\{0,1,\ldots,N\}\cup\mathcal{D}$).
Section \ref{sec:summary} is for a summary and comments.
Data for the twelve orthogonal polynomials are presented in
Appendix\,\ref{app:data}. Data for various other quantities and several
formulas are also presented.
In Appendix\,\ref{app:n->di} we discuss the ``$n\to d_i$ limit.''

\section{Darboux Transformations}
\label{sec:DT}

In this section, after recapitulating the finite type rdQM systems
\cite{os12,os24}, we discuss the multi-step Darboux transformations with seed
solutions $\phi_n(x)$ \cite{os22}, especially their algebraic aspects.

\subsection{Original systems}
\label{sec:orgsys}

We consider the finite type rdQM systems as a starting point,
whose eigenvectors are described by the following twelve orthogonal polynomials
\cite{os12,os24}: Hahn (H), Krawtchouk (K), Racah (R), dual Hahn (dH),
dual quantum $q$-Krawtchouk (dq$q$K), $q$-Hahn ($q$H), $q$-Krawtchouk ($q$K),
quantum $q$-Krawtchouk (q$q$K), affine $q$-Krawtchouk (a$q$K), $q$-Racah ($q$R),
dual $q$-Hahn (d$q$H) and dual $q$-Krawtchouk (d$q$K).
The data of these polynomials are presented in Appendix\,\ref{app:poly}.

Let $N$ be a positive integer.
The Hamiltonian $\mathcal{H}$ of a finite type rdQM is a real symmetric
(tridiagonal in this case) matrix of order $N+1$,
\begin{align}
  \mathcal{H}=&(\mathcal{H}_{x,y})_{x,y=0,1,\ldots,N},\quad
  \mathcal{H}_{x,x}=B(x)+D(x),\quad\mathcal{H}_{x,y}=0\ \ (|x-y|>1),\nonumber\\
  &\mathcal{H}_{x,x+1}=-\sqrt{B(x)D(x+1)},\quad
  \mathcal{H}_{x,x-1}=-\sqrt{B(x-1)D(x)}.
  \label{orgH}
\end{align}
Here the potential functions $B(x)$ and $D(x)$ are positive but vanish at the
boundary,
\begin{align}
  &B(x)>0\ \ (x=0,1,\ldots,N-1),\quad D(x)>0\ \ (x=1,2,\ldots,N),\nonumber\\
  &B(N)=0,\quad D(0)=0.
  \label{B,D>0}
\end{align}
We write a matrix \eqref{orgH} as
\begin{align}
  \mathcal{H}&=-\sqrt{B(x)D(x+1)}\,e^{\partial}
  -\sqrt{B(x-1)D(x)}\,e^{-\partial}+B(x)+D(x)\nonumber\\
  &=-\sqrt{B(x)}\,e^{\partial}\sqrt{D(x)}
  -\sqrt{D(x)}\,e^{-\partial}\sqrt{B(x)}+B(x)+D(x),
  \label{Hdef}
\end{align}
where $e^{\pm\partial}$ is a matrix whose $(x,y)$-element is $\delta_{x\pm1,y}$
($\Rightarrow$ $(e^{\partial})^{\dagger}=e^{-\partial}$),
and $A(x)$ means a diagonal matrix $A(x)=\text{diag}(A(0),A(1),\ldots,A(N))$.
Note that $e^{\pm\partial}e^{\mp\partial}\neq1$ due to the effect of boundaries,
\begin{equation}
  e^{\partial}e^{-\partial}=\text{diag}(1,1,\ldots,1,0),\quad
  e^{-\partial}e^{\partial}=\text{diag}(0,1,1,\ldots,1).
\end{equation}
We have $B(x)e^{\partial}e^{-\partial}=B(x)$,
$D(x)e^{-\partial}e^{\partial}=D(x)$, etc.
The Schr\"odinger equation of rdQM is a matrix eigenvalue problem,
\begin{equation}
  \mathcal{H}\phi_n(x)=\mathcal{E}_n\phi_n(x)\ \ (n=0,1,\ldots,N),
  \label{Hphin}
\end{equation}
where $\mathcal{H}\phi_n(x)\stackrel{\text{def}}{=}\sum_{y=0}^N\mathcal{H}_{x,y}\phi_n(y)$.
The Hamiltonian \eqref{Hdef} can be written in a factorized form,
\begin{equation}
  \mathcal{H}=\mathcal{A}^{\dagger}\mathcal{A},\quad
  \mathcal{A}\stackrel{\text{def}}{=}\sqrt{B(x)}-e^{\partial}\sqrt{D(x)},\quad
  \mathcal{A}^{\dagger}\stackrel{\text{def}}{=}\sqrt{B(x)}-\sqrt{D(x)}\,e^{-\partial}.
  \label{H=AdA}
\end{equation}
The tridiagonality \eqref{orgH} and factorization \eqref{H=AdA} imply
$0=\mathcal{E}_0<\mathcal{E}_1<\cdots<\mathcal{E}_N$.
The ground state eigenvector $\phi_0(x)$ is characterized by
$\mathcal{A}\phi_0(x)=0$,
\begin{equation}
  \sqrt{B(x)}\,\phi_0(x)=\sqrt{D(x+1)}\,\phi_0(x+1),
  \label{phi0eq}
\end{equation}
and given by
\begin{equation}
  \phi_0(x)\stackrel{\text{def}}{=}\sqrt{\prod_{y=0}^{x-1}\frac{B(y)}{D(y+1)}},
  \label{phi0def}
\end{equation}
which satisfies the normalization $\phi_0(x)=1$ by the convention
$\prod_{j=n}^{n-1}*\stackrel{\text{def}}{=} 1$.
For the twelve systems under consideration, the eigenvectors have the following
form,
\begin{equation}
  \phi_n(x)=\phi_0(x)\check{P}_n(x),\quad
  \check{P}_n(x)\stackrel{\text{def}}{=} P_n\bigl(\eta(x)\bigr),
  \label{phin=phi0cPn}
\end{equation}
where $P_n(\eta(x))$ is a polynomial of degree $n$ in the sinusoidal
coordinate $\eta(x)$ \cite{os7,os12}. We take the normalization as
\begin{equation}
  \check{P}_n(0)=P_n(0)=1.
  \label{cPn0=1}
\end{equation}
The similarity transformed Hamiltonian $\widetilde{\mathcal{H}}$ is defined by
\begin{align}
  &\widetilde{\mathcal{H}}\stackrel{\text{def}}{=}\phi_0(x)^{-1}\circ\mathcal{H}\circ\phi_0(x)
  =B(x)(1-e^{\partial})+D(x)(1-e^{-\partial})
  \label{Htdef}\\
  &\Bigl(\Rightarrow\widetilde{\mathcal{H}}_{x,x+1}=-B(x),
  \ \widetilde{\mathcal{H}}_{x,x-1}=-D(x),
  \ \widetilde{\mathcal{H}}_{x,x}=B(x)+D(x),
  \ \widetilde{\mathcal{H}}_{x,y}=0\ (|x-y|>1)\Bigr),
  \nonumber
\end{align}
and its eigenvalue problem is solved by the polynomial $\check{P}_n(x)$,
\begin{equation}
  \widetilde{\mathcal{H}}\check{P}_n(x)=\mathcal{E}_n\check{P}_n(x)
  \ \ (n=0,1,\ldots,N).
  \label{HtcPn}
\end{equation}
The orthogonality relations of $\check{P}_n(x)$ are
\begin{equation}
  \sum_{x=0}^N\phi_0(x)^2\check{P}_n(x)\check{P}_m(x)
  =\frac{\delta_{nm}}{d_n^2}\ \ (n,m=0,1,\ldots,N).
  \label{orthocPn}
\end{equation}

We have five families of the sinusoidal coordinates $\eta(x)$ \cite{os12},
\begin{equation}
  \begin{array}{rll}
  \text{(\romannumeral1)}:&\eta(x)=x
  &:\text{H,\,K}\\[3pt]
  \text{(\romannumeral2)}:&\eta(x)=x(x+d)
  &:\text{R},\,\text{dH}(d=a+b-1)\\[3pt]
  \text{(\romannumeral3)}:&\eta(x)=1-q^x
  &:\text{dq$q$K}\\[3pt]
  \text{(\romannumeral4)}:&\eta(x)=q^{-x}-1
  &:\text{$q$H,\,$q$K,\,q$q$K,\,a$q$K}\\[3pt]
  \text{(\romannumeral5)}:&\eta(x)=(q^{-x}-1)(1-dq^x)
  &:\text{$q$R},\,\text{d$q$H}(d=abq^{-1}),\,\text{d$q$K}(d=-p).
  \end{array}
  \label{etadef}
\end{equation}
We also have five families of the energy eigenvalues $\mathcal{E}_n$ \cite{os12},
\begin{equation}
  \begin{array}{rll}
  \text{(\romannumeral1)}':&\mathcal{E}_n=n
  &:\text{dH,\,K}\\[3pt]
  \text{(\romannumeral2)}':&\mathcal{E}_n=n(n+\tilde{d})
  &:\text{R},\,\text{H}(\tilde{d}=a+b-1)\\[3pt]
  \text{(\romannumeral3)}':&\mathcal{E}_n=1-q^n
  &:\text{q$q$K}\\[3pt]
  \text{(\romannumeral4)}':&\mathcal{E}_n=q^{-n}-1
  &:\text{d$q$H,\,d$q$K,\,dq$q$K,\,a$q$K}\\[3pt]
  \text{(\romannumeral5)}':&\mathcal{E}_n=(q^{-n}-1)(1-\tilde{d}q^n)
  &:\text{$q$R},\,\text{$q$H}(\tilde{d}=abq^{-1}),\,\text{$q$K}(\tilde{d}=-p).
  \end{array}
  \label{Endef}
\end{equation}
Note that $\eta(0)=\mathcal{E}_0=0$.
The constants $\rho$ and $\kappa$ are defined as follows,
\begin{equation}
  \rho\stackrel{\text{def}}{=}\left\{
  \begin{array}{ll}
  1&:\text{(\romannumeral1)},\,\text{(\romannumeral2)}\\
  q&:\text{(\romannumeral3)}\\
  q^{-1}&:\text{(\romannumeral4)},\,\text{(\romannumeral5)}
  \end{array}\right.,\quad
  \kappa\stackrel{\text{def}}{=}\left\{
  \begin{array}{ll}
  1&:\text{(\romannumeral1)}',\,\text{(\romannumeral2)}'\\
  q&:\text{(\romannumeral3)}'\\
  q^{-1}&:\text{(\romannumeral4)}',\,\text{(\romannumeral5)}'
  \end{array}\right..
  \label{defrhokappa}
\end{equation}
The rdQM systems have a set of parameters
$\boldsymbol{\lambda}=(\lambda_1,\lambda_2,\ldots)$ including the parameter $N$,
and various quantities depend on
$\boldsymbol{\lambda}$. Their dependence is expressed like, $f=f(\boldsymbol{\lambda})$,
$f(x)=f(x;\boldsymbol{\lambda})$.
The parameter $q$ is $0<q<1$ and $q^{\boldsymbol{\lambda}}$ stands for
$q^{(\lambda_1,\lambda_2,\ldots)}=(q^{\lambda_1},q^{\lambda_2},\ldots)$.
We omit writing $q$-dependence and sometimes omit writing
$\boldsymbol{\lambda}$-dependence, when it does not cause confusion.

\subsection{Darboux transformations}
\label{sec:Darb}

\subsubsection{difference equations}
\label{sec:sabuneq}

The matrix eigenvalue problem \eqref{HtcPn} is written in components as,
\begin{equation}
  B(x)\bigl(\check{P}_n(x)-\check{P}_n(x+1)\bigr)
  +D(x)\bigl(\check{P}_n(x)-\check{P}_n(x-1)\bigr)
  =\mathcal{E}_n\check{P}_n(x).
  \label{sabuneqcPn}
\end{equation}
Eq.\eqref{HtcPn} means that \eqref{sabuneqcPn} holds for $n=0,1,\ldots,N$ and
$x=0,1,\ldots,N$. However, we remark that this difference equation
\eqref{sabuneqcPn} holds for $x\in\mathbb{R}$. Moreover \eqref{sabuneqcPn}
holds for $n\in\mathbb{Z}_{\geq 0}$ (exactly speaking, we need replace
$\check{P}_n(x)$ with $\check{P}^{\text{monic}}_n(x)$, see
\S\,\ref{sec:newmiop}).
We also remark that the positive integer parameter $N$ can be extended to a real
value in \eqref{sabuneqcPn}.
Thus the difference equation \eqref{sabuneqcPn}
(with the replacement $\check{P}_n(x)\to\check{P}^{\text{monic}}_n(x)$) holds
for $x\in\mathbb{R}$, $n\in\mathbb{Z}_{\geq0}$ and $N\in\mathbb{R}$.
This is an important point for the Darboux transformations in this subsection
and constructions of new multi-indexed orthogonal polynomials in
\S\,\ref{sec:newmiop}.
The component form of \eqref{Hphin} is
\begin{equation}
  -\sqrt{B(x)D(x+1)}\,\phi_n(x+1)-\sqrt{B(x-1)D(x)}\,\phi_n(x-1)
  +\bigl(B(x)+D(x)\bigr)\phi_n(x)
  =\mathcal{E}_n\phi_n(x).
  \label{sabuneqphin}
\end{equation}
Since the explicit forms of $\phi_0(x)^2$ in Appendix\,\ref{app:poly} are
expressed in terms of $\alpha^x$,
$(\alpha)_x=(\alpha)_{\infty}/(\alpha+x)_{\infty}$ and
$(\alpha\,;q)_x=(\alpha\,;q)_{\infty}/(\alpha q^x\,;q)_{\infty}$,
we can consider $\phi_0(x)^2$ for $x\in\mathbb{R}$ (if needed, we shift $N$ to
a non-integer value).
Let us illustrate this situation using $q$R case as an example.
We rewrite $\phi_0(x)^2$ in Appendix\,\ref{app:qR} as follows:
\begin{align*}
  \phi_0(x)^2
  &=\frac{(a,b,c,d\,;q)_x}{(a^{-1}dq,b^{-1}dq,c^{-1}dq,q\,;q)_x\tilde{d}^x}
  \frac{1-dq^{2x}}{1-d}\nonumber\\
  &=\frac{(a,b,c,d\,;q)_{\infty}}{(aq^x,bq^x,cq^x,dq^x\,;q)_{\infty}}
  \frac{(a^{-1}dq^{x+1},b^{-1}dq^{x+1},c^{-1}dq^{x+1},q^{x+1}\,;q)_{\infty}}
  {(a^{-1}dq,b^{-1}dq,c^{-1}dq,q\,;q)_{\infty}\tilde{d}^x}
  \frac{1-dq^{2x}}{1-d}.
\end{align*}
This is defined for generic values of $x\in\mathbb{R}$ and satisfies
$B(x)\phi_0(x)^2=D(x+1)\phi_0(x)^2$ ($x\in\mathbb{R}$).
However, the factor $(a\,;q)_{\infty}$ vanishes because of $a=q^{-N}$.
To avoid this, we shift $N$ slightly from an integer value.
For $\phi_0(x)^2$ at $x\in\mathbb{Z}$, after simplifying
$(a\,;q)_{\infty}/(aq^x\,;q)_{\infty}$,
we shift $N$ back to an integer value.
For $N\in\mathbb{Z}_{>0}$ and $x\in\mathbb{Z}$, $\phi_0(x;\boldsymbol{\lambda})^2$ is
non-vanishing only for $x=0,1,\ldots,N$.
If we ignore the positivity of the square root argument, the difference
equation \eqref{sabuneqphin} (with the replacement
$\phi_n(x)\to\phi^{\text{monic}}_n(x)=\phi_0(x)\check{P}^{\text{monic}}_n(x)$)
also holds for
$x\in\mathbb{R}$, $n\in\mathbb{Z}_{\geq0}$ and $N\in\mathbb{R}$.
Since these difference equations \eqref{sabuneqcPn} and \eqref{sabuneqphin} are
algebraic relations, they hold for any parameter range of $\boldsymbol{\lambda}$
(unless they are ill-defined).
The orthogonality relations \eqref{orthocPn} with positive weight restrict the
parameter range of $\boldsymbol{\lambda}$ (see the parameter range for the
positivity in Appendix\,\ref{app:poly}).
If we do not require the positivity of the weight factor, the relations
\eqref{orthocPn} themselves are valid for any parameter range of $\boldsymbol{\lambda}$
(unless they are ill-defined), because the relations \eqref{orthocPn} are
finite sums and algebraic relations.

To write the difference equations \eqref{sabuneqcPn} and \eqref{sabuneqphin}
for $x\in\mathbb{R}$ compactly, let us introduce the shift operators
$e^{\pm\hat{\partial}}$ acting on functions of $x\in\mathbb{R}$ as
$e^{\pm\hat{\partial}}f(x)\stackrel{\text{def}}{=} f(x\pm1)$ ($x\in\mathbb{R}$).
They are related as $(e^{\hat{\partial}})^{\dagger}=e^{-\hat{\partial}}$ and
satisfy $e^{\pm\hat{\partial}}e^{\mp\hat{\partial}}=1$.
By replacing $e^{\pm\partial}$ with $e^{\pm\hat{\partial}}$ in \eqref{Hdef} and
\eqref{Htdef}, we define the following operators acting on functions of
$x\in\mathbb{R}$,
\begin{align}
  \mathcal{H}^{\text{op}}&\stackrel{\text{def}}{=}-\sqrt{B(x)}\,e^{\hat{\partial}}\sqrt{D(x)}
  -\sqrt{D(x)}\,e^{-\hat{\partial}}\sqrt{B(x)}+B(x)+D(x),
  \label{Hdef2}\\
  \widetilde{\mathcal{H}}^{\text{op}}
  &\stackrel{\text{def}}{=}\phi_0(x)^{-1}\circ\mathcal{H}^{\text{op}}\circ\phi_0(x)
  =B(x)(1-e^{\hat{\partial}})+D(x)(1-e^{-\hat{\partial}}).
  \label{Htdef2}
\end{align}
Then the difference equations \eqref{sabuneqphin} and \eqref{sabuneqcPn} for
$x\in\mathbb{R}$ are expressed as
\begin{equation}
  \mathcal{H}^{\text{op}}\phi_n(x)=\mathcal{E}_n\phi_n(x),\quad
  \widetilde{\mathcal{H}}^{\text{op}}\check{P}_n(x)=\mathcal{E}_n\check{P}_n(x)
  \quad(n\in\mathbb{Z}_{\geq0}),
  \label{Hopphin}
\end{equation}
where we ignore the positivity of the square root argument, and $\phi_n(x)$ and
$\check{P}_n(x)$ should be replaced with $\phi^{\text{monic}}_n(x)$ and
$\check{P}^{\text{monic}}_n(x)$ for $n>N$, respectively.

\subsubsection{Darboux transformations}
\label{sec:Dar}

To describe the multi-step Darboux transformations, the Casorati determinant
(Casoratian) is needed. The Casoratian for $n$ functions $f_j(x)$ is defined as
\begin{equation}
  \text{W}_{\text{C}}[f_1,f_2,\ldots,f_n](x)
  \stackrel{\text{def}}{=}\det\Bigl(f_k(x+j-1)\Bigr)_{1\leq j,k\leq n},
\end{equation}
(for $n=0$, we set $\text{W}_{\text{C}}[\cdot](x)=1$) and the following
properties are used,
\begin{align}
  &\text{W}_{\text{C}}[gf_1,gf_2,\ldots,gf_n](x)
  =\prod_{k=0}^{n-1}g(x+k)\cdot\text{W}_{\text{C}}[f_1,f_2,\ldots,f_n](x),
  \label{WCformula1}\\
  &\text{W}_{\text{C}}\bigl[\text{W}_{\text{C}}[f_1,f_2,\ldots,f_n,g],
  \text{W}_{\text{C}}[f_1,f_2,\ldots,f_n,h]\,\bigr](x)\nonumber\\
  &=\text{W}_{\text{C}}[f_1,f_2,\ldots,f_n](x+1)\,
  \text{W}_{\text{C}}[f_1,f_2,\ldots,f_n,g,h](x)
  \quad(n\geq 0).
  \label{WCformula2}
\end{align}
See \cite{wcid} for further properties of the Casoratian.

The Darboux transformations for rdQM with seed solutions $\phi_n(x)$ are
studied in \cite{os22}.
Since we are interested in their algebraic aspect here, we consider the
Darboux transformations for $\mathcal{H}^{\text{op}}$ rather than $\mathcal{H}$.
The algebraic calculations are exactly the same in both cases.
In this subsection, we ignore the positivity of the square root argument and
adopt the rule $\sqrt{A^2}=A$ (instead of $\sqrt{A^2}=|A|$ for $A\in\mathbb{R}$).
Although $\phi_n(x)$ should be replaced with
$\phi^{\text{monic}}_n(x)$ for $n>N$ case, we write it as
$\phi_n(x)$ in this subsection for simplicity of presentation.

We consider $M$-step Darboux transformations with the seed solutions
$\phi_{d_1}(x),\phi_{d_2}(x),\ldots,$ $\phi_{d_M}(x)$
satisfying \eqref{Hopphin}.
Let us denote a set of labels of seed solutions as $\mathcal{D}$
(exactly speaking an ordered set)
\begin{equation}
  \mathcal{D}=\{d_1,d_2,\ldots,d_M\}\ \ (d_j\in\mathbb{Z}_{\geq 0}
  \text{ : mutually distinct}).
  \label{Ddef}
\end{equation}
(Although this notation $d_j$ conflicts with the notation of the normalization
constant $d_n$ in \eqref{orthocPn}, we think this does not cause any confusion
because the latter appears as $\delta_{nm}/d_n^2$.)
For later use, let us define $\ell_{\mathcal{D}}$,
$\ell^{\text{KA}}_{\mathcal{D}}$ and $\mathcal{D}^{[i]}$ for
$i\geq-\min\mathcal{D}$,
\begin{align}
  &\ell_{\mathcal{D}}\stackrel{\text{def}}{=}\sum_{j=1}^Md_j-\frac12M(M-1),\quad
  \ell^{\text{KA}}_{\mathcal{D}}\stackrel{\text{def}}{=}\sum_{j=1}^Md_j-\frac12M(M+1),
  \label{ellDdef}\\
  &\mathcal{D}^{[i]}\stackrel{\text{def}}{=}\{d_1+i,d_2+i\ldots,d_M+i\}.
  \label{D[i]def}
\end{align}
The Hamiltonian \eqref{Hdef2} is factorized as
\begin{equation}
  \mathcal{H}^{\text{op}}
  =\mathcal{A}^{\text{op}\,\dagger}\mathcal{A}^{\text{op}},\quad
  \mathcal{A}^{\text{op}}\stackrel{\text{def}}{=}\sqrt{B(x)}-e^{\hat{\partial}}\sqrt{D(x)},\quad
  \mathcal{A}^{\text{op}\,\dagger}\stackrel{\text{def}}{=}
  \sqrt{B(x)}-\sqrt{D(x)}\,e^{-\hat{\partial}},
\end{equation}
and the difference equation and its solution are
\begin{equation}
  \mathcal{H}^{\text{op}}\phi_n(x)
  =\mathcal{E}_n\phi_n(x)\ \ (n\in\mathbb{Z}_{\geq0}).
  \label{Hopphin2}
\end{equation}
Let $s=0$ be the quantity at this starting point, and define the quantity at
the $s$-th step as follows \cite{os22}:
\begin{align}
  &\mathcal{H}^{\text{op}}_{d_1\ldots d_s}
  \stackrel{\text{def}}{=}\hat{\mathcal{A}}^{\text{op}}_{d_1\ldots d_s}
  \hat{\mathcal{A}}_{d_1\ldots d_s}^{\text{op}\,\dagger}+\mathcal{E}_{d_s},\\
  &\hat{\mathcal{A}}^{\text{op}}_{d_1\ldots d_s}
  \stackrel{\text{def}}{=}\sqrt{\hat{B}_{d_1\ldots d_s}(x)}
  -e^{\hat{\partial}}\sqrt{\hat{D}_{d_1\ldots d_s}(x)},
  \ \ \hat{\mathcal{A}}_{d_1\ldots d_s}^{\text{op}\,\dagger}
  \stackrel{\text{def}}{=}\sqrt{\hat{B}_{d_1\ldots d_s}(x)}
  -\sqrt{\hat{D}_{d_1\ldots d_s}(x)}\,e^{-\hat{\partial}},\\
  &\hat{B}_{d_1\ldots d_s}(x)\stackrel{\text{def}}{=}\left\{
  \begin{array}{ll}
  {\displaystyle\sqrt{\hat{B}_{d_1\ldots d_{s-1}}(x+1)
  \hat{D}_{d_1\ldots d_{s-1}}(x+1)}\,
  \frac{\phi_{d_1\ldots d_s}(x+1)}{\phi_{d_1\ldots d_s}(x)}}&(s\geq1)\\
  B(x-1)&(s=0)
  \end{array}\right.,\\
  &\hat{D}_{d_1\ldots d_s}(x)\stackrel{\text{def}}{=}\left\{
  \begin{array}{ll}
  {\displaystyle\sqrt{\hat{B}_{d_1\ldots d_{s-1}}(x)
  \hat{D}_{d_1\ldots d_{s-1}}(x)}\,
  \frac{\phi_{d_1\ldots d_s}(x-1)}{\phi_{d_1\ldots d_s}(x)}}&(s\geq1)\\
  D(x)&(s=0)
  \end{array}\right.,\\
  &\phi_{d_1\ldots d_s\,n}(x)\stackrel{\text{def}}{=}\hat{\mathcal{A}}^{\text{op}}_{d_1\ldots d_s}
  \phi_{d_1\ldots d_{s-1}\,n}(x).
\end{align}
Note that $\phi_{d_1\ldots d_s\,n}(x)$ is defined for $n\in\mathbb{Z}_{\geq 0}$,
but we have $\phi_{d_1\ldots d_s\,n}(x)=0$ for $n\in\{d_1,\ldots,d_s\}$.
Then we obtain
\begin{align}
  &\mathcal{H}^{\text{op}}_{d_1\ldots d_s}\phi_{d_1\ldots d_s\,n}(x)
  =\mathcal{E}_n\,\phi_{d_1\ldots d_s\,n}(x),\\
  &\phi_{d_1\ldots d_s\,n}(x)
  =(-1)^s\sqrt{\prod_{k=1}^s\hat{B}_{d_1\ldots d_k}(x)}\,
  \frac{\text{W}_{\text{C}}[\phi_{d_1},\ldots,\phi_{d_s},\phi_n](x)}
  {\text{W}_{\text{C}}[\phi_{d_1},\ldots,\phi_{d_s}](x+1)}\nonumber\\
  &\phantom{\phi_{d_1\ldots d_s\,n}(x)}
  =(-1)^s\sqrt{\prod_{k=1}^s\hat{D}_{d_1\ldots d_k}(x+s+1-k)}\,
  \frac{\text{W}_{\text{C}}[\phi_{d_1},\ldots,\phi_{d_s},\phi_n](x)}
  {\text{W}_{\text{C}}[\phi_{d_1},\ldots,\phi_{d_s}](x)},\\
  &\phi_{d_1\ldots d_{s-1}\,n}(x)
  =\frac{\mathcal{A}_{d_1\ldots d_s}^{\text{op}\,\dagger}}
  {\mathcal{E}_n-\mathcal{E}_{d_s}}
  \phi_{d_1\ldots d_s\,n}(x)\ \ (s\geq 1),\\
  &\hat{B}_{d_1\ldots d_s}(x)=\hat{D}_{d_1\ldots d_s}(x+1)
  \Bigl(\frac{\phi_{d_1\ldots d_s}(x+1)}{\phi_{d_1\ldots d_s}(x)}\Bigr)^2,\\
  &\hat{B}_{d_1\ldots d_s}(x)\hat{D}_{d_1\ldots d_s}(x+1)
  =\hat{B}_{d_1\ldots d_{s-1}}(x+1)\hat{D}_{d_1\ldots d_{s-1}}(x+1)
  \ \ (s\geq1),\\
  &\hat{B}_{d_1\ldots d_s}(x)+\hat{D}_{d_1\ldots d_s}(x)+\mathcal{E}_{d_s}
  =\left\{\begin{array}{ll}
  \hat{B}_{d_1\ldots d_{s-1}}(x)+\hat{D}_{d_1\ldots d_{s-1}}(x+1)
  +\mathcal{E}_{d_{s-1}}&(s\geq2)\\[3pt]
  B(x)+D(x)&(s=1)
  \end{array}\right.,\\
  &\mathcal{H}^{\text{op}}_{d_1\ldots d_s}=
  \hat{\mathcal{A}}_{d_1\ldots d_{s+1}}^{\text{op}\,\dagger}
  \hat{\mathcal{A}}^{\text{op}}_{d_1\ldots d_{s+1}}+\mathcal{E}_{d_{s+1}},\\
  &\hat{B}_{d_1\ldots d_s}(x)=\sqrt{B(x+s-1)D(x+s)}\,
  \frac{w_{s-1}(x)}{w_{s-1}(x+1)}\frac{w_s(x+1)}{w_s(x)}\ \ (s\geq1),\\
  &\hat{D}_{d_1\ldots d_s}(x)=\sqrt{B(x-1)D(x)}\,
  \frac{w_{s-1}(x+1)}{w_{s-1}(x)}\frac{w_s(x-1)}{w_s(x)}\ \ (s\geq1),\\
  &\prod_{k=1}^s\hat{B}_{d_1\ldots d_k}(x)
  =\sqrt{\prod_{k=1}^sB(x+k-1)D(x+k)}\,
  \frac{w_s(x+1)}{w_s(x)},\\
  &\prod_{k=1}^s\hat{D}_{d_1\ldots d_k}(x+s+1-k)
  =\sqrt{\prod_{k=1}^sB(x+k-1)D(x+k)}\,
  \frac{w_s(x)}{w_s(x+1)},\\
  &\phi_{d_1\ldots d_s\,n}(x)
  =(-1)^s\Bigl(\prod_{k=1}^sB(x+k-1)D(x+k)\Bigr)^{\frac14}\,
  \frac{\text{W}_{\text{C}}[\phi_{d_1},\ldots,\phi_{d_s},\phi_n](x)}
  {\sqrt{w_s(x)w_s(x+1)}},
\end{align}
etc. Here $w_s(x)$ is
\begin{equation}
  w_s(x)=\text{W}_{\text{C}}[\phi_{d_1},\ldots,\phi_{d_s}](x).
\end{equation}

For $s=M$, we write $\mathcal{H}^{\text{op}}_{d_1\ldots d_M}
=\mathcal{H}^{\text{op}}_{\mathcal{D}}$,
$\hat{\mathcal{A}}^{\text{op}}_{d_1\ldots d_M}
=\hat{\mathcal{A}}^{\text{op}}_{\mathcal{D}}$,
$\hat{B}_{d_1\ldots d_M}(x)=\hat{B}_{\mathcal{D}}(x)$,
$\hat{D}_{d_1\ldots d_M}(x)=\hat{D}_{\mathcal{D}}(x)$,
$\phi_{d_1\ldots d_M\,n}(x)=\phi_{\mathcal{D}\,n}(x)$,
then we have
\begin{equation}
  \mathcal{H}^{\text{op}}_{\mathcal{D}}
  =\hat{\mathcal{A}}^{\text{op}}_{\mathcal{D}}
  \hat{\mathcal{A}}_{\mathcal{D}}^{\text{op}\,\dagger}
  +\mathcal{E}_{d_M},\quad
  \mathcal{H}^{\text{op}}_{\mathcal{D}}\phi_{\mathcal{D}\,n}(x)
  =\mathcal{E}_n\phi_{\mathcal{D}\,n}(x).
\end{equation}
By expressing this Hamiltonian in the standard form, we have
\begin{align}
  &\mathcal{H}^{\text{op}}_{\mathcal{D}}
  =\mathcal{A}_{\mathcal{D}}^{\text{op}\,\dagger}
  \mathcal{A}^{\text{op}}_{\mathcal{D}}+\mathcal{E}_{\mu},\quad
  \mu\stackrel{\text{def}}{=}\min(\mathbb{Z}_{\geq0}\backslash\mathcal{D}),
  \label{HopD}\\
  &\mathcal{A}^{\text{op}}_{\mathcal{D}}\stackrel{\text{def}}{=}\sqrt{B_{\mathcal{D}}(x)}
  -e^{\hat{\partial}}\sqrt{D_{\mathcal{D}}(x)},\quad
  \mathcal{A}_{\mathcal{D}}^{\text{op}\,\dagger}\stackrel{\text{def}}{=}\sqrt{B_{\mathcal{D}}(x)}
  -\sqrt{D_{\mathcal{D}}(x)}\,e^{-\hat{\partial}},\\
  &B_{\mathcal{D}}(x)\stackrel{\text{def}}{=}
  \sqrt{\hat{B}_{\mathcal{D}}(x+1)\hat{D}_{\mathcal{D}}(x+1)}\,
  \frac{\phi_{\mathcal{D}\,\mu}(x+1)}{\phi_{\mathcal{D}\,\mu}(x)},
  \label{BD}\\
  &D_{\mathcal{D}}(x)\stackrel{\text{def}}{=}
  \sqrt{\hat{B}_{\mathcal{D}}(x)\hat{D}_{\mathcal{D}}(x)}\,
  \frac{\phi_{\mathcal{D}\,\mu}(x-1)}{\phi_{\mathcal{D}\,\mu}(x)},
  \label{DD}\\
  &\phi_{\mathcal{D}\,n}(x)
  =(-1)^M\Bigl(\prod_{k=1}^MB(x+k-1)D(x+k)\Bigr)^{\frac14}\,
  \frac{\text{W}_{\text{C}}[\phi_{d_1},\ldots,\phi_{d_M},\phi_n](x)}
  {\sqrt{w_M(x)w_M(x+1)}},
  \label{phiDn}
\end{align}
and
\begin{align}
  &B_{\mathcal{D}}(x)+D_{\mathcal{D}}(x)+\mathcal{E}_{\mu}
  =\hat{B}_{\mathcal{D}}(x)+\hat{D}_{\mathcal{D}}(x+1)+\mathcal{E}_{d_M},\\
  &B_{\mathcal{D}}(x)=\sqrt{B(x+M)D(x+M+1)}\,
  \frac{w_M(x)}{w_M(x+1)}
  \frac{\text{W}_{\text{C}}[\phi_{d_1},\ldots,\phi_{d_M},\phi_{\mu}](x+1)}
  {\text{W}_{\text{C}}[\phi_{d_1},\ldots,\phi_{d_M},\phi_{\mu}](x)},\\
  &D_{\mathcal{D}}(x)=\sqrt{B(x-1)D(x)}\,
  \frac{w_M(x+1)}{w_M(x)}
  \frac{\text{W}_{\text{C}}[\phi_{d_1},\ldots,\phi_{d_M},\phi_{\mu}](x-1)}
  {\text{W}_{\text{C}}[\phi_{d_1},\ldots,\phi_{d_M},\phi_{\mu}](x)}.
\end{align}

We remark that the results so far have used \eqref{Hopphin2}, but not
\eqref{phin=phi0cPn}. The information \eqref{phin=phi0cPn} will be used in
\S\,\ref{sec:KAmiop} and \S\,\ref{sec:newmiop}.
In the following, we will denote $\phi_{\mathcal{D}\,n}(x)$ in \eqref{phiDn}
as $\phi_{\mathcal{D}\,n}^{\text{gen}}(x)$.

\subsubsection{orthogonality relations}
\label{sec:orthorel}

In \S\,\ref{sec:Dar} we consider a continuous variable $x\in\mathbb{R}$, and
the results are valid for a discrete variable $x\in\mathbb{Z}$.
Let us consider the case $x\in\mathbb{Z}$. This means the embedding of the
finite system in \S\,\ref{sec:orgsys} into the infinite system.
The matrices $e^{\pm\partial}$ satisfy $e^{\pm\partial}e^{\mp\partial}=1$.
The inner product of vectors $f$ and $g$ is defined by
$(f,g)=\sum_{x\in\mathbb{Z}}f(x)g(x)$.
Then the formal calculation gives
\begin{align}
  (\phi_{d_1\ldots d_s\,n},\phi_{d_1\ldots d_s\,m})
  &=(\hat{\mathcal{A}}_{d_1\ldots d_s}\phi_{d_1\ldots d_{s-1}\,n},
  \hat{\mathcal{A}}_{d_1\ldots d_s}\phi_{d_1\ldots d_{s-1}\,m})\nonumber\\
  &=(\hat{\mathcal{A}}_{d_1\ldots d_s}^{\dagger}
  \hat{\mathcal{A}}_{d_1\ldots d_s}
  \phi_{d_1\ldots d_{s-1}\,n},\phi_{d_1\ldots d_{s-1}\,m})\nonumber\\
  &=\bigl((\mathcal{H}_{d_1\ldots d_{s-1}}-\mathcal{E}_{d_s})
  \phi_{d_1\ldots d_{s-1}\,n},\phi_{d_1\ldots d_{s-1}\,m}\bigr)\nonumber\\
  &=(\mathcal{E}_n-\mathcal{E}_{d_s})
  (\phi_{d_1\ldots d_{s-1}\,n},\phi_{d_1\ldots d_{s-1}\,m}\bigr)
  \ \ (s\geq 1),
\end{align}
and this gives
\begin{equation}
  (\phi_{d_1\ldots d_s\,n},\phi_{d_1\ldots d_s\,m})
  =(\phi_n,\phi_m)\prod_{j=1}^s(\mathcal{E}_n-\mathcal{E}_{d_j}).
\end{equation}
If everything goes well after $M$-steps, we obtain
\begin{equation}
  \sum_{x\in\mathbb{Z}}\phi^{\text{gen}}_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})
  \phi^{\text{gen}}_{\mathcal{D}\,m}(x;\boldsymbol{\lambda})
  =\frac{\delta_{nm}}{d_n(\boldsymbol{\lambda})^2}\prod_{j=1}^M
  \bigl(\mathcal{E}_n(\boldsymbol{\lambda})-\mathcal{E}_{d_j}(\boldsymbol{\lambda})\bigr).
  \label{sumphigen}
\end{equation}
As discussed in \S\,\ref{sec:sabuneq}, $\phi_0(x;\boldsymbol{\lambda})^2$ for
$N\in\mathbb{Z}_{>0}$ and $x\in\mathbb{Z}$ is non-vanishing only for
$x=0,1,\ldots,N$.
The function $\phi^{\text{gen}}_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})$ contains
the function ``$\phi_0(x)$'' as a factor.
For the original system ($M=0$ case), the ``$\phi_0(x)$ factor'' is
$\phi_0(x;\boldsymbol{\lambda})$, and the sum in \eqref{sumphigen} is reduced to
$\sum_{x=0}^N$.
In \S\,\ref{sec:KAmiop} and \S\,\ref{sec:newmiop} we will see the following
situations.
For the systems in \S\,\ref{sec:KAmiop}, the ``$\phi_0(x)$ factor'' is
$\phi_0(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})$, whose parameter $N$ is $N-M$, and
the sum in \eqref{sumphigen} is reduced to $\sum_{x=0}^{N-M}$.
For the systems in \S\,\ref{sec:newmiop}, the ``$\phi_0(x)$ factor'' is
$\phi_0(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})$, whose parameter $N$ is $N+M$,
and the sum in \eqref{sumphigen} is reduced to $\sum_{x=-M}^N$.

\section{Krein-Adler Type Multi-Indexed Orthogonal Polynomials From
State-Deleting Darboux Transformations}
\label{sec:KAmiop}

The Darboux transformations with seed solutions $\phi_n(x)$ ($n\leq N$) are
studied in \cite{os22} and we call the resulting multi-indexed polynomials
Krein-Adler type (KA-type) multi-indexed polynomials.
Explicit forms of their difference equations, orthogonal relations etc. are
not given in \cite{os22}, and we present them here.

The multi-index set we will consider is
\begin{equation}
  \mathcal{D}=\{d_1,d_2,\ldots,d_M\}\ \ (0\leq d_j\leq N
  \text{ : mutually distinct}),
  \label{DdefKA}
\end{equation}
and $M$ should be $M\leq N$.
By using the information \eqref{phin=phi0cPn}, the eigenfunctions
$\phi_{\mathcal{D}\,n}^{\text{gen}}(x)$ are expressed in terms of $\phi_0$
and multi-indexed polynomials \cite{os22}.

The denominator polynomial $\Xi^{\text{KA}}_{\mathcal{D}}(\eta)$ and the
multi-indexed polynomial $P^{\text{KA}}_{\mathcal{D},n}(\eta)$ are defined by
(see Appendix\,\ref{app:poly} for $\boldsymbol{\delta}$ and Appendix\,\ref{app:varphiM}
for $\varphi_M(x)$)
\begin{align}
  \check{\Xi}^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})
  &\stackrel{\text{def}}{=}\Xi^{\text{KA}}_{\mathcal{D}}
  \bigl(\eta(x;\boldsymbol{\lambda}+(M-1)\boldsymbol{\delta});\boldsymbol{\lambda}\bigr)\nonumber\\
  &\stackrel{\text{def}}{=}\mathcal{C}^{\text{KA}}_{\mathcal{D}}(\boldsymbol{\lambda})^{-1}
  \varphi_M(x;\boldsymbol{\lambda})^{-1}
  \text{W}_{\text{C}}[\check{P}_{d_1},\ldots,\check{P}_{d_M}](x;\boldsymbol{\lambda}),
  \label{cXiDKAdef}\\
  \check{P}^{\text{KA}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})
  &\stackrel{\text{def}}{=} P^{\text{KA}}_{\mathcal{D},n}\bigl(\eta(x;\boldsymbol{\lambda}+M\boldsymbol{\delta});
  \boldsymbol{\lambda}\bigr)\quad(n=0,1,\ldots,N)\nonumber\\
  &\stackrel{\text{def}}{=}\mathcal{C}^{\text{KA}}_{\mathcal{D},n}(\boldsymbol{\lambda})^{-1}
  \varphi_{M+1}(x;\boldsymbol{\lambda})^{-1}
  \text{W}_{\text{C}}[\check{P}_{d_1},\ldots,\check{P}_{d_M},\check{P}_n]
  (x;\boldsymbol{\lambda}).
  \label{cPDnKAdef}
\end{align}
We remark that $\check{P}^{\text{KA}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})=0$
for $n\in\mathcal{D}$, and $\check{P}^{\text{KA}}_{\mathcal{D},n}(x)$ can be
defined for $n\in\mathbb{Z}_{\geq0}$ by replacing $\check{P}_n$ with
$\check{P}^{\text{monic}}_n$ for $n>N$, see \S\,\ref{sec:newmiop}.
They are
\begin{align}
  \check{\Xi}^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})&:
  \text{a polynomial of degree
  $\ell^{\text{KA}}_{\mathcal{D}}+M=\ell_{\mathcal{D}}$ in
  $\eta(x;\boldsymbol{\lambda}+(M-1)\boldsymbol{\delta})$},\nonumber\\
  \check{P}^{\text{KA}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})&:
  \text{a polynomial of degree $\ell^{\text{KA}}_{\mathcal{D}}+n$ in
  $\eta(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})$},
  \label{cPDnKAdeg}
\end{align}
where $\ell^{\text{KA}}_{\mathcal{D}}$ and $\ell_{\mathcal{D}}$ are given by
\eqref{ellDdef}.
The constants $\mathcal{C}^{\text{KA}}_{\mathcal{D}}(\boldsymbol{\lambda})$ and
$\mathcal{C}^{\text{KA}}_{\mathcal{D},n}(\boldsymbol{\lambda})$ are
determined by the following normalization conditions,
\begin{equation}
  \check{\Xi}^{\text{KA}}_{\mathcal{D}}(0;\boldsymbol{\lambda})
  =\check{P}^{\text{KA}}_{\mathcal{D},n}(0;\boldsymbol{\lambda})=1.
  \label{cPKADn(0)=1}
\end{equation}
In contrast to the case-(1) multi-indexed polynomials in \cite{os26},
the denominator polynomial and the multi-indexed polynomial are essentially
the same,
\begin{equation}
  \check{P}^{\text{KA}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})
  =\check{\Xi}^{\text{KA}}_{\mathcal{D}'}(x;\boldsymbol{\lambda}),\quad
  \mathcal{D}'=\{d_1,d_2,\ldots,d_M,n\}.
\end{equation}
Thus we have
\begin{equation}
  \mathcal{C}^{\text{KA}}_{\mathcal{D},n}(\boldsymbol{\lambda})
  =\mathcal{C}^{\text{KA}}_{\mathcal{D}'}(\boldsymbol{\lambda}),\quad
  \mathcal{D}'=\{d_1,d_2,\ldots,d_M,n\}.
  \label{cCXiP}
\end{equation}
The constant $\mathcal{C}^{\text{KA}}_{\mathcal{D}}(\boldsymbol{\lambda})$ is given by
((A.16) of \cite{os22})
\begin{equation}
  \mathcal{C}^{\text{KA}}_{\mathcal{D}}(\boldsymbol{\lambda})
  =(-1)^{\binom{M}{2}}\kappa^{-\binom{M}{3}}
  \prod_{1\leq j<k\leq M}
  \frac{\mathcal{E}_{d_k}(\boldsymbol{\lambda})-\mathcal{E}_{d_j}(\boldsymbol{\lambda})}
  {B(0;\boldsymbol{\lambda}+(j-1)\boldsymbol{\delta})},
\end{equation}
and \eqref{cCXiP} gives the constant
$\mathcal{C}^{\text{KA}}_{\mathcal{D},n}(\boldsymbol{\lambda})$,
\begin{equation}
  \mathcal{C}^{\text{KA}}_{\mathcal{D},n}(\boldsymbol{\lambda})
  =\mathcal{C}^{\text{KA}}_{\mathcal{D}}(\boldsymbol{\lambda})
  (-1)^M\kappa^{-\binom{M}{2}}\prod_{j=1}^M
  \frac{\mathcal{E}_n(\boldsymbol{\lambda})-\mathcal{E}_{d_j}(\boldsymbol{\lambda})}
  {B(0;\boldsymbol{\lambda}+(j-1)\boldsymbol{\delta})}.
  \label{CKADn=}
\end{equation}

Let us calculate $\phi_{\mathcal{D}\,n}^{\text{gen}}(x)$ \eqref{phiDn}.
{}From the property \eqref{WCformula1} and the definitions
\eqref{cXiDKAdef}--\eqref{cPDnKAdef}, we have
\begin{align}
  &w_M(x)=\text{W}_{\text{C}}[\phi_{d_1},\ldots,\phi_{d_M}](x)
  =\prod_{k=1}^M\phi_0(x+k-1;\boldsymbol{\lambda})\cdot
  \mathcal{C}^{\text{KA}}_{\mathcal{D}}(\boldsymbol{\lambda})
  \varphi_M(x;\boldsymbol{\lambda})
  \check{\Xi}^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda}),\\
  &\text{W}_{\text{C}}[\phi_{d_1},\ldots,\phi_{d_M},\phi_n](x)
  =\prod_{k=1}^{M+1}\phi_0(x+k-1;\boldsymbol{\lambda})\cdot
  \mathcal{C}^{\text{KA}}_{\mathcal{D},n}(\boldsymbol{\lambda})
  \varphi_{M+1}(x;\boldsymbol{\lambda})
  \check{P}^{\text{KA}}_{\mathcal{D},n}(x;\boldsymbol{\lambda}).
\end{align}
By using these and \eqref{phi0eq}, \eqref{varphiM+1M} and \eqref{B(x+j),D(x)},
we obtain
\begin{align}
  \phi_{\mathcal{D}\,n}^{\text{gen}}(x;\boldsymbol{\lambda})
  &=\kappa^{-\frac12\binom{M}{2}}\prod_{j=1}^M
  \frac{\mathcal{E}_n(\boldsymbol{\lambda})-\mathcal{E}_{d_j}(\boldsymbol{\lambda})}
  {\sqrt{B(0;\boldsymbol{\lambda}+(j-1)\boldsymbol{\delta})}}
  \times\phi^{\text{KA}}_{\mathcal{D}\,n}(x;\boldsymbol{\lambda}),
  \label{phiDngenKA}\\
  \phi^{\text{KA}}_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})
  &\stackrel{\text{def}}{=}\psi^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})
  \check{P}^{\text{KA}}_{\mathcal{D},n}(x;\boldsymbol{\lambda}),
  \label{phiDnKA}\\
  \psi^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \frac{\phi_0(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})}
  {\sqrt{\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})
  \check{\Xi}^{\text{KA}}_{\mathcal{D}}(x+1;\boldsymbol{\lambda})}}.
  \label{psiDKA}
\end{align}
Next let us calculate the potential functions \eqref{BD}--\eqref{DD}.
By using \eqref{varphiM+1M} and \eqref{B(x+j),D(x)}, we obtain
\begin{align}
  B_{\mathcal{D}}(x;\boldsymbol{\lambda})&=\kappa^MB(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})
  \frac{\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})}
  {\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x+1;\boldsymbol{\lambda})}
  \frac{\check{P}^{\text{KA}}_{\mathcal{D},\mu}(x+1;\boldsymbol{\lambda})}
  {\check{P}^{\text{KA}}_{\mathcal{D},\mu}(x;\boldsymbol{\lambda})},
  \label{BDKA}\\
  D_{\mathcal{D}}(x;\boldsymbol{\lambda})&=\kappa^MD(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})
  \frac{\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x+1;\boldsymbol{\lambda})}
  {\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})}
  \frac{\check{P}^{\text{KA}}_{\mathcal{D},\mu}(x-1;\boldsymbol{\lambda})}
  {\check{P}^{\text{KA}}_{\mathcal{D},\mu}(x;\boldsymbol{\lambda})}.
  \label{DDKA}
\end{align}
The deformed Hamiltonian \eqref{HopD} is similarity transformed to
\begin{align}
  \widetilde{\mathcal{H}}^{\text{KA\,op}}_{\mathcal{D}}(\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \psi^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})^{-1}\circ
  \mathcal{H}^{\text{op}}_{\mathcal{D}}(\boldsymbol{\lambda})\circ
  \psi^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})\nonumber\\
  &=\kappa^MB(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})
  \frac{\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})}
  {\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x+1;\boldsymbol{\lambda})}\Bigl(
  \frac{\check{P}^{\text{KA}}_{\mathcal{D},\mu}(x+1;\boldsymbol{\lambda})}
  {\check{P}^{\text{KA}}_{\mathcal{D},\mu}(x;\boldsymbol{\lambda})}
  -e^{\hat{\partial}}\Bigr)\nonumber\\
  &\quad+\kappa^MD(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})
  \frac{\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x+1;\boldsymbol{\lambda})}
  {\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})}\Bigl(
  \frac{\check{P}^{\text{KA}}_{\mathcal{D},\mu}(x-1;\boldsymbol{\lambda})}
  {\check{P}^{\text{KA}}_{\mathcal{D},\mu}(x;\boldsymbol{\lambda})}
  -e^{-\hat{\partial}}\Bigr)
  +\mathcal{E}_{\mu}(\boldsymbol{\lambda}).
  \label{HtDopdef}
\end{align}
{}From the result in \S\,\ref{sec:Dar}, the multi-indexed polynomials
$\check{P}^{\text{KA}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ satisfy the difference
equation,
\begin{equation}
  \widetilde{\mathcal{H}}^{\text{KA\,op}}_{\mathcal{D}}(\boldsymbol{\lambda})
  \check{P}^{\text{KA}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})
  =\mathcal{E}_n(\boldsymbol{\lambda})
  \check{P}^{\text{KA}}_{\mathcal{D},n}(x;\boldsymbol{\lambda}),
  \label{HtDopcPDnKA}
\end{equation}
for $x\in\mathbb{R}$.
We remark that this difference equation holds for any $\boldsymbol{\lambda}$.
In the expressions \eqref{psiDKA}--\eqref{DDKA}, they contain $\phi_0$, $B$
and $D$ with the parameter $\boldsymbol{\lambda}+M\boldsymbol{\delta}$, which means that the
parameter $N$ is $N-M$.
Therefore the matrix
$\widetilde{\mathcal{H}}^{\text{KA}}_{\mathcal{D}}(\boldsymbol{\lambda})$ is of order
$N-M+1$,
\begin{align}
  \widetilde{\mathcal{H}}^{\text{KA}}_{\mathcal{D}}(\boldsymbol{\lambda})&=
  \bigl(\widetilde{\mathcal{H}}^{\text{KA}}_{\mathcal{D}\,x,y}
  (\boldsymbol{\lambda})\bigr)_{x,y=0,1,\ldots,N-M},\nonumber\\
  \widetilde{\mathcal{H}}^{\text{KA}}_{\mathcal{D}}(\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \widetilde{\mathcal{H}}^{\text{KA\,op}}_{\mathcal{D}}(\boldsymbol{\lambda})
  \bigm|_{\hat{\partial}\to\partial}.
  \label{HtDKAdef}
\end{align}
Since $B(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})$ and $D(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})$ 
satisfy the boundary conditions $B(N-M;\boldsymbol{\lambda}+M\boldsymbol{\delta})=0$ at $x=N-M$
and $D(0;\boldsymbol{\lambda}+M\boldsymbol{\delta})=0$ at $x=0$, the difference equation
\eqref{HtDopcPDnKA} implies the following.
The matrix eigenvalue problem for $\widetilde{\mathcal{H}}_{\mathcal{D}}$ is
solved by the multi-indexed polynomials
$\check{P}^{\text{KA}}_{\mathcal{D},n}(x)$,
\begin{equation}
  \widetilde{\mathcal{H}}^{\text{KA}}_{\mathcal{D}}(\boldsymbol{\lambda})
  \check{P}^{\text{KA}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})
  =\mathcal{E}_n(\boldsymbol{\lambda})
  \check{P}^{\text{KA}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})
  \ \ (n\in\{0,1,\ldots,N\}\backslash\mathcal{D}).
\end{equation}

The orthogonality relations for $\check{P}^{\text{KA}}_{\mathcal{D},n}(x)$
can be read from \eqref{sumphigen}.
{}From \eqref{phiDngenKA}--\eqref{psiDKA}, 
$\phi_{\mathcal{D}\,n}^{\text{gen}}(x;\boldsymbol{\lambda})$ contains
$\phi_0(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})$, which does not vanish
at $x=0,1,\ldots,N-M$ for $x\in\mathbb{Z}$.
So the sum in \eqref{sumphigen} is reduced to $\sum_{x=0}^{N-M}$.
We obtain the orthogonality relations,
\begin{align}
  &\quad\sum_{x=0}^{N-M}\frac{\phi_0(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})^2}
  {\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})
  \check{\Xi}^{\text{KA}}_{\mathcal{D}}(x+1;\boldsymbol{\lambda})}
  \check{P}^{\text{KA}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})
  \check{P}^{\text{KA}}_{\mathcal{D},m}(x;\boldsymbol{\lambda})\nonumber\\
  &=\frac{\delta_{nm}}{d^{\text{KA}}_{\mathcal{D},n}(\boldsymbol{\lambda})^2}
  \ \ (n,m\in\{0,1,\ldots,N\}\backslash\mathcal{D}),
  \label{orthoKA}
\end{align}
where the normalization constants
$d^{\text{KA}}_{\mathcal{D},n}(\boldsymbol{\lambda})^2$ are obtained from
\eqref{sumphigen} and \eqref{phiDngenKA}--\eqref{psiDKA} as follows,
\begin{equation}
  d^{\text{KA}}_{\mathcal{D},n}(\boldsymbol{\lambda})^2
  =d_n(\boldsymbol{\lambda})^2\kappa^{-\binom{M}{2}}
  \prod_{j=1}^M\frac{\mathcal{E}_n(\boldsymbol{\lambda})-\mathcal{E}_{d_j}(\boldsymbol{\lambda})}
  {B(0;\boldsymbol{\lambda}+(j-1)\boldsymbol{\delta})}.
\end{equation}
In order to call the relations \eqref{orthoKA} truly orthogonality relations,
the positivity of the weight factor is necessary.
The positivity of $\phi_0(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})^2$ can be easily
achieved by choosing parameters $\boldsymbol{\lambda}$ so that
$B(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})$ and $D(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})$ are
positive (except for the boundaries). The positivity of
$\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})
\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x+1;\boldsymbol{\lambda})$
(for $x=0,1,\ldots,N-M$) can be achieved by imposing the following
Krein-Adler condition on $\mathcal{D}$ \cite{os22},
\begin{equation}
  \prod_{j=1}^M(m-d_j)\geq 0\ \ (\forall m\in\mathbb{Z}_{\geq0}).
  \label{KA}
\end{equation}
We remark that this condition \eqref{KA} is rewritten as
\begin{equation}
  e_j-e_{j-1}\equiv 1\ \ (\text{mod }2)\ \ (j=2,3,\ldots),
  \label{KA2}
\end{equation}
where $\mathbb{Z}_{\geq0}\backslash\mathcal{D}=\{e_1,e_2,\ldots\}$
with $e_1<e_2<\cdots$.
In this situation, the Hamiltonian of the deformed system
\begin{align}
  \mathcal{H}^{\text{KA}}_{\mathcal{D}}(\boldsymbol{\lambda})&=
  \bigl(\mathcal{H}^{\text{KA}}_{\mathcal{D}\,x,y}
  (\boldsymbol{\lambda})\bigr)_{x,y=0,1,\ldots,N-M},\nonumber\\
  \mathcal{H}^{\text{KA}}_{\mathcal{D}}(\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \psi^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})\circ
  \widetilde{\mathcal{H}}^{\text{KA}}_{\mathcal{D}}(\boldsymbol{\lambda})\circ
  \psi^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})^{-1}\nonumber\\
  &=-\sqrt{B_{\mathcal{D}}(x;\boldsymbol{\lambda})D_{\mathcal{D}}(x+1;\boldsymbol{\lambda})}
  \,e^{\partial}
  -\sqrt{B_{\mathcal{D}}(x-1;\boldsymbol{\lambda})D_{\mathcal{D}}(x;\boldsymbol{\lambda})}
  \,e^{-\partial}\nonumber\\
  &\quad+B_{\mathcal{D}}(x;\boldsymbol{\lambda})+D_{\mathcal{D}}(x;\boldsymbol{\lambda})
  +\mathcal{E}_{\mu}(\boldsymbol{\lambda})
\end{align}
is a real symmetric matrix, and the matrix eigenvalue problem for
$\mathcal{H}^{\text{KA}}_{\mathcal{D}}$ is solved by
$\phi^{\text{KA}}_{\mathcal{D}\,n}(x)$ \eqref{phiDnKA},
\begin{equation}
  \mathcal{H}^{\text{KA}}_{\mathcal{D}}(\boldsymbol{\lambda})
  \phi^{\text{KA}}_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})
  =\mathcal{E}_n(\boldsymbol{\lambda})
  \phi^{\text{KA}}_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})
  \ \ (n\in\{0,1,\ldots,N\}\backslash\mathcal{D}).
\end{equation}
However, we remark that the relations \eqref{orthoKA} hold for any
$\boldsymbol{\lambda}$ and $\mathcal{D}$, because the expressions are algebraic and
the sum in \eqref{orthoKA} is a finite sum.

We have shown the state-deleting property of the Darboux transformations by
using the embedding of the finite system into the infinite system.
On the other hand, in \cite{os22}, the state-deleting property was shown without
using such a trick.
It was shown that each step of the Darboux transformation reduces the size of
the matrix (Hamiltonian) by one.

The original systems in \S\,\ref{sec:orgsys} have shape invariance.
As its consequence, the forward and backward shift relations hold.
Here we consider the forward shift relation,
\begin{equation}
  \mathcal{F}^{\text{op}}(\boldsymbol{\lambda})=B(0;\boldsymbol{\lambda})\varphi(x)^{-1}
  (1-e^{\hat{\partial}}),\quad
  \mathcal{F}^{\text{op}}(\boldsymbol{\lambda})\check{P}_n(x;\boldsymbol{\lambda})
  =\mathcal{E}_n(\boldsymbol{\lambda})\check{P}_{n-1}(x;\boldsymbol{\lambda}+\boldsymbol{\delta}).
  \label{FopcP=}
\end{equation}
By using this and properties of determinant and the notation \eqref{D[i]def},
we can show the following for $d_j\geq1$,
\begin{equation}
  \check{P}^{\text{KA}}_{\mathcal{D},0}(x;\boldsymbol{\lambda})
  =\check{\Xi}^{\text{KA}}_{\mathcal{D}^{[-1]}}(x;\boldsymbol{\lambda}+\boldsymbol{\delta}).
  \label{cPKAD0}
\end{equation}
A slightly different but similar relation exists for the case-(1) multi-indexed
polynomials in \cite{os26}.
Let us consider the case $d_j\geq1$.
Since this condition means $\mu=0$, \eqref{HtDopdef} is expressed as
\begin{align}
  \widetilde{\mathcal{H}}^{\text{KA\,op}}_{\mathcal{D}}(\boldsymbol{\lambda})
  &=\kappa^MB(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})
  \frac{\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})}
  {\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x+1;\boldsymbol{\lambda})}\Bigl(
  \frac{\check{\Xi}^{\text{KA}}_{\mathcal{D}^{[-1]}}
  (x+1;\boldsymbol{\lambda}+\boldsymbol{\delta})}
  {\check{\Xi}^{\text{KA}}_{\mathcal{D}^{[-1]}}(x;\boldsymbol{\lambda}+\boldsymbol{\delta})}
  -e^{\hat{\partial}}\Bigr)\nonumber\\
  &\quad+\kappa^MD(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})
  \frac{\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x+1;\boldsymbol{\lambda})}
  {\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})}\Bigl(
  \frac{\check{\Xi}^{\text{KA}}_{\mathcal{D}^{[-1]}}
  (x-1;\boldsymbol{\lambda}+\boldsymbol{\delta})}
  {\check{\Xi}^{\text{KA}}_{\mathcal{D}^{[-1]}}(x;\boldsymbol{\lambda}+\boldsymbol{\delta})}
  -e^{-\hat{\partial}}\Bigr).
  \label{HtDop0}
\end{align}
By replacing $\hat{\partial}$ with $\partial$, an expression of
$\widetilde{\mathcal{H}}^{\text{KA}}_{\mathcal{D}}(\boldsymbol{\lambda})$ is obtained.
This resembles the similarity transformed Hamiltonian for the case-(1)
multi-indexed polynomials in \cite{os26}.

\section{New Multi-Indexed Orthogonal Polynomials From State-Adding
Darboux Transformations}
\label{sec:newmiop}

In this section we consider the Darboux transformations with seed solutions
$\phi_n(x)$ ($n>N$) \cite{os40} and present new multi-indexed orthogonal
polynomials $\check{Q}_{\mathcal{D}',n}(x)$.

The multi-index set we will consider is
\begin{align}
  \mathcal{D}&=\{d_1,d_2,\ldots,d_M\}\ \ (d_j>N
  \text{ : mutually distinct}),\\
  \label{Ddefnew}
  &\quad d_j=N+1+m_j,\ \ m_j\in\mathbb{Z}_{\geq0},
\end{align}
and we assume $M\leq N$.
In this section we use the following notations,
\begin{align}
  \mathcal{D}'&=\mathcal{D}^{[-N-1]}=\{m_1,m_2,\ldots,m_M\},
  \label{D'def}\\
  \boldsymbol{\lambda}'&=\boldsymbol{\lambda}+(N+1)(\boldsymbol{\delta}+\boldsymbol{\bar{\delta}}),
  \label{lambda'def}
\end{align}
where $\boldsymbol{\bar{\delta}}$ is given in Appendix\,\ref{app:poly}.
We remark that $N$ is one element of $\boldsymbol{\lambda}$.

The twelve polynomials $\check{P}_n(x)$ in Appendix\,\ref{app:poly} with the
normalization \eqref{cPn0=1} are ill-defined for $n>N$, because there is a
factor $1/(-N)_k$ or $1/(q^{-N}\,;q)_k$ in the sum $\sum_{k=0}^n$ of the
(basic) hypergeometric series expansion.
To avoid this, let us consider the monic polynomial,
\begin{equation}
  \check{P}^{\text{monic}}_n(x;\boldsymbol{\lambda})
  =c_n(\boldsymbol{\lambda})^{-1}\check{P}_n(x;\boldsymbol{\lambda}),
\end{equation}
where $c_n(\boldsymbol{\lambda})$ is the coefficient of the highest degree term,
\begin{equation}
  \check{P}_n(x;\boldsymbol{\lambda})
  =c_n(\boldsymbol{\lambda})\eta(x;\boldsymbol{\lambda})^n+(\text{lower degree terms}).
  \label{cndef}
\end{equation}
The monic polynomial $\check{P}^{\text{monic}}_n(x;\boldsymbol{\lambda})$ is
well-defined for $n\in\mathbb{Z}_{\geq0}$.
Explicit forms of $c_n(\boldsymbol{\lambda})$ are given in Appendix\,\ref{app:poly} and
the universal expression of $c_n(\boldsymbol{\lambda})$ is ((A.14) in \cite{os22})
\begin{equation}
  c_n(\boldsymbol{\lambda})=(-1)^n\kappa^{-\binom{n}{2}}
  \prod_{j=1}^n\frac{\mathcal{E}_n(\boldsymbol{\lambda})
  -\mathcal{E}_{j-1}(\boldsymbol{\lambda})}
  {\eta(j;\boldsymbol{\lambda})B(0;\boldsymbol{\lambda}+(j-1)\boldsymbol{\delta})}.
  \label{cnuniv}
\end{equation}
For the monic polynomial, the orthogonality relations \eqref{orthocPn} become
\begin{equation}
  \sum_{x=0}^N\phi_0(x;\boldsymbol{\lambda})^2
  \check{P}^{\text{monic}}_n(x;\boldsymbol{\lambda})
  \check{P}^{\text{monic}}_m(x;\boldsymbol{\lambda})
  =\frac{\delta_{nm}}{d^{\text{monic}}_n(\boldsymbol{\lambda})^2},\quad
  d^{\text{monic}}_n(\boldsymbol{\lambda})
  =c_n(\boldsymbol{\lambda})d_n(\boldsymbol{\lambda}).
  \label{orthocPnmonic}
\end{equation}

The polynomials 
$\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})$ \eqref{cXiDKAdef} and
$\check{P}^{\text{KA}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ \eqref{cPDnKAdef} with
the normalization conditions \eqref{cPKADn(0)=1} are defined for $d_j,n\leq N$.
By replacing $\check{P}_{d_j}$ and $\check{P}_n$ with
$\check{P}^{\text{monic}}_{d_j}$ and $\check{P}^{\text{monic}}_n$, their
definitions are extended to $d_j,n\in\mathbb{Z}_{\geq0}$.
The monic version of these polynomials are given by
\begin{align}
  \check{\Xi}^{\text{KA\,monic}}_{\mathcal{D}}(x;\boldsymbol{\lambda})
  &=c^{\eta}_{\mathcal{D}}(\boldsymbol{\lambda})^{-1}
  \varphi_M(x;\boldsymbol{\lambda})^{-1}
  \text{W}_{\text{C}}[\check{P}^{\text{monic}}_{d_1},\ldots,
  \check{P}^{\text{monic}}_{d_M}](x;\boldsymbol{\lambda}),
  \label{cXiDKAmonic}\\
  \check{P}^{\text{KA\,monic}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})
  &=c^{\eta}_{\mathcal{D},n}(\boldsymbol{\lambda})^{-1}
  \varphi_{M+1}(x;\boldsymbol{\lambda})^{-1}
  \text{W}_{\text{C}}[\check{P}^{\text{monic}}_{d_1},\ldots,
  \check{P}^{\text{monic}}_{d_M},
  \check{P}^{\text{monic}}_n](x;\boldsymbol{\lambda}),
  \label{cPDnKAmonic}
\end{align}
where the constants $c^{\eta}_{\mathcal{D}}(\boldsymbol{\lambda})$ and
$c^{\eta}_{\mathcal{D},n}(\boldsymbol{\lambda})$ are given by \eqref{cetaD} and
\eqref{cetaDn} respectively.
These monic polynomials and their non-monic versions are related as
\begin{align}
  c^{\eta}_{\mathcal{D}}(\boldsymbol{\lambda})
  \check{\Xi}^{\text{KA\,monic}}_{\mathcal{D}}(x;\boldsymbol{\lambda})
  &=\mathcal{C}^{\text{KA}}_{\mathcal{D}}(\boldsymbol{\lambda})
  \frac{\check{\Xi}^{\text{KA}}_{\mathcal{D}}(x;\boldsymbol{\lambda})}
  {\prod_{j=1}^Mc_{d_j}(\boldsymbol{\lambda})},
  \label{cXiDKAmonicnonmonic}\\
  c^{\eta}_{\mathcal{D},n}(\boldsymbol{\lambda})
  \check{P}^{\text{KA\,monic}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})
  &=\mathcal{C}^{\text{KA}}_{\mathcal{D},n}(\boldsymbol{\lambda})
  \frac{\check{P}^{\text{KA}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})}
  {\prod_{j=1}^Mc_{d_j}(\boldsymbol{\lambda})\cdot c_n(\boldsymbol{\lambda})}.
  \label{cPDnKAmonicnonmonic}
\end{align}
This $\check{P}^{\text{KA\,monic}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$
satisfies the difference equation \eqref{HtDopcPDnKA} by replacing
$\check{\Xi}^{\text{KA}}_{\mathcal{D}}$ and
$\check{P}^{\text{KA}}_{\mathcal{D},n}$ with
$\check{\Xi}^{\text{KA\,monic}}_{\mathcal{D}}$ and
$\check{P}^{\text{KA\,monic}}_{\mathcal{D},n}$ in \eqref{HtDopdef}.

For $n>N$, by setting $n=N+1+m$ ($m\in\mathbb{Z}_{\geq 0}$),
the monic polynomial $\check{P}_n(x;\boldsymbol{\lambda})$ has a factorization property
\cite{os40} (see \eqref{eta(x-N-1)} with $M=1$),
\begin{equation}
  \check{P}^{\text{monic}}_n(x;\boldsymbol{\lambda})
  =\Lambda(x;\boldsymbol{\lambda})\rho^{(N+1)m}
  \check{P}^{\text{monic}}_m(x-N-1;\boldsymbol{\lambda}'),
  \label{fac_cPn}
\end{equation}
where $\Lambda(x;\boldsymbol{\lambda})$ is defined by
\begin{equation}
  \Lambda(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=}\prod_{k=0}^N\bigl(
  \eta(x;\boldsymbol{\lambda})-\eta(k;\boldsymbol{\lambda})\bigr),
  \label{Lambdadef}
\end{equation}
and its properties are given in Appendix\,\ref{app:Lam}.
{}From \eqref{cXiDKAmonic} and \eqref{fac_cPn}, we have
\begin{align}
  &\quad\check{\Xi}^{\text{KA\,monic}}_{\mathcal{D}}(x;\boldsymbol{\lambda})\nonumber\\
  &=c^{\eta}_{\mathcal{D}}(\boldsymbol{\lambda})^{-1}
  \varphi_M(x;\boldsymbol{\lambda})^{-1}
  \prod_{j=1}^M\Lambda(x+j-1;\boldsymbol{\lambda})\cdot
  \prod_{k=1}^M\rho^{(N+1)m_k}\nonumber\\
  &\quad\times
  \text{W}_{\text{C}}[\check{P}^{\text{monic}}_{m_1},\ldots,
  \check{P}^{\text{monic}}_{m_M}](x-N-1;\boldsymbol{\lambda}'),\nonumber\\
  &=\frac{c^{\eta}_{\mathcal{D}'}(\boldsymbol{\lambda}')}
  {c^{\eta}_{\mathcal{D}}(\boldsymbol{\lambda})}
  \frac{\varphi_M(0;\boldsymbol{\lambda}')}{\varphi_M(N+1;\boldsymbol{\lambda})}
  \prod_{j=1}^M\Lambda(x+j-1;\boldsymbol{\lambda})\cdot
  \prod_{k=1}^M\rho^{(N+1)m_k}\cdot
  \check{\Xi}^{\text{KA\,monic}}_{\mathcal{D}'}(x-N-1;\boldsymbol{\lambda}'),\!
  \label{cXiDKAmonic2}
\end{align}
where \eqref{WCformula1} and \eqref{varphiM/varphiM} are used.
So $\check{\Xi}^{\text{KA\,monic}}_{\mathcal{D}}(x;\boldsymbol{\lambda})$ is
divisible by $\prod_{j=1}^M\Lambda(x+j-1;\boldsymbol{\lambda})$, which is a polynomial
in $\eta(x;\boldsymbol{\lambda}+(M-1)\boldsymbol{\delta})$, see \eqref{Lam_poly}
(with $M\to M-1$).
The relation \eqref{cPKAD0} implies that
$\check{P}^{\text{KA\,monic}}_{\mathcal{D},0}(x;\boldsymbol{\lambda})$ is
divisible by $\prod_{j=1}^M\Lambda(x+j-1;\boldsymbol{\lambda}+\boldsymbol{\delta})$,
which is a polynomial in $\eta(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})$ \eqref{Lam_poly2}.
It is expected that
$\check{P}^{\text{KA\,monic}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ is also
divisible by $\prod_{j=1}^M\Lambda(x+j-1;\boldsymbol{\lambda}+\boldsymbol{\delta})$, and
this is indeed the case.

Based on $\check{P}^{\text{KA\,monic}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$,
we define new multi-indexed polynomials
$\check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})$ for
$n\in\mathbb{Z}_{\geq0}$.
First, let us define $\check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})$
for $n\in\mathbb{Z}_{\geq0}\backslash\mathcal{D}$ as follows
(see \eqref{Lam_poly2}),
\begin{align}
  \check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})
  &\stackrel{\text{def}}{=}\rho^{\binom{M}{2}N}
  \prod_{j=1}^M\Lambda(x+j-1;\boldsymbol{\lambda}+\boldsymbol{\delta})^{-1}\cdot
  \check{P}^{\text{KA\,monic}}_{\mathcal{D},n}(x;\boldsymbol{\lambda})
  \label{cQmonicdef}\\
  &=c^{\eta}_{\mathcal{D},n}(\boldsymbol{\lambda})^{-1}\rho^{\binom{M}{2}N}
  \prod_{j=1}^M\Lambda(x+j-1;\boldsymbol{\lambda}+\boldsymbol{\delta})^{-1}\nonumber\\
  &\quad\times\varphi_{M+1}(x;\boldsymbol{\lambda})^{-1}
  \text{W}_{\text{C}}[\check{P}^{\text{monic}}_{d_1},\ldots,
  \check{P}^{\text{monic}}_{d_M},\check{P}^{\text{monic}}_n](x;\boldsymbol{\lambda}).
  \label{cQmonic}
\end{align}
For $n>N$ and $n\not\in\mathcal{D}$ (we set $n=N+1+m$), by using
\eqref{fac_cPn} and \eqref{WCformula1}, we have
\begin{align}
  &\quad\check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})\nonumber\\
  &=c^{\eta}_{\mathcal{D},n}(\boldsymbol{\lambda})^{-1}\rho^{\binom{M}{2}N}
  \prod_{k=1}^M\rho^{(N+1)m_k}\cdot\rho^{(N+1)m}
  \frac{\prod_{j=1}^{M+1}\Lambda(x+j-1;\boldsymbol{\lambda})}
  {\prod_{j=1}^M\Lambda(x+j-1;\boldsymbol{\lambda}+\boldsymbol{\delta})}
  \frac{\varphi_{M+1}(x-N-1;\boldsymbol{\lambda}')}{\varphi_{M+1}(x;\boldsymbol{\lambda})}\nonumber\\
  &\quad\times
  \varphi_{M+1}(x-N-1;\boldsymbol{\lambda}')^{-1}
  \text{W}_{\text{C}}[\check{P}^{\text{monic}}_{m_1},\ldots,
  \check{P}^{\text{monic}}_{m_M},\check{P}^{\text{monic}}_m]
  (x-N-1;\boldsymbol{\lambda}'),
  \label{cQmonic_n=N+1+m}
\end{align}
which vanishes for $x\in\{-M,-M+1,\ldots,N\}$ by \eqref{prodLam/Lam} and
\eqref{varphiM/varphiM}.
Next, let us define $\check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})$
for $n\in\mathcal{D}$, based on the expression \eqref{cQmonic}.
For $n=d_i\in\mathcal{D}$, the Casoratian $\text{W}_{\text{C}}[\cdots]$ in
\eqref{cQmonic} vanishes, but
$c^{\eta}_{\mathcal{D},n}(\boldsymbol{\lambda})$ \eqref{cetaDn} also vanishes.
So, by taking certain appropriate ``$n\to d_i$ limit'', we may obtain a finite
quantity.
Following the prescription explained in Appendix\,\ref{app:n->di}, we define
$\check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})$ with
$n=d_i\in\mathcal{D}$ as follows,
\begin{align}
  \check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})
  &\stackrel{\text{def}}{=}\rho^{-n}
  c^{\prime\,\eta}_{\mathcal{D},n}(\boldsymbol{\lambda})^{-1}\rho^{\binom{M}{2}N}
  \prod_{j=1}^M\Lambda(x+j-1;\boldsymbol{\lambda}+\boldsymbol{\delta})^{-1}
  \label{cQmonic2}\\
  &\quad\times\varphi_{M+1}(x;\boldsymbol{\lambda})^{-1}
  \left|
  \begin{array}{cccc}
  \check{P}^{\text{monic}}_{d_1}(x_1;\boldsymbol{\lambda})&\cdots
  &\check{P}^{\text{monic}}_{d_M}(x_1;\boldsymbol{\lambda})
  &R_1(x)\\
  \check{P}^{\text{monic}}_{d_1}(x_2;\boldsymbol{\lambda})&\cdots
  &\check{P}^{\text{monic}}_{d_M}(x_2;\boldsymbol{\lambda})
  &R_2(x)\\
  \vdots&\cdots&\vdots&\vdots\\
  \check{P}^{\text{monic}}_{d_1}(x_{M+1};\boldsymbol{\lambda})&\cdots
  &\check{P}^{\text{monic}}_{d_M}(x_{M+1};\boldsymbol{\lambda})
  &R_{M+1}(x)
  \end{array}\right|,
  \nonumber
\end{align}
where $x_j=x+j-1$ and $R_j(x)$ are given by
\begin{align}
  &\quad R_j(x)=R_{M,j,n}(x;\boldsymbol{\lambda})
  \label{Rjdef}\\
  &\stackrel{\text{def}}{=}\rho^{(N+1)m_i}Y_{M,j}(x;\boldsymbol{\lambda})
  \check{P}^{\text{monic}}_{m_i}(x_j-N-1;\boldsymbol{\lambda}')
  +\check{\tilde{P}}^{(1)}_n(x_j;\boldsymbol{\lambda})
  +\rho^{(N+1)m_i}\Lambda(x_j;\boldsymbol{\lambda})
  \check{\tilde{P}}^{(2)}_{m_i}(x_j;\boldsymbol{\lambda}),
  \nonumber
\end{align}
and $Y_{M,j}(x;\boldsymbol{\lambda})$, $\check{\tilde{P}}^{(1)}_n(x;\boldsymbol{\lambda})$,
$\check{\tilde{P}}^{(2)}_m(x;\boldsymbol{\lambda})$ and
$c^{\prime\,\eta}_{\mathcal{D},n}(\boldsymbol{\lambda})$ are defined by \eqref{YMj},
\eqref{ctP1def}, \eqref{ctP2def} and \eqref{cetaDn'}, respectively.
We remark that the contribution to
$\check{Q}^{\text{monic}}_{\mathcal{D}',n}(x)$ \eqref{cQmonic2} coming from
the third term of $R_j(x)$ \eqref{Rjdef} (the term containing
$\check{\tilde{P}}^{(2)}_{m_i}(x)$ : contribution (c) in
Appendix\,\ref{app:(c)cPm}) has the following form,
\begin{equation*}
  \Lambda(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})\times
  \bigl(\text{a polynomial in $\eta(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})$}\bigr),
\end{equation*}
which vanishes for $x\in\{-M,-M+1,\ldots,N\}$.
So, when dealing with $\check{Q}^{\text{monic}}_{\mathcal{D}',d_i}(x)$ 
for $x\in\{-M,-M+1,\ldots,N\}$, we can ignore the third term of $R_j(x)$.
These $\check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})$
($n\in\mathbb{Z}_{\geq0}$) are
\begin{align}
  &\check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda}):
  \text{a monic polynomial of degree $\ell_{\mathcal{D}'}+n$ in
  $\eta(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})$},\nonumber\\
  &\check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})
  \stackrel{\text{def}}{=} Q^{\text{monic}}_{\mathcal{D}',n}\bigl(
  \eta(x;\boldsymbol{\lambda}+M\boldsymbol{\delta});\boldsymbol{\lambda}\bigr),
  \label{cQD'nmonic}
\end{align}
where $\ell_{\mathcal{D}'}$ is given by \eqref{ellDdef}.
We define the non-monic version $\check{Q}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})$
as follows,
\begin{equation}
  \check{Q}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})
  \stackrel{\text{def}}{=}\frac{\check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})}
  {\check{Q}^{\text{monic}}_{\mathcal{D}',n}(-M;\boldsymbol{\lambda})}
  \ \ \bigl(\Rightarrow\check{Q}_{\mathcal{D}',n}(-M;\boldsymbol{\lambda})=1\bigr).
\end{equation}

Let us calculate $\phi_{\mathcal{D}\,n}^{\text{gen}}(x)$ \eqref{phiDn} for
$n\in\mathbb{Z}_{\geq0}\backslash\mathcal{D}$.
By using \eqref{phiDngenKA},
\eqref{cXiDKAmonicnonmonic}--\eqref{cPDnKAmonicnonmonic},
\eqref{cXiDKAmonic2}, \eqref{cQmonicdef}, \eqref{phi0Lam} and \eqref{CKADn=},
we obtain
\begin{align}
  \phi_{\mathcal{D}\,n}^{\text{gen}}(x;\boldsymbol{\lambda})
  &=(-1)^Mc_n(\boldsymbol{\lambda})\kappa^{\frac12\binom{M}{2}}
  \rho^{-(2N+1)\binom{M}{2}}
  \prod_{j=1}^M\sqrt{B(0;\boldsymbol{\lambda}+(j-1)\boldsymbol{\delta})}\nonumber\\
  &\quad\times
  \frac{c^{\eta}_{\mathcal{D},n}(\boldsymbol{\lambda})}
  {c^{\eta}_{\mathcal{D}'}(\boldsymbol{\lambda}')}
  \frac{\Lambda_M(0;\boldsymbol{\lambda})^{\frac12}}
  {\phi_0(M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})}
  \frac{\varphi_M(N+1;\boldsymbol{\lambda})}{\varphi_M(0;\boldsymbol{\lambda}')}
  \times\phi^{Q\,\text{monic}}_{\mathcal{D}'\,n}(x;\boldsymbol{\lambda}),
  \label{phiDngenQ}\\
  \phi^{Q\,\text{monic}}_{\mathcal{D}'\,n}(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \psi^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda})
  \check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda}),
  \label{phiDnQ}\\
  \psi^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \frac{\phi_0(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})}
  {\sqrt{\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda})
  \check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x+1;\boldsymbol{\lambda})}}.
  \label{psiDQ}
\end{align}
Here $\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda})$ is given by
\begin{align}
  &\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda})
  \stackrel{\text{def}}{=}\rho^{(N+1)\ell_{\mathcal{D}'}}
  \check{\Xi}^{\text{KA\,monic}}_{\mathcal{D}'}(x-N-1;\boldsymbol{\lambda}')\nonumber\\
  &\qquad:\text{a monic polynomial of degree $\ell_{\mathcal{D}'}$ in
  $\eta(x;\boldsymbol{\lambda}+(M-1)\boldsymbol{\delta})$},
\end{align}
which is shown by \eqref{cPDnKAdeg} and
\begin{equation}
  \eta(x-N-1;\boldsymbol{\lambda}'+(M-1)\boldsymbol{\delta})
  =\rho^{-N-1}\bigl(\eta(x;\boldsymbol{\lambda}+(M-1)\boldsymbol{\delta})
  -\eta(N+1;\boldsymbol{\lambda}+(M-1)\boldsymbol{\delta})\bigr).
  \label{eta(x-N-1)}
\end{equation}
Since the parameter $N$ for $\boldsymbol{\lambda}+\boldsymbol{\delta}$ is $N-1$, we have
\begin{equation}
  \check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda}+\boldsymbol{\delta})
  =\rho^{N\ell_{\mathcal{D}'}}
  \check{\Xi}^{\text{KA\,monic}}_{\mathcal{D}'}
  \bigl(x-N;\boldsymbol{\lambda}+\boldsymbol{\delta}+N(\boldsymbol{\delta}+\boldsymbol{\bar{\delta}})\bigr)
  =\rho^{N\ell_{\mathcal{D}'}}
  \check{\Xi}^{\text{KA\,monic}}_{\mathcal{D}'}
  (x-N;\boldsymbol{\lambda}'-\boldsymbol{\bar{\delta}}).
\end{equation}
We define the non-monic version $\check{\Xi}^Q_{\mathcal{D}'}(x;\boldsymbol{\lambda})$
as follows,
\begin{equation}
  \check{\Xi}^Q_{\mathcal{D}'}(x;\boldsymbol{\lambda})
  \stackrel{\text{def}}{=}\frac{\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda})}
  {\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(-M;\boldsymbol{\lambda})}
  \ \ \bigl(\Rightarrow\check{\Xi}^Q_{\mathcal{D}'}(-M;\boldsymbol{\lambda})=1\bigr).
\end{equation}
Next let us calculate the potential functions \eqref{BD}--\eqref{DD}.
Recall \eqref{BDKA}--\eqref{DDKA}, $\mu=0$ now and \eqref{cPKAD0}.
{}From \eqref{cXiDKAmonic2} with the replacements
$\mathcal{D}\to\mathcal{D}^{[-1]}$ and
$\boldsymbol{\lambda}\to\boldsymbol{\lambda}+\boldsymbol{\delta}$, we have
\begin{align}
  &\quad\check{\Xi}^{\text{KA}\,\text{monic}}_{\mathcal{D}^{[-1]}}
  (x;\boldsymbol{\lambda}+\boldsymbol{\delta})\nonumber\\
  &=(\text{const})\times\prod_{j=1}^M\Lambda(x+j-1;\boldsymbol{\lambda}+\boldsymbol{\delta})
  \cdot\check{\Xi}^{\text{KA}\,\text{monic}}_{(\mathcal{D}^{[-1]})^{[-N]}}
  \bigl(x-N;\boldsymbol{\lambda}+\boldsymbol{\delta}+N(\boldsymbol{\delta}+\boldsymbol{\bar{\delta}})\bigr)\nonumber\\
  &=(\text{const})\times\prod_{j=1}^M\Lambda(x+j-1;\boldsymbol{\lambda}+\boldsymbol{\delta})
  \cdot\check{\Xi}^{\text{KA}\,\text{monic}}_{\mathcal{D}^{[-N-1]}}
  (x-N;\boldsymbol{\lambda}'-\boldsymbol{\bar{\delta}})\nonumber\\
  &=(\text{const})\times\prod_{j=1}^M\Lambda(x+j-1;\boldsymbol{\lambda}+\boldsymbol{\delta})
  \cdot\rho^{-N\ell_{\mathcal{D}'}}\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}
  (x;\boldsymbol{\lambda}+\boldsymbol{\delta}).
  \label{cXiD[-1]KAmonic}
\end{align}
By using \eqref{cXiDKAmonicnonmonic}--\eqref{cPDnKAmonicnonmonic},
\eqref{cPKAD0}, \eqref{cXiDKAmonic2}, \eqref{cXiD[-1]KAmonic} and
\eqref{B(x+M),D(x+M)}, \eqref{BDKA}--\eqref{DDKA} become
\begin{align}
  B_{\mathcal{D}}(x;\boldsymbol{\lambda})&=B(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})
  \frac{\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda})}
  {\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x+1;\boldsymbol{\lambda})}
  \frac{\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}
  (x+1;\boldsymbol{\lambda}+\boldsymbol{\delta})}
  {\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}
  (x;\boldsymbol{\lambda}+\boldsymbol{\delta})},
  \label{BDnew}\\
  D_{\mathcal{D}}(x;\boldsymbol{\lambda})&=D(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})
  \frac{\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x+1;\boldsymbol{\lambda})}
  {\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda})}
  \frac{\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}
  (x-1;\boldsymbol{\lambda}+\boldsymbol{\delta})}
  {\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}
  (x;\boldsymbol{\lambda}+\boldsymbol{\delta})}.
  \label{DDnew}
\end{align}
These potential functions resemble those of the case-(1) multi-indexed
polynomials in \cite{os26}.
The deformed Hamiltonian \eqref{HopD} is similarity transformed to
\begin{align}
  &\quad\widetilde{\mathcal{H}}^{Q\,\text{op}}_{\mathcal{D}'}(\boldsymbol{\lambda})
  \stackrel{\text{def}}{=}\psi^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda})^{-1}\circ
  \mathcal{H}^{\text{op}}_{\mathcal{D}}(\boldsymbol{\lambda})\circ
  \psi^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda})\nonumber\\
  &=B(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})
  \frac{\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda})}
  {\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x+1;\boldsymbol{\lambda})}\Bigl(
  \frac{\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}
  (x+1;\boldsymbol{\lambda}+\boldsymbol{\delta})}
  {\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}
  (x;\boldsymbol{\lambda}+\boldsymbol{\delta})}
  -e^{\hat{\partial}}\Bigr)\nonumber\\
  &\quad+D(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})
  \frac{\check{\Xi}^{Q,\,\text{monic}}_{\mathcal{D}'}(x+1;\boldsymbol{\lambda})}
  {\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda})}\Bigl(
  \frac{\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}
  (x-1;\boldsymbol{\lambda}+\boldsymbol{\delta})}
  {\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}
  (x;\boldsymbol{\lambda}+\boldsymbol{\delta})}
  -e^{-\hat{\partial}}\Bigr).
  \label{HtDQop}
\end{align}
{}From the result in \S\,\ref{sec:Dar}, the multi-indexed polynomials
$\check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})$ satisfy the
difference equation,
\begin{equation}
  \widetilde{\mathcal{H}}^{Q\,\text{op}}_{\mathcal{D}'}(\boldsymbol{\lambda})
  \check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})
  =\mathcal{E}_n(\boldsymbol{\lambda})
  \check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda}),
  \label{HtDQopcQDn}
\end{equation}
for $x\in\mathbb{R}$ and $n\in\mathbb{Z}_{\geq0}$.
The results \eqref{HtDQopcQDn} for $n\in\mathcal{D}$ are obtained from those
for $n\in\mathbb{Z}_{\geq0}\backslash\mathcal{D}$ by taking the
``$n\to d_i$ limit.''
Exactly speaking, the results \eqref{HtDQopcQDn} for $n\in\mathcal{D}$ may not
be proven, because there is ambiguity in the ``$n\to d_i$ limit'' and the
prescription given in Appendix\,\ref{app:n->di} is one way of taking that limit.
However, we can verify \eqref{HtDQopcQDn} for $n\in\mathbb{Z}_{\geq0}$ by
direct calculation (using Mathematica) for small $N$, $M$, $d_j$ and $n$.
We remark that this difference equation holds for any $\boldsymbol{\lambda}$.
In the expressions \eqref{psiDQ} and \eqref{BDnew}--\eqref{DDnew}, they contain
$\phi_0$, $B$ and $D$ with the coordinate $x+M$ and the parameter
$\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}}$, which means that the parameter $N$ is $N+M$.
Therefore the matrix $\widetilde{\mathcal{H}}^Q_{\mathcal{D}'}(\boldsymbol{\lambda})$ is
of order $N+M+1$,
\begin{align}
  \widetilde{\mathcal{H}}^Q_{\mathcal{D}'}(\boldsymbol{\lambda})&=
  \bigl(\widetilde{\mathcal{H}}^Q_{\mathcal{D}'\,x,y}
  (\boldsymbol{\lambda})\bigr)_{x,y=-M,-M+1,\ldots,N},\nonumber\\
  \widetilde{\mathcal{H}}^Q_{\mathcal{D}'}(\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \widetilde{\mathcal{H}}^{Q\,\text{op}}_{\mathcal{D}'}(\boldsymbol{\lambda})
  \bigm|_{\hat{\partial}\to\partial}.
  \label{HtDQdef}
\end{align}
Since $B(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})$ and
$D(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})$ satisfy the boundary conditions
$B(N+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})=0$ at $x=N$ and
$D(0;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})=0$ at $x=-M$, the difference equation
\eqref{HtDQopcQDn} implies the following.
The matrix eigenvalue problem for $\widetilde{\mathcal{H}}^Q_{\mathcal{D}'}$ is
solved by the multi-indexed polynomials
$\check{Q}^{\text{monic}}_{\mathcal{D}',n}(x)$,
\begin{equation}
  \widetilde{\mathcal{H}}^Q_{\mathcal{D}'}(\boldsymbol{\lambda})
  \check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})
  =\mathcal{E}_n(\boldsymbol{\lambda})
  \check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})
  \ \ (n\in\{0,1,\ldots,N\}\cup\mathcal{D}).
\end{equation}
For $n>N$ and $n\not\in\mathcal{D}$, this equation also holds as $0=0$,
see \eqref{cQmonic_n=N+1+m} and its comment.

The orthogonality relations for $\check{Q}^{\text{monic}}_{\mathcal{D}',n}(x)$
can be read from \eqref{sumphigen}.
{}From \eqref{phiDngenQ}--\eqref{psiDQ},
$\phi_{\mathcal{D}\,n}^{\text{gen}}(x;\boldsymbol{\lambda})$ contains
$\phi_0(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})$, which does not vanish
at $x=-M,-M+1,\ldots,N$ for $x\in\mathbb{Z}$.
So the sum in \eqref{sumphigen} is reduced to $\sum_{x=-M}^N$.
We obtain the orthogonality relations,
\begin{align}
  &\quad\sum_{x=-M}^N\frac{\phi_0(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})^2}
  {\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda})
  \check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x+1;\boldsymbol{\lambda})}
  \check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})
  \check{Q}^{\text{monic}}_{\mathcal{D}',m}(x;\boldsymbol{\lambda})\nonumber\\
  &=\frac{\delta_{nm}}{d^{Q\,\text{monic}}_{\mathcal{D}',n}(\boldsymbol{\lambda})^2}
  \ \ (n,m\in\{0,1,\ldots,N\}\cup\mathcal{D}).
  \label{orthoQ}
\end{align}
Here the normalization constants
$d^{Q\,\text{monic}}_{\mathcal{D}',n}(\boldsymbol{\lambda})^2$ for
$n\in\{0,1,\ldots,N\}$ are obtained from \eqref{sumphigen} and
\eqref{phiDngenQ}--\eqref{psiDQ} as follows,
\begin{align}
  d^{Q\,\text{monic}}_{\mathcal{D}',n}(\boldsymbol{\lambda})^2
  &=d^{\text{monic}}_n(\boldsymbol{\lambda})^2\,
  \kappa^{\binom{M}{2}}\rho^{-2(2N+1)\binom{M}{2}}
  \prod_{j=1}^M\frac{B(0;\boldsymbol{\lambda}+(j-1)\boldsymbol{\delta})}
  {\mathcal{E}_n(\boldsymbol{\lambda})-\mathcal{E}_{d_j}(\boldsymbol{\lambda})}\nonumber\\
  &\quad\times
  \frac{c^{\eta}_{\mathcal{D},n}(\boldsymbol{\lambda})^2}
  {c^{\eta}_{\mathcal{D}'}(\boldsymbol{\lambda}')^2}
  \frac{\Lambda_M(0;\boldsymbol{\lambda})}
  {\phi_0(M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})^2}
  \frac{\varphi_M(N+1;\boldsymbol{\lambda})^2}{\varphi_M(0;\boldsymbol{\lambda}')^2}.
  \label{dQmonicD'n}
\end{align}
For $n=d_i\in\mathcal{D}$, we take the ``$n\to d_i$ limit'' of
\eqref{dQmonicD'n}.
By using \eqref{cetaDn=-A}, \eqref{En-Edi=-A} and \eqref{dnmonic=-A},
we obtain $d^{Q\,\text{monic}}_{\mathcal{D}',n}(\boldsymbol{\lambda})^2$ with
$n=d_i\in\mathcal{D}$,
\begin{align}
  d^{Q\,\text{monic}}_{\mathcal{D}',n}(\boldsymbol{\lambda})^2
  &=d^{\,\prime\,\text{monic}}_n(\boldsymbol{\lambda})^2\,
  \kappa^{\binom{M}{2}}\rho^{-2(2N+1)\binom{M}{2}}
  \frac{\prod_{j=1}^MB(0;\boldsymbol{\lambda}+(j-1)\boldsymbol{\delta})}
  {\prod_{\substack{j=1\\ j\neq i}}^M
  (\mathcal{E}_n(\boldsymbol{\lambda})-\mathcal{E}_{d_j}(\boldsymbol{\lambda}))}\nonumber\\
  &\quad\times
  \frac{c^{\prime\,\eta}_{\mathcal{D},n}(\boldsymbol{\lambda})^2}
  {c^{\eta}_{\mathcal{D}'}(\boldsymbol{\lambda}')^2}
  \frac{\Lambda_M(0;\boldsymbol{\lambda})}
  {\phi_0(M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})^2}
  \frac{\varphi_M(N+1;\boldsymbol{\lambda})^2}{\varphi_M(0;\boldsymbol{\lambda}')^2}\nonumber\\
  &\quad\times\rho^{2n}\kappa^{-n}\times\left\{
  \begin{array}{ll}
  1&:\text{(\romannumeral1)}',\,\text{(\romannumeral3)}',\,
  \text{(\romannumeral4)}'\\
  (2n+\tilde{d}\,)^{-1}&:\text{(\romannumeral2)}'\\[2pt]
  (1-\tilde{d}q^{2n})^{-1}&:\text{(\romannumeral5)}'
  \end{array}\right..
  \label{dQminicD'di}
\end{align}
In order to call the relations \eqref{orthoQ} truly orthogonality relations,
the positivity of the weight factor is necessary.
We will discuss this problem in the next subsection.
If the positivity is satisfied, the Hamiltonian of the deformed system
\begin{align}
  \mathcal{H}^Q_{\mathcal{D}'}(\boldsymbol{\lambda})&=
  \bigl(\mathcal{H}^Q_{\mathcal{D}'\,x,y}
  (\boldsymbol{\lambda})\bigr)_{x,y=-M,-M+1,\ldots,N},\nonumber\\
  \mathcal{H}^Q_{\mathcal{D}'}(\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \psi^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda})\circ
  \widetilde{\mathcal{H}}^Q_{\mathcal{D}'}(\boldsymbol{\lambda})\circ
  \psi^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda})^{-1}\nonumber\\
  &=-\sqrt{B_{\mathcal{D}}(x;\boldsymbol{\lambda})D_{\mathcal{D}}(x+1;\boldsymbol{\lambda})}
  \,e^{\partial}
  -\sqrt{B_{\mathcal{D}}(x-1;\boldsymbol{\lambda})D_{\mathcal{D}}(x;\boldsymbol{\lambda})}
  \,e^{-\partial}\nonumber\\
  &\quad+B_{\mathcal{D}}(x;\boldsymbol{\lambda})+D_{\mathcal{D}}(x;\boldsymbol{\lambda})
\end{align}
is a real symmetric matrix, and the matrix eigenvalue problem for
$\mathcal{H}^Q_{\mathcal{D}'}$ is solved by
$\phi^{Q\,\text{monic}}_{\mathcal{D}'\,n}(x)$ \eqref{phiDnQ},
\begin{equation}
  \mathcal{H}^Q_{\mathcal{D}'}(\boldsymbol{\lambda})
  \phi^{Q\,\text{monic}}_{\mathcal{D}'\,n}(x;\boldsymbol{\lambda})
  =\mathcal{E}_n(\boldsymbol{\lambda})
  \phi^{Q\,\text{monic}}_{\mathcal{D}'\,n}(x;\boldsymbol{\lambda})
  \ \ (n\in\{0,1,\ldots,N\}\cup\mathcal{D}).
\end{equation}
However, we remark that the relations \eqref{orthoQ} hold for any $\boldsymbol{\lambda}$
and $\mathcal{D}$, because the expressions are algebraic and the sum in
\eqref{orthoQ} is a finite sum.
We can verify \eqref{orthoQ} by direct calculation (using Mathematica) for
small $N$, $M$, $d_j$ and $n$.

We comment on the special case $\mathcal{D}=\{N+1,N+2,\ldots,N+M\}$ studied in
\cite{os40}. In this case we have $\mathcal{D}'=\{0,1,\ldots,M-1\}$ and the
denominator polynomial $\check{\Xi}^{\text{KA}}_{\mathcal{D}'}(x)$ becomes
a constant, $\check{\Xi}^{\text{KA}}_{\mathcal{D}'}(x)=1$ \cite{os22},
namely $\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x)=1$.
So $\widetilde{\mathcal{H}}^{Q\,\text{op}}_{\mathcal{D}'}(\boldsymbol{\lambda})$
\eqref{HtDQop} becomes
\begin{equation}
  \widetilde{\mathcal{H}}^{Q\,\text{op}}_{\mathcal{D}'}(\boldsymbol{\lambda})
  =B(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})(1-e^{\hat{\partial}})
  +D(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})(1-e^{-\hat{\partial}}).
\end{equation}
This and \eqref{HtDQopcQDn} imply
\begin{equation}
  \check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})
  =\rho^{-nM}\check{P}^{\text{monic}}_n(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})
  \ \ (n=0,1,\ldots,N+M),
  \label{cQD'n=cPn}
\end{equation}
which is shown by \eqref{cQD'nmonic} and
\begin{equation}
  \eta(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})
  =\rho^M\bigl(\eta(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})
  -\eta(-M;\boldsymbol{\lambda}+M\boldsymbol{\delta})\bigr).
\end{equation}
We remark that verification of this relation \eqref{cQD'n=cPn} based on
\eqref{cQmonic} and \eqref{cQmonic2} is non-trivial.
In \cite{os40}, where
$\check{Q}^{\text{monic}}_{\mathcal{D}',n}(x;\boldsymbol{\lambda})$ is not given,
the fact that the eigenpolynomials are given by
$\check{P}^{\text{monic}}_n(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})$
is shown based on direct calculations of Casoratians.

\subsection{Positivity of the weight}
\label{sec:positivity}

Let us consider the condition for the positivity of the weight factor in
\eqref{orthoQ}.
The positivity of $\phi_0(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})^2$ can be easily
achieved by choosing parameters $\boldsymbol{\lambda}$ so that
$B(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})$ and
$D(x+M;\boldsymbol{\lambda}-M\boldsymbol{\bar{\delta}})$ are positive (except for the
boundaries). From the positivity conditions for $B(x;\boldsymbol{\lambda})$ and
$D(x;\boldsymbol{\lambda})$ given in Appendix\,\ref{app:poly}, we have
\begin{align}
  \text{H}:&\ \ a,b>0\nonumber\\[-5pt]
  \text{K}:&\ \ 0<p<1\nonumber\\[-5pt]
  \text{R}:&\ \ b>N+d,\ \ d>M,\ \ 0<c<1+d-M,\ \ (d\neq M+1)\nonumber\\[-5pt]
  \text{dH}:&\ \ a>0,\ \ b>M,\ \ (a+b\neq M+1,M+2)\nonumber\\[-5pt]
  \text{dq$q$K}:&\ \ p>q^{-N}\nonumber\\[-5pt]
  \text{$q$H}:&\ \ 0<a,b<1\nonumber\\[-5pt]
  \text{$q$K}:&\ \ p>0
  \label{rangeM}\\[-5pt]
  \text{q$q$K}:&\ \ p>q^{-N-M}\nonumber\\[-5pt]
  \text{a$q$K}:&\ \ 0<p<q^{-1}\nonumber\\[-5pt]
  \text{$q$R}:&\ \ 0<b<dq^N,\ \ d<q^M,\ \ dq^{1-M}<c<1,
  \ \ (d\neq q^{M+1})\nonumber\\[-5pt]
  \text{d$q$H}:&\ \ 0<a<1,\ \ 0<bq^{-M}<1,\ \ (ab\neq q^{M+1},q^{M+2})\nonumber\\[-5pt]
  \text{d$q$K}:&\ \ p>0.
  \nonumber
\end{align}
We divide these into two classes,
\begin{equation}
  \text{(a)}:\text{H, K, dq$q$K, $q$H, $q$K, a$q$K, d$q$K},\quad
  \text{(b)}:\text{R, dH, q$q$K, $q$R, d$q$H}.
  \label{classab}
\end{equation}
The parameter ranges \eqref{rangeM} are independent of $M$ for class (a) and
dependent on $M$ for class (b).
The positivity of
$\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x;\boldsymbol{\lambda})
\check{\Xi}^{Q\,\text{monic}}_{\mathcal{D}'}(x+1;\boldsymbol{\lambda})$
(for $x=-M,-M+1,\ldots,N$) is highly non-trivial.
We try this problem by numerical calculations.
In the rest of this subsection we assume $d_1<d_2<\cdots<d_M$ and
the parameters $\boldsymbol{\lambda}$ satisfying \eqref{rangeM}.

Based on numerical calculations, we have observed the followings.
For class (a), the condition for the positivity is given by
\begin{align}
  &d_j-d_{j-1}\equiv 1\ \ (\text{mod }2)
  \ \ (j=1,2,\ldots,M\ ;\ d_0=N),\nonumber\\
  &\Bigl(\Leftrightarrow
  m_{2k-1}:\text{even}\ \bigl(k=1,2,\ldots,[\tfrac{M+1}{2}]\bigr),
  \ m_{2k}:\text{odd}\ \bigl(k=1,2,\ldots,[\tfrac{M}{2}]\bigr)\Bigr).
  \label{djsaodd}
\end{align}
For class (b), the situation is more complicated.
For $\mathcal{D}$ satisfying \eqref{djsaodd}, the positivity is satisfied by
further restricting the range of parameters \eqref{rangeM}.
We have not yet found a definite range, but the following ranges seem to work
well: q$q$K: large $p$, dH: large $b$, d$q$H: small $b$, R: large $d$,
$q$R: small $d$.
Even for $\mathcal{D}$ not satisfying \eqref{djsaodd}, the positivity
may be satisfied.
For example, in the case of $M=1$, if the parameters are well chosen, there
are no restrictions on $d_1$. For general $M$, we have not yet found the
conditions for $\mathcal{D}$.

\section{Summary and Comments}
\label{sec:summary}

The Hamiltonian of a finite type rdQM is a real symmetric matrix of order $N+1$
(the coordinate $x\in\{0,1,\ldots,N\}$) and the Schr\"odinger equation is a
matrix eigenvalue problem, whose eigenvectors are $\phi_n(x)$ with
$n\in\{0,1,\ldots,N\}$. This eigenvector $\phi_n(x)$ is extended to the
function $\phi_n(x)$ ($x\in\mathbb{R}$, $n\in\mathbb{Z}_{\geq 0}$) (for $n>N$,
$\phi_n(x)$ is replaced with $\phi^{\text{monic}}_n(x)$), which satisfy the
difference equation.
Based on such rdQM systems described by the twelve orthogonal polynomials,
we have considered their deformation by the multi-step Darboux transformations
with seed solutions $\phi_n(x)$.
The seed solution $\phi_n(x)$ with $n>N$ corresponds to the overshoot
eigenfunction $\tilde{\phi}^{\text{os}}_n(x)$ in oQM systems with a finite
number of eigenstates \cite{os28}, and it is a ``zero norm'' eigenvector because
it vanishes at $x\in\{0,1,\ldots,N\}$ \cite{os40}.
For seed solutions $\phi_n(x)$ with $n\leq N$
($n\in\mathcal{D}=\{d_1,d_2,\ldots,d_M\}$), the Darboux transformations are
state-deleting and the deformed Hamiltonian has order $N-M+1$
($x\in\{0,1,\ldots,N-M\}$) and the eigenvectors are described by the
Krein-Alder type multi-indexed orthogonal polynomials
$\check{P}^{\text{KA}}_{\mathcal{D},n}(x)$
($n\in\{0,1,\ldots,N\}\backslash\mathcal{D}$) \cite{os22}.
For seed solutions $\phi_n(x)$ with $n>N$ ($n\in\mathcal{D}$), the Darboux
transformations are state-adding and the deformed Hamiltonian has order $N+M+1$
($x\in\{-M,-M+1,\ldots,N\}$) and the eigenvectors are described by new
multi-indexed orthogonal polynomials $\check{Q}_{\mathcal{D}',n}(x)$
($n\in\{0,1,\ldots,N\}\cup\mathcal{D}$).
Explicit forms of the difference equation and the orthogonality relations for
$\check{P}^{\text{KA}}_{\mathcal{D},n}(x)$ are new results and given in
\S\,\ref{sec:KAmiop}.
New multi-indexed orthogonal polynomials $\check{Q}_{\mathcal{D}',n}(x)$ are
main results of this paper. Their definitions, the difference equation and the
orthogonality relations are given in \S\,\ref{sec:newmiop}.
The energy eigenvalues of $M$ added states are $\mathcal{E}_n$
($n\in\mathcal{D}$). This is in contrast to the case of state-adding Darboux
transformations with the pseudo virtual states as seed solutions
\cite{casoidrdqm}, in which the energy eigenvalues of $M$ added states are
$\mathcal{E}_{-n-1}$ ($n\in\mathcal{D}$).
The positivity condition of the weight factor is discussed in
\S\,\ref{sec:positivity}. Partial results are obtained, but they are still
unsatisfactory. It is an important problem to clarify the positivity condition.
In \S\,\ref{sec:newmiop} the condition $M\leq N$ is assumed, and it is an
interesting problem to consider the case $M>N$.

We have shown the state-adding property of the Darboux transformations with
seed solutions $\phi_n(x)$ ($n>N$) by using the embedding of the finite system
into the infinite system.
In Appendix B of \cite{casoidrdqm}, where the state-adding Darboux
transformations with the pseudo virtual states
$\tilde{\phi}^{\text{pv}}_{\text{v}}(x)$ as seed solutions are discussed,
the state-adding property was shown without using such a trick.
It was shown that each step of the Darboux transformation increases the size of
the matrix (Hamiltonian) by one.
We think that such a proof is possible in the present case as well.
We hope that we will be able to report on this subject in detail elsewhere.

The multi-indexed orthogonal polynomials do not satisfy the three term
recurrence relations, which characterize the ordinary orthogonal polynomials
\cite{ismail}. They satisfy the recurrence relations with more terms
(\cite{d14,d17,rrmiop5,mtv22} for rdQM).
It is an interesting problem to study the recurrence relations for the new
multi-indexed orthogonal polynomials $\check{Q}_{\mathcal{D}',n}(x)$.

\section*{Acknowledgements}

I thank Ryu Sasaki for discussions in the early stage of this work.
This work is supported by JSPS KAKENHI Grant Number JP19K03667.

\bigskip
