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

The exceptional and multi-indexed orthogonal polynomials
$\{P_{\mathcal{D},n}(\eta)|n\in\mathbb{Z}_{\geq 0}\}$
\cite{gkm08_1}--\cite{os35}
are new type of orthogonal polynomials which
form a complete set of orthogonal basis in an appropriate Hilbert space
in spite of missing degrees.
Instead of the three term recurrence relations which characterize the
ordinary orthogonal polynomials
\cite{szego_1,szego_2,szego_3}, they satisfy the recurrence
relations with more terms \cite{stz10}--\cite{gkkm15}, and the constraints
by Bochner's theorem and its generalizations
\cite{bochner_1,bochner_2,szego_1,szego_2,szego_3} are avoided.
They satisfy second order differential or difference equations and
quantum mechanical formulation has played an important role in constructing
these new orthogonal polynomials.
The quantum mechanical systems described by the classical orthogonal
polynomials can be iso-spectrally (or almost iso-spectrally) deformed by
applying the multi-step Darboux transformations
\cite{darb_1}--\cite{os22}
and the main part of the eigenfunctions of the deformed systems are
multi-indexed orthogonal polynomials.
We distinguish the following two cases;
the set of missing degrees $\mathcal{I}=\mathbb{Z}_{\geq 0}\backslash
\{\text{deg}\,P_{\mathcal{D},n}(\eta)|n\in\mathbb{Z}_{\geq 0}\}$ is
case (1): $\mathcal{I}=\{0,1,\ldots,\ell-1\}$, or
case (2): $\mathcal{I}\neq\{0,1,\ldots,\ell-1\}$, where $\ell$ is a positive
integer. The situation of case (1) is called stable in \cite{gkm11}.
Case (1) is obtained by taking virtual states as the seed solutions of
the Darboux transformation \cite{os25,os27,os26,os35} and case (2) is
obtained by employing eigenstates or pseudo virtual states
\cite{q12b}--\cite{os30}.
In the quantum mechanical formulation, the multi-indexed orthogonal
polynomials appear as polynomials in the sinusoidal coordinate $\eta(x)$
\cite{os7_1,os7_2},
$P_{\mathcal{D},n}\bigl(\eta(x)\bigr)$, where $x$ is the
coordinate of the quantum system.

The eigenstates of the deformed system are expressed in terms of
determinants: the Wronskians for the ordinary quantum mechanics (in which
the Schr\"odinger equation is second order differential equation) and the
Casoratians for the discrete quantum mechanics (in which the Schr\"odinger
equations are second order difference equations).
In our previous paper \cite{os37}, we have presented simplified determinant
expressions of the case (1) multi-indexed Laguerre and Jacobi polynomials
in the ordinary quantum mechanics. We rewrite the Wronskians
by using the identities of the original Laguerre and Jacobi polynomials and
two simplified expressions A and B are presented. The identities used in
case A are essentially forward shift relations. For case B, the Schr\"odinger
equation is used and it is regarded as a combination of the forward and
backward shift relations, which are the consequences of the shape invariance.
Thus both A and B expressions are derived by using the properties of the
Wronskians and the shape invariance (see \S\,\ref{sec:summary} for details).
This idea is also applicable to the discrete quantum mechanics.
In this paper we consider the case (1) multi-indexed orthogonal polynomials
appearing in the discrete quantum mechanics:
the Meixner, little $q$-Laguerre, little $q$-Jacobi, Racah, $q$-Racah, Wilson
and Askey-Wilson types.
The first five appear in the discrete quantum mechanics with real shifts
(rdQM) and the last two belong to the discrete quantum mechanics with pure
imaginary shifts (idQM).
We rewrite the Casoratians by using their properties and the identities of
the original polynomials, which are the consequences of the shape invariance.
Various new determinant expressions of the multi-indexed orthogonal
polynomials are obtained.
Corresponding to the cases A and B in \cite{os37}, we present two typical
expressions explicitly.
The matrix elements of the Casoratians are functions of $x$ at various
points, $x+j-1$ for rdQM and $x+i(\frac{n+1}{2}-j)\gamma$ for idQM.
In some of new determinant expressions, the matrix elements are functions
of the same point $x$. Therefore the multi-indexed orthogonal polynomials
are expressed in terms of $\eta$ only, namely without explicit $x$-dependence.
For the Racah and $q$-Racah cases, the corresponding matrix elements for
the multi-indexed polynomials do have explicit $x$-dependence because of
the parameter ($\boldsymbol{\lambda}$) dependence of the sinusoidal coordinates,
$\eta(x;\boldsymbol{\lambda})$.
These new equivalent expressions show explicitly the constituents of the
multi-indexed orthogonal polynomials and they are helpful for deeper
understanding.

This paper is organized as follows.
We derive the new equivalent expressions for the multi-indexed Meixner,
little $q$-Laguerre, little $q$-Jacobi, Racah and $q$-Racah polynomials in
section \ref{sec:rdQM}. Those for the multi-indexed Wilson and Askey-Wilson
polynomials are presented in section \ref{sec:idQM}.
The discussions in these two sections go parallelly.
In \S\,\ref{sec:rdQM_org} (\S\,\ref{sec:idQM_org}) the discrete quantum
mechanics with real (pure imaginary) shifts is recapitulated and the data
for the original systems are presented.
The virtual state vectors (wavefunctions) and related functions are
introduced in \S\,\ref{sec:rdQM_vs} (\S\,\ref{sec:idQM_vs}).
The original systems have shape invariance which entails the well-known
forward and backward shift relations of the eigenpolynomials $P_n$.
The forward and backward shift relations for virtual state polynomials
$\xi_{\text{v}}$ are presented in \S\,\ref{sec:rdQM_si}
(\S\,\ref{sec:idQM_si}).
In \S\,\ref{sec:rdQM_miop} (\S\,\ref{sec:idQM_miop}) the definitions of the
multi-indexed orthogonal polynomials given in \cite{os26,os35} (\cite{os27})
are recapitulated.
They are expressed in terms of Casoratians.
Sections \ref{sec:rdQM_exp} and \ref{sec:idQM_exp} are the main contents
of this paper.
The Casoratians are rewritten by using their properties and the identities
stemming from the shape invariance, and new determinant expressions of the
multi-indexed orthogonal polynomials are obtained.
We present two typical examples.
Section \ref{sec:summary} is for a summary and comments.

\section{Multi-indexed Orthogonal Polynomials in rdQM}
\label{sec:rdQM}

In this section we derive various equivalent determinant expressions for
the multi-indexed orthogonal polynomials in the framework of the discrete
quantum mechanics with real shifts.
They are the multi-indexed Meixner, little $q$-Laguerre, little $q$-Jacobi,
Racah and $q$-Racah polynomials.

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

The Hamiltonian of the discrete quantum mechanics with real shifts (rdQM)
\cite{os12,os34} $\mathcal{H}=(\mathcal{H}_{x,y})$ is a tri-diagonal real
symmetric (Jacobi) matrix and its rows and columns are indexed by integers
$x$ and $y$, which take values in $\{0,1,\ldots,N\}$ (finite) or
$\mathbb{Z}_{\geq 0}$ (semi-infinite),
\begin{equation}
  \mathcal{H}_{x,y}\stackrel{\text{def}}{=}
  -\sqrt{B(x)D(x+1)}\,\delta_{x+1,y}-\sqrt{B(x-1)D(x)}\,\delta_{x-1,y}
  +\bigl(B(x)+D(x)\bigr)\delta_{x,y},
  \label{rdQM:Hxy}
\end{equation}
where the potential functions $B(x)$ and $D(x)$ are real and positive
but vanish at the boundary, $D(0)=0$ and $B(N)=0$ for finite cases.
This Hamiltonian can be expressed in a factorized form:
\begin{align}
  &\mathcal{H}=\mathcal{A}^{\dagger}\mathcal{A},\qquad
  \mathcal{A}=(\mathcal{A}_{x,y}),
  \ \ \mathcal{A}^{\dagger}=((\mathcal{A}^{\dagger})_{x,y})
  =(\mathcal{A}_{y,x}),
  \label{factor}\\
  &\mathcal{A}_{x,y}\stackrel{\text{def}}{=}
  \sqrt{B(x)}\,\delta_{x,y}-\sqrt{D(x+1)}\,\delta_{x+1,y},\quad
  (\mathcal{A}^{\dagger})_{x,y}=
  \sqrt{B(x)}\,\delta_{x,y}-\sqrt{D(x)}\,\delta_{x-1,y}.
\end{align}
For simplicity in notation, we write $\mathcal{H}$, $\mathcal{A}$ and
$\mathcal{A}^{\dagger}$ as follows:
\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\\
  &\phantom{\mathcal{H}}=-\sqrt{B(x)}\,e^{\partial}\sqrt{D(x)}
  -\sqrt{D(x)}\,e^{-\partial}\sqrt{B(x)}+B(x)+D(x),
  \label{rdQM:H}\\
  &\mathcal{A}=\sqrt{B(x)}-e^{\partial}\sqrt{D(x)},\quad
  \mathcal{A}^{\dagger}=\sqrt{B(x)}-\sqrt{D(x)}\,e^{-\partial},
  \label{rdQM:A,Ad}
\end{align}
where matrices $e^{\pm\partial}$ are
\begin{equation}
  e^{\pm\partial}=((e^{\pm\partial})_{x,y}),\quad
  (e^{\pm\partial})_{x,y}\stackrel{\text{def}}{=}\delta_{x\pm 1,y},\quad
  (e^{\partial})^{\dagger}=e^{-\partial},
  \label{partdef}
\end{equation}
and we suppress the unit matrix $\boldsymbol{1}=(\delta_{x,y})$:
$\bigl(B(x)+D(x)\bigr)\boldsymbol{1}$ in \eqref{rdQM:H},
$\sqrt{B(x)}\,\boldsymbol{1}$ in \eqref{rdQM:A,Ad}.
(The notation $f(x)Ag(x)$, where $f(x)$ and $g(x)$ are functions of $x$ and
$A$ is a matrix $A=(A_{x,y})$, stands for a matrix whose $(x,y)$-element
is $f(x)A_{x,y}g(y)$.)
Note that the matrices $e^{\partial}$ and $e^{-\partial}$ are not inverse
to each other: $e^{\pm\partial}e^{\mp\partial}\neq\boldsymbol{1}$ for finite systems
and $e^{-\partial}e^{\partial}\neq\boldsymbol{1}$ for semi-infinite systems.
The Schr\"odinger equation is the eigenvalue problem for the hermitian
matrix $\mathcal{H}$,
\begin{equation}
  \mathcal{H}\phi_n(x)=\mathcal{E}_n\phi_n(x)\quad
  (n=0,1,2,\ldots),\quad
  0=\mathcal{E}_0<\mathcal{E}_1<\mathcal{E}_2<\cdots,
  \label{Scheq}
\end{equation}
($n=0,1,\ldots,N$ for a finite case).

We consider rdQM described by the Meixner (M), little $q$-Laguerre (l$q$L),
little $q$-Jacobi (l$q$J), Racah (R) and $q$-Racah ($q$R) polynomials.
The first three are semi-infinite systems and the last two are finite systems.
Eigenvectors have the following factorized 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).
\end{equation}
Here the ground state eigenvector $\phi_0(x)$, which is characterized by
$\mathcal{A}\phi_0(x)=0$, is chosen as
\begin{equation}
  \phi_0(x)=\prod_{y=0}^{x-1}\sqrt{\frac{B(y)}{D(y+1)}}\,>0.
  \label{rdQM:phi0}
\end{equation}
We use the convention: $\prod\limits_{k=n}^{n-1}*=1$, which means the
normalization $\phi_0(0)=1$.
The other function $P_n(\eta)$ is a degree $n$ polynomial in $\eta$, and
$\eta(x)$ is the sinusoidal coordinate satisfying the boundary
condition $\eta(0)=0$.
We adopt the universal normalization condition \cite{os12,os34} as
\begin{equation}
  P_n(0)=1\ \bigl(\Leftrightarrow \check{P}_n(0)=1\bigr).
  \label{Pzero}
\end{equation}
Various quantities depend on a set of parameters
$\boldsymbol{\lambda}=(\lambda_1,\lambda_2,\ldots)$ and their dependence is
expressed like,
$\mathcal{H}=\mathcal{H}(\boldsymbol{\lambda})$,
$\mathcal{A}=\mathcal{A}(\boldsymbol{\lambda})$,
$\mathcal{E}_n=\mathcal{E}_n(\boldsymbol{\lambda})$,
$B(x)=B(x;\boldsymbol{\lambda})$,
$\phi_n(x)=\phi_n(x;\boldsymbol{\lambda})$,
$\check{P}_n(x)=\check{P}_n(x;\boldsymbol{\lambda})
=P_n\bigl(\eta(x;\boldsymbol{\lambda});\boldsymbol{\lambda}\bigr)$, etc.
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)$.
The symbols $(a)_n$ and $(a;q)_n$ are ($q$-)shifted factorials \cite{kls}
and $[x]$ denotes the greatest integer not exceeding $x$.

Parameters of the systems are
\begin{align}
  \text{M}:\ \ &\boldsymbol{\lambda}=(\beta,c),\quad
  \boldsymbol{\delta}=(1,0),\quad \kappa=1,\quad \beta>0,\quad 0<c<1,\\
  \text{l$q$L}:\ \ &q^{\boldsymbol{\lambda}}=a,\quad
  \boldsymbol{\delta}=1,\quad \kappa=q^{-1},\quad 0<a<q^{-1},\\
  \text{l$q$J}:\ \ &q^{\boldsymbol{\lambda}}=(a,b),\quad
  \boldsymbol{\delta}=(1,1),\quad \kappa=q^{-1},\quad 0<a<q^{-1},\quad b<q^{-1},\\
  \text{R}:\ \ &\boldsymbol{\lambda\,}=(a,b,c,d),\quad \boldsymbol{\delta}=(1,1,1,1),
  \quad\kappa=1,\\
  \text{$q$R}:\ \ &q^{\boldsymbol{\lambda}}=(a,b,c,d),\quad \boldsymbol{\delta}=(1,1,1,1),
  \quad\kappa=q^{-1},
\end{align}
and we adopt the following choice of the parameter ranges for R and $q$R:
\begin{align}
  \text{R}:\ \ &a=-N,\quad 0<d<a+b,\quad 0<c<1+d,\\
  \text{$q$R}:\ \ &a=q^{-N},\quad 0<ab<d<1,\quad qd<c<1.
\end{align}
We list the fundamental data \cite{os12}.\\
$\circ$ {\bf Meixner:}\\
We rescale the overall normalization of the Hamiltonian in \cite{os12}:
$(\mathcal{H}\text{ in \cite{os12}})\times(1-c)\to\mathcal{H}$.
\begin{align}
  &B(x;\boldsymbol{\lambda})=c(x+\beta),\quad
  D(x)=x,
  \label{MBD}\\
  &\mathcal{E}_n(\boldsymbol{\lambda})=(1-c)n,\quad\eta(x)=x,\quad
  \varphi(x)=1,
  \label{Meta}\\
  &\check{P}_n(x;\boldsymbol{\lambda})
  ={}_2F_1\Bigl(
  \genfrac{}{}{0pt}{}{-n,\,-x}{\beta}\Bigm|1-c^{-1}\Bigr)
  =M_n(x;\beta,c),
  \label{MP}\\
  &\phi_0(x;\boldsymbol{\lambda})=\sqrt{\frac{(\beta)_x\,c^x}{x!}},
  \label{Mphi0}
\end{align}
where $M_n(x;\beta,c)$ is the Meixner polynomial \cite{kls}.\\
$\circ$ {\bf little $q$-Laguerre and little $q$-Jacobi:}
\begin{align}
  &B(x;\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  aq^{-x}&:\text{l$q$L}\\
  a(q^{-x}-bq)&:\text{l$q$J}\\
  \end{array}\right.,\quad
  D(x)=q^{-x}-1,\\[2pt]
  &\mathcal{E}_n(\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  q^{-n}-1&:\text{l$q$L}\\
  (q^{-n}-1)(1-abq^{n+1})&:\text{l$q$J}
  \end{array}\right.,\quad
  \eta(x)=1-q^x,\quad\varphi(x)=q^x,
  \label{lqJeta}\\[2pt]
  &\check{P}_n(x;\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  {\displaystyle
  {}_2\phi_0\Bigl(
  \genfrac{}{}{0pt}{}{q^{-n},\,q^{-x}}{-}\Bigm|q;\frac{q^x}{a}\Bigr)
  =(a^{-1}q^{-n}\,;q)_n\,{}_2\phi_1\Bigl(
  \genfrac{}{}{0pt}{}{q^{-n},\,0}{aq}\Bigm|q;q^{x+1}\Bigr)}&\\[4pt]
  {\displaystyle
  {}_3\phi_1\Bigl(
  \genfrac{}{}{0pt}{}{q^{-n},\,abq^{n+1},\,q^{-x}}{bq}\Bigm|
  q\,;\frac{q^x}{a}\Bigr)
  =\frac{(a^{-1}q^{-n};q)_n}{(bq;q)_n}\,
  {}_2\phi_1\Bigl(
  \genfrac{}{}{0pt}{}{q^{-n},\,abq^{n+1}}{aq}\Bigm|q\,;q^{x+1}\Bigr)}&
  \end{array}\right.\nonumber\\
  &\phantom{\check{P}_n(x;\boldsymbol{\lambda})}
  =\left\{
  \begin{array}{ll}
  (a^{-1}q^{-n}\,;q)_n\,p_n\bigl(1-\eta(x);a|q\bigr)&:\text{l$q$L}\\[4pt]
  (bq;q)_n^{-1}(a^{-1}q^{-n};q)_n\,
  p_n\bigl(1-\eta(x);a,b|q\bigr)&:\text{l$q$J}
  \end{array}\right.,\\[2pt]
  &\phi_0(x;\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  \sqrt{(q;q)_x^{-1}(aq)^x}&:\text{l$q$L}\\[4pt]
  \sqrt{(bq;q)_x(q;q)_x^{-1}(aq)^x}&:\text{l$q$J}
  \end{array}\right.,
\end{align}
where $p_n(y;a|q)$ and $p_n(y;a,b|q)$ are the little $q$-Laguerre and little
$q$-Jacobi polynomials in the conventional definition \cite{kls},
respectively.\\
$\circ$ {\bf Racah and $q$-Racah:}
\begin{align}
  &B(x;\boldsymbol{\lambda})=
  \left\{
  \begin{array}{ll}
  {\displaystyle
  -\frac{(x+a)(x+b)(x+c)(x+d)}{(2x+d)(2x+1+d)}}&:\text{R}\\[8pt]
  {\displaystyle-\frac{(1-aq^x)(1-bq^x)(1-cq^x)(1-dq^x)}
  {(1-dq^{2x})(1-dq^{2x+1})}}&:\text{$q$R}
  \end{array}\right.,
  \label{Bform}\\
  &D(x;\boldsymbol{\lambda})=
  \left\{
  \begin{array}{ll}
  {\displaystyle
  -\frac{(x+d-a)(x+d-b)(x+d-c)x}{(2x-1+d)(2x+d)}}&:\text{R}\\[8pt]
  {\displaystyle-\tilde{d}\,
  \frac{(1-a^{-1}dq^x)(1-b^{-1}dq^x)(1-c^{-1}dq^x)(1-q^x)}
  {(1-dq^{2x-1})(1-dq^{2x})}}&:\text{$q$R}
  \end{array}\right.,
  \label{Dform}\\
  &\mathcal{E}_n(\boldsymbol{\lambda})=
  \left\{
  \begin{array}{ll}
  n(n+\tilde{d})&:\text{R}\\[2pt]
  (q^{-n}-1)(1-\tilde{d}q^n)&:\text{$q$R}
  \end{array}\right.\!,\quad
  \eta(x;\boldsymbol{\lambda})=
  \left\{
  \begin{array}{ll}
  x(x+d)&:\text{R}\\[2pt]
  (q^{-x}-1)(1-dq^x)&:\text{$q$R}
  \end{array}\right.,
  \label{etadefs}\\
  &\varphi(x;\boldsymbol{\lambda})=
  \left\{
  \begin{array}{ll}
  {\displaystyle\frac{2x+d+1}{d+1}}&:\text{R}\\[6pt]
  {\displaystyle\frac{q^{-x}-dq^{x+1}}{1-dq}}&:\text{$q$R}
  \end{array}\right.,\quad
  \tilde{d}\stackrel{\text{def}}{=}
  \left\{
  \begin{array}{ll}
  a+b+c-d-1&:\text{R}\\[2pt]
  abcd^{-1}q^{-1}&:\text{$q$R}
  \end{array}\right.,
  \label{vphidef}\\
  &\check{P}_n(x;\boldsymbol{\lambda})
  =P_n\bigl(\eta(x;\boldsymbol{\lambda});\boldsymbol{\lambda}\bigr)
  =\left\{
  \begin{array}{ll}
  {\displaystyle
  {}_4F_3\Bigl(
  \genfrac{}{}{0pt}{}{-n,\,n+\tilde{d},\,-x,\,x+d}
  {a,\,b,\,c}\Bigm|1\Bigr)}\\[8pt]
  {\displaystyle
  {}_4\phi_3\Bigl(
  \genfrac{}{}{0pt}{}{q^{-n},\,\tilde{d}q^n,\,q^{-x},\,dq^x}
  {a,\,b,\,c}\Bigm|q\,;q\Bigr)}
  \end{array}\right.\nonumber\\
  &\phantom{\check{P}_n(x;\boldsymbol{\lambda})
  =P_n\bigl(\eta(x;\boldsymbol{\lambda});\boldsymbol{\lambda}\bigr)}
  =\left\{
  \begin{array}{ll}
  {\displaystyle
  R_n\bigl(\eta(x;\boldsymbol{\lambda});a-1,\tilde{d}-a,c-1,d-c\bigr)}
  &:\text{R}\\[2pt]
  {\displaystyle
  R_n\bigl(1+d+\eta(x;\boldsymbol{\lambda});
  aq^{-1},\tilde{d}a^{-1},cq^{-1},dc^{-1}|q\bigr)}&:\text{$q$R}
  \end{array}\right.,
  \label{Pn=R,qR}\\[4pt]
  &\phi_0(x;\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  {\displaystyle
  \sqrt{\frac{(a,b,c,d)_x}{(d-a+1,d-b+1,d-c+1,1)_x}\,
  \frac{2x+d}{d}}}&:\text{R}\\[12pt]
  {\displaystyle
  \sqrt{\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}}}&:\text{$q$R}
  \end{array}\right.,
\end{align}
where $R_n\bigl(x(x+\gamma+\delta+1);\alpha,\beta,\gamma,\delta\bigr)$ and
$R_n(q^{-x}+\gamma\delta q^{x+1};\alpha,\beta,\gamma,\delta|q)$ are the Racah
and $q$-Racah polynomials \cite{kls}, respectively.
It should be emphasized that for these two polynomials, the sinusoidal
coordinates $\eta(x)$ \eqref{etadefs} and the auxiliary function
$\varphi(x)$ \eqref{vphidef} depend on the parameters $\eta(x;\boldsymbol{\lambda})$,
$\varphi(x;\boldsymbol{\lambda})$ in distinction with the other polynomials in this
section.

Note that the following relations are valid for all the systems discussed
in this section,
\begin{align}
  &\varphi(x;\boldsymbol{\lambda})=
  \frac{\eta(x+1;\boldsymbol{\lambda})-\eta(x;\boldsymbol{\lambda})}{\eta(1;\boldsymbol{\lambda})},\\
  &\varphi(x;\boldsymbol{\lambda})=\sqrt{\frac{B(0;\boldsymbol{\lambda})}{B(x;\boldsymbol{\lambda})}}
  \,\frac{\phi_0(x;\boldsymbol{\lambda}+\boldsymbol{\delta})}{\phi_0(x;\boldsymbol{\lambda})},\quad
  \varphi(x;\boldsymbol{\lambda})=\sqrt{\frac{B(0;\boldsymbol{\lambda})}{D(x+1;\boldsymbol{\lambda})}}
  \,\frac{\phi_0(x;\boldsymbol{\lambda}+\boldsymbol{\delta})}{\phi_0(x+1;\boldsymbol{\lambda})},
  \label{OS12(4.12,13)}\\
  &\frac{B(x;\boldsymbol{\lambda}+\boldsymbol{\delta})}{B(x+1;\boldsymbol{\lambda})}
  =\kappa^{-1}\frac{\varphi(x+1;\boldsymbol{\lambda})}{\varphi(x;\boldsymbol{\lambda})},\quad
  \frac{D(x;\boldsymbol{\lambda}+\boldsymbol{\delta})}{D(x;\boldsymbol{\lambda})}
  =\kappa^{-1}\frac{\varphi(x-1;\boldsymbol{\lambda})}{\varphi(x;\boldsymbol{\lambda})}.
  \label{OS12(5.81etc)}
\end{align}

\subsection{Virtual states}
\label{sec:rdQM_vs}

The virtual state vectors are obtained by using the discrete symmetries
of the Hamiltonian \cite{os26,os35}.
We define the twist operation $\mathfrak{t}$,
\begin{align}
  &\mathfrak{t}(\boldsymbol{\lambda})\stackrel{\text{def}}{=}\left\{
  \begin{array}{ll}
  (\lambda_1,\lambda_2^{-1})&:\text{M}\\
  -\lambda_1&:\text{l$q$L}\\
  (-\lambda_1,\lambda_2)&:\text{l$q$J}\\
  (\lambda_4-\lambda_1+1,\lambda_4-\lambda_2+1,\lambda_3,\lambda_4)
  &:\text{R,\,$q$R}
  \end{array}\right.,\nonumber\\[2pt]
  &\text{namely}\ \ \left\{
  \begin{array}{ll}
  \mathfrak{t}(\boldsymbol{\lambda})=(\beta,c^{-1})&:\text{M}\\
  q^{\mathfrak{t}(\boldsymbol{\lambda})}=a^{-1}&:\text{l$q$L}\\
  q^{\mathfrak{t}(\boldsymbol{\lambda})}=(a^{-1},b)&:\text{l$q$J}\\
  \mathfrak{t}(\boldsymbol{\lambda})=(d-a+1,d-b+1,c,d)&:\text{R}\\
  q^{\mathfrak{t}(\boldsymbol{\lambda})}=(a^{-1}dq,b^{-1}dq,c,d)&:\text{$q$R}
  \end{array}\right.,
  \label{rdQM:twist}
\end{align}
which is an involution $\mathfrak{t}^2=\text{id}$ and satisfies
\begin{equation}
  \mathfrak{t}(\boldsymbol{\lambda})+u\boldsymbol{\delta}
  =\mathfrak{t}(\boldsymbol{\lambda}+u\tilde{\boldsymbol{\delta}})
  \ \ (\forall u\in\mathbb{R}),\quad
  \tilde{\boldsymbol{\delta}}\stackrel{\text{def}}{=}\left\{
  \begin{array}{ll}
  (1,0)&:\text{M}\\
  -1&:\text{l$q$L}\\
  (-1,1)&:\text{l$q$J}\\
  (0,0,1,1)&:\text{R,\,$q$R}
  \end{array}\right..
\end{equation}
Note that this twist operation does not affect $\eta(x;\boldsymbol{\lambda})$ and
$\varphi(x;\boldsymbol{\lambda})$,
\begin{equation}
  \eta\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})\bigr)=\eta(x;\boldsymbol{\lambda}),\quad
  \varphi\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})\bigr)=\varphi(x;\boldsymbol{\lambda}).
\end{equation}
We define two functions $B'(x;\boldsymbol{\lambda})$ and $D'(x;\boldsymbol{\lambda})$,
\begin{equation}
  B'(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=} B\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})\bigr),\quad
  D'(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=} D\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})\bigr).
  \label{B'D'}
\end{equation}
By using the twist operation, the virtual state vectors
$\tilde{\phi}_{\text{v}}(x;\boldsymbol{\lambda})$ are introduced:
\begin{align}
  \tilde{\phi}_{\text{v}}(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \phi_{\text{v}}\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})\bigr)
  =\tilde{\phi}_0(x;\boldsymbol{\lambda})
  \check{\xi}_{\text{v}}(x;\boldsymbol{\lambda}),\quad
  \tilde{\phi}_0(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \phi_0\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})\bigr),\nonumber\\
  \check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \xi_{\text{v}}\bigl(\eta(x;\boldsymbol{\lambda});\boldsymbol{\lambda}\bigr)\stackrel{\text{def}}{=}
  \check{P}_{\text{v}}\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})\bigr)
  =P_{\text{v}}\bigl(\eta(x;\boldsymbol{\lambda});
  \mathfrak{t}(\boldsymbol{\lambda})\bigr).
  \label{rdQM:phitv}
\end{align}
The virtual state polynomials $\xi_{\text{v}}(\eta;\boldsymbol{\lambda})$
are polynomials of degree $\text{v}$ in $\eta$.
The virtual state vectors satisfy the Schr\"{o}dinger equation
$\mathcal{H}(\boldsymbol{\lambda})\tilde{\phi}_{\text{v}}(x;\boldsymbol{\lambda})
=\tilde{\mathcal{E}}_{\text{v}}(\boldsymbol{\lambda})
\tilde{\phi}_{\text{v}}(x;\boldsymbol{\lambda})$
(except for the upper boundary $x=N$ for finite cases)
with the virtual energies $\tilde{\mathcal{E}}_{\text{v}}(\boldsymbol{\lambda})$,
\begin{equation}
  \tilde{\mathcal{E}}_{\text{v}}(\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  -(1-c)(\text{v}+\beta)&:\text{M}\\
  -(1-aq^{-\text{v}})&:\text{l$q$L}\\
  -(1-aq^{-\text{v}})(1-bq^{\text{v}+1})&:\text{l$q$J}\\
  -(c+\text{v})(a+b-d-\text{v}-1)&:\text{R}\\
  -(1-cq^{\text{v}})(1-abd^{-1}q^{-\text{v}-1})&:\text{$q$R}
  \end{array}\right..
  \label{rdQM:Etv}
\end{equation}
Other data $\alpha(\boldsymbol{\lambda})$ and $\alpha'(\boldsymbol{\lambda})$
(see \cite{os26,os35} for details) are
\begin{equation}
  \alpha(\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  c&:\text{M}\\
  a&:\text{l$q$L,\,l$q$J}\\
  1&:\text{R}\\
  abd^{-1}q^{-1}&:\text{$q$R}
  \end{array}\right.,\quad
  \alpha^{\prime}(\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  -(1-c)\beta&:\text{M}\\
  -(1-a)&:\text{l$q$L}\\
  -(1-a)(1-bq)&:\text{l$q$J}\\
  -c(a+b-d-1)&:\text{R}\\
  -(1-c)(1-abd^{-1}q^{-1})&:\text{$q$R}
  \end{array}\right..
\end{equation}
In order to obtain well-defined quantum systems, we have to restrict the
degree $\text{v}$ and parameters $\boldsymbol{\lambda}$ so that the conditions,
$\tilde{\mathcal{E}}_{\text{v}}(\boldsymbol{\lambda})<0$, $\alpha(\boldsymbol{\lambda})>0$,
$\alpha'(\boldsymbol{\lambda})<0$, no zeros of $\xi_{\text{v}}(\eta;\boldsymbol{\lambda})$
in the domain etc., should be satisfied, see \cite{os26,os35}.
In this paper, however, we consider algebraic properties only and we
do not bother about the ranges of $\text{v}$ and $\boldsymbol{\lambda}$.

The function $\nu(x;\boldsymbol{\lambda})$ is defined as the ratio of $\phi_0$ and
$\tilde{\phi}_0$,
\begin{equation}
  \nu(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \frac{\phi_0(x;\boldsymbol{\lambda})}{\tilde{\phi}_0(x;\boldsymbol{\lambda})}.
  \label{rdQM:nu}
\end{equation}
Since the potential functions $B(x;\boldsymbol{\lambda})$ and $D(x;\boldsymbol{\lambda})$
are rational functions of $x$ or $q^x$, this function $\nu(x;\boldsymbol{\lambda})$
can be analytically continued into a meromorphic function of $x$ or $q^x$
through the functional relations:
\begin{equation}
  \nu(x+1;\boldsymbol{\lambda})=\frac{B(x;\boldsymbol{\lambda})}
  {\alpha(\boldsymbol{\lambda})B'(x;\boldsymbol{\lambda})}\nu(x;\boldsymbol{\lambda}),\quad
  \nu(x-1;\boldsymbol{\lambda})=\frac{D(x;\boldsymbol{\lambda})}
  {\alpha(\boldsymbol{\lambda})D'(x;\boldsymbol{\lambda})}\nu(x;\boldsymbol{\lambda}).
  \label{nurel}
\end{equation}
Explicitly it is
\begin{equation}
  \nu(x;\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  c^x&:\text{M}\\
  a^x&:\text{l$q$L,\,l$q$J}\\[2pt]
  {\displaystyle
  \frac{\Gamma(1-a)\Gamma(x+b)\Gamma(d-a+1)\Gamma(b-d-x)}
  {\Gamma(1-a-x)\Gamma(b)\Gamma(x+d-a+1)\Gamma(b-d)}}
  &:\text{R}\\[8pt]
  {\displaystyle
  \frac{(a^{-1}q^{1-x},b,a^{-1}dq^{x+1},bd^{-1};q)_{\infty}}
  {(a^{-1}q,bq^x,a^{-1}dq,bd^{-1}q^{-x};q)_{\infty}}}
  &:\text{$q$R}
  \end{array}\right..
\end{equation}
For non-negative integer $x$, it reduces to
\begin{equation}
  \nu(x;\boldsymbol{\lambda})=\frac{(1-a-x,b)_x}{(d-a+1,b-d-x)_x}
  \ \ \text{: R},\quad
  \nu(x;\boldsymbol{\lambda})
  =\frac{(a^{-1}q^{1-x},b;q)_x}{(a^{-1}dq,bd^{-1}q^{-x};q)_x}
  \ \ \text{: $q$R}.
\end{equation}
The functions $r_j(x;\boldsymbol{\lambda},M)$ are defined as the ratio of $\nu$'s,
\begin{equation}
  r_j(x+j-1;\boldsymbol{\lambda},M)\stackrel{\text{def}}{=}
  \frac{\nu(x+j-1;\boldsymbol{\lambda})}
  {\nu\bigl(x;\boldsymbol{\lambda}+M\tilde{\boldsymbol{\delta}}\bigr)}
  \ \ (1\leq j\leq M+1),
  \label{rdQM:rj}
\end{equation}
whose explicit forms are
\begin{equation}
  r_j(x+j-1;\boldsymbol{\lambda},M)=\left\{
  \begin{array}{ll}
  c^{j-1}&:\text{M}\\
  a^{j-1}q^{Mx}&:\text{l$q$L,\,l$q$J}\\[2pt]
  {\displaystyle
  \frac{(x+a,x+b)_{j-1}(x+d-a+j,x+d-b+j)_{M+1-j}}
  {(d-a+1,d-b+1)_M}}&:\text{R}\\[10pt]
  {\displaystyle
  \frac{(aq^x,bq^x;q)_{j-1}(a^{-1}dq^{x+j},b^{-1}dq^{x+j};q)_{M+1-j}}
  {(abd^{-1}q^{-1})^{j-1}q^{Mx}(a^{-1}dq,b^{-1}dq;q)_M}}&:\text{$q$R}
  \end{array}\right..
\end{equation}
The auxiliary function $\varphi_M(x)$ \cite{os22} is defined by
\begin{align}
  \varphi_M(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}\prod_{1\leq j<k\leq M}
  \frac{\eta(x+k-1;\boldsymbol{\lambda})-\eta(x+j-1;\boldsymbol{\lambda})}
  {\eta(k-j;\boldsymbol{\lambda})}\nonumber\\
  &=\prod_{1\leq j<k\leq M}
  \varphi\bigl(x+j-1;\boldsymbol{\lambda}+(k-j-1)\boldsymbol{\delta}\bigr),
  \label{rdQM:varphiM}
\end{align}
and $\varphi_0(x)=\varphi_1(x)=1$.
For M, l$q$L and l$q$J, its explicit form is
$\varphi_M(x)=\kappa^{-\frac12M(M-1)x-\frac16M(M-1)(M-2)}$.

In the following we adopt the convention,
\begin{equation}
  \check{P}_n(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=} 0\ \ (n\in\mathbb{Z}_{<0}),\quad
  \check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=} 0
  \ \ (\text{v}\in\mathbb{Z}_{<0}).
  \label{Pn<0=0}
\end{equation}

\subsection{Shape invariance}
\label{sec:rdQM_si}

The original systems in \S\,\ref{sec:rdQM_org} are shape invariant
\cite{os12} and they satisfy the relation,
\begin{equation}
  \mathcal{A}(\boldsymbol{\lambda})\mathcal{A}(\boldsymbol{\lambda})^{\dagger}
  =\kappa\mathcal{A}(\boldsymbol{\lambda}+\boldsymbol{\delta})^{\dagger}
  \mathcal{A}(\boldsymbol{\lambda}+\boldsymbol{\delta})+\mathcal{E}_1(\boldsymbol{\lambda}),
\end{equation}
which is a sufficient condition for exact solvability.
As a consequence of the shape invariance, the action of the operators
$\mathcal{A}(\boldsymbol{\lambda})$ and $\mathcal{A}(\boldsymbol{\lambda})^{\dagger}$
on the eigenvectors is
\begin{equation}
  \mathcal{A}(\boldsymbol{\lambda})\phi_n(x;\boldsymbol{\lambda})
  =\frac{\mathcal{E}_n(\boldsymbol{\lambda})}{\sqrt{B(0;\boldsymbol{\lambda})}}\,
  \phi_{n-1}(x;\boldsymbol{\lambda}+\boldsymbol{\delta}),\quad
  \mathcal{A}(\boldsymbol{\lambda})^{\dagger}
  \phi_{n-1}(x;\boldsymbol{\lambda}+\boldsymbol{\delta})
  =\sqrt{B(0;\boldsymbol{\lambda})}\,\phi_n(x;\boldsymbol{\lambda}).
  \label{rdQM:Aphi=}
\end{equation}
These relations \eqref{rdQM:Aphi=} are equivalent to the forward and
backward shift relations of the orthogonal polynomials $P_n$ \cite{kls},
\begin{equation}
  \mathcal{F}(\boldsymbol{\lambda})\check{P}_n(x;\boldsymbol{\lambda})
  =\mathcal{E}_n(\boldsymbol{\lambda})\check{P}_{n-1}(x;\boldsymbol{\lambda}+\boldsymbol{\delta}),
  \quad
  \mathcal{B}(\boldsymbol{\lambda})\check{P}_{n-1}(x;\boldsymbol{\lambda}+\boldsymbol{\delta})
  =\check{P}_n(x;\boldsymbol{\lambda}),
  \label{rdQM:FP=}
\end{equation}
where the forward and backward shift operators $\mathcal{F}(\boldsymbol{\lambda})$
and $\mathcal{B}(\boldsymbol{\lambda})$ are defined by
\begin{align}
  \mathcal{F}(\boldsymbol{\lambda})
  &\stackrel{\text{def}}{=}\sqrt{B(0;\boldsymbol{\lambda})}\,\phi_0(x;\boldsymbol{\lambda}+\boldsymbol{\delta})^{-1}
  \circ\mathcal{A}(\boldsymbol{\lambda})\circ\phi_0(x;\boldsymbol{\lambda})
  =B(0;\boldsymbol{\lambda})\varphi(x;\boldsymbol{\lambda})^{-1}(1-e^{\partial}),
  \label{rdQM:Fdef}\\
  \mathcal{B}(\boldsymbol{\lambda})
  &\stackrel{\text{def}}{=}\frac{1}{\sqrt{B(0;\boldsymbol{\lambda})}}\,
  \phi_0(x;\boldsymbol{\lambda})^{-1}\circ\mathcal{A}(\boldsymbol{\lambda})^{\dagger}
  \circ\phi_0(x;\boldsymbol{\lambda}+\boldsymbol{\delta})\nonumber\\
  &=B(0;\boldsymbol{\lambda})^{-1}
  \bigl(B(x;\boldsymbol{\lambda})-D(x;\boldsymbol{\lambda})e^{-\partial}\bigr)
  \varphi(x;\boldsymbol{\lambda}).
  \label{rdQM:Bdef}
\end{align}
(For a finite system, the first equation in \eqref{rdQM:FP=} holds for
$x\leq N-1$ as a matrix and vector equation.)
The similarity transformed Hamiltonian $\widetilde{\mathcal{H}}(\boldsymbol{\lambda})$
acting on the polynomial eigenvectors is square root free. It is defined by
\begin{align}
  &\widetilde{\mathcal{H}}(\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \phi_0(x;\boldsymbol{\lambda})^{-1}\circ\mathcal{H}(\boldsymbol{\lambda})
  \circ\phi_0(x;\boldsymbol{\lambda})
  =\mathcal{B}(\boldsymbol{\lambda})\mathcal{F}(\boldsymbol{\lambda})\nonumber\\
  &\phantom{\widetilde{\mathcal{H}}_{\ell}(\boldsymbol{\lambda})}
  =B(x;\boldsymbol{\lambda})(1-e^{\partial})
  +D(x;\boldsymbol{\lambda})(1-e^{-\partial}),\\
  &\widetilde{\mathcal{H}}(\boldsymbol{\lambda})\check{P}_n(x;\boldsymbol{\lambda})
  =\mathcal{E}_n(\boldsymbol{\lambda})\check{P}_n(x;\boldsymbol{\lambda}).
  \label{rdQM:HtP=EP}
\end{align}
Note that the relations \eqref{rdQM:FP=} and \eqref{rdQM:HtP=EP} (after
writing down in components) are valid for any $x\in\mathbb{R}$,
because they are `polynomial' equations.

The action of $\mathcal{F}(\boldsymbol{\lambda})$ and $\mathcal{B}(\boldsymbol{\lambda})$
on $(\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda})$
$\stackrel{\text{def}}{=}\nu(x;\boldsymbol{\lambda})^{-1}\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})$ is
\begin{equation}
  \mathcal{F}(\boldsymbol{\lambda})(\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda})
  =\tilde{\mathcal{E}}_{\text{v}}(\boldsymbol{\lambda})
  (\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda}+\boldsymbol{\delta}),
  \ \ \mathcal{B}(\boldsymbol{\lambda})
  (\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda}+\boldsymbol{\delta})
  =(\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda}).
  \label{rdQM:Fnuxi=}
\end{equation}
(For a finite system, the first equation holds for $x\leq N-1$ as a matrix
and vector equation.)
Explicitly they are equivalent to the following forward and backward shift
relations of the virtual state polynomial $\check{\xi}_{\text{v}}$ :
\begin{align}
  &A\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})
  -B\check{\xi}_{\text{v}}(x+1;\boldsymbol{\lambda})
  =C\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda}+\boldsymbol{\delta}),
  \label{rdQM:xiforward}\\
  &A=\left\{
  \begin{array}{ll}
  1&\!\!:\text{M,\,l$q$L,\,l$q$J}\\
  (x+a)(x+b)&\!\!:\text{R}\\[2pt]
  (1-aq^x)(1-bq^x)&\!\!:\text{$q$R}
  \end{array}\right.\!\!\!\!,
  \ \ B=\left\{
  \begin{array}{ll}
  c^{-1}&\!\!:\text{M}\\
  a^{-1}&\!\!:\text{l$q$L,\,l$q$J}\\
  (x+d-a+1)(x+d-b+1)&\!\!:\text{R}\\[2pt]
  \frac{ab}{dq}(1-a^{-1}dq^{x+1})(1-b^{-1}dq^{x+1})&\!\!:\text{$q$R}
  \end{array}\right.\!\!,\nonumber\\
  &C=\left\{
  \begin{array}{ll}
  -(c\beta)^{-1}(1-c)(\text{v}+\beta)&\!\!:\text{M}\\[2pt]
  -a^{-1}(1-aq^{-\text{v}})&\!\!:\text{l$q$L}\\[2pt]
  -a^{-1}(1-bq)^{-1}(1-aq^{-\text{v}})(1-bq^{\text{v+1}})
  &\!\!:\text{l$q$J}\\[2pt]
  c^{-1}(c+\text{v})(a+b-d-1-\text{v})(2x+d+1)&\!\!:\text{R}\\[2pt]
  (1-c)^{-1}(1-cq^{\text{v}})(1-abd^{-1}q^{-\text{v}-1})(1-dq^{2x+1})
  &\!\!:\text{$q$R}
  \end{array}\right.\!\!,
  \nonumber
\end{align}
and
\begin{align}
  &A\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda}+\boldsymbol{\delta})
  -B\check{\xi}_{\text{v}}(x-1;\boldsymbol{\lambda}+\boldsymbol{\delta})
  =C\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda}),
  \label{rdQM:xibackward}\\
  &A=\left\{
  \begin{array}{ll}
  x+\beta&\!\!:\text{M}\\
  1&\!\!:\text{l$q$L}\\
  1-bq^{x+1}&\!\!:\text{l$q$J}\\
  (x+c)(x+d)&\!\!:\text{R}\\
  (1-cq^x)(1-dq^x)&\!\!:\text{$q$R}
  \end{array}\right.\!\!,
  \ \ B=\left\{
  \begin{array}{ll}
  x&\!\!:\text{M}\\
  1-q^x&\!\!:\text{l$q$L,\,l$q$J}\\
  (x+d-c)x&\!\!:\text{R}\\
  c(1-c^{-1}dq^x)(1-q^x)&\!\!:\text{$q$R}
  \end{array}\right.\!\!,\nonumber\\
  &C=\left\{
  \begin{array}{ll}
  \beta&\!\!:\text{M}\\
  q^x&\!\!:\text{l$q$L}\\
  (1-bq)q^x&\!\!:\text{l$q$J}\\
  c(2x+d)&\!\!:\text{R}\\
  (1-c)(1-dq^{2x})&\!\!:\text{$q$R}
  \end{array}\right.\!\!.
  \nonumber
\end{align}
Note that the relations \eqref{rdQM:Fnuxi=} (after writing down in components)
and \eqref{rdQM:xiforward}--\eqref{rdQM:xibackward} are valid for any
$x\in\mathbb{R}$, because $\check{\xi}_{\text{v}}(x)$ is a `polynomial' and
$\nu(x)$ is a meromorphic function.

\subsection{Multi-indexed orthogonal polynomials}
\label{sec:rdQM_miop}

The original systems in \S\,\ref{sec:rdQM_org} are iso-spectrally deformed
by applying multiple Darboux transformations with virtual state vectors
$\{\tilde{\phi}_{\text{v}}(x)\}$ as seed solutions.
The virtual state vectors are labeled by the degree $\text{v}$ of the
polynomial part $\xi_{\text{v}}$.
We take $M$ virtual state vectors specified by the multi-index set
$\mathcal{D}=\{d_1,\ldots,d_M\}$ (ordered set).
The deformed Hamiltonian is denoted as
$\mathcal{H}_{\mathcal{D}}(\boldsymbol{\lambda})$.
The general formula for the eigenfunctions
$\{\phi^{\text{gen}}_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})\}$ of the deformed
system $\mathcal{H}_{\mathcal{D}}(\boldsymbol{\lambda})
\phi^{\text{gen}}_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})=\mathcal{E}_n(\boldsymbol{\lambda})
\phi^{\text{gen}}_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})$ is \cite{os26,os35}
\begin{align}
  \phi^{\text{gen}}_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})&=
  \frac{\mathcal{S}_{\mathcal{D}}(\boldsymbol{\lambda})
  \sqrt{\prod_{j=1}^M\alpha(\boldsymbol{\lambda})B'(x+j-1;\boldsymbol{\lambda})}
  \,\tilde{\phi}_0(x;\boldsymbol{\lambda})\,
  \text{W}_{\text{C}}[\check{\xi}_{d_1},\ldots,\check{\xi}_{d_M},
  \nu\check{P}_n](x;\boldsymbol{\lambda})}
  {\sqrt{\text{W}_{\text{C}}[\check{\xi}_{d_1},\ldots,\check{\xi}_{d_M}]
  (x;\boldsymbol{\lambda})\,
  \text{W}_{\text{C}}[\check{\xi}_{d_1},\ldots,\check{\xi}_{d_M}]
  (x+1;\boldsymbol{\lambda})}},\nonumber\\
  \mathcal{S}_{\mathcal{D}}(\boldsymbol{\lambda})&\stackrel{\text{def}}{=}(-1)^M
  \prod_{1\leq i<j\leq M}\text{sgn}\bigl(
  \tilde{\mathcal{E}}_{d_i}(\boldsymbol{\lambda})
  -\tilde{\mathcal{E}}_{d_j}(\boldsymbol{\lambda})\bigr),
\end{align}
where $\text{W}_{\text{C}}[f_1,\ldots,f_n](x)$ is the Casoratian
\eqref{rdQM:Wdef}.

The multi-indexed orthogonal polynomials
$\{\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})\}$ are the main parts
of the eigenfunctions $\{\phi_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})\}$ of the
deformed system $\mathcal{H}_{\mathcal{D}}(\boldsymbol{\lambda})
\phi_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})=\mathcal{E}_n(\boldsymbol{\lambda})
\phi_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})$ \cite{os26,os35}:
\begin{align}
  \phi_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})
  &\stackrel{\text{def}}{=}\psi_{\mathcal{D}}(x;\boldsymbol{\lambda})
  \check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})
  \propto\phi^{\text{gen}}_{\mathcal{D}\,n}(x;\boldsymbol{\lambda}),\quad
  \phi_{\mathcal{D}\,n}(0;\boldsymbol{\lambda})=1,\\
  \psi_{\mathcal{D}}(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \sqrt{\check{\Xi}_{\mathcal{D}}(1;\boldsymbol{\lambda})}\,
  \frac{\phi_0(x;\boldsymbol{\lambda}+M\tilde{\boldsymbol{\delta}})}
  {\sqrt{\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})\,
  \check{\Xi}_{\mathcal{D}}(x+1;\boldsymbol{\lambda})}},\quad
  \psi_{\mathcal{D}}(0;\boldsymbol{\lambda})=1,\\
  \check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \text{W}_{\text{C}}[\check{\xi}_{d_1},\ldots,\check{\xi}_{d_M}]
  (x;\boldsymbol{\lambda})
  \times\mathcal{C}_{\mathcal{D}}(\boldsymbol{\lambda})^{-1}
  \varphi_M(x;\boldsymbol{\lambda})^{-1},
  \label{rdQM:cXiD}\\
  \check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \text{W}[\check{\xi}_{d_1},\ldots,\check{\xi}_{d_M},\nu\check{P}_n]
  (x;\boldsymbol{\lambda})
  \times\mathcal{C}_{\mathcal{D},n}(\boldsymbol{\lambda})^{-1}
  \varphi_{M+1}(x;\boldsymbol{\lambda})^{-1}
  \nu(x;\boldsymbol{\lambda}+M\tilde{\boldsymbol{\delta}})^{-1},
  \label{rdQM:cPDn}\\
  \check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  P_{\mathcal{D},n}\bigl(\eta(x;\boldsymbol{\lambda}+M\boldsymbol{\delta});
  \boldsymbol{\lambda}\bigr),\quad
  \check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \Xi_{\mathcal{D}}\bigl(\eta(x;\boldsymbol{\lambda}+(M-1)\boldsymbol{\delta});
  \boldsymbol{\lambda}\bigr).
  \label{rdQM:PDn}
\end{align}
Here the universal normalization $\check{\Xi}_{\mathcal{D}}(0;\boldsymbol{\lambda})=1$
and $\check{P}_{\mathcal{D},n}(0;\boldsymbol{\lambda})=1$ are adopted,
which determines the constants $\mathcal{C}_{\mathcal{D}}(\boldsymbol{\lambda})$
and $\mathcal{C}_{\mathcal{D},n}(\boldsymbol{\lambda})$
(convention: $\prod\limits_{1\leq j<k\leq M}\!\!\!\!\!*=1$ for $M=1$),
\begin{align}
  \mathcal{C}_{\mathcal{D}}(\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \frac{1}{\varphi_M(0;\boldsymbol{\lambda})}
  \prod_{1\leq j<k\leq M}
  \frac{\tilde{\mathcal{E}}_{d_j}(\boldsymbol{\lambda})
  -\tilde{\mathcal{E}}_{d_k}(\boldsymbol{\lambda})}
  {\alpha(\boldsymbol{\lambda})B'(j-1;\boldsymbol{\lambda})},
  \label{rdQM:CD}\\
  \mathcal{C}_{\mathcal{D},n}(\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  (-1)^M\mathcal{C}_{\mathcal{D}}(\boldsymbol{\lambda})
  \tilde{d}_{\mathcal{D},n}(\boldsymbol{\lambda})^2,\quad
  \tilde{d}_{\mathcal{D},n}(\boldsymbol{\lambda})^2\stackrel{\text{def}}{=}
  \frac{\varphi_M(0;\boldsymbol{\lambda})}{\varphi_{M+1}(0;\boldsymbol{\lambda})}
  \prod_{j=1}^M\frac{\mathcal{E}_n(\boldsymbol{\lambda})
  -\tilde{\mathcal{E}}_{d_j}(\boldsymbol{\lambda})}
  {\alpha(\boldsymbol{\lambda})B'(j-1;\boldsymbol{\lambda})}.
  \label{rdQM:CDn}
\end{align}
The denominator polynomial $\Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})$ and the
multi-indexed orthogonal polynomial $P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})$
\eqref{rdQM:PDn} are polynomials in $\eta$ and their degrees are generically
$\ell_{\mathcal{D}}$ and $\ell_{\mathcal{D}}+n$, respectively.
Here $\ell_{\mathcal{D}}$ is
$\ell_{\mathcal{D}}\stackrel{\text{def}}{=}\sum_{j=1}^{M}d_j-\tfrac12M(M-1)$.
Note that
$\check{P}_{\mathcal{D},0}(x;\boldsymbol{\lambda})
=\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda}+\boldsymbol{\delta})$.

In \cite{os26}, the multi-indexed orthogonal polynomials
$\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ \eqref{rdQM:cPDn} for R and $q$R
are expressed as
\begin{align}
  \check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})
  &=\mathcal{C}_{\mathcal{D},n}(\boldsymbol{\lambda})^{-1}
  \varphi_{M+1}(x;\boldsymbol{\lambda})^{-1}\nonumber\\
  &\quad\times\left|
  \begin{array}{cccc}
  \check{\xi}_{d_1}(x_1)&\cdots&\check{\xi}_{d_M}(x_1)
  &r_1(x_1)\check{P}_n(x_1)\\
  \check{\xi}_{d_1}(x_2)&\cdots&\check{\xi}_{d_M}(x_2)
  &r_2(x_2)\check{P}_n(x_2)\\
  \vdots&\cdots&\vdots&\vdots\\
  \check{\xi}_{d_1}(x_{M+1})&\cdots&\check{\xi}_{d_M}(x_{M+1})
  &r_{M+1}(x_{M+1})\check{P}_n(x_{M+1})\\
  \end{array}\right|,
  \label{cPDn}
\end{align}
where $x_j\stackrel{\text{def}}{=} x+j-1$ and $r_j(x)=r_j(x;\boldsymbol{\lambda},M)$
($1\leq j\leq M+1$) are given by \eqref{rdQM:rj}.
This expression is also valid for M, l$q$L and l$q$J.

In the next subsection we will rewrite \eqref{rdQM:cXiD}--\eqref{rdQM:cPDn}
by using the identities implied by the shape invariance and the properties
of the Casoratians.

\subsection{New determinant expressions}
\label{sec:rdQM_exp}

\subsubsection{Casoratian}
\label{sec:rdQM_Caso}

The Casorati determinant of a set of $n$ functions $\{f_j(x)\}$ is defined by
\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},
  \label{rdQM:Wdef}
\end{equation}
(for $n=0$, we set $\text{W}_{\text{C}}[\cdot](x)=1$), which satisfies
\begin{equation}
  \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}
\end{equation}

In the rest of this section we consider the Casoratians for $\check{P}_n(x)$,
$\check{\xi}_{\text{v}}(x)$ and $\nu(x)$.
Since they are polynomials or meromorphic functions, we can realize the
shift operator as an exponential of the differential operator,
$e^{a\frac{d}{dx}}f(x)=f(x+a)$ ($a\in\mathbb{R}$).
Contrary to the matrices $e^{\partial}$ and $e^{-\partial}$, operators
$e^{\frac{d}{dx}}$ and $e^{-\frac{d}{dx}}$ are inverse to each other.

The following determinant formula holds for any smooth functions $q_j(x)$ and
$r_j(x)$ ($j=1,2,\ldots$):
\begin{equation}
  \det\Bigl(\laprod{l=1}{j-1}\hat{D}_l\cdot f_k(x)\Bigr)_{1\leq j,k\leq n}
  =\prod_{l=1}^n\prod_{m=1}^{n-l}r_m(x+l-1)\cdot
  \text{W}_{\text{C}}[f_1,\ldots,f_n](x),
  \label{rdQM:WCid}
\end{equation}
where operators $\hat{D}_j$ ($j=1,2,\ldots$) are
\begin{equation}
  \hat{D}_j=q_j(x)+r_j(x)e^{\frac{d}{dx}},
\end{equation}
and the ordered product is
\begin{equation}
  \laprod{j=1}{n}a_j\stackrel{\text{def}}{=} a_n\cdots a_2a_1.
\end{equation}
This formula is shown by using the properties of the determinants
(row properties) and induction in $n$.

We define operators $\hat{\mathcal{F}}'(\boldsymbol{\lambda})$ and
$\hat{\mathcal{B}}'(\boldsymbol{\lambda})$ as
$\mathcal{F}\bigl(\mathfrak{t}(\boldsymbol{\lambda})\bigr)$ and
$\mathcal{B}\bigl(\mathfrak{t}(\boldsymbol{\lambda})\bigr)$ with the replacement
$e^{\pm\partial}\mapsto e^{\pm\frac{d}{dx}}$,
\begin{align}
  \hat{\mathcal{F}}'(\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  B'(0;\boldsymbol{\lambda})\varphi(x;\boldsymbol{\lambda})^{-1}(1-e^{\frac{d}{dx}}),\nonumber\\
  \hat{\mathcal{B}}'(\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  B'(0;\boldsymbol{\lambda})^{-1}
  \bigl(B'(x;\boldsymbol{\lambda})-D'(x;\boldsymbol{\lambda})e^{-\frac{d}{dx}}\bigr)
  \varphi(x;\boldsymbol{\lambda}).
\end{align}
Then the relations \eqref{rdQM:FP=} and \eqref{rdQM:Fnuxi=} with the
replacement $\boldsymbol{\lambda}\mapsto\mathfrak{t}(\boldsymbol{\lambda})$ and
exchange $n\leftrightarrow\text{v}$ give
\begin{alignat}{2}
  \hat{\mathcal{F}}'(\boldsymbol{\lambda})\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})
  &=\mathcal{E}'_{\text{v}}(\boldsymbol{\lambda})
  \check{\xi}_{\text{v}-1}(x;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}}),&\quad
  \hat{\mathcal{B}}'(\boldsymbol{\lambda})
  \check{\xi}_{\text{v}-1}(x;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})
  &=\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda}),
  \label{rdQM:Fh'xiv=}\\
  \hat{\mathcal{F}}'(\boldsymbol{\lambda})(\nu\check{P}_n)(x;\boldsymbol{\lambda})
  &=\tilde{\mathcal{E}}'_n(\boldsymbol{\lambda})
  (\nu\check{P}_n)(x;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}}),&\quad
  \hat{\mathcal{B}}'(\boldsymbol{\lambda})
  (\nu\check{P}_n)(x;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})
  &=(\nu\check{P}_n)(x;\boldsymbol{\lambda}),
  \label{rdQM:Fh'nuPn=}
\end{alignat}
where $(\nu\check{P}_n)(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=}\nu(x;\boldsymbol{\lambda})
\check{P}_n(x;\boldsymbol{\lambda})$, and $\mathcal{E}'_{\text{v}}(\boldsymbol{\lambda})$ and
$\tilde{\mathcal{E}}'_n(\boldsymbol{\lambda})$ are
\begin{equation}
  \mathcal{E}'_{\text{v}}(\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \mathcal{E}_{\text{v}}\bigl(\mathfrak{t}(\boldsymbol{\lambda})\bigr),\quad
  \tilde{\mathcal{E}}'_n(\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \tilde{\mathcal{E}}_n\bigl(\mathfrak{t}(\boldsymbol{\lambda})\bigr).
\end{equation}
These relations \eqref{rdQM:Fh'xiv=}--\eqref{rdQM:Fh'nuPn=} are valid for
any $x\in\mathbb{R}$.

By taking $\hat{\mathcal{F}}'$ or $\hat{\mathcal{B}}'$ as $\hat{D}_j$
and using shape invariance properties
\eqref{rdQM:Fh'xiv=}--\eqref{rdQM:Fh'nuPn=}, the Casoratians
$\text{W}_{\text{C}}[\check{\xi}_{d_1},\ldots,\check{\xi}_{d_M}]
(x;\boldsymbol{\lambda})$ and $\text{W}_{\text{C}}[\check{\xi}_{d_1},\ldots,
\check{\xi}_{d_M},\nu\check{P}_n](x;\boldsymbol{\lambda})$
in \eqref{rdQM:cXiD}--\eqref{rdQM:cPDn} can be rewritten in various ways.
We consider two typical cases:
\begin{align}
  \text{A}:\ \,&
  \hat{D}_j=\hat{\mathcal{F}}'\bigl(\boldsymbol{\lambda}
  +(j-1)\tilde{\boldsymbol{\delta}}\bigr)\ \ (j=1,2,\ldots),
  \label{rdQM:A}\\
  \text{B}:\ \,&
  \hat{D}_{2l-1}=e^{(l-1)\frac{d}{dx}}\circ\hat{\mathcal{F}}'(\boldsymbol{\lambda})
  \circ e^{-(l-1)\frac{d}{dx}},
  \ \hat{D}_{2l}=e^{l\frac{d}{dx}}\circ\hat{\mathcal{B}}'(\boldsymbol{\lambda})
  \circ e^{-(l-1)\frac{d}{dx}}
  \ \ (l=1,2,\ldots),\!
  \label{rdQM:B}
\end{align}
which correspond to the cases A and B in \cite{os37}.

\subsubsection{case A}
\label{sec:rdQM_caseA}

Firstly we consider case A \eqref{rdQM:A}.
The functions $q_j(x)$ and $r_j(x)$ are
\begin{equation}
  -q_j(x)=r_j(x)=-B'\bigl(0;\boldsymbol{\lambda}+(j-1)\tilde{\boldsymbol{\delta}}\bigr)
  \varphi\bigl(0;\boldsymbol{\lambda}+(j-1)\boldsymbol{\delta}\bigr)^{-1}.
\end{equation}
The shape invariance properties \eqref{rdQM:Fh'xiv=}--\eqref{rdQM:Fh'nuPn=}
give ($k=0,1,\ldots$)
\begin{align}
  \laprod{l=1}{k}\hat{D}_l\cdot\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})
  &=\prod_{m=0}^{k-1}\mathcal{E}'_{\text{v}-m}
  (\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})\cdot
  \check{\xi}_{\text{v}-k}(x;\boldsymbol{\lambda}+k\tilde{\boldsymbol{\delta}}),\nonumber\\[-2pt]
  \laprod{l=1}{k}\hat{D}_l\cdot(\nu\check{P}_n)(x;\boldsymbol{\lambda})
  &=\prod_{m=0}^{k-1}\tilde{\mathcal{E}}'_n
  (\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})\cdot
  (\nu\check{P}_n)(x;\boldsymbol{\lambda}+k\tilde{\boldsymbol{\delta}}).
  \label{rdQM:siA}
\end{align}

We rewrite \eqref{rdQM:cXiD}--\eqref{rdQM:cPDn} in the following way.
(\romannumeral1) By using \eqref{rdQM:WCid},
rewrite the Casoratians in \eqref{rdQM:cXiD}--\eqref{rdQM:cPDn} as
determinants $\det(a_{j,k})_{1\leq j,k\leq M'}$ ($M'=M$ for
$\check{\Xi}_{\mathcal{D}}$ and $M'=M+1$ for $\check{P}_{\mathcal{D},n}$).
(\romannumeral2) Rewrite each matrix element $a_{j,k}$ by \eqref{rdQM:siA}.
For $\check{P}_{\mathcal{D},n}$, we do the following (\romannumeral3) and
(\romannumeral4).
(\romannumeral3) Divide the $(M+1)$-th column of $\det(a_{j,k})$ by
$\nu(x;\boldsymbol{\lambda}+M\tilde{\boldsymbol{\delta}})$, and multiply the determinant
by $\nu(x;\boldsymbol{\lambda}+M\tilde{\boldsymbol{\delta}})$.
(\romannumeral4) Rewrite
$\nu\bigl(x;\boldsymbol{\lambda}+(j-1)\tilde{\boldsymbol{\delta}}\bigr)/
\nu(x;\boldsymbol{\lambda}+M\tilde{\boldsymbol{\delta}})$
in the $(M+1)$-th column of $\det(a_{j,k})$ by
\begin{equation}
  \frac{\nu\bigl(x;\boldsymbol{\lambda}+(j-1)\tilde{\boldsymbol{\delta}}\bigr)}
  {\nu(x;\boldsymbol{\lambda}+M\tilde{\boldsymbol{\delta}})}=\left\{
  \begin{array}{ll}
  1&:\text{M}\\
  q^{(M+1-j)x}=\bigl(1-\eta(x)\bigr)^{M+1-j}&:\text{l$q$L,\,l$q$J}\\[2pt]
  {\displaystyle
  \frac{(x+d-a+j,x+d-b+j)_{M+1-j}}{(d-a+j,d-b+j)_{M+1-j}}}&:\text{R}\\[10pt]
  {\displaystyle
  q^{-(M+1-j)x}\frac{(a^{-1}dq^{x+j},b^{-1}dq^{x+j};q)_{M+1-j}}
  {(a^{-1}dq^j,b^{-1}dq^j;q)_{M+1-j}}}&:\text{$q$R}
  \end{array}\right..
  \label{rdQM:A:nu/nu}
\end{equation}
Then $\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})$ and
$\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ are expressed in terms of
$\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda}')$,
$\check{P}_n(x;\boldsymbol{\lambda}')$, $\varphi(x';\boldsymbol{\lambda}')$ and
\eqref{rdQM:A:nu/nu}.
Straightforward calculation shows that the factors
$\varphi(x';\boldsymbol{\lambda}')$ are canceled out.
Thus, for M, l$q$L and l$q$J cases,
$\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})$ and
$\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ are expressed in terms of
$\eta(x)$, namely
$\Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})$ and
$P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})$ are expressed without explicit
$x$-dependence (which is trivial for M, because $\eta(x)=x$).
Their final forms are as follows.\\
$\circ$ {\bf Meixner:}
The denominator polynomial $\Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})$ is
\begin{equation}
  \Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})
  =\frac{(-c)^{\frac12M(M-1)}}{\mathcal{C}_{\mathcal{D}}(\boldsymbol{\lambda})
  \prod_{m=1}^{M-1}(\beta)_m}\det(a_{j,k})_{1\leq j,k\leq M},
\end{equation}
where $a_{j,k}$ are
\begin{equation}
  a_{j,k}=\prod_{m=0}^{j-2}\mathcal{E}'_{d_k-m}
  (\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})\cdot
  \xi_{d_k+1-j}\bigl(\eta;\boldsymbol{\lambda}+(j-1)\tilde{\boldsymbol{\delta}}\bigr)
  \ \ (j,k=1,\ldots,M),
  \label{rdQM:A:M:ajk}
\end{equation}
and the multi-indexed orthogonal polynomial
$P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})$ is
\begin{equation}
  P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})
  =\frac{(-c)^{\frac12M(M+1)}}{\mathcal{C}_{\mathcal{D},n}(\boldsymbol{\lambda})
  \prod_{m=1}^M(\beta)_m}\det(a_{j,k})_{1\leq j,k\leq M+1},
\end{equation}
where $a_{j,k}$ are
\begin{align}
  a_{j,k}&=\text{eq.\,\eqref{rdQM:A:M:ajk}}
  \ \ (j=1,\ldots,M+1;k=1,\ldots,M),\nonumber\\
  a_{j,M+1}&=\prod_{m=0}^{j-2}\tilde{\mathcal{E}}'_n
  (\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})\cdot
  P_n\bigl(\eta;\boldsymbol{\lambda}+(j-1)\tilde{\boldsymbol{\delta}}\bigr)
  \ \ (j=1,\ldots,M+1).
\end{align}
$\circ$ {\bf little $q$-Laguerre and little $q$-Jacobi:}\\
The denominator polynomial $\Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})$ is
\begin{equation}
  \Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})
  =\frac{(-a)^{\frac12M(M-1)}}
  {q^{\frac16M(M-1)(M-2)}\mathcal{C}_{\mathcal{D}}(\boldsymbol{\lambda})}
  \det(a_{j,k})_{1\leq j,k\leq M}\times\left\{
  \begin{array}{ll}
  1&:\text{l$q$L}\\
  \prod_{m=1}^{M-1}(bq;q)_m^{-1}&:\text{l$q$J}
  \end{array}\right.,
\end{equation}
where $a_{j,k}$ are given by \eqref{rdQM:A:M:ajk},
and the multi-indexed orthogonal polynomial
$P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})$ is
\begin{equation}
  P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})
  =\frac{(-a)^{\frac12M(M+1)}}
  {q^{\frac16(M+1)M(M-1)}\mathcal{C}_{\mathcal{D},n}(\boldsymbol{\lambda})}
  \det(a_{j,k})_{1\leq j,k\leq M+1}\times\left\{
  \begin{array}{ll}
  1&:\text{l$q$L}\\
  \prod_{m=1}^M(bq;q)_m^{-1}&:\text{l$q$J}
  \end{array}\right.,
\end{equation}
where $a_{j,k}$ are
\begin{align}
  a_{j,k}&=\text{eq.\,\eqref{rdQM:A:M:ajk}}
  \ \ (j=1,\ldots,M+1;k=1,\ldots,M),\nonumber\\
  a_{j,M+1}&=\prod_{m=0}^{j-2}\tilde{\mathcal{E}}'_n
  (\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})\cdot(1-\eta)^{M+1-j}
  P_n\bigl(\eta;\boldsymbol{\lambda}+(j-1)\tilde{\boldsymbol{\delta}}\bigr)
  \ \ (j=1,\ldots,M+1).
\end{align}
$\circ$ {\bf Racah and $q$-Racah:}\\
The denominator polynomial $\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})$ is
\begin{equation}
  \check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})
  =\frac{(-1)^{\frac12M(M-1)}}
  {\mathcal{C}_{\mathcal{D}}(\boldsymbol{\lambda})
  \prod_{m=0}^{M-1}B'(0;\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})^{M-1-m}}
  \det(a_{j,k})_{1\leq j,k\leq M},
\end{equation}
where $a_{j,k}$ are
\begin{equation}
  a_{j,k}=\prod_{m=0}^{j-2}\mathcal{E}'_{d_k-m}
  (\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})\cdot
  \check{\xi}_{d_k+1-j}\bigl(x;\boldsymbol{\lambda}+(j-1)\tilde{\boldsymbol{\delta}}\bigr)
  \ \ (j,k=1,\ldots,M),
  \label{rdQM:A:ajk}
\end{equation}
and the multi-indexed orthogonal polynomial
$\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ is
\begin{equation}
  \check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})
  =\frac{(-1)^{\frac12M(M+1)}}
  {\mathcal{C}_{\mathcal{D},n}(\boldsymbol{\lambda})
  \prod_{m=0}^MB'(0;\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})^{M-m}}
  \det(a_{j,k})_{1\leq j,k\leq M+1},
\end{equation}
where $a_{j,k}$ are
\begin{align}
  a_{j,k}&=\text{eq.\,\eqref{rdQM:A:ajk}}
  \ \ (j=1,\ldots,M+1;k=1,\ldots,M),\nonumber\\
  a_{j,M+1}&=\prod_{m=0}^{j-2}\tilde{\mathcal{E}}'_n
  (\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})\cdot
  P_n\bigl(\eta;\boldsymbol{\lambda}+(j-1)\tilde{\boldsymbol{\delta}}\bigr)\nonumber\\
  &\quad\times\left\{
  \begin{array}{ll}
  {\displaystyle
  \frac{(x+d-a+j,x+d-b+j)_{M+1-j}}{(d-a+j,d-b+j)_{M+1-j}}}&:\text{R}\\[10pt]
  {\displaystyle
  q^{-(M+1-j)x}\frac{(a^{-1}dq^{x+j},b^{-1}dq^{x+j};q)_{M+1-j}}
  {(a^{-1}dq^j,b^{-1}dq^j;q)_{M+1-j}}}&:\text{$q$R}
  \end{array}\right.
  \ \ (j=1,\ldots,M+1).
\end{align}

\subsubsection{case B}
\label{sec:rdQM_caseB}

Secondly we consider case B \eqref{rdQM:B}.
The functions $q_j(x)$ and $r_j(x)$ are ($l=1,2,\ldots$)
\begin{align}
  -q_{2l-1}(x)&=r_{2l-1}(x)=-B'(0;\boldsymbol{\lambda})
  \varphi(x+l-1;\boldsymbol{\lambda})^{-1},\nonumber\\
  q_{2l}(x)&=-B'(0;\boldsymbol{\lambda})^{-1}D'(x+l;\boldsymbol{\lambda})
  \varphi(x+l-1;\boldsymbol{\lambda}),\\
  r_{2l}(x)&=B'(0;\boldsymbol{\lambda})^{-1}B'(x+l;\boldsymbol{\lambda})
  \varphi(x+l;\boldsymbol{\lambda}).
  \nonumber
\end{align}
The shape invariance properties \eqref{rdQM:Fh'xiv=}--\eqref{rdQM:Fh'nuPn=}
give ($l=1,2,\ldots$)
\begin{align}
  \laprod{m=1}{2l-2}\hat{D}_m\cdot\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})
  &=\mathcal{E}'_{\text{v}}(\boldsymbol{\lambda})^{l-1}
  \check{\xi}_{\text{v}}(x+l-1;\boldsymbol{\lambda}),\nonumber\\[-2pt]
  \laprod{m=1}{2l-1}\hat{D}_m\cdot\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})
  &=\mathcal{E}'_{\text{v}}(\boldsymbol{\lambda})^l
  \check{\xi}_{\text{v}-1}(x+l-1;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}}),\nonumber\\[-2pt]
  \laprod{m=1}{2l-2}\hat{D}_m\cdot(\nu\check{P}_n)(x;\boldsymbol{\lambda})
  &=\tilde{\mathcal{E}}'_n(\boldsymbol{\lambda})^{l-1}
  (\nu\check{P}_n)(x+l-1;\boldsymbol{\lambda}),
  \label{rdQM:siB}\\
  \laprod{m=1}{2l-1}\hat{D}_m\cdot(\nu\check{P}_n)(x;\boldsymbol{\lambda})
  &=\tilde{\mathcal{E}}'_n(\boldsymbol{\lambda})^l
  (\nu\check{P}_n)(x+l-1;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}}).
  \nonumber
\end{align}

Like case A, we rewrite \eqref{rdQM:cXiD}--\eqref{rdQM:cPDn} by the steps
(\romannumeral1)--(\romannumeral4).
In (\romannumeral4) we use
\begin{align}
  &\frac{\nu(x+l-1;\boldsymbol{\lambda})}
  {\nu(x;\boldsymbol{\lambda}+M\tilde{\boldsymbol{\delta}})}=\left\{
  \begin{array}{ll}
  c^{l-1}&:\text{M}\\
  a^{l-1}q^{Mx}&:\text{l$q$L,\,l$q$J}\\[2pt]
  {\displaystyle
  \frac{(x+d-a+l,x+d-b+l)_{M+1-l}(x+a,x+b)_{l-1}}
  {(d-a+1,d-b+1)_M}}&:\text{R}\\[10pt]
  {\displaystyle
  \frac{(a^{-1}dq^{x+l},b^{-1}dq^{x+l};q)_{M+1-l}(aq^x,bq^x;q)_{l-1}}
  {(abd^{-1}q^{-1})^{l-1}q^{Mx}
  (a^{-1}dq,b^{-1}dq;q)_M}}&:\text{$q$R}
  \end{array}\right.,\nonumber\\
  &\frac{\nu(x+l-1;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})}
  {\nu(x;\boldsymbol{\lambda}+M\tilde{\boldsymbol{\delta}})}
  \label{rdQM:B:nu/nu}\\
  &=\left\{
  \begin{array}{ll}
  c^{l-1}&:\text{M}\\
  (aq^{-1})^{l-1}q^{(M-1)x}&:\text{l$q$L,\,l$q$J}\\[2pt]
  {\displaystyle
  \frac{(x+d-a+l+1,x+d-b+l+1)_{M-l}(x+a,x+b)_{l-1}}
  {(d-a+2,d-b+2)_{M-1}}}&:\text{R}\\[10pt]
  {\displaystyle
  \frac{(a^{-1}dq^{x+l+1},b^{-1}dq^{x+l+1};q)_{M-l}(aq^x,bq^x;q)_{l-1}}
  {(abd^{-1}q^{-2})^{l-1}q^{(M-1)x}
  (a^{-1}dq^2,b^{-1}dq^2;q)_{M-1}}}&:\text{$q$R}
  \end{array}\right..
  \nonumber
\end{align}
Contrary to case A, the determinants contain
$\check{\xi}_{\text{v}}(x';\boldsymbol{\lambda}')$ and
$\check{P}_n(x';\boldsymbol{\lambda}')$ at various points $x'$,
and the factors $\varphi(x';\boldsymbol{\lambda}')$ are not canceled out.
Their final forms are as follows.
The denominator polynomial $\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})$ is
\begin{align}
  \check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})
  &=\frac{(-1)^{\frac12(1+(-1)^M)\frac{M}{2}}}
  {\mathcal{C}_{\mathcal{D}}(\boldsymbol{\lambda})B'(0;\boldsymbol{\lambda})^{[\frac{M}{2}]}}
  \prod_{m=0}^{[\frac{M-2}{2}]}\varphi(x+m;\boldsymbol{\lambda})\cdot
  \varphi_M(x;\boldsymbol{\lambda})^{-1}\nonumber\\
  &\quad\times\biggl(
  \prod_{m=0}^{[\frac{M-2}{2}]}\prod_{l=1}^{M-2-2m}
  B'(x+l+m;\boldsymbol{\lambda})\biggr)^{-1}
  \det(a_{j,k})_{1\leq j,k\leq M},
\end{align}
where $a_{j,k}$ are
\begin{align}
  a_{2l-1,k}&=\mathcal{E}'_{d_k}(\boldsymbol{\lambda})^{l-1}
  \check{\xi}_{d_k}(x+l-1;\boldsymbol{\lambda})
  \ \ (l=1,\ldots,[\tfrac{M+1}{2}];k=1,\ldots,M),\nonumber\\
  a_{2l,k}&=\mathcal{E}'_{d_k}(\boldsymbol{\lambda})^l
  \check{\xi}_{d_k-1}(x+l-1;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})
  \ \ (l=1,\ldots,[\tfrac{M}{2}];k=1,\ldots,M),
\end{align}
and the multi-indexed orthogonal polynomial
$\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ is
\begin{align}
  \check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})
  &=\frac{(-1)^{\frac12(1-(-1)^M)\frac{M+1}{2}}}
  {\mathcal{C}_{\mathcal{D},n}(\boldsymbol{\lambda})
  B'(0;\boldsymbol{\lambda})^{[\frac{M+1}{2}]}}
  \prod_{m=0}^{[\frac{M-1}{2}]}\varphi(x+m;\boldsymbol{\lambda})\cdot
  \varphi_{M+1}(x;\boldsymbol{\lambda})^{-1}\nonumber\\
  &\quad\times\biggl(
  \prod_{m=0}^{[\frac{M-1}{2}]}\prod_{l=1}^{M-1-2m}
  B'(x+l+m;\boldsymbol{\lambda})\biggr)^{-1}
  \det(a_{j,k})_{1\leq j,k\leq M+1},
\end{align}
where $a_{j,k}$ are
\begin{align}
  a_{2l-1,k}&=\mathcal{E}'_{d_k}(\boldsymbol{\lambda})^{l-1}
  \check{\xi}_{d_k}(x+l-1;\boldsymbol{\lambda})
  \ \ (l=1,\ldots,[\tfrac{M+2}{2}];k=1,\ldots,M),\nonumber\\
  a_{2l,k}&=\mathcal{E}'_{d_k}(\boldsymbol{\lambda})^l
  \check{\xi}_{d_k-1}(x+l-1;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})
  \ \ (l=1,\ldots,[\tfrac{M+1}{2}];k=1,\ldots,M),\nonumber\\
  a_{2l-1,M+1}&=\tilde{\mathcal{E}}'_n(\boldsymbol{\lambda})^{l-1}
  \check{P}_n(x+l-1;\boldsymbol{\lambda})
  \qquad(l=1,\ldots,[\tfrac{M+2}{2}])\nonumber\\
  &\quad\times\left\{
  \begin{array}{ll}
  c^{l-1}&:\text{M}\\
  a^{l-1}q^{Mx}&:\text{l$q$L,\,l$q$J}\\[2pt]
  {\displaystyle
  \frac{(x+d-a+l,x+d-b+l)_{M+1-l}(x+a,x+b)_{l-1}}
  {(d-a+1,d-b+1)_M}}&:\text{R}\\[10pt]
  {\displaystyle
  \frac{(a^{-1}dq^{x+l},b^{-1}dq^{x+l};q)_{M+1-l}(aq^x,bq^x;q)_{l-1}}
  {(abd^{-1}q^{-1})^{l-1}q^{Mx}
  (a^{-1}dq,b^{-1}dq;q)_M}}&:\text{$q$R}
  \end{array}\right.,\nonumber\\
  a_{2l,M+1}&=\tilde{\mathcal{E}}'_n(\boldsymbol{\lambda})^l
  \check{P}_n(x+l-1;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})
  \qquad(l=1,\ldots,[\tfrac{M+1}{2}])\\
  &\quad\times\left\{
  \begin{array}{ll}
  c^{l-1}&:\text{M}\\
  (aq^{-1})^{l-1}q^{(M-1)x}&:\text{l$q$L,\,l$q$J}\\[2pt]
  {\displaystyle
  \frac{(x+d-a+l+1,x+d-b+l+1)_{M-l}(x+a,x+b)_{l-1}}
  {(d-a+2,d-b+2)_{M-1}}}&:\text{R}\\[10pt]
  {\displaystyle
  \frac{(a^{-1}dq^{x+l+1},b^{-1}dq^{x+l+1};q)_{M-l}(aq^x,bq^x;q)_{l-1}}
  {(abd^{-1}q^{-2})^{l-1}q^{(M-1)x}
  (a^{-1}dq^2,b^{-1}dq^2;q)_{M-1}}}&:\text{$q$R}
  \end{array}\right..
  \nonumber
\end{align}

\section{Multi-indexed Orthogonal Polynomials in idQM}
\label{sec:idQM}

In this section we derive various equivalent expressions for the
multi-indexed Wilson and Askey-Wilson polynomials in the framework of
the discrete quantum mechanics with pure imaginary shifts.

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

The Hamiltonian of the discrete quantum mechanics with pure imaginary shifts
(idQM) \cite{os13_1,os13_2},
whose dynamical variables are the real coordinate
$x$ ($x_1\leq x\leq x_2$) and the conjugate momentum $p=-i\frac{d}{dx}$, is
\begin{align}
  &\mathcal{H}\stackrel{\text{def}}{=}\sqrt{V(x)}\,e^{\gamma p}\sqrt{V^*(x)}
  +\!\sqrt{V^*(x)}\,e^{-\gamma p}\sqrt{V(x)}
  -V(x)-V^*(x)=\mathcal{A}^{\dagger}\mathcal{A},
  \label{idQM:H}\\
  &\mathcal{A}\stackrel{\text{def}}{=} i\bigl(e^{\frac{\gamma}{2}p}\sqrt{V^*(x)}
  -e^{-\frac{\gamma}{2}p}\sqrt{V(x)}\,\bigr),\quad
  \mathcal{A}^{\dagger}\stackrel{\text{def}}{=} -i\bigl(\sqrt{V(x)}\,e^{\frac{\gamma}{2}p}
  -\sqrt{V^*(x)}\,e^{-\frac{\gamma}{2}p}\bigr).
\end{align}
Here the potential function $V(x)$ is an analytic function of $x$ and
$\gamma$ is a real constant.
The $*$-operation on an analytic function $f(x)=\sum_na_nx^n$
($a_n\in\mathbb{C}$) is defined by $f^*(x)=\sum_na_n^*x^n$, in which
$a_n^*$ is the complex conjugation of $a_n$.
Since the momentum operator appears in exponentiated forms,
the Schr\"{o}dinger equation,
\begin{equation}
  \mathcal{H}\phi_n(x)=\mathcal{E}_n\phi_n(x)
  \ \ (n=0,1,2,\ldots),\quad
  0=\mathcal{E}_0<\mathcal{E}_1<\mathcal{E}_2<\cdots,
\end{equation}
is an analytic difference equation with pure imaginary shifts.

We consider idQM described by the Wilson and Askey-Wilson polynomials.
Various parameters are
\begin{alignat}{2}
  \text{W}:\ \ &x_1=0,\ x_2=\infty,\ \gamma=1,
  &\ \ \boldsymbol{\lambda}=(a_1,a_2,a_3,a_4),
  &\ \ \boldsymbol{\delta}=(\tfrac12,\tfrac12,\tfrac12,\tfrac12),
  \ \ \kappa=1,\\
  \text{AW}:\ \ &x_1=0,\ x_2=\pi,\ \gamma=\log q,
  &\ \ q^{\boldsymbol{\lambda}}=(a_1,a_2,a_3,a_4),
  &\ \ \boldsymbol{\delta}=(\tfrac12,\tfrac12,\tfrac12,\tfrac12),
  \ \ \kappa=q^{-1},
\end{alignat}
where $0<q<1$.
The parameters are restricted by
\begin{equation}
  \{a_1^*,a_2^*,a_3^*,a_4^*\}=\{a_1,a_2,a_3,a_4\}\ \ (\text{as a set});\quad
  \text{W}:\ \text{Re}\,a_i>0,\quad
  \text{AW}:\ |a_i|<1.
  \label{rangeorg}
\end{equation}
Here are the fundamental data \cite{os13_1}:
\begin{align}
  &V(x;\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  \bigl(2ix(2ix+1)\bigr)^{-1}\prod_{j=1}^4(a_j+ix)&:\text{W}\\[2pt]
  \bigl((1-e^{2ix})(1-qe^{2ix})\bigr)^{-1}\prod_{j=1}^4(1-a_je^{ix})
  &:\text{AW}
  \end{array}\right.,
  \label{Vform}\\
  &\eta(x)=\left\{
  \begin{array}{ll}
  x^2&:\text{W}\\
  \cos x&:\text{AW}
  \end{array}\right.,\quad
  \varphi(x)=\left\{
  \begin{array}{ll}
  2x&:\text{W}\\
  2\sin x&:\text{AW}
  \end{array}\right.,
  \label{etadef}\\
  &\mathcal{E}_n(\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  n(n+b_1-1)&:\text{W}\\[2pt]
  (q^{-n}-1)(1-b_4q^{n-1})&:\text{AW}
  \end{array}\right.,\quad
  \begin{array}{l}
  b_1\stackrel{\text{def}}{=} a_1+a_2+a_3+a_4,\\
  b_4\stackrel{\text{def}}{=} a_1a_2a_3a_4,
  \end{array}\\
  &\phi_n(x;\boldsymbol{\lambda})
  =\phi_0(x;\boldsymbol{\lambda})\check{P}_n(x;\boldsymbol{\lambda}),
  \label{factphin}\\
  &\check{P}_n(x;\boldsymbol{\lambda})=P_n\bigl(\eta(x);\boldsymbol{\lambda}\bigr)
  =\left\{\begin{array}{ll}
  W_n\bigl(\eta(x);a_1,a_2,a_3,a_4\bigr)&:\text{W}\\[2pt]
  p_n\bigl(\eta(x);a_1,a_2,a_3,a_4|q\bigr)&:\text{AW}
  \end{array}\right.\nonumber\\
  &\phantom{\check{P}_n(x;\boldsymbol{\lambda})}=\left\{\begin{array}{ll}
  {\displaystyle
  (a_1+a_2,a_1+a_3,a_1+a_4)_n}\\[2pt]
  {\displaystyle
  \quad\times
  {}_4F_3\Bigl(
  \genfrac{}{}{0pt}{}{-n,\,n+b_1-1,\,a_1+ix,\,a_1-ix}
  {a_1+a_2,\,a_1+a_3,\,a_1+a_4}\Bigm|1\Bigr)
  }&:\text{W}\\[8pt]
  {\displaystyle
  a_1^{-n}(a_1a_2,a_1a_3,a_1a_4;q)_n}\\[2pt]
  {\displaystyle
  \quad\times
  {}_4\phi_3\Bigl(\genfrac{}{}{0pt}{}{q^{-n},\,b_4q^{n-1},\,
  a_1e^{ix},\,a_1e^{-ix}}{a_1a_2,\,a_1a_3,\,a_1a_4}\!\!\Bigm|\!q\,;q\Bigr)
  }&:\text{AW}
  \end{array}\right.
  \label{Pn=W,AW}\\[2pt]
  &\phi_0(x;\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  \sqrt{\bigl(\Gamma(2ix)\Gamma(-2ix)\bigr)^{-1}\prod_{j=1}^4
  \Gamma(a_j+ix)\Gamma(a_j-ix)}&:\text{W}\\[4pt]
  \sqrt{(e^{2ix},e^{-2ix};q)_{\infty}
  \prod_{j=1}^4(a_je^{ix},a_je^{-ix};q)_{\infty}^{-1}}
  &:\text{AW}
  \end{array}\right.,
  \label{phi0=W,AW}
\end{align}
where $W_n(\eta;a_1,a_2,a_3,a_4)$ and $p_n(\eta;a_1,a_2,a_3,a_4|q)$ are
the Wilson and Askey-Wilson polynomials \cite{kls}, respectively.
Note that
\begin{align}
  &\phi_0(x;\boldsymbol{\lambda}+\boldsymbol{\delta})
  =\varphi(x)\sqrt{V(x+i\tfrac{\gamma}{2};\boldsymbol{\lambda})}\,
  \phi_0(x+i\tfrac{\gamma}{2};\boldsymbol{\lambda}),
  \label{phi0(l+d)}\\
  &V(x;\boldsymbol{\lambda}+\boldsymbol{\delta})
  =\kappa^{-1}\frac{\varphi(x-i\gamma)}{\varphi(x)}
  V(x-i\tfrac{\gamma}{2};\boldsymbol{\lambda}).
  \label{varphiprop3}
\end{align}

\subsection{Virtual states}
\label{sec:idQM_vs}

The virtual state wavefunctions are obtained by using the discrete symmetries
of the Hamiltonian \cite{os27}.
In the following we restrict the parameters as follows:
\begin{equation}
  a_1,a_2\in\mathbb{R}\ \ \text{or}\ \ a_2^*=a_1\ ;\quad
  a_3,a_4\in\mathbb{R}\ \ \text{or}\ \ a_4^*=a_3.
\end{equation}
We define two twist operations $\mathfrak{t}$, type $\text{\Romannumeral{1}}$ and type $\text{\Romannumeral{2}}$,
\begin{equation}
  \mathfrak{t}^{\text{\Romannumeral{1}}}(\boldsymbol{\lambda})
  \stackrel{\text{def}}{=}(1-\lambda_1,1-\lambda_2,\lambda_3,\lambda_4),\quad
  \mathfrak{t}^{\text{\Romannumeral{2}}}(\boldsymbol{\lambda})
  \stackrel{\text{def}}{=}(\lambda_1,\lambda_2,1-\lambda_3,1-\lambda_4),
  \label{idQM:twist}
\end{equation}
which are involutions $\mathfrak{t}^2=\text{id}$ and satisfy
\begin{equation}
  \mathfrak{t}(\boldsymbol{\lambda})+u\boldsymbol{\delta}
  =\mathfrak{t}(\boldsymbol{\lambda}+u\tilde{\boldsymbol{\delta}})
  \ \ (\forall u\in\mathbb{R}),\quad
  \tilde{\boldsymbol{\delta}}^{\text{\Romannumeral{1}}}\stackrel{\text{def}}{=}(-\tfrac12,-\tfrac12,\tfrac12,\tfrac12),
  \ \ \tilde{\boldsymbol{\delta}}^{\text{\Romannumeral{2}}}\stackrel{\text{def}}{=}(\tfrac12,\tfrac12,-\tfrac12,-\tfrac12).
\end{equation}
By using these two types of twist operations, two types of virtual state
wavefunctions, $\tilde{\phi}_{\text{v}}^{\text{\Romannumeral{1}}}(x;\boldsymbol{\lambda})$ and
$\tilde{\phi}_{\text{v}}^{\text{\Romannumeral{2}}}(x;\boldsymbol{\lambda})$, are introduced:
\begin{align}
  \tilde{\phi}_{\text{v}}(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \phi_{\text{v}}\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})\bigr)
  =\tilde{\phi}_0(x;\boldsymbol{\lambda})
  \check{\xi}_{\text{v}}(x;\boldsymbol{\lambda}),\quad
  \tilde{\phi}_0(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \phi_0\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})\bigr),\nonumber\\
  \check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \xi_{\text{v}}\bigl(\eta(x);\boldsymbol{\lambda}\bigr)\stackrel{\text{def}}{=}
  \check{P}_{\text{v}}\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})\bigr)
  =P_{\text{v}}\bigl(\eta(x);\mathfrak{t}(\boldsymbol{\lambda})\bigr).
  \label{idQM:phitv}
\end{align}
The virtual state polynomials $\xi_{\text{v}}(\eta;\boldsymbol{\lambda})$
(two types $\xi^{\text{\Romannumeral{1}}}_{\text{v}}$ and $\xi^{\text{\Romannumeral{2}}}_{\text{v}}$)
are polynomials of degree $\text{v}$ in $\eta$.
The virtual state wavefunctions satisfy the Schr\"{o}dinger equation
$\mathcal{H}(\boldsymbol{\lambda})\tilde{\phi}_{\text{v}}(x;\boldsymbol{\lambda})
=\tilde{\mathcal{E}}_{\text{v}}(\boldsymbol{\lambda})
\tilde{\phi}_{\text{v}}(x;\boldsymbol{\lambda})$
with the virtual energies $\tilde{\mathcal{E}}_{\text{v}}(\boldsymbol{\lambda})$,
\begin{align}
  &\tilde{\mathcal{E}}^{\text{\Romannumeral{1}}}_{\text{v}}(\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  -(a_1+a_2-\text{v}-1)(a_3+a_4+\text{v})&:\text{W}\\
  -(1-a_1a_2q^{-\text{v}-1})(1-a_3a_4q^{\text{v}})&:\text{AW}
  \end{array}\right.,\nonumber\\
  &\tilde{\mathcal{E}}^{\text{\Romannumeral{2}}}_{\text{v}}(\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  -(a_3+a_4-\text{v}-1)(a_1+a_2+\text{v})&:\text{W}\\
  -(1-a_3a_4q^{-\text{v}-1})(1-a_1a_2q^{\text{v}})&:\text{AW}
  \end{array}\right..
  \label{idQM:Etv}
\end{align}
Other data $\alpha(\boldsymbol{\lambda})$ and $\alpha'(\boldsymbol{\lambda})$
(see \cite{os27} for details) are
\begin{align}
  &\alpha^{\text{\Romannumeral{1}}}(\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  1&:\text{W}\\
  a_1a_2q^{-1}&:\text{AW}
  \end{array}\right.,\quad
  \alpha^{\prime\,\text{\Romannumeral{1}}}(\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  -(a_1+a_2-1)(a_3+a_4)&:\text{W}\\
  -(1-a_1a_2q^{-1})(1-a_3a_4)&:\text{AW}
  \end{array}\right.,\nonumber\\
  &\alpha^{\text{\Romannumeral{2}}}(\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  1&:\text{W}\\
  a_3a_4q^{-1}&:\text{AW}
  \end{array}\right.,\quad
  \alpha^{\prime\,\text{\Romannumeral{2}}}(\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  -(a_3+a_4-1)(a_1+a_2)&:\text{W}\\
  -(1-a_3a_4q^{-1})(1-a_1a_2)&:\text{AW}
  \end{array}\right..
\end{align}
In order to obtain well-defined quantum systems, we have to restrict the
degree $\text{v}$ and the parameters $a_i$ so that the conditions,
$\tilde{\mathcal{E}}_{\text{v}}(\boldsymbol{\lambda})<0$, $\alpha(\boldsymbol{\lambda})>0$,
$\alpha'(\boldsymbol{\lambda})<0$, no zeros of $\xi_{\text{v}}(\eta;\boldsymbol{\lambda})$
in the domain etc., are satisfied, see \cite{os27}.
In this paper, however, we consider algebraic properties only and we
do not bother about the ranges of $\text{v}$ and $a_i$.

The functions $\nu(x;\boldsymbol{\lambda})$ (two types $\nu^{\text{\Romannumeral{1}}}$ and $\nu^{\text{\Romannumeral{2}}}$)
are defined as the ratio of $\phi_0$ and $\tilde{\phi}_0$,
\begin{equation}
  \nu(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \frac{\phi_0(x;\boldsymbol{\lambda})}{\tilde{\phi}_0(x;\boldsymbol{\lambda})}.
  \label{idQM:nu}
\end{equation}
The functions $r_j(x^{(M)}_j;\boldsymbol{\lambda},M)$ (two types $r^{\text{\Romannumeral{1}}}_j$ and
$r^{\text{\Romannumeral{2}}}_j$) are defined as the ratio of $\nu$'s,
\begin{equation}
  r_j(x^{(M)}_j;\boldsymbol{\lambda},M)\stackrel{\text{def}}{=}
  \frac{\nu(x^{(M)}_j;\boldsymbol{\lambda})}
  {\nu\bigl(x;\boldsymbol{\lambda}+(M-1)\tilde{\boldsymbol{\delta}}\bigr)}
  \ \ (1\leq j\leq M),
  \label{idQM:rj}
\end{equation}
where $x^{(M)}_j\stackrel{\text{def}}{=} x+i(\frac{M+1}{2}-j)\gamma$.
Their explicit forms are
\begin{align}
  r^{\text{\Romannumeral{1}}}_j(x^{(M)}_j;\boldsymbol{\lambda},M)
  &=\alpha^{\text{\Romannumeral{1}}}\bigl(\boldsymbol{\lambda}
  +(M-1)\tilde{\boldsymbol{\delta}}^{\text{\Romannumeral{1}}}\bigr)^{-\frac12(M-1)}
  \kappa^{\frac12(M-1)^2-(j-1)(M-j)}\\
  &\quad\times\left\{
  \begin{array}{ll}
  {\displaystyle
  \prod_{k=1,2}(a_k-\tfrac{M-1}{2}+ix)_{j-1}(a_k-\tfrac{M-1}{2}-ix)_{M-j}
  }&:\text{W}\\
  {\displaystyle
  e^{ix(M+1-2j)}\prod_{k=1,2}(a_kq^{-\frac{M-1}{2}}e^{ix};q)_{j-1}
  (a_kq^{-\frac{M-1}{2}}e^{-ix};q)_{M-j}
  }&:\text{AW}
  \end{array}\right.,\nonumber\\
  r^{\text{\Romannumeral{2}}}_j(x^{(M)}_j;\boldsymbol{\lambda},M)
  &=\alpha^{\text{\Romannumeral{2}}}\bigl(\boldsymbol{\lambda}
  +(M-1)\tilde{\boldsymbol{\delta}}^{\text{\Romannumeral{2}}}\bigr)^{-\frac12(M-1)}
  \kappa^{\frac12(M-1)^2-(j-1)(M-j)}\\
  &\quad\times\left\{
  \begin{array}{ll}
  {\displaystyle
  \prod_{k=3,4}(a_k-\tfrac{M-1}{2}+ix)_{j-1}(a_k-\tfrac{M-1}{2}-ix)_{M-j}
  }&:\text{W}\\
  {\displaystyle
  e^{ix(M+1-2j)}\prod_{k=3,4}(a_kq^{-\frac{M-1}{2}}e^{ix};q)_{j-1}
  (a_kq^{-\frac{M-1}{2}}e^{-ix};q)_{M-j}
  }&:\text{AW}
  \end{array}\right..\nonumber
\end{align}
The auxiliary function $\varphi_M(x)$ \cite{gos} is defined by
\begin{align}
  \varphi_M(x)&\stackrel{\text{def}}{=}
  \varphi(x)^{[\frac{M}{2}]}\prod_{k=1}^{M-2}
  \bigl(\varphi(x-i\tfrac{k}{2}\gamma)\varphi(x+i\tfrac{k}{2}\gamma)
  \bigr)^{[\frac{M-k}{2}]}\nonumber\\
  &=\prod_{1\leq j<k\leq M}
  \frac{\eta(x^{(M)}_j)-\eta(x^{(M)}_k)}
  {\varphi(i\frac{j}{2}\gamma)}
  \times\left\{
  \begin{array}{ll}
  1&:\text{W}\\
  (-2)^{\frac12M(M-1)}&:\text{AW}
  \end{array}\right.,
  \label{idQM:varphiM}
\end{align}
and $\varphi_0(x)=\varphi_1(x)=1$.

For $-M\leq k\leq M$ and $k\equiv M$ ($\text{mod}\,2$),
we can show that
\begin{equation}
  \frac{\nu(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})}{\nu(x;\boldsymbol{\lambda}+k\boldsymbol{\delta})}
  =\prod_{l=0}^{\frac{M-k}{2}-1}
  \check{U}\bigl(x;\boldsymbol{\lambda}+(k+2l)\boldsymbol{\delta}\bigr),
  \label{nuMk}
\end{equation}
where the functions $\check{U}(x;\boldsymbol{\lambda})$ (two types $\check{U}^{\text{\Romannumeral{1}}}$
and $\check{U}^{\text{\Romannumeral{2}}}$) are
\begin{align}
  \check{U}^{\text{\Romannumeral{1}}}(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}\left\{
  \begin{array}{ll}
  \prod_{k=1}^2(a_k+ix)(a_k-ix)&:\text{W}\\[2pt]
  (a_1a_2)^{-1}\prod_{k=1}^2(1-a_ke^{ix})(1-a_ke^{-ix})&:\text{AW}
  \end{array}\right.,\nonumber\\
  \check{U}^{\text{\Romannumeral{2}}}(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}\left\{
  \begin{array}{ll}
  \prod_{k=3}^4(a_k+ix)(a_k-ix)&:\text{W}\\[2pt]
  (a_3a_4)^{-1}\prod_{k=3}^4(1-a_ke^{ix})(1-a_ke^{-ix})&:\text{AW}
  \end{array}\right..
\end{align}
In fact, these functions $\check{U}(x;\boldsymbol{\lambda})$ are polynomials
$U(\eta;\boldsymbol{\lambda})$ in $\eta$ as defined by
\begin{align}
  &U\bigl(\eta(x);\boldsymbol{\lambda}\bigr)\stackrel{\text{def}}{=}\check{U}(x;\boldsymbol{\lambda}),\\
  &U^{\text{\Romannumeral{1}}}(\eta;\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  \prod_{k=1}^2(a_k^2+\eta)&:\text{W}\\[2pt]
  \prod_{k=1}^2(a_k+a_k^{-1}-2\eta)&:\text{AW}
  \end{array}\right.\!\!,
  \ \ U^{\text{\Romannumeral{2}}}(\eta;\boldsymbol{\lambda})=\left\{
  \begin{array}{ll}
  \prod_{k=3}^4(a_k^2+\eta)&:\text{W}\\[2pt]
  \prod_{k=3}^4(a_k+a_k^{-1}-2\eta)&:\text{AW}
  \end{array}\right..
  \nonumber
\end{align}
Since $\nu(x;\boldsymbol{\lambda})$ satisfies
$\nu(x;\boldsymbol{\lambda}+k\tilde{\boldsymbol{\delta}})=\nu(x;\boldsymbol{\lambda}-k\boldsymbol{\delta})$
($k\in\mathbb{Z}$), \eqref{nuMk} gives
\begin{equation}
  \frac{\nu(x;\boldsymbol{\lambda}+M\tilde{\boldsymbol{\delta}})}
  {\nu(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})}
  =\prod_{l=0}^{M-1}
  \check{U}\bigl(x;\boldsymbol{\lambda}+(-M+2l)\boldsymbol{\delta}\bigr)^{-1}.
  \label{nuMtM}
\end{equation}

In the following we adopt the convention,
\begin{equation}
  \check{P}_n(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=} 0\ \ (n\in\mathbb{Z}_{<0}),\quad
  \check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=} 0
  \ \ (\text{v}\in\mathbb{Z}_{<0}).
\end{equation}

\subsection{Shape invariance}
\label{sec:idQM_si}

The original systems in \S\,\ref{sec:idQM_org} are shape invariant
\cite{os13_1,os13_2}
and they satisfy the relation
\begin{equation}
  \mathcal{A}(\boldsymbol{\lambda})\mathcal{A}(\boldsymbol{\lambda})^{\dagger}
  =\kappa\mathcal{A}(\boldsymbol{\lambda}+\boldsymbol{\delta})^{\dagger}
  \mathcal{A}(\boldsymbol{\lambda}+\boldsymbol{\delta})+\mathcal{E}_1(\boldsymbol{\lambda}),
\end{equation}
which is a sufficient condition for exact solvability.
As a consequence of the shape invariance, the action of the operators
$\mathcal{A}(\boldsymbol{\lambda})$ and $\mathcal{A}(\boldsymbol{\lambda})^{\dagger}$
on the eigenfunctions is
\begin{equation}
  \mathcal{A}(\boldsymbol{\lambda})\phi_n(x;\boldsymbol{\lambda})
  =f_n(\boldsymbol{\lambda})\phi_{n-1}(x;\boldsymbol{\lambda}+\boldsymbol{\delta}),\quad
  \mathcal{A}(\boldsymbol{\lambda})^{\dagger}
  \phi_{n-1}(x;\boldsymbol{\lambda}+\boldsymbol{\delta})
  =b_{n-1}(\boldsymbol{\lambda})\phi_n(x;\boldsymbol{\lambda}),
  \label{idQM:Aphi=}
\end{equation}
where the factors of the energy eigenvalue, $f_n(\boldsymbol{\lambda})$ and
$b_{n-1}(\boldsymbol{\lambda})$,
$\mathcal{E}_n(\boldsymbol{\lambda})=f_n(\boldsymbol{\lambda})b_{n-1}(\boldsymbol{\lambda})$,
are given by
\begin{equation}
  f_n(\boldsymbol{\lambda})\stackrel{\text{def}}{=}\left\{
  \begin{array}{ll}
  -n(n+b_1-1)&:\text{W}\\
  q^{\frac{n}{2}}(q^{-n}-1)(1-b_4q^{n-1})&:\text{AW}
  \end{array}\right.,
  \quad
  b_{n-1}(\boldsymbol{\lambda})\stackrel{\text{def}}{=}\left\{
  \begin{array}{ll}
  -1&:\text{W}\\
  q^{-\frac{n}{2}}&:\text{AW}
  \end{array}\right..
\end{equation}
The relations \eqref{idQM:Aphi=} are equivalent to the forward and
backward shift relations of the orthogonal polynomials $P_n$ \cite{kls},
\begin{equation}
  \mathcal{F}(\boldsymbol{\lambda})\check{P}_n(x;\boldsymbol{\lambda})
  =f_n(\boldsymbol{\lambda})\check{P}_{n-1}(x;\boldsymbol{\lambda}+\boldsymbol{\delta}),\quad
  \mathcal{B}(\boldsymbol{\lambda})\check{P}_{n-1}(x;\boldsymbol{\lambda}+\boldsymbol{\delta})
  =b_{n-1}(\boldsymbol{\lambda})\check{P}_n(x;\boldsymbol{\lambda}),
  \label{idQM:FP=}
\end{equation}
where the forward and backward shift operators $\mathcal{F}(\boldsymbol{\lambda})$
and $\mathcal{B}(\boldsymbol{\lambda})$ are defined by
\begin{align}
  &\mathcal{F}(\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \phi_0(x;\boldsymbol{\lambda}+\boldsymbol{\delta})^{-1}\circ
  \mathcal{A}(\boldsymbol{\lambda})\circ\phi_0(x;\boldsymbol{\lambda})
  =i\varphi(x)^{-1}(e^{\frac{\gamma}{2}p}-e^{-\frac{\gamma}{2}p}),
  \label{idQM:Fdef}\\
  &\mathcal{B}(\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \phi_0(x;\boldsymbol{\lambda})^{-1}\circ
  \mathcal{A}(\boldsymbol{\lambda})^{\dagger}
  \circ\phi_0(x;\boldsymbol{\lambda}+\boldsymbol{\delta})
  =-i\bigl(V(x;\boldsymbol{\lambda})e^{\frac{\gamma}{2}p}
  -V^*(x;\boldsymbol{\lambda})e^{-\frac{\gamma}{2}p}\bigr)\varphi(x).
  \label{idQM:Bdef}
\end{align}
Note that the forward shift operator $\mathcal{F}(\boldsymbol{\lambda})$ is parameter
independent.
The second order difference operator $\widetilde{\mathcal{H}}(\boldsymbol{\lambda})$
acting on the polynomial eigenfunctions is square root free. It is defined by
\begin{align}
  &\widetilde{\mathcal{H}}(\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \phi_0(x;\boldsymbol{\lambda})^{-1}\circ\mathcal{H}(\boldsymbol{\lambda})
  \circ\phi_0(x;\boldsymbol{\lambda})
  =\mathcal{B}(\boldsymbol{\lambda})\mathcal{F}(\boldsymbol{\lambda})\nonumber\\
  &\phantom{\widetilde{\mathcal{H}}_{\ell}(\boldsymbol{\lambda})}
  =V(x;\boldsymbol{\lambda})(e^{\gamma p}-1)
  +V^*(x;\boldsymbol{\lambda})(e^{-\gamma p}-1),\\
  &\widetilde{\mathcal{H}}(\boldsymbol{\lambda})\check{P}_n(x;\boldsymbol{\lambda})
  =\mathcal{E}_n(\boldsymbol{\lambda})\check{P}_n(x;\boldsymbol{\lambda}).
  \label{idQM:HtP=EP}
\end{align}

The action of the operators
$\mathcal{A}(\boldsymbol{\lambda})$ and $\mathcal{A}(\boldsymbol{\lambda})^{\dagger}$
on the virtual state wavefunctions $\tilde{\phi}_{\text{v}}(x;\boldsymbol{\lambda})$
($\tilde{\phi}^{\text{\Romannumeral{1}}}_{\text{v}}$ and $\tilde{\phi}^{\text{\Romannumeral{2}}}_{\text{v}}$) is
\begin{equation}
  \mathcal{A}(\boldsymbol{\lambda})\tilde{\phi}_{\text{v}}(x;\boldsymbol{\lambda})
  =\tilde{f}_{\text{v}}(\boldsymbol{\lambda})
  \tilde{\phi}_{\text{v}}(x;\boldsymbol{\lambda}+\boldsymbol{\delta}),
  \ \ \mathcal{A}(\boldsymbol{\lambda})^{\dagger}
  \tilde{\phi}_{\text{v}}(x;\boldsymbol{\lambda}+\boldsymbol{\delta})
  =\tilde{b}_{\text{v}}(\boldsymbol{\lambda})
  \tilde{\phi}_{\text{v}}(x;\boldsymbol{\lambda}),
  \label{idQM:Atphi=}
\end{equation}
where the factors of the virtual energy eigenvalue,
$\tilde{f}_{\text{v}}(\boldsymbol{\lambda})$ and $\tilde{b}_{\text{v}}(\boldsymbol{\lambda})$,
$\tilde{\mathcal{E}}_{\text{v}}(\boldsymbol{\lambda})
=\tilde{f}_{\text{v}}(\boldsymbol{\lambda})\tilde{b}_{\text{v}}(\boldsymbol{\lambda})$,
are given by
\begin{align}
  \tilde{f}^{\text{\Romannumeral{1}}}_{\text{v}}(\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \alpha^{\text{\Romannumeral{1}}}(\boldsymbol{\lambda})^{-\frac12}\times\left\{
  \begin{array}{ll}
  a_1+a_2-\text{v}-1&:\text{W}\\
  -q^{\frac{\text{v}}{2}}(1-a_1a_2q^{-\text{v}-1})&:\text{AW}
  \end{array}\right.,\nonumber\\
  \tilde{b}^{\text{\Romannumeral{1}}}_{\text{v}}(\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \alpha^{\text{\Romannumeral{1}}}(\boldsymbol{\lambda})^{\frac12}\times\left\{
  \begin{array}{ll}
  -(a_3+a_4+\text{v})&:\text{W}\\
  q^{-\frac{\text{v}}{2}}(1-a_3a_4q^{\text{v}})&:\text{AW}
  \end{array}\right.,\nonumber\\
  \tilde{f}^{\text{\Romannumeral{2}}}_{\text{v}}(\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \alpha^{\text{\Romannumeral{2}}}(\boldsymbol{\lambda})^{-\frac12}\times\left\{
  \begin{array}{ll}
  a_3+a_4-\text{v}-1&:\text{W}\\
  -q^{\frac{\text{v}}{2}}(1-a_3a_4q^{-\text{v}-1})&:\text{AW}
  \end{array}\right.,\\
  \tilde{b}^{\text{\Romannumeral{2}}}_{\text{v}}(\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \alpha^{\text{\Romannumeral{2}}}(\boldsymbol{\lambda})^{\frac12}\times\left\{
  \begin{array}{ll}
  -(a_1+a_2+\text{v})&:\text{W}\\
  q^{-\frac{\text{v}}{2}}(1-a_1a_2q^{\text{v}})&:\text{AW}
  \end{array}\right..
  \nonumber
\end{align}
The relations \eqref{idQM:Atphi=} are equivalent to the forward and
backward shift relations of
$(\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda})$
$\stackrel{\text{def}}{=}\nu(x;\boldsymbol{\lambda})^{-1}\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})$
($\nu^{\text{\Romannumeral{1}}\,-1}\check{\xi}^{\text{\Romannumeral{1}}}_{\text{v}}$ and
$\nu^{\text{\Romannumeral{2}}\,-1}\check{\xi}^{\text{\Romannumeral{2}}}_{\text{v}}$),
\begin{equation}
  \mathcal{F}(\boldsymbol{\lambda})(\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda})
  =\tilde{f}_{\text{v}}(\boldsymbol{\lambda})
  (\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda}+\boldsymbol{\delta}),
  \ \ \mathcal{B}(\boldsymbol{\lambda})
  (\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda}+\boldsymbol{\delta})
  =\tilde{b}_{\text{v}}(\boldsymbol{\lambda})
  (\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda}).
  \label{idQM:Fnuxi=}
\end{equation}
Explicitly they give square root free relations: forward and backward shift
relations of the virtual polynomials $\xi_{\text{v}}$ ($\xi^{\text{\Romannumeral{1}}}_{\text{v}}$
and $\xi^{\text{\Romannumeral{2}}}_{\text{v}}$) \cite{os27},
\begin{align}
  i\varphi(x)^{-1}\bigl(
  v_1^*(x;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})e^{\frac{\gamma}{2}p}
  -v_1(x;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})e^{-\frac{\gamma}{2}p}\bigr)
  \check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})
  &=\alpha(\boldsymbol{\lambda})^{\frac12}\tilde{f}_{\text{v}}(\boldsymbol{\lambda})
  \check{\xi}_{\text{v}}(x;\boldsymbol{\lambda}+\boldsymbol{\delta}),\nonumber\\
  i\varphi(x)^{-1}\bigl(
  v_2(x;\boldsymbol{\lambda})e^{\frac{\gamma}{2}p}
  -v_2^*(x;\boldsymbol{\lambda})e^{-\frac{\gamma}{2}p}\bigr)
  \check{\xi}_{\text{v}}(x;\boldsymbol{\lambda}+\boldsymbol{\delta})
  &=\alpha(\boldsymbol{\lambda})^{-\frac12}\tilde{b}_{\text{v}}(\boldsymbol{\lambda})
  \check{\xi}_{\text{v}}(x;\boldsymbol{\lambda}),
  \label{idQM:Fxi}
\end{align}
where the functions $v_1(x;\boldsymbol{\lambda})$ and $v_2(x;\boldsymbol{\lambda})$ are
\begin{equation}
  v_1(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=}\left\{
  \begin{array}{ll}
  \!\prod_{j=1}^2(a_j+ix)&\!\!:\text{W}\\[2pt]
  \!e^{-ix}\prod_{j=1}^2(1-a_je^{ix})&\!\!:\text{AW}
  \end{array}\right.\!\!\!,
  \ \ v_2(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=}\left\{
  \begin{array}{ll}
  \!\prod_{j=3}^4(a_j+ix)&\!\!:\text{W}\\[2pt]
  \!e^{-ix}\prod_{j=3}^4(1-a_je^{ix})&\!\!:\text{AW}
  \end{array}\right.\!\!\!.\!\!\!
  \label{v1v2}
\end{equation}

\subsection{Multi-indexed orthogonal polynomials}
\label{sec:idQM_miop}

The original systems in \S\,\ref{sec:idQM_org} are iso-spectrally deformed
by applying multiple Darboux transformations with virtual state wavefunctions
$\{\tilde{\phi}_{\text{v}}(x)\}$ as seed solutions.
The virtual states wavefunctions are labeled by ($\text{v}$, $\text{t}$):
$\text{v}$ is the degree of the polynomial part $\xi_{\text{v}}$ and
$\text{t}$ stands for type $\text{\Romannumeral{1}}$ or $\text{\Romannumeral{2}}$.
We take $M$ virtual state wavefunctions specified by the multi-index set
$\mathcal{D}=\{(d_1,\text{t}_1),\ldots,(d_M,\text{t}_M)\}$
(or shortly $\mathcal{D}=\{d_1,\ldots,d_M\}$) (ordered set),
$M=M_{\text{\Romannumeral{1}}}+M_{\text{\Romannumeral{2}}}$, $M_{\text{\Romannumeral{1}}}\stackrel{\text{def}}{=}\#\{d_j|(d_j,\text{\Romannumeral{1}})\in\mathcal{D}\}$,
$M_{\text{\Romannumeral{2}}}\stackrel{\text{def}}{=}\#\{d_j|(d_j,\text{\Romannumeral{2}})\in\mathcal{D}\}$.
The deformed Hamiltonian is denoted as
$\mathcal{H}_{\mathcal{D}}(\boldsymbol{\lambda})$.

The multi-indexed orthogonal polynomials
$\{\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})\}$ are the main parts
of the eigenfunctions $\{\phi_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})\}$ of the
deformed system $\mathcal{H}_{\mathcal{D}}(\boldsymbol{\lambda})
\phi_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})=\mathcal{E}_n(\boldsymbol{\lambda})
\phi_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})$ \cite{os27}:
\begin{align}
  \phi_{\mathcal{D}\,n}(x;\boldsymbol{\lambda})&=\left(
  \frac{\sqrt{\prod_{j=0}^{M-1}
  V\bigl(x+i(\frac{M}{2}-j)\gamma;\boldsymbol{\lambda}\bigr)
  V^*\bigl(x-i(\frac{M}{2}-j)\gamma;\boldsymbol{\lambda}\bigr)}}
  {\text{W}_{\gamma}[\tilde{\phi}_{d_1},\ldots,\tilde{\phi}_{d_M}]
  (x-i\frac{\gamma}{2};\boldsymbol{\lambda})
  \text{W}_{\gamma}[\tilde{\phi}_{d_1},\ldots,\tilde{\phi}_{d_M}]
  (x+i\frac{\gamma}{2};\boldsymbol{\lambda})}\right)^{\frac12}\nonumber\\
  &\quad\times
  \text{W}_{\gamma}[\tilde{\phi}_{d_1},\ldots,\tilde{\phi}_{d_M},\phi_n]
  (x;\boldsymbol{\lambda})\nonumber\\
  &=\alpha^{\text{\Romannumeral{1}}}(\boldsymbol{\lambda}^{[M_{\text{\Romannumeral{1}}},M_{\text{\Romannumeral{2}}}]})^{\frac12M_{\text{\Romannumeral{1}}}}
  \alpha^{\text{\Romannumeral{2}}}(\boldsymbol{\lambda}^{[M_{\text{\Romannumeral{1}}},M_{\text{\Romannumeral{2}}}]})^{\frac12M_{\text{\Romannumeral{2}}}}
  \kappa^{-\frac14M_{\text{\Romannumeral{1}}}(M_{\text{\Romannumeral{1}}}+1)-\frac14M_{\text{\Romannumeral{2}}}(M_{\text{\Romannumeral{2}}}+1)
  +\frac52M_{\text{\Romannumeral{1}}}M_{\text{\Romannumeral{2}}}}\nonumber\\
  &\quad\times
  \psi_{\mathcal{D}}(x;\boldsymbol{\lambda})
  \check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda}),
  \label{idQM:phiDn}\\
  \psi_{\mathcal{D}}(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  \frac{\phi_0(x;\boldsymbol{\lambda}^{[M_{\text{\Romannumeral{1}}},M_{\text{\Romannumeral{2}}}]})}
  {\sqrt{\check{\Xi}_{\mathcal{D}}(x-i\frac{\gamma}{2};\boldsymbol{\lambda})
  \check{\Xi}_{\mathcal{D}}(x+i\frac{\gamma}{2};\boldsymbol{\lambda})}},\quad
  \boldsymbol{\lambda}^{[M_{\text{\Romannumeral{1}}},M_{\text{\Romannumeral{2}}}]}\stackrel{\text{def}}{=}
  \boldsymbol{\lambda}+M_{\text{\Romannumeral{1}}}\tilde{\boldsymbol{\delta}}^{\text{\Romannumeral{1}}}+M_{\text{\Romannumeral{2}}}\tilde{\boldsymbol{\delta}}^{\text{\Romannumeral{2}}},
  \label{idQM:psiD}\\
  \check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  P_{\mathcal{D},n}\bigl(\eta(x);\boldsymbol{\lambda}\bigr),\quad
  \check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \Xi_{\mathcal{D}}\bigl(\eta(x);\boldsymbol{\lambda}\bigr).
  \label{idQM:PDn}
\end{align}
Here $\text{W}_{\gamma}[f_1,\ldots,f_n](x)$ is the Casoratian
\eqref{idQM:Wdef}, and $\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})$ and
$\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ are the main polynomial part of
the Casoratians
$\text{W}_{\gamma}[\tilde{\phi}_{d_1},\ldots,\tilde{\phi}_{d_M}]
(x;\boldsymbol{\lambda})$ and
$\text{W}_{\gamma}[\tilde{\phi}_{d_1},\ldots,\tilde{\phi}_{d_M},\phi_n]
(x;\boldsymbol{\lambda})$, respectively.
The denominator polynomial $\Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})$ and the
multi-indexed orthogonal polynomial $P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})$
are polynomials in $\eta$ and their degrees are generically
$\ell_{\mathcal{D}}$ and $\ell_{\mathcal{D}}+n$, respectively.
Here $\ell_{\mathcal{D}}$ is
$\ell_{\mathcal{D}}\stackrel{\text{def}}{=}\sum_{j=1}^{M}d_j-\tfrac12M(M-1)+2M_{\text{\Romannumeral{1}}}M_{\text{\Romannumeral{2}}}$.

In \cite{os27}, the denominator polynomials and the multi-indexed orthogonal
polynomials for the Wilson and Askey-Wilson types are expressed as
\begin{align}
  &\quad\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})\nonumber\\
  &\stackrel{\text{def}}{=} i^{\frac12M(M-1)}\left|
  \begin{array}{lll}
  \vec{X}^{(M)}_{d_1}&\cdots&\vec{X}^{(M)}_{d_M}
  \end{array}\right|
  \times\varphi_M(x)^{-1}\kappa^{-\frac12M_{\text{\Romannumeral{1}}}M_{\text{\Romannumeral{2}}}(M-2)}
  \label{idQM:cXiDdef}\\
  &\quad\times\biggl(
  \prod_{l=0}^{M_{\text{\Romannumeral{1}}}-1}\check{U}^{\text{\Romannumeral{2}}}\bigl(x;\boldsymbol{\lambda}
  +(1-M+2l)\boldsymbol{\delta}\bigr)^{M_{\text{\Romannumeral{1}}}-1-l}\cdot
  \prod_{l=0}^{M_{\text{\Romannumeral{2}}}-1}\check{U}^{\text{\Romannumeral{1}}}\bigl(x;\boldsymbol{\lambda}
  +(1-M+2l)\boldsymbol{\delta}\bigr)^{M_{\text{\Romannumeral{2}}}-1-l}\biggr)^{-1},\nonumber\\
  &\quad\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})\nonumber\\
  &\stackrel{\text{def}}{=} i^{\frac12M(M+1)}\left|
  \begin{array}{llll}
  \vec{X}^{(M+1)}_{d_1}&\cdots&\vec{X}^{(M+1)}_{d_M}&\vec{Z}^{(M+1)}_n
  \end{array}\right|
  \times\varphi_{M+1}(x)^{-1}\kappa^{-\frac12M_{\text{\Romannumeral{1}}}M_{\text{\Romannumeral{2}}}(M+2)}
  \label{idQM:cPDndef}\\
  &\quad\times\biggl(
  \prod_{l=0}^{M_{\text{\Romannumeral{1}}}-1}\check{U}^{\text{\Romannumeral{2}}}\bigl(x;\boldsymbol{\lambda}
  +(-M+2l)\boldsymbol{\delta}\bigr)^{M_{\text{\Romannumeral{1}}}-l}\cdot
  \prod_{l=0}^{M_{\text{\Romannumeral{2}}}-1}\check{U}^{\text{\Romannumeral{1}}}\bigl(x;\boldsymbol{\lambda}
  +(-M+2l)\boldsymbol{\delta}\bigr)^{M_{\text{\Romannumeral{2}}}-l}\biggr)^{-1},
  \nonumber
\end{align}
where $\vec{X}^{(M)}_{\text{v}}$ ($\vec{X}^{(M)\,\text{\Romannumeral{1}}}_{\text{v}}$ or
$\vec{X}^{(M)\,\text{\Romannumeral{2}}}_{\text{v}}$) and $\vec{Z}^{(M)}_n$ are
\begin{align}
  &\bigl(\vec{X}^{(M)\,\text{t}}_{\text{v}}\bigr)_j
  =r_j^{\bar{\text{t}}}(x^{(M)}_j;\boldsymbol{\lambda},M)
  \check{\xi}^{\text{t}}_{\text{v}}(x^{(M)}_j;\boldsymbol{\lambda}),\quad
  \bar{\text{t}}\stackrel{\text{def}}{=}\left\{
  \begin{array}{ll}
  \text{\Romannumeral{2}}&:\text{t}=\text{\Romannumeral{1}}\\
  \text{\Romannumeral{1}}&:\text{t}=\text{\Romannumeral{2}}
  \end{array}\right.,\nonumber\\
  &\bigl(\vec{Z}^{(M)}_n\bigr)_j
  =r^{\text{\Romannumeral{2}}}_j(x^{(M)}_j;\boldsymbol{\lambda},M)r^{\text{\Romannumeral{1}}}_j(x^{(M)}_j;\boldsymbol{\lambda},M)
  \check{P}_n(x^{(M)}_j;\boldsymbol{\lambda})
  \ \ (1\leq j\leq M).
\end{align}
By using the properties of the determinant, \eqref{dWformula1},
\eqref{idQM:rj}, \eqref{nuMk} and \eqref{nuMtM}, we can rewrite these as
follows:
\begin{align}
  &\quad\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})\nonumber\\
  &=\text{W}_{\gamma}[\nu^{-1}\check{\xi}_{d_1},\ldots,
  \nu^{-1}\check{\xi}_{d_M}](x;\boldsymbol{\lambda})
  \times\varphi_M(x)^{-1}\kappa^{-\frac12M_{\text{\Romannumeral{1}}}M_{\text{\Romannumeral{2}}}(M-2)}\nonumber\\
  &\quad\times
  \nu^{\text{\Romannumeral{1}}}\bigl(x;\boldsymbol{\lambda}+(M-1)\boldsymbol{\delta}\bigr)^{M_{\text{\Romannumeral{1}}}}
  \nu^{\text{\Romannumeral{2}}}\bigl(x;\boldsymbol{\lambda}+(M-1)\boldsymbol{\delta}\bigr)^{M_{\text{\Romannumeral{2}}}}
  \label{idQM:cXiD2}\\
  &\quad\times\biggl(
  \prod_{m=1}^{M_{\text{\Romannumeral{1}}}}\check{U}^{\text{\Romannumeral{1}}}\bigl(x;\boldsymbol{\lambda}
  +(M-1-2m)\boldsymbol{\delta}\bigr)^{M_{\text{\Romannumeral{1}}}-m}\cdot
  \prod_{m=1}^{M_{\text{\Romannumeral{2}}}}\check{U}^{\text{\Romannumeral{2}}}\bigl(x;\boldsymbol{\lambda}
  +(M-1-2m)\boldsymbol{\delta}\bigr)^{M_{\text{\Romannumeral{2}}}-m}\biggr)^{-1},\nonumber\\
  &\quad\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})\nonumber\\
  &=\text{W}_{\gamma}[\nu^{-1}\check{\xi}_{d_1},\ldots,
  \nu^{-1}\check{\xi}_{d_M},\check{P}_n](x;\boldsymbol{\lambda})
  \times\varphi_{M+1}(x)^{-1}\kappa^{-\frac12M_{\text{\Romannumeral{1}}}M_{\text{\Romannumeral{2}}}(M+2)}\nonumber\\
  &\quad\times
  \nu^{\text{\Romannumeral{1}}}(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})^{M_{\text{\Romannumeral{1}}}}
  \nu^{\text{\Romannumeral{2}}}(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})^{M_{\text{\Romannumeral{2}}}}
  \label{idQM:cPDn2}\\
  &\quad\times\biggl(
  \prod_{m=1}^{M_{\text{\Romannumeral{1}}}}\check{U}^{\text{\Romannumeral{1}}}\bigl(x;\boldsymbol{\lambda}
  +(M-2m)\boldsymbol{\delta}\bigr)^{M_{\text{\Romannumeral{1}}}-m}\cdot
  \prod_{m=1}^{M_{\text{\Romannumeral{2}}}}\check{U}^{\text{\Romannumeral{2}}}\bigl(x;\boldsymbol{\lambda}
  +(M-2m)\boldsymbol{\delta}\bigr)^{M_{\text{\Romannumeral{2}}}-m}\biggr)^{-1}.
  \nonumber
\end{align}
For the special situations when the index set $\mathcal{D}$ consists type
$\text{\Romannumeral{1}}$ indices only (or type $\text{\Romannumeral{2}}$ indices only),
the above Casoratians can be simplified as
\begin{align*}
  \text{W}_{\gamma}[\nu^{-1}\check{\xi}_{d_1},\ldots,
  \nu^{-1}\check{\xi}_{d_M}](x;\boldsymbol{\lambda})
  &=\prod_{j=1}^M\nu(x^{(M)}_j;\boldsymbol{\lambda})^{-1}\cdot
  \text{W}_{\gamma}[\check{\xi}_{d_1},\ldots,\check{\xi}_{d_M}]
  (x;\boldsymbol{\lambda}),\\
  \text{W}_{\gamma}[\nu^{-1}\check{\xi}_{d_1},\ldots,
  \nu^{-1}\check{\xi}_{d_M},\check{P}_n](x;\boldsymbol{\lambda})
  &=\prod_{j=1}^{M+1}\nu(x^{(M+1)}_j;\boldsymbol{\lambda})^{-1}\cdot
  \text{W}_{\gamma}[\check{\xi}_{d_1},\ldots,\check{\xi}_{d_M},
  \nu\check{P}_n](x;\boldsymbol{\lambda}),
\end{align*}
by \eqref{dWformula1}.
However, when $\mathcal{D}$ contains both $\text{\Romannumeral{1}}$ and $\text{\Romannumeral{2}}$ types,
such rewriting is impossible due to
$\nu^{\text{\Romannumeral{1}}}(x;\boldsymbol{\lambda})\neq\nu^{\text{\Romannumeral{2}}}(x;\boldsymbol{\lambda})$.

In the next subsection we will rewrite \eqref{idQM:cXiD2}--\eqref{idQM:cPDn2}
by using the identities obtained from the shape invariance and the
Casoratian.

\subsection{New determinant expressions}
\label{sec:idQM_exp}

\subsubsection{Casoratian}
\label{sec:idQM_Caso}

The Casorati determinant of a set of $n$ functions $\{f_j(x)\}$ is defined by
\begin{equation}
  \text{W}_{\gamma}[f_1,f_2,\ldots,f_n](x)
  \stackrel{\text{def}}{=} i^{\frac12n(n-1)}
  \det\Bigl(f_k\bigl(x^{(n)}_j\bigr)\Bigr)_{1\leq j,k\leq n},\quad
  x_j^{(n)}\stackrel{\text{def}}{=} x+i(\tfrac{n+1}{2}-j)\gamma,
  \label{idQM:Wdef}
\end{equation}
(for $n=0$, we set $\text{W}_{\gamma}[\cdot](x)=1$),
which satisfies
\begin{equation}
  \text{W}_{\gamma}[gf_1,gf_2,\ldots,gf_n]
  =\prod_{j=1}^ng\bigl(x^{(n)}_j\bigr)\cdot
  \text{W}_{\gamma}[f_1,f_2,\ldots,f_n](x).
  \label{dWformula1}
\end{equation}

The following determinant formula holds for any analytic functions $q_j(x)$,
$r_j(x)$, $q'_j(x)$ and $r'_j(x)$ ($j=1,2,\ldots$):
\begin{align}
  &\quad i^{\frac12n(n-1)}\det(a_{j,k})_{1\leq j,k\leq n}\nonumber\\
  &=\prod_{m=1}^{[\frac{n}{2}]}D_{n+1-2m}(x)\cdot
  \prod_{m=0}^{n-2}\prod_{l=0}^{[\frac{m-1}{2}]}
  r_{n-1-m}\bigl(x+i(\tfrac{m}{2}-l)\gamma\bigr)
  q_{n-1-m}\bigl(x-i(\tfrac{m}{2}-l)\gamma\bigr)\nonumber\\
  &\quad\times\text{W}_{\gamma}[f_1,\ldots,f_n](x),
  \label{idQM:Wgid}\\
  &\left\{
  \begin{array}{ll}
  {\displaystyle
  a_{j,k}=\laprod{l=1}{n+1-2j}\hat{D}_l\cdot f_k(x)}
  &(j=1,\ldots,[\tfrac{n+1}{2}])\\
  {\displaystyle
  a_{j,k}=\hat{D}'_{2j-n-1}\laprod{l=1}{2j-n-2}\hat{D}_l\cdot f_k(x)}
  &(j=[\tfrac{n+3}{2}],\ldots,n)
  \end{array}\right.,\ \ (k=1,\ldots,n),
  \label{idQM:Wgid,ajk}
\end{align}
where operators $\hat{D}_j$, $\hat{D}'_j$ and functions $D_j(x)$
($j=1,2,\ldots$) are
\begin{align}
  &\hat{D}_j=q_j(x)e^{\frac{\gamma}{2}p}+r_j(x)e^{-\frac{\gamma}{2}p},\quad
  \hat{D}'_j=q'_j(x)e^{\frac{\gamma}{2}p}+r'_j(x)e^{-\frac{\gamma}{2}p},\\
  &D_j(x)=r_j(x)q'_j(x)-r'_j(x)q_j(x).
  \label{Dj}
\end{align}
This formula can be proven by using the properties of the determinant
(row properties) and induction in $n$ (even $n$ $\to$ odd $n+1$,
odd $n$ $\to$ even $n+1$).

By taking $\mathcal{F}$ or $\mathcal{B}$ as $\hat{D}_j$ and $\hat{D}'_j$
and using shape invariance properties \eqref{idQM:FP=} and \eqref{idQM:Fnuxi=},
the Casoratians $\text{W}_{\gamma}[\nu^{-1}\check{\xi}_{d_1},\ldots,
\nu^{-1}\check{\xi}_{d_M}](x;\boldsymbol{\lambda})$
and $\text{W}_{\gamma}[\nu^{-1}\check{\xi}_{d_1},\ldots,
\nu^{-1}\check{\xi}_{d_M},\check{P}_n](x;\boldsymbol{\lambda})$
in \eqref{idQM:cXiD2}--\eqref{idQM:cPDn2}
can be rewritten in various ways.
We consider two typical cases:
\begin{align}
  \text{A}:\ \ &
  \hat{D}_j=\mathcal{F}\bigl(\boldsymbol{\lambda}+(j-1)\boldsymbol{\delta}\bigr),
  \ \ \hat{D}'_j=\mathcal{B}\bigl(\boldsymbol{\lambda}+(j-2)\boldsymbol{\delta}\bigr)
  \ \ (j=1,2,\ldots),
  \label{idQM:A}\\
  \text{B}:\ \ &
  \hat{D}_{2l-1}=\mathcal{F}(\boldsymbol{\lambda}),
  \ \hat{D}_{2l}=\mathcal{B}(\boldsymbol{\lambda}),
  \ \hat{D}'_{2l-1}=\mathcal{B}(\boldsymbol{\lambda}-\boldsymbol{\delta}),
  \ \hat{D}'_{2l}=\mathcal{F}(\boldsymbol{\lambda}+\boldsymbol{\delta})
  \ \ (l=1,2,\ldots),
  \label{idQM:B}
\end{align}
which correspond to the cases A and B in \cite{os37}.
For the function $D_j(x)$ \eqref{Dj}, we define a function
$\check{S}(x;\boldsymbol{\lambda})$ and a polynomial $S(\eta;\boldsymbol{\lambda})$:
\begin{align}
  \check{S}(x;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}
  -i\bigl(V(x;\boldsymbol{\lambda})\varphi(x-i\tfrac{\gamma}{2})
  -V^*(x;\boldsymbol{\lambda})\varphi(x+i\tfrac{\gamma}{2})\bigr),\nonumber\\
  S(\eta;\boldsymbol{\lambda})&\stackrel{\text{def}}{=}\left\{
  \begin{array}{ll}
  b_1\eta-b_3&:\text{W}\\
  q^{-\frac12}\bigl(2(1-b_4)\eta+b_3-b_1\bigr)&:\text{AW}
  \end{array}\right.,\quad
  \check{S}(x;\boldsymbol{\lambda})=S\bigl(\eta(x);\boldsymbol{\lambda}\bigr),
\end{align}
where $b_3$ is $b_3\stackrel{\text{def}}{=} a_1a_2a_3+a_1a_2a_4+a_1a_3a_4+a_2a_3a_4$.

\subsubsection{case A}
\label{sec:idQM_caseA}

Firstly we consider case A \eqref{idQM:A}.
The functions $q_j(x)$, $r_j(x)$, $q'_j(x)$, $r'_j(x)$ and $D_j(x)$ are
\begin{align}
  q_j(x)&=i\varphi(x)^{-1},\quad r_j(x)=-i\varphi(x)^{-1},\nonumber\\
  q'_j(x)&=-iV\bigl(x;\boldsymbol{\lambda}+(j-2)\boldsymbol{\delta}\bigr)
  \varphi(x-i\tfrac{\gamma}{2}),\quad
  r'_j(x)=iV^*\bigl(x;\boldsymbol{\lambda}+(j-2)\boldsymbol{\delta}\bigr)
  \varphi(x+i\tfrac{\gamma}{2}),\\
  D_j(x)&=-i\varphi(x)^{-1}
  \check{S}\bigl(x;\boldsymbol{\lambda}+(j-2)\boldsymbol{\delta}\bigr).
  \nonumber
\end{align}
The shape invariance properties \eqref{idQM:FP=} and \eqref{idQM:Fnuxi=}
give ($k=0,1,\ldots$)
\begin{align}
  &\laprod{l=1}{k}\hat{D}_l\cdot\check{P}_n(x;\boldsymbol{\lambda})
  =\prod_{m=0}^{k-1}f_{n-m}(\boldsymbol{\lambda}+m\boldsymbol{\delta})\cdot
  \check{P}_{n-k}(x;\boldsymbol{\lambda}+k\boldsymbol{\delta}),\nonumber\\[-2pt]
  &\laprod{l=1}{k}\hat{D}_l\cdot
  (\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda})
  =\prod_{m=0}^{k-1}\tilde{f}_{\text{v}}(\boldsymbol{\lambda}+m\boldsymbol{\delta})\cdot
  (\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda}+k\boldsymbol{\delta}),
  \label{idQM:siA}\\[-2pt]
  &\hat{D}'_{k+1}\laprod{l=1}{k}\hat{D}_l\cdot\check{P}_n(x;\boldsymbol{\lambda})
  =b_{n-k}\bigl(\boldsymbol{\lambda}+(k-1)\boldsymbol{\delta}\bigr)
  \prod_{m=0}^{k-1}f_{n-m}(\boldsymbol{\lambda}+m\boldsymbol{\delta})\cdot
  \check{P}_{n+1-k}\bigl(x;\boldsymbol{\lambda}+(k-1)\boldsymbol{\delta}\bigr),\nonumber\\[-2pt]
  &\hat{D}'_{k+1}\laprod{l=1}{k}\hat{D}_l\cdot
  (\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda})
  =\tilde{b}_{\text{v}}\bigl(\boldsymbol{\lambda}+(k-1)\boldsymbol{\delta}\bigr)
  \prod_{m=0}^{k-1}\tilde{f}_{\text{v}}(\boldsymbol{\lambda}+m\boldsymbol{\delta})\cdot
  (\nu^{-1}\check{\xi}_{\text{v}})\bigl(x;\boldsymbol{\lambda}+(k-1)\boldsymbol{\delta}\bigr).
  \nonumber
\end{align}

We rewrite \eqref{idQM:cXiD2}--\eqref{idQM:cPDn2} in the following way.
(\romannumeral1) By using \eqref{idQM:Wgid}--\eqref{idQM:Wgid,ajk},
rewrite the Casoratians in \eqref{idQM:cXiD2}--\eqref{idQM:cPDn2} as
determinants $\det(a_{j,k})_{1\leq j,k\leq M'}$ ($M'=M$ for
$\check{\Xi}_{\mathcal{D}}$ and $M'=M+1$ for $\check{P}_{\mathcal{D},n}$).
(\romannumeral2) Rewrite each matrix element $a_{j,k}$ by \eqref{idQM:siA}.
(\romannumeral3) Multiply the $k$-th column ($1\leq k\leq M$) of
$\det(a_{j,k})$ by $\nu\bigl(x;\boldsymbol{\lambda}+(M'-1)\boldsymbol{\delta}\bigr)$
($\nu^{\text{\Romannumeral{1}}}$ or $\nu^{\text{\Romannumeral{2}}}$ depending on $d_k$), and divide the determinant by
$\nu^{\text{\Romannumeral{1}}}\bigl(x;\boldsymbol{\lambda}+(M'-1)\boldsymbol{\delta}\bigr)^{M_{\text{\Romannumeral{1}}}}
\nu^{\text{\Romannumeral{2}}}\bigl(x;\boldsymbol{\lambda}+(M'-1)\boldsymbol{\delta}\bigr)^{M_{\text{\Romannumeral{2}}}}$.
(\romannumeral4) Rewrite
$\nu\bigl(x;\boldsymbol{\lambda}+(M'-1)\boldsymbol{\delta}\bigr)/
\nu(x;\boldsymbol{\lambda}+l\boldsymbol{\delta})$ in $a_{j,k}$ by \eqref{nuMk}.
Then $\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})$ and
$\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ are expressed in terms of
$\check{\xi}_{d_k}(x;\boldsymbol{\lambda}')$,
$\check{P}_{n'}(x;\boldsymbol{\lambda}')$, $\check{U}(x;\boldsymbol{\lambda}')$,
$\check{S}(x;\boldsymbol{\lambda}')$ and $\varphi(x')$.
Straightforward calculation shows that the factors $\varphi(x')$ are
canceled out.
Thus $\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})$ and
$\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ are expressed in terms of
$\eta(x)$, namely
$\Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})$ and
$P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})$ are expressed without explicit
$x$-dependence.
Their final forms are as follows.
The denominator polynomial $\Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})$ is
\begin{align}
  &\quad\Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})\nonumber\\
  &=(-1)^{\frac16(M-1)M(M+1)}\kappa^{-\frac12M_{\text{\Romannumeral{1}}}M_{\text{\Romannumeral{2}}}(M-2)}
  \det(a_{j,k})_{1\leq j,k\leq M}\nonumber\\
  &\quad\times\biggl(
  \prod_{l=0}^{M_{\text{\Romannumeral{1}}}-1}U^{\text{\Romannumeral{1}}}\bigl(\eta;\boldsymbol{\lambda}
  +(M_{\text{\Romannumeral{2}}}-M_{\text{\Romannumeral{1}}}-1+2l)\boldsymbol{\delta}\bigr)^l\cdot
  \prod_{l=0}^{M_{\text{\Romannumeral{2}}}-1}U^{\text{\Romannumeral{2}}}\bigl(\eta;\boldsymbol{\lambda}
  +(M_{\text{\Romannumeral{1}}}-M_{\text{\Romannumeral{2}}}-1+2l)\boldsymbol{\delta}\bigr)^l\nonumber\\
  &\qquad\times\prod_{m=1}^{[\frac{M}{2}]}
  S\bigl(\eta;\boldsymbol{\lambda}+(M-1-2m)\boldsymbol{\delta}\bigr)\biggr)^{-1},
  \label{idQM:A:XiD}
\end{align}
where $a_{j,k}$ are
\begin{align}
  a_{j,k}&=\prod_{m=0}^{M-2j}\tilde{f}_{d_k}(\boldsymbol{\lambda}+m\boldsymbol{\delta})\cdot
  \prod_{l=1}^{j-1}U\bigl(\eta;\boldsymbol{\lambda}+(M-1-2l)\boldsymbol{\delta}\bigr)\cdot
  \xi_{d_k}\bigl(\eta;\boldsymbol{\lambda}+(M+1-2j)\boldsymbol{\delta}\bigr)\nonumber\\
  &\qquad(j=1,\ldots,[\tfrac{M+1}{2}];k=1,\ldots,M),\nonumber\\
  a_{j,k}&=\tilde{b}_{d_k}\bigl(\boldsymbol{\lambda}+(2j-M-3)\boldsymbol{\delta}\bigr)
  \prod_{m=0}^{2j-M-3}\tilde{f}_{d_k}(\boldsymbol{\lambda}+m\boldsymbol{\delta})\cdot
  \prod_{l=1}^{M+1-j}U\bigl(\eta;\boldsymbol{\lambda}+(M-1-2l)\boldsymbol{\delta}\bigr)\nonumber\\
  &\quad\times
  \xi_{d_k}\bigl(\eta;\boldsymbol{\lambda}+(2j-M-3)\boldsymbol{\delta}\bigr)
  \ \ (j=[\tfrac{M+3}{2}],\ldots,M;k=1,\ldots,M).
\end{align}
The multi-indexed orthogonal polynomial
$P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})$ is
\begin{align}
  &\quad P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})\nonumber\\
  &=(-1)^{\frac16M(M+1)(M+2)}\kappa^{-\frac12M_{\text{\Romannumeral{1}}}M_{\text{\Romannumeral{2}}}(M+2)}
  \det(a_{j,k})_{1\leq j,k\leq M+1}\nonumber\\
  &\quad\times\biggl(
  \prod_{l=0}^{M_{\text{\Romannumeral{1}}}-1}U^{\text{\Romannumeral{1}}}\bigl(\eta;\boldsymbol{\lambda}
  +(M_{\text{\Romannumeral{2}}}-M_{\text{\Romannumeral{1}}}+2l)\boldsymbol{\delta}\bigr)^l\cdot
  \prod_{l=0}^{M_{\text{\Romannumeral{2}}}-1}U^{\text{\Romannumeral{2}}}\bigl(\eta;\boldsymbol{\lambda}
  +(M_{\text{\Romannumeral{1}}}-M_{\text{\Romannumeral{2}}}+2l)\boldsymbol{\delta}\bigr)^l\nonumber\\
  &\qquad\times\prod_{m=1}^{[\frac{M+1}{2}]}
  S\bigl(\eta;\boldsymbol{\lambda}+(M-2m)\boldsymbol{\delta}\bigr)\biggr)^{-1},
  \label{idQM:A:PDn}
\end{align}
where $a_{j,k}$ are
\begin{align}
  a_{j,k}&=\prod_{m=0}^{M+1-2j}\tilde{f}_{d_k}(\boldsymbol{\lambda}+m\boldsymbol{\delta})\cdot
  \prod_{l=1}^{j-1}U\bigl(\eta;\boldsymbol{\lambda}+(M-2l)\boldsymbol{\delta}\bigr)\cdot
  \xi_{d_k}\bigl(\eta;\boldsymbol{\lambda}+(M+2-2j)\boldsymbol{\delta}\bigr)\nonumber\\
  &\qquad(j=1,\ldots,[\tfrac{M+2}{2}];k=1,\ldots,M),\nonumber\\
  a_{j,k}&=\tilde{b}_{d_k}\bigl(\boldsymbol{\lambda}+(2j-M-4)\boldsymbol{\delta}\bigr)
  \prod_{m=0}^{2j-M-4}\tilde{f}_{d_k}(\boldsymbol{\lambda}+m\boldsymbol{\delta})\cdot
  \prod_{l=1}^{M+2-j}U\bigl(\eta;\boldsymbol{\lambda}+(M-2l)\boldsymbol{\delta}\bigr)\nonumber\\
  &\quad\times
  \xi_{d_k}\bigl(\eta;\boldsymbol{\lambda}+(2j-M-4)\boldsymbol{\delta}\bigr)
  \ \ (j=[\tfrac{M+4}{2}],\ldots,M+1;k=1,\ldots,M),\nonumber\\
  a_{j,M+1}&=\prod_{m=0}^{M+1-2j}f_{n-m}(\boldsymbol{\lambda}+m\boldsymbol{\delta})\cdot
  P_{n-M-2+2j}\bigl(\eta;\boldsymbol{\lambda}+(M+2-2j)\boldsymbol{\delta}\bigr)
  \ \ (j=1,\ldots,[\tfrac{M+2}{2}]),\nonumber\\
  a_{j,M+1}&=b_{n-2j+M+3}\bigl(\boldsymbol{\lambda}+(2j-M-4)\boldsymbol{\delta}\bigr)
  \prod_{m=0}^{2j-M-4}f_{n-m}(\boldsymbol{\lambda}+m\boldsymbol{\delta})\\
  &\quad\times
  P_{n-2j+M+4}\bigl(\eta;\boldsymbol{\lambda}+(2j-M-4)\boldsymbol{\delta}\bigr)
  \ \ (j=[\tfrac{M+4}{2}],\ldots,M+1).
  \nonumber
\end{align}

\subsubsection{case B}
\label{sec:idQM_caseB}

Secondly we consider case B \eqref{idQM:B}.
The functions $q_j(x)$, $r_j(x)$, $q'_j(x)$, $r'_j(x)$ and $D_j(x)$ are
($l=1,2,\ldots$)
\begin{align}
  q_{2l-1}(x)&=i\varphi(x)^{-1},\quad r_{2l-1}(x)=-i\varphi(x)^{-1},\nonumber\\
  q_{2l}(x)&=-iV(x;\boldsymbol{\lambda})\varphi(x-i\tfrac{\gamma}{2}),\quad
  r_{2l}(x)=iV^*(x;\boldsymbol{\lambda})\varphi(x+i\tfrac{\gamma}{2}),\nonumber\\
  q'_{2l-1}(x)&=-iV(x;\boldsymbol{\lambda}-\boldsymbol{\delta})
  \varphi(x-i\tfrac{\gamma}{2}),\quad
  r'_{2l-1}=iV^*(x;\boldsymbol{\lambda}-\boldsymbol{\delta})
  \varphi(x+i\tfrac{\gamma}{2}),\\
  q'_{2l}(x)&=i\varphi(x)^{-1},\quad r'_{2l}(x)=-i\varphi(x)^{-1},\nonumber\\
  D_{2l-1}(x)&=-i\varphi(x)^{-1}\check{S}(x;\boldsymbol{\lambda}-\boldsymbol{\delta}),\quad
  D_{2l}(x)=i\varphi(x)^{-1}\check{S}(x;\boldsymbol{\lambda}).
  \nonumber
\end{align}
The shape invariance properties \eqref{idQM:FP=} and \eqref{idQM:Fnuxi=}
give ($l=1,2,\ldots$)
\begin{align}
  &\laprod{m=1}{2l-2}\hat{D}_m\cdot\check{P}_n(x;\boldsymbol{\lambda})
  =\mathcal{E}_n(\boldsymbol{\lambda})^{l-1}
  \check{P}_n(x;\boldsymbol{\lambda}),\nonumber\\[-2pt]
  &\laprod{m=1}{2l-1}\hat{D}_m\cdot\check{P}_n(x;\boldsymbol{\lambda})
  =f_n(\boldsymbol{\lambda})\mathcal{E}_n(\boldsymbol{\lambda})^{l-1}
  \check{P}_{n-1}(x;\boldsymbol{\lambda}+\boldsymbol{\delta}),\nonumber\\[-2pt]
  &\laprod{m=1}{2l-2}\hat{D}_m\cdot
  (\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda})
  =\tilde{\mathcal{E}}_{\text{v}}(\boldsymbol{\lambda})^{l-1}
  (\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda}),\nonumber\\[-2pt]
  &\laprod{m=1}{2l-1}\hat{D}_m\cdot
  (\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda})
  =\tilde{f}_{\text{v}}(\boldsymbol{\lambda})
  \tilde{\mathcal{E}}_{\text{v}}(\boldsymbol{\lambda})^{l-1}
  (\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda}+\boldsymbol{\delta}),\\[-2pt]
  &\hat{D}'_{2l-1}\laprod{m=1}{2l-2}\hat{D}_m\cdot\check{P}_n(x;\boldsymbol{\lambda})
  =b_n(\boldsymbol{\lambda}-\boldsymbol{\delta})\mathcal{E}_n(\boldsymbol{\lambda})^{l-1}
  \check{P}_{n+1}(x;\boldsymbol{\lambda}-\boldsymbol{\delta}),\nonumber\\[-2pt]
  &\hat{D}'_{2l}\laprod{m=1}{2l-1}\hat{D}_m\cdot\check{P}_n(x;\boldsymbol{\lambda})
  =f_{n-1}(\boldsymbol{\lambda}+\boldsymbol{\delta})f_n(\boldsymbol{\lambda})
  \mathcal{E}_n(\boldsymbol{\lambda})^{l-1}
  \check{P}_{n-2}(x;\boldsymbol{\lambda}+2\boldsymbol{\delta}),\nonumber\\[-2pt]
  &\hat{D}'_{2l-1}\laprod{m=1}{2l-2}\hat{D}_m\cdot
  (\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda})
  =\tilde{b}_{\text{v}}(\boldsymbol{\lambda}-\boldsymbol{\delta})
  \tilde{\mathcal{E}}_{\text{v}}(\boldsymbol{\lambda})^{l-1}
  (\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda}-\boldsymbol{\delta}),\nonumber\\[-2pt]
  &\hat{D}'_{2l}\laprod{m=1}{2l-1}\hat{D}_m\cdot
  (\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda})
  =\tilde{f}_{\text{v}}(\boldsymbol{\lambda}+\boldsymbol{\delta})
  \tilde{f}_{\text{v}}(\boldsymbol{\lambda})
  \tilde{\mathcal{E}}_{\text{v}}(\boldsymbol{\lambda})^{l-1}
  (\nu^{-1}\check{\xi}_{\text{v}})(x;\boldsymbol{\lambda}+2\boldsymbol{\delta}).
  \nonumber
\end{align}

Like case A, we rewrite \eqref{idQM:cXiD2}--\eqref{idQM:cPDn2} by the steps
(\romannumeral1)--(\romannumeral4).
By using \eqref{varphiprop3} and
\begin{equation}
  V(x;\boldsymbol{\lambda})V^*(x;\boldsymbol{\lambda})
  =\alpha^{\text{\Romannumeral{1}}}(\boldsymbol{\lambda})\alpha^{\text{\Romannumeral{2}}}(\boldsymbol{\lambda})\kappa^{-1}
  \frac{\check{U}^{\text{\Romannumeral{1}}}(x;\boldsymbol{\lambda})\check{U}^{\text{\Romannumeral{2}}}(x;\boldsymbol{\lambda})}
  {\varphi(x-i\frac{\gamma}{2})\varphi(x)^2\varphi(x+i\frac{\gamma}{2})}
  \label{VV*},
\end{equation}
$\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})$ and
$\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ are expressed in terms of
$\check{\xi}_{d_k}(x;\boldsymbol{\lambda}')$,
$\check{P}_{n'}(x;\boldsymbol{\lambda}')$, $\check{U}(x;\boldsymbol{\lambda}')$,
$\check{S}(x;\boldsymbol{\lambda}')$ and $\varphi(x')$.
Straightforward calculation shows that the factors $\varphi(x')$ are
canceled out.
Thus $\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})$ and
$\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ are expressed in terms of
$\eta(x)$, namely
$\Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})$ and
$P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})$ are expressed without explicit
$x$-dependence.
Their final forms are as follows.
The denominator polynomial $\Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})$ is
\begin{align}
  &\quad\Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})\nonumber\\
  &=(-1)^{[\frac{M}{2}]}\kappa^{-\frac12M_{\text{\Romannumeral{1}}}M_{\text{\Romannumeral{2}}}(M-2)}
  \bigl(\alpha^{\text{\Romannumeral{1}}}(\boldsymbol{\lambda})\alpha^{\text{\Romannumeral{2}}}(\boldsymbol{\lambda})\kappa^{-1}
  \bigr)^{-\frac12[\frac{M-2}{2}][\frac{M}{2}]}
  \det(a_{j,k})_{1\leq j,k\leq M}\nonumber\\
  &\quad\times\biggl(
  \prod_{m=1}^{M_{\text{\Romannumeral{1}}}}U^{\text{\Romannumeral{1}}}\bigl(\eta;\boldsymbol{\lambda}
  +(M-1-2m)\boldsymbol{\delta}\bigr)^{M_{\text{\Romannumeral{1}}}-m}\cdot
  \prod_{m=1}^{M_{\text{\Romannumeral{2}}}}U^{\text{\Romannumeral{2}}}\bigl(\eta;\boldsymbol{\lambda}
  +(M-1-2m)\boldsymbol{\delta}\bigr)^{M_{\text{\Romannumeral{2}}}-m}\nonumber\\
  &\qquad\times
  \prod_{m=0}^{[\frac{M-2}{2}]}U^{\text{\Romannumeral{1}}}\bigl(\eta;\boldsymbol{\lambda}
  +(M-1-2m)\boldsymbol{\delta}\bigr)^m\cdot
  \prod_{m=0}^{[\frac{M-2}{2}]}U^{\text{\Romannumeral{2}}}\bigl(\eta;\boldsymbol{\lambda}
  +(M-1-2m)\boldsymbol{\delta}\bigr)^m\nonumber\\
  &\qquad\times
  S\bigl(\eta;\boldsymbol{\lambda}-\tfrac12(1+(-1)^M)\boldsymbol{\delta}
  \bigr)^{[\frac{M}{2}]}\biggr)^{-1},
  \label{idQM:B:XiD}
\end{align}
where $a_{j,k}$ are
\begin{align}
  a_{j,k}&=\tilde{\mathcal{E}}_{d_k}(\boldsymbol{\lambda})^{\frac{M+1}{2}-j}
  \prod_{l=0}^{\frac{M-3}{2}}U(\eta;\boldsymbol{\lambda}+2l\boldsymbol{\delta})\cdot
  \xi_{d_k}(\eta;\boldsymbol{\lambda})
  \ \ (j=1,\ldots,\tfrac{M+1}{2};k=1,\ldots,M),\nonumber\\
  a_{j,k}&=\tilde{f}_{d_k}(\boldsymbol{\lambda}+\boldsymbol{\delta})
  \tilde{f}_{d_k}(\boldsymbol{\lambda})
  \tilde{\mathcal{E}}_{d_k}(\boldsymbol{\lambda})^{j-\frac{M+3}{2}}
  \prod_{l=1}^{\frac{M-3}{2}}U(\eta;\boldsymbol{\lambda}+2l\boldsymbol{\delta})\cdot
  \xi_{d_k}(\eta;\boldsymbol{\lambda}+2\boldsymbol{\delta})\\
  &\qquad(j=\tfrac{M+3}{2},\ldots,M;k=1,\ldots,M),
  \nonumber
\end{align}
for odd $M$ and
\begin{align}
  a_{j,k}&=\tilde{f}_{d_k}(\boldsymbol{\lambda})
  \tilde{\mathcal{E}}_{d_k}(\boldsymbol{\lambda})^{\frac{M}{2}-j}
  \prod_{l=1}^{\frac{M-2}{2}}
  U\bigl(\eta;\boldsymbol{\lambda}+(2l-1)\boldsymbol{\delta}\bigr)\cdot
  \xi_{d_k}(\eta;\boldsymbol{\lambda}+\boldsymbol{\delta}),\nonumber\\
  &\qquad(j=1,\ldots,\tfrac{M}{2};k=1,\ldots,M),\\
  a_{j,k}&=\tilde{b}_{d_k}(\boldsymbol{\lambda}-\boldsymbol{\delta})
  \tilde{\mathcal{E}}_{d_k}(\boldsymbol{\lambda})^{j-\frac{M+2}{2}}
  \prod_{l=0}^{\frac{M-2}{2}}
  U\bigl(\eta;\boldsymbol{\lambda}+(2l-1)\boldsymbol{\delta}\bigr)\cdot
  \xi_{d_k}(\eta;\boldsymbol{\lambda}-\boldsymbol{\delta})
  \ \ (j=\tfrac{M+2}{2},\ldots,M),
  \nonumber
\end{align}
for even $M$.
The multi-indexed orthogonal polynomial
$P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})$ is
\begin{align}
  &\quad P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})\nonumber\\
  &=(-1)^{[\frac{M+1}{2}]}\kappa^{-\frac12M_{\text{\Romannumeral{1}}}M_{\text{\Romannumeral{2}}}(M+2)}
  \bigl(\alpha^{\text{\Romannumeral{1}}}(\boldsymbol{\lambda})\alpha^{\text{\Romannumeral{2}}}(\boldsymbol{\lambda})\kappa^{-1}
  \bigr)^{-\frac12[\frac{M-1}{2}][\frac{M+1}{2}]}
  \det(a_{j,k})_{1\leq j,k\leq M+1}\nonumber\\
  &\quad\times\biggl(
  \prod_{m=1}^{M_{\text{\Romannumeral{1}}}}U^{\text{\Romannumeral{1}}}\bigl(\eta;\boldsymbol{\lambda}
  +(M-2m)\boldsymbol{\delta}\bigr)^{M_{\text{\Romannumeral{1}}}-m}\cdot
  \prod_{m=1}^{M_{\text{\Romannumeral{2}}}}U^{\text{\Romannumeral{2}}}\bigl(\eta;\boldsymbol{\lambda}
  +(M-2m)\boldsymbol{\delta}\bigr)^{M_{\text{\Romannumeral{2}}}-m}\nonumber\\
  &\qquad\times
  \prod_{m=0}^{[\frac{M-1}{2}]}U^{\text{\Romannumeral{1}}}\bigl(\eta;\boldsymbol{\lambda}
  +(M-2m)\boldsymbol{\delta}\bigr)^m\cdot
  \prod_{m=0}^{[\frac{M-1}{2}]}U^{\text{\Romannumeral{2}}}\bigl(\eta;\boldsymbol{\lambda}
  +(M-2m)\boldsymbol{\delta}\bigr)^m\nonumber\\
  &\qquad\times
  S\bigl(\eta;\boldsymbol{\lambda}-\tfrac12(1-(-1)^M)\boldsymbol{\delta}
  \bigr)^{[\frac{M+1}{2}]}\biggr)^{-1},
  \label{idQM:B:PDn}
\end{align}
where $a_{j,k}$ are
\begin{align}
  a_{j,k}&=\tilde{\mathcal{E}}_{d_k}(\boldsymbol{\lambda})^{\frac{M+2}{2}-j}
  \prod_{l=0}^{\frac{M}{2}-1}U(\eta;\boldsymbol{\lambda}+2l\boldsymbol{\delta})\cdot
  \xi_{d_k}(\eta;\boldsymbol{\lambda})
  \ \ (j=1,\ldots,\tfrac{M+2}{2};k=1,\ldots,M),\nonumber\\
  a_{j,k}&=\tilde{f}_{d_k}(\boldsymbol{\lambda}+\boldsymbol{\delta})
  \tilde{f}_{d_k}(\boldsymbol{\lambda})
  \tilde{\mathcal{E}}_{d_k}(\boldsymbol{\lambda})^{j-\frac{M+4}{2}}
  \prod_{l=1}^{\frac{M}{2}-1}U(\eta;\boldsymbol{\lambda}+2l\boldsymbol{\delta})\cdot
  \xi_{d_k}(\eta;\boldsymbol{\lambda}+2\boldsymbol{\delta})\nonumber\\
  &\qquad(j=\tfrac{M+4}{2},\ldots,M+1;k=1,\ldots,M),\\
  a_{j,M+1}&=\mathcal{E}_n(\boldsymbol{\lambda})^{\frac{M+2}{2}-j}
  P_n(\eta;\boldsymbol{\lambda})
  \ \ (j=1,\ldots,\tfrac{M+2}{2}),\nonumber\\
  a_{j,M+1}&=f_{n-1}(\boldsymbol{\lambda}+\boldsymbol{\delta})f_n(\boldsymbol{\lambda})
  \mathcal{E}_n(\boldsymbol{\lambda})^{j-\frac{M+4}{2}}
  P_{n-2}(\eta;\boldsymbol{\lambda}+2\boldsymbol{\delta})
  \ \ (j=\tfrac{M+4}{2},\ldots,M+1),
  \nonumber
\end{align}
for even $M$ and
\begin{align}
  a_{j,k}&=\tilde{f}_{d_k}(\boldsymbol{\lambda})
  \tilde{\mathcal{E}}_{d_k}(\boldsymbol{\lambda})^{\frac{M+1}{2}-j}
  \prod_{l=1}^{\frac{M-1}{2}}
  U\bigl(\eta;\boldsymbol{\lambda}+(2l-1)\boldsymbol{\delta}\bigr)\cdot
  \xi_{d_k}(\eta;\boldsymbol{\lambda}+\boldsymbol{\delta}),\nonumber\\
  &\qquad(j=1,\ldots,\tfrac{M+1}{2};k=1,\ldots,M),\nonumber\\
  a_{j,k}&=\tilde{b}_{d_k}(\boldsymbol{\lambda}-\boldsymbol{\delta})
  \tilde{\mathcal{E}}_{d_k}(\boldsymbol{\lambda})^{j-\frac{M+3}{2}}
  \prod_{l=0}^{\frac{M-1}{2}}
  U\bigl(\eta;\boldsymbol{\lambda}+(2l-1)\boldsymbol{\delta}\bigr)\cdot
  \xi_{d_k}(\eta;\boldsymbol{\lambda}-\boldsymbol{\delta})\nonumber\\
  &\qquad(j=\tfrac{M+3}{2},\ldots,M+1;k=1,\ldots,M),\\
  a_{j,M+1}&=f_n(\boldsymbol{\lambda})\mathcal{E}_n(\boldsymbol{\lambda})^{\frac{M+1}{2}-j}
  P_{n-1}(\eta;\boldsymbol{\lambda}+\boldsymbol{\delta})
  \ \ (j=1,\ldots,\tfrac{M+1}{2}),\nonumber\\
  a_{j,M+1}&=b_n(\boldsymbol{\lambda}-\boldsymbol{\delta})
  \mathcal{E}_n(\boldsymbol{\lambda})^{j-\frac{M+3}{2}}
  P_{n+1}(\eta;\boldsymbol{\lambda}-\boldsymbol{\delta})
  \ \ (j=\tfrac{M+3}{2},\ldots,M+1),
  \nonumber
\end{align}
for odd $M$.

\subsubsection{single type}
\label{sec:idQM_single}

As remarked after \eqref{idQM:cPDn2}, when the index set $\mathcal{D}$
consists of type $\text{\Romannumeral{1}}$ indices only (or type $\text{\Romannumeral{2}}$ indices only),
simplifications occur.
We consider such a single type index set here.
As shown in \cite{equiv_miop,t13}, the multi-indexed polynomials with both
type indices can always be recovered from those with single type indices
$\mathcal{D}'$ (type $\text{\Romannumeral{1}}$ only or type $\text{\Romannumeral{2}}$ only) with shifted parameters
$\boldsymbol{\lambda}'$, namely
\begin{equation}
  P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})
  =\frac{c_{\mathcal{D},n}^{P}(\boldsymbol{\lambda})}
  {c_{\mathcal{D}',n}^{P}(\boldsymbol{\lambda}')}
  P_{\mathcal{D}',n}(\eta;\boldsymbol{\lambda}'),
  \label{PDn=PD'n}
\end{equation}
where $c_{\mathcal{D},n}^{P}(\boldsymbol{\lambda})$ is the coefficient of the
highest degree term of $P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})$ \cite{os27},
and $\mathcal{D}'$ and $\boldsymbol{\lambda}'$ are determined by $\mathcal{D}$ and
$\boldsymbol{\lambda}$ \cite{equiv_miop}.

Let us define $\mathcal{F}'(\boldsymbol{\lambda})$,
$\mathcal{B}'(\boldsymbol{\lambda})$, $f'_n(\boldsymbol{\lambda})$,
$\mathcal{E}'_n(\boldsymbol{\lambda})$, etc. as follows:
\begin{align}
  &\mathcal{F}'(\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \mathcal{F}\bigl(\mathfrak{t}(\boldsymbol{\lambda})\bigr),\quad
  \mathcal{B}'(\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \mathcal{B}\bigl(\mathfrak{t}(\boldsymbol{\lambda})\bigr),\nonumber\\
  &f'_n(\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  f_n\bigl(\mathfrak{t}(\boldsymbol{\lambda})\bigr),\quad
  b'_n(\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  b_n\bigl(\mathfrak{t}(\boldsymbol{\lambda})\bigr),\quad
  \tilde{f}'_{\text{v}}(\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \tilde{f}_{\text{v}}\bigl(\mathfrak{t}(\boldsymbol{\lambda})\bigr),\quad
  \tilde{b}'_{\text{v}}(\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \tilde{b}_{\text{v}}\bigl(\mathfrak{t}(\boldsymbol{\lambda})\bigr),\\
  &\mathcal{E}'_n(\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \mathcal{E}_n\bigl(\mathfrak{t}(\boldsymbol{\lambda})\bigr)
  =f'_n(\boldsymbol{\lambda})b'_{n-1}(\boldsymbol{\lambda}),\quad
  \tilde{\mathcal{E}}'_{\text{v}}(\boldsymbol{\lambda})\stackrel{\text{def}}{=}
  \tilde{\mathcal{E}}_{\text{v}}\bigl(\mathfrak{t}(\boldsymbol{\lambda})\bigr)
  =\tilde{f}'_{\text{v}}(\boldsymbol{\lambda})\tilde{b}'_{\text{v}}(\boldsymbol{\lambda}).
  \nonumber
\end{align}
For $\mathcal{F}'(\boldsymbol{\lambda})$, $\mathcal{B}'(\boldsymbol{\lambda})$,
$f'_n(\boldsymbol{\lambda})$, $b'_n(\boldsymbol{\lambda})$ and $\mathcal{E}'_n(\boldsymbol{\lambda})$,
the type of $\mathfrak{t}$ is $\text{\Romannumeral{1}}$ or $\text{\Romannumeral{2}}$,
e.g., $\mathcal{F}^{\prime\,\text{\Romannumeral{1}}}(\boldsymbol{\lambda})
=\mathcal{F}\bigl(\mathfrak{t}^{\text{\Romannumeral{1}}}(\boldsymbol{\lambda})\bigr)$ and
$\mathcal{F}^{\prime\,\text{\Romannumeral{2}}}(\boldsymbol{\lambda})
=\mathcal{F}\bigl(\mathfrak{t}^{\text{\Romannumeral{2}}}(\boldsymbol{\lambda})\bigr)$.
For $\tilde{f}'_{\text{v}}(\boldsymbol{\lambda})$,
$\tilde{b}'_{\text{v}}(\boldsymbol{\lambda})$ and
$\tilde{\mathcal{E}}'_{\text{v}}(\boldsymbol{\lambda})$,
the type of $\mathfrak{t}$ is matched to that of them,
e.g., $\tilde{f}^{\text{\Romannumeral{1}}\,\prime}_{\text{v}}(\boldsymbol{\lambda})
=\tilde{f}^{\text{\Romannumeral{1}}}_{\text{v}}\bigl(\mathfrak{t}^{\text{\Romannumeral{1}}}(\boldsymbol{\lambda})\bigr)$ and
$\tilde{f}^{\text{\Romannumeral{2}}\,\prime}_{\text{v}}(\boldsymbol{\lambda})
=\tilde{f}^{\text{\Romannumeral{2}}}_{\text{v}}\bigl(\mathfrak{t}^{\text{\Romannumeral{2}}}(\boldsymbol{\lambda})\bigr)$.
Then the relations \eqref{idQM:FP=} and \eqref{idQM:Fnuxi=} with the
replacement $\boldsymbol{\lambda}\mapsto\mathfrak{t}(\boldsymbol{\lambda})$
(the type of $\mathfrak{t}$ is $\text{\Romannumeral{1}}$ or $\text{\Romannumeral{2}}$ for \eqref{idQM:FP=},
and matched to $\nu^{-1}\check{\xi}_{\text{v}}$ for \eqref{idQM:Fnuxi=})
and exchange $n\leftrightarrow\text{v}$ give
\begin{alignat}{2}
  \hat{\mathcal{F}}'(\boldsymbol{\lambda})\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})
  &=f'_{\text{v}}(\boldsymbol{\lambda})
  \check{\xi}_{\text{v}-1}(x;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}}),
  &\ \ \hat{\mathcal{B}}'(\boldsymbol{\lambda})
  \check{\xi}_{\text{v}-1}(x;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})
  &=b'_{\text{v}-1}(\boldsymbol{\lambda})\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda}),
  \label{idQM:F'xiv=}\\
  \!\!\hat{\mathcal{F}}'(\boldsymbol{\lambda})(\nu\check{P}_n)(x;\boldsymbol{\lambda})
  &=\tilde{f}'_n(\boldsymbol{\lambda})
  (\nu\check{P}_n)(x;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}}),
  &\ \ \hat{\mathcal{B}}'(\boldsymbol{\lambda})
  (\nu\check{P}_n)(x;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})
  &=\tilde{b}'_n(\boldsymbol{\lambda})(\nu\check{P}_n)(x;\boldsymbol{\lambda}),\!\!\!
  \label{idQM:F'nuPn=}
\end{alignat}
where $(\nu\check{P}_n)(x;\boldsymbol{\lambda})\stackrel{\text{def}}{=}\nu(x;\boldsymbol{\lambda})
\check{P}_n(x;\boldsymbol{\lambda})$.

In the following we consider the index set $\mathcal{D}$ with single type
twists. The type of unspecified $\mathfrak{t}$, $\nu(x;\boldsymbol{\lambda})$,
$\check{U}(x;\boldsymbol{\lambda})$, $\mathcal{F}'(\boldsymbol{\lambda})$,
$\mathcal{B}'(\boldsymbol{\lambda})$, $f'_n(\boldsymbol{\lambda})$,
$\mathcal{E}'_n(\boldsymbol{\lambda})$, $\tilde{\boldsymbol{\delta}}$ etc.\ is same as
the type of the indices of $\mathcal{D}$.
The expressions \eqref{idQM:cXiD2}--\eqref{idQM:cPDn2} are simplified to
\begin{align}
  \check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})
  &=\text{W}_{\gamma}[\check{\xi}_{d_1},\ldots,
  \check{\xi}_{d_M}](x;\boldsymbol{\lambda})
  \times\varphi_M(x)^{-1}
  \prod_{j=1}^M\frac{\nu\bigl(x;\boldsymbol{\lambda}+(M-1)\boldsymbol{\delta}\bigr)}
  {\nu(x^{(M)}_j;\boldsymbol{\lambda})}\nonumber\\
  &\quad\times\biggl(
  \prod_{m=1}^M\check{U}\bigl(x;\boldsymbol{\lambda}
  +(M-1-2m)\boldsymbol{\delta}\bigr)^{M-m}\biggr)^{-1},
  \label{idQM:cXiD3}\\
  \check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})
  &=\text{W}_{\gamma}[\check{\xi}_{d_1},\ldots,
  \check{\xi}_{d_M},\nu\check{P}_n](x;\boldsymbol{\lambda})
  \times\varphi_{M+1}(x)^{-1}
  \prod_{j=1}^{M+1}\frac{\nu(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})}
  {\nu(x^{(M+1)}_j;\boldsymbol{\lambda})}\nonumber\\
  &\quad\times
  \nu(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})^{-1}
  \biggl(\prod_{m=1}^M\check{U}\bigl(x;\boldsymbol{\lambda}
  +(M-2m)\boldsymbol{\delta}\bigr)^{M-m}\biggr)^{-1}.
  \label{idQM:cPDn3}
\end{align}
By taking $\mathcal{F}'$ or $\mathcal{B}'$ as $\hat{D}_j$ and $\hat{D}'_j$
and using shape invariance properties
\eqref{idQM:F'xiv=}--\eqref{idQM:F'nuPn=},
the Casoratians $\text{W}_{\gamma}[\check{\xi}_{d_1},\ldots,
\check{\xi}_{d_M}](x;\boldsymbol{\lambda})$
and $\text{W}_{\gamma}[\check{\xi}_{d_1},\ldots,
\check{\xi}_{d_M},\nu\check{P}_n](x;\boldsymbol{\lambda})$
in \eqref{idQM:cXiD3}--\eqref{idQM:cPDn3}
can be rewritten in various ways.
We consider typical two cases:
\begin{align}
  \!\text{A}:\ \,&
  \hat{D}_j=\mathcal{F}'\bigl(\boldsymbol{\lambda}+(j-1)\tilde{\boldsymbol{\delta}}\bigr),
  \ \ \hat{D}'_j=\mathcal{B}'\bigl(\boldsymbol{\lambda}+(j-2)\tilde{\boldsymbol{\delta}}\bigr)
  \ \ (j=1,2,\ldots),
  \label{idQM:singleA}\\
  \!\text{B}:\ \,&
  \hat{D}_{2l-1}=\mathcal{F}'(\boldsymbol{\lambda}),
  \,\hat{D}_{2l}=\mathcal{B}'(\boldsymbol{\lambda}),
  \,\hat{D}'_{2l-1}=\mathcal{B}'(\boldsymbol{\lambda}-\tilde{\boldsymbol{\delta}}),
  \,\hat{D}'_{2l}=\mathcal{F}'(\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})
  \ \ (l=1,2,\ldots),\!
  \label{idQM:singleB}
\end{align}
which correspond to the cases A and B in \cite{os37}.

\medskip\noindent
{\bf case A:}\\
Firstly we consider case A \eqref{idQM:singleA}.
The functions $q_j(x)$, $r_j(x)$, $q'_j(x)$, $r'_j(x)$ and $D_j(x)$ are
\begin{align}
  q_j(x)&=i\varphi(x)^{-1},\quad r_j(x)=-i\varphi(x)^{-1},\nonumber\\
  q'_j(x)&=-iV\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})+(j-2)\boldsymbol{\delta}\bigr)
  \varphi(x-i\tfrac{\gamma}{2}),\quad
  r'_j(x)=iV^*\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})+(j-2)\boldsymbol{\delta}\bigr)
  \varphi(x+i\tfrac{\gamma}{2}),\nonumber\\
  D_j(x)&=-i\varphi(x)^{-1}
  \check{S}\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})+(j-2)\boldsymbol{\delta}\bigr).
\end{align}
The shape invariance properties \eqref{idQM:F'xiv=}--\eqref{idQM:F'nuPn=}
give ($k=0,1,\ldots$)
\begin{align}
  &\laprod{l=1}{k}\hat{D}_l\cdot\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})
  =\prod_{m=0}^{k-1}f'_{\text{v}-m}(\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})\cdot
  \check{\xi}_{\text{v}-k}(x;\boldsymbol{\lambda}+k\tilde{\boldsymbol{\delta}}),\nonumber\\[-2pt]
  &\laprod{l=1}{k}\hat{D}_l\cdot
  (\nu\check{P}_n)(x;\boldsymbol{\lambda})
  =\prod_{m=0}^{k-1}\tilde{f}'_n(\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})\cdot
  (\nu\check{P}_n)(x;\boldsymbol{\lambda}+k\tilde{\boldsymbol{\delta}}),
  \label{idQM:singlesiA}\\[-2pt]
  &\hat{D}'_{k+1}\laprod{l=1}{k}\hat{D}_l\cdot
  \check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})
  =b'_{\text{v}-k}\bigl(\boldsymbol{\lambda}+(k-1)\tilde{\boldsymbol{\delta}}\bigr)
  \prod_{m=0}^{k-1}f'_{\text{v}-m}(\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})\cdot
  \check{\xi}_{\text{v}+1-k}\bigl(x;\boldsymbol{\lambda}
  +(k-1)\tilde{\boldsymbol{\delta}}\bigr),\nonumber\\[-2pt]
  &\hat{D}'_{k+1}\laprod{l=1}{k}\hat{D}_l\cdot
  (\nu\check{P}_n)(x;\boldsymbol{\lambda})
  =\tilde{b}'_n\bigl(\boldsymbol{\lambda}+(k-1)\tilde{\boldsymbol{\delta}}\bigr)
  \prod_{m=0}^{k-1}\tilde{f}'_n(\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})\cdot
  (\nu\check{P}_n)\bigl(x;\boldsymbol{\lambda}+(k-1)\tilde{\boldsymbol{\delta}}\bigr).
  \nonumber
\end{align}

We rewrite \eqref{idQM:cXiD3}--\eqref{idQM:cPDn3} in the following way.
(\romannumeral1) By using \eqref{idQM:Wgid}--\eqref{idQM:Wgid,ajk},
rewrite the Casoratians in \eqref{idQM:cXiD3}--\eqref{idQM:cPDn3} as
determinants $\det(a_{j,k})_{1\leq j,k\leq M'}$ ($M'=M$ for
$\check{\Xi}_{\mathcal{D}}$ and $M'=M+1$ for $\check{P}_{\mathcal{D},n}$).
(\romannumeral2) Rewrite each matrix element $a_{j,k}$ by
\eqref{idQM:singlesiA}.
For $\check{P}_{\mathcal{D},n}$, we do the following
(\romannumeral3)--(\romannumeral5).
(\romannumeral3) Divide the $(M+1)$-th column of $\det(a_{j,k})$ by
$\nu(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})$, and multiply the determinant
by $\nu(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})$.
(\romannumeral4) Rewrite
$\nu(x;\boldsymbol{\lambda}')/\nu(x;\boldsymbol{\lambda}+M\boldsymbol{\delta})$
in the $(M+1)$-th column of $\det(a_{j,k})$ by \eqref{nuMk}.
(\romannumeral5) Multiply the $(M+1)$-th column of $\det(a_{j,k})$ by
$\prod_{m=1}^M\check{U}\bigl(x;\boldsymbol{\lambda}+(M-2m)\boldsymbol{\delta}\bigl)$,
and divide the determinant by
$\prod_{m=1}^M\check{U}\bigl(x;\boldsymbol{\lambda}+(M-2m)\boldsymbol{\delta}\bigl)$.
Then $\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})$ and
$\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ are expressed in terms of
$\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda}')$,
$\check{P}_n(x;\boldsymbol{\lambda}')$, $\check{U}(x;\boldsymbol{\lambda}')$,
$\check{S}(x;\boldsymbol{\lambda}')$ and $\varphi(x')$.
Straightforward calculation shows that the factors $\varphi(x')$ are
canceled out.
Thus $\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})$ and
$\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ are expressed in terms of
$\eta(x)$, namely
$\Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})$ and
$P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})$ are expressed without explicit
$x$-dependence.
Their final forms are as follows.
The denominator polynomial $\Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})$ is
\begin{equation}
  \Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})
  =(-1)^{\frac16(M-1)M(M+1)}
  \biggl(\prod_{m=1}^{[\frac{M}{2}]}
  S\bigl(\eta;\mathfrak{t}(\boldsymbol{\lambda})+(M-1-2m)\boldsymbol{\delta}\bigr)\biggr)^{-1}
  \det(a_{j,k})_{1\leq j,k\leq M},
  \label{idQM:singleA:XiD}
\end{equation}
where $a_{j,k}$ are
\begin{align}
  a_{j,k}&=\prod_{m=0}^{M-2j}f'_{d_k-m}(\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})
  \cdot\xi_{d_k-M-1+2j}\bigl(\eta;\boldsymbol{\lambda}
  +(M+1-2j)\tilde{\boldsymbol{\delta}}\bigr)\nonumber\\
  &\qquad(j=1,\ldots,[\tfrac{M+1}{2}];k=1,\ldots,M),\nonumber\\
  a_{j,k}&=b'_{d_k-2j+M+2}\bigl(\boldsymbol{\lambda}+(2j-M-3)\tilde{\boldsymbol{\delta}}\bigr)
  \prod_{m=0}^{2j-M-3}f'_{d_k-m}(\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})\\
  &\quad\times
  \xi_{d_k-2j+M+3}\bigl(\eta;\boldsymbol{\lambda}+(2j-M-3)\tilde{\boldsymbol{\delta}}\bigr)
  \ \ (j=[\tfrac{M+3}{2}],\ldots,M;k=1,\ldots,M).
  \nonumber
\end{align}
The multi-indexed orthogonal polynomial
$P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})$ is
\begin{equation}
  P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})
  =(-1)^{\frac16M(M+1)(M+2)}
  \biggl(\prod_{m=1}^{[\frac{M+1}{2}]}
  S\bigl(\eta;\mathfrak{t}(\boldsymbol{\lambda})+(M-2m)\boldsymbol{\delta}\bigr)\biggr)^{-1}
  \det(a_{j,k})_{1\leq j,k\leq M+1},
  \label{idQM:singleA:PDn}
\end{equation}
where $a_{j,k}$ are
\begin{align}
  a_{j,k}&=\prod_{m=0}^{M+1-2j}f'_{d_k-m}(\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})
  \cdot
  \xi_{d_k-M-2+2j}\bigl(\eta;\boldsymbol{\lambda}+(M+2-2j)\tilde{\boldsymbol{\delta}}\bigr)\nonumber\\
  &\qquad(j=1,\ldots,[\tfrac{M+2}{2}];k=1,\ldots,M),\nonumber\\
  a_{j,k}&=b'_{d_k-2j+M+3}\bigl(\boldsymbol{\lambda}+(2j-M-4)\tilde{\boldsymbol{\delta}}\bigr)
  \prod_{m=0}^{2j-M-4}f'_{d_k-m}(\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})\nonumber\\
  &\quad\times
  \xi_{d_k-2j+M+4}\bigl(\eta;\boldsymbol{\lambda}+(2j-M-4)\tilde{\boldsymbol{\delta}}\bigr)
  \ \ (j=[\tfrac{M+4}{2}],\ldots,M+1;k=1,\ldots,M),\nonumber\\
  a_{j,M+1}&=\prod_{m=0}^{M+1-2j}
  \tilde{f}'_n(\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})\cdot
  \prod_{l=M+2-j}^MU\bigl(\eta;\boldsymbol{\lambda}+(M-2l)\boldsymbol{\delta}\bigr)\cdot
  P_n\bigl(\eta;\boldsymbol{\lambda}+(M+2-2j)\tilde{\boldsymbol{\delta}}\bigr)\nonumber\\
  &\qquad(j=1,\ldots,[\tfrac{M+2}{2}]),\\
  a_{j,M+1}&=\tilde{b}'_n\bigl(\boldsymbol{\lambda}+(2j-M-4)\tilde{\boldsymbol{\delta}}\bigr)
  \prod_{m=0}^{2j-M-4}\tilde{f}'_n(\boldsymbol{\lambda}+m\tilde{\boldsymbol{\delta}})\cdot
  \prod_{l=j-1}^MU\bigl(\eta;\boldsymbol{\lambda}+(M-2l)\boldsymbol{\delta}\bigr)\nonumber\\
  &\quad\times
  P_n\bigl(\eta;\boldsymbol{\lambda}+(2j-M-4)\tilde{\boldsymbol{\delta}}\bigr)
  \ \ (j=[\tfrac{M+4}{2}],\ldots,M+1).
  \nonumber
\end{align}

\medskip\noindent
{\bf case B:}\\
Secondly we consider case B \eqref{idQM:singleB}.
The functions $q_j(x)$, $r_j(x)$, $q'_j(x)$, $r'_j(x)$ and $D_j(x)$ are
($l=1,2,\ldots$)
\begin{align}
  q_{2l-1}(x)&=i\varphi(x)^{-1},\quad r_{2l-1}(x)=-i\varphi(x)^{-1},\nonumber\\
  q_{2l}(x)&=-iV\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})\bigr)
  \varphi(x-i\tfrac{\gamma}{2}),\quad
  r_{2l}(x)=iV^*\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})\bigr)
  \varphi(x+i\tfrac{\gamma}{2}),\nonumber\\
  q'_{2l-1}(x)&=-iV\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})-\boldsymbol{\delta}\bigr)
  \varphi(x-i\tfrac{\gamma}{2}),\quad
  r'_{2l-1}=iV^*\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})-\boldsymbol{\delta}\bigr)
  \varphi(x+i\tfrac{\gamma}{2}),\\
  q'_{2l}(x)&=i\varphi(x)^{-1},\quad r'_{2l}(x)=-i\varphi(x)^{-1},\nonumber\\
  D_{2l-1}(x)&=-i\varphi(x)^{-1}
  \check{S}\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})-\boldsymbol{\delta}\bigr),\quad
  D_{2l}(x)=i\varphi(x)^{-1}
  \check{S}\bigl(x;\mathfrak{t}(\boldsymbol{\lambda})\bigr).
  \nonumber
\end{align}
The shape invariance properties \eqref{idQM:F'xiv=}--\eqref{idQM:F'nuPn=}
give ($l=1,2,\ldots$)
\begin{align}
  &\laprod{m=1}{2l-2}\hat{D}_m\cdot\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})
  =\mathcal{E}'_{\text{v}}(\boldsymbol{\lambda})^{l-1}
  \check{\xi}_{\text{v}}(x;\boldsymbol{\lambda}),\nonumber\\[-2pt]
  &\laprod{m=1}{2l-1}\hat{D}_m\cdot\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})
  =f'_{\text{v}}(\boldsymbol{\lambda})\mathcal{E}'_{\text{v}}(\boldsymbol{\lambda})^{l-1}
  \check{\xi}_{\text{v}-1}(x;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}}),\nonumber\\[-2pt]
  &\laprod{m=1}{2l-2}\hat{D}_m\cdot
  (\nu\check{P}_n)(x;\boldsymbol{\lambda})
  =\tilde{\mathcal{E}}'_n(\boldsymbol{\lambda})^{l-1}
  (\nu\check{P}_n)(x;\boldsymbol{\lambda}),\nonumber\\[-2pt]
  &\laprod{m=1}{2l-1}\hat{D}_m\cdot
  (\nu\check{P}_n)(x;\boldsymbol{\lambda})
  =\tilde{f}'_n(\boldsymbol{\lambda})
  \tilde{\mathcal{E}}'_n(\boldsymbol{\lambda})^{l-1}
  (\nu\check{P}_n)(x;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}}),\\[-2pt]
  &\hat{D}'_{2l-1}\laprod{m=1}{2l-2}\hat{D}_m\cdot
  \check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})
  =b'_{\text{v}}(\boldsymbol{\lambda}-\tilde{\boldsymbol{\delta}})
  \mathcal{E}'_{\text{v}}(\boldsymbol{\lambda})^{l-1}
  \check{\xi}_{\text{v}+1}(x;\boldsymbol{\lambda}-\tilde{\boldsymbol{\delta}}),\nonumber\\[-2pt]
  &\hat{D}'_{2l}\laprod{m=1}{2l-1}\hat{D}_m\cdot
  \check{\xi}_{\text{v}}(x;\boldsymbol{\lambda})
  =f'_{\text{v}-1}(\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})f'_{\text{v}}(\boldsymbol{\lambda})
  \mathcal{E}'_{\text{v}}(\boldsymbol{\lambda})^{l-1}
  \check{\xi}_{\text{v}-2}(x;\boldsymbol{\lambda}+2\tilde{\boldsymbol{\delta}}),\nonumber\\[-2pt]
  &\hat{D}'_{2l-1}\laprod{m=1}{2l-2}\hat{D}_m\cdot
  (\nu\check{P}_n)(x;\boldsymbol{\lambda})
  =\tilde{b}'_n(\boldsymbol{\lambda}-\tilde{\boldsymbol{\delta}})
  \tilde{\mathcal{E}}'_n(\boldsymbol{\lambda})^{l-1}
  (\nu\check{P}_n)(x;\boldsymbol{\lambda}-\tilde{\boldsymbol{\delta}}),\nonumber\\[-2pt]
  &\hat{D}'_{2l}\laprod{m=1}{2l-1}\hat{D}_m\cdot
  (\nu\check{P}_n)(x;\boldsymbol{\lambda})
  =\tilde{f}'_n(\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})
  \tilde{f}'_n(\boldsymbol{\lambda})
  \tilde{\mathcal{E}}'_n(\boldsymbol{\lambda})^{l-1}
  (\nu\check{P}_n)(x;\boldsymbol{\lambda}+2\tilde{\boldsymbol{\delta}}).
  \nonumber
\end{align}

Like case A, we rewrite \eqref{idQM:cXiD3}--\eqref{idQM:cPDn3} by the steps
(\romannumeral1)--(\romannumeral5).
In (\romannumeral5), multiply and divide by
$\prod_{l=1}^{[\frac{M+2}{2}]}
\check{U}\bigl(x;\boldsymbol{\lambda}+(M-2l)\boldsymbol{\delta}\bigr)$.
By using \eqref{varphiprop3} and \eqref{VV*},
$\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})$ and
$\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ are expressed in terms of
$\check{\xi}_{\text{v}}(x;\boldsymbol{\lambda}')$,
$\check{P}_n(x;\boldsymbol{\lambda}')$, $\check{U}(x;\boldsymbol{\lambda}')$,
$\check{S}(x;\boldsymbol{\lambda}')$ and $\varphi(x')$.
Straightforward calculation shows that the factors $\varphi(x')$ are
canceled out.
Thus $\check{\Xi}_{\mathcal{D}}(x;\boldsymbol{\lambda})$ and
$\check{P}_{\mathcal{D},n}(x;\boldsymbol{\lambda})$ are expressed in terms of
$\eta(x)$, namely
$\Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})$ and
$P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})$ are expressed without explicit
$x$-dependence.
Their final forms are as follows.
The denominator polynomial $\Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})$ is
\begin{align}
  \Xi_{\mathcal{D}}(\eta;\boldsymbol{\lambda})
  &=(-1)^{[\frac{M}{2}]}
  \bigl(\alpha^{\text{t}}(\boldsymbol{\lambda})
  \alpha^{\bar{\text{t}}}(\boldsymbol{\lambda})^{-1}\kappa
  \bigr)^{\frac12[\frac{M-2}{2}][\frac{M}{2}]}
  \det(a_{j,k})_{1\leq j,k\leq M}\nonumber\\
  &\quad\times\biggl(
  \prod_{m=[\frac{M+3}{2}]}^MU\bigl(\eta;\boldsymbol{\lambda}
  +(M-1-2m)\boldsymbol{\delta}\bigr)^{M-m}\cdot
  \prod_{m=0}^{[\frac{M-2}{2}]}U^{\bar{\text{t}}}
  \bigl(\eta;\boldsymbol{\lambda}+(M-1-2m)\boldsymbol{\delta}\bigr)^m\nonumber\\
  &\qquad\times
  S\bigl(\eta;\mathfrak{t}(\boldsymbol{\lambda})-\tfrac12(1+(-1)^M)\boldsymbol{\delta}
  \bigr)^{[\frac{M}{2}]}\biggr)^{-1},
  \label{idQM:singleB:XiD}
\end{align}
where $\text{t}$ is the type of the indices of $\mathcal{D}$ and $a_{j,k}$ are
\begin{align}
  a_{j,k}&=\mathcal{E}'_{d_k}(\boldsymbol{\lambda})^{\frac{M+1}{2}-j}
  \xi_{d_k}(\eta;\boldsymbol{\lambda})
  \ \ (j=1,\ldots,\tfrac{M+1}{2};k=1,\ldots,M),\\
  a_{j,k}&=f'_{d_k-1}(\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})
  f'_{d_k}(\boldsymbol{\lambda})
  \mathcal{E}'_{d_k}(\boldsymbol{\lambda})^{j-\frac{M+3}{2}}
  \xi_{d_k-2}(\eta;\boldsymbol{\lambda}+2\tilde{\boldsymbol{\delta}})
  \ \ (j=\tfrac{M+3}{2},\ldots,M;k=1,\ldots,M),
  \nonumber
\end{align}
for odd $M$ and
\begin{align}
  a_{j,k}&=f'_{d_k}(\boldsymbol{\lambda})
  \mathcal{E}'_{d_k}(\boldsymbol{\lambda})^{\frac{M}{2}-j}
  \xi_{d_k-1}(\eta;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})
  \ \ (j=1,\ldots,\tfrac{M}{2};k=1,\ldots,M),\nonumber\\
  a_{j,k}&=b'_{d_k}(\boldsymbol{\lambda}-\tilde{\boldsymbol{\delta}})
  \mathcal{E}'_{d_k}(\boldsymbol{\lambda})^{j-\frac{M+2}{2}}
  \xi_{d_k+1}(\eta;\boldsymbol{\lambda}-\tilde{\boldsymbol{\delta}})
  \ \ (j=\tfrac{M+2}{2},\ldots,M;k=1,\ldots,M),
\end{align}
for even $M$.
The multi-indexed orthogonal polynomial
$P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})$ is
\begin{align}
  P_{\mathcal{D},n}(\eta;\boldsymbol{\lambda})
  &=(-1)^{[\frac{M+1}{2}]}
  \bigl(\alpha^{\text{t}}(\boldsymbol{\lambda})
  \alpha^{\bar{\text{t}}}(\boldsymbol{\lambda})^{-1}\kappa
  \bigr)^{\frac12[\frac{M-1}{2}][\frac{M+1}{2}]}
  \det(a_{j,k})_{1\leq j,k\leq M+1}\nonumber\\
  &\quad\times\biggl(
  \prod_{m=[\frac{M+4}{2}]}^MU\bigl(\eta;\boldsymbol{\lambda}
  +(M-2m)\boldsymbol{\delta}\bigr)^{M-m}\cdot
  \prod_{m=0}^{[\frac{M-1}{2}]}U^{\bar{\text{t}}}
  \bigl(\eta;\boldsymbol{\lambda}
  +(M-2m)\boldsymbol{\delta}\bigr)^m\nonumber\\
  &\qquad\times
  S\bigl(\eta;\mathfrak{t}(\boldsymbol{\lambda})-\tfrac12(1-(-1)^M)\boldsymbol{\delta}
  \bigr)^{[\frac{M+1}{2}]}\biggr)^{-1},
  \label{idQM:singleB:PDn}
\end{align}
where $\text{t}$ is the type of the indices of $\mathcal{D}$ and $a_{j,k}$ are
\begin{align}
  a_{j,k}&=\mathcal{E}'_{d_k}(\boldsymbol{\lambda})^{\frac{M+2}{2}-j}
  \xi_{d_k}(\eta;\boldsymbol{\lambda})
  \ \ (j=1,\ldots,\tfrac{M+2}{2};k=1,\ldots,M),\nonumber\\
  a_{j,k}&=f'_{d_k-1}(\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})
  f'_{d_k}(\boldsymbol{\lambda})
  \mathcal{E}'_{d_k}(\boldsymbol{\lambda})^{j-\frac{M+4}{2}}
  \xi_{d_k-2}(\eta;\boldsymbol{\lambda}+2\tilde{\boldsymbol{\delta}})\nonumber\\
  &\qquad(j=\tfrac{M+4}{2},\ldots,M+1;k=1,\ldots,M),\nonumber\\
  a_{j,M+1}&=\tilde{\mathcal{E}}'_n(\boldsymbol{\lambda})^{\frac{M+2}{2}-j}
  U(\eta;\boldsymbol{\lambda}-2\boldsymbol{\delta})P_n(\eta;\boldsymbol{\lambda})
  \ \ (j=1,\ldots,\tfrac{M+2}{2}),\\
  a_{j,M+1}&=\tilde{f}'_n(\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})
  \tilde{f}'_n(\boldsymbol{\lambda})
  \tilde{\mathcal{E}}'_n(\boldsymbol{\lambda})^{j-\frac{M+4}{2}}
  P_n(\eta;\boldsymbol{\lambda}+2\tilde{\boldsymbol{\delta}})
  \ \ (j=\tfrac{M+4}{2},\ldots,M+1),
  \nonumber
\end{align}
for even $M$ and
\begin{align}
  a_{j,k}&=f'_{d_k}(\boldsymbol{\lambda})
  \mathcal{E}'_{d_k}(\boldsymbol{\lambda})^{\frac{M+1}{2}-j}
  \xi_{d_k-1}(\eta;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})
  \ \ (j=1,\ldots,\tfrac{M+1}{2};k=1,\ldots,M),\nonumber\\
  a_{j,k}&=b'_{d_k}(\boldsymbol{\lambda}-\tilde{\boldsymbol{\delta}})
  \mathcal{E}'_{d_k}(\boldsymbol{\lambda})^{j-\frac{M+3}{2}}
  \xi_{d_k+1}(\eta;\boldsymbol{\lambda}-\tilde{\boldsymbol{\delta}})
  \ \ (j=\tfrac{M+3}{2},\ldots,M+1;k=1,\ldots,M),\nonumber\\
  a_{j,M+1}&=\tilde{f}'_n(\boldsymbol{\lambda})
  \tilde{\mathcal{E}}'_n(\boldsymbol{\lambda})^{\frac{M+1}{2}-j}
  P_n(\eta;\boldsymbol{\lambda}+\tilde{\boldsymbol{\delta}})
  \ \ (j=1,\ldots,\tfrac{M+1}{2}),\\
  a_{j,M+1}&=\tilde{b}'_n(\boldsymbol{\lambda}-\tilde{\boldsymbol{\delta}})
  \tilde{\mathcal{E}}'_n(\boldsymbol{\lambda})^{j-\frac{M+3}{2}}
  U(\eta;\boldsymbol{\lambda}-\boldsymbol{\delta})
  P_n(\eta;\boldsymbol{\lambda}-\tilde{\boldsymbol{\delta}})
  \ \ (j=\tfrac{M+3}{2},\ldots,M+1),
  \nonumber
\end{align}
for odd $M$.

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

The multi-indexed orthogonal polynomials (the Meixner, little $q$-Laguerre,
little $q$-Jacobi, Racah, $q$-Racah, Wilson and Askey-Wilson types)
introduced in the framework of the discrete quantum mechanics are expressed
in terms of Casoratians.
Various new determinant expressions of the multi-indexed orthogonal
polynomials are derived by using the properties of the Casoratians and
shape invariance.
Two typical cases are presented explicitly.
For the Meixner, little $q$-Laguerre, little $q$-Jacobi, Wilson and
Askey-Wilson cases, the new expressions for the multi-indexed polynomials
$P_{\mathcal{D},n}(\eta)$ do not have explicit dependence on $x$,
which is the coordinate of the quantum system.
We hope that these new expressions are helpful for deeper understanding of
the multi-indexed orthogonal polynomials.

In \cite{os37}, simplified expressions of the multi-indexed Laguerre
and Jacobi polynomials are derived.
From the viewpoint adopted in this paper, the calculation in \cite{os37}
is regarded as follows.
The eigenfunctions of the deformed systems, namely the multi-indexed orthogonal
polynomials are expressed in terms of the Wronskian.
The Wronskian of a set of $n$ functions $\{f_j(x)\}$ is defined by
\begin{equation}
  \text{W}[f_1,f_2,\ldots,f_n](x)\stackrel{\text{def}}{=}
  \det\Bigl(\frac{d^{j-1}f_k(x)}{dx^{j-1}}\Bigr)_{1\leq j,k\leq n},
\end{equation}
(for $n=0$, we set $\text{W}\,[\cdot](x)=1$), and the following formula
holds for any smooth functions $q_j(x)$ ($j=1,2,\ldots$):
\begin{equation}
  \det\Bigl(\laprod{l=1}{j-1}\hat{D}_l\cdot f_k(x)
  \Bigr)_{1\leq j,k\leq n}
  =\text{W}[f_1,\ldots,f_n](x),
  \ \ \hat{D}_j=\frac{d}{dx}-q_j(x).
  \label{QM:Wid}
\end{equation}
The original systems have shape invariance
$\mathcal{A}(\boldsymbol{\lambda})\mathcal{A}(\boldsymbol{\lambda})^{\dagger}
=\mathcal{A}(\boldsymbol{\lambda}+\boldsymbol{\delta})^{\dagger}
\mathcal{A}(\boldsymbol{\lambda}+\boldsymbol{\delta})
+\mathcal{E}_1(\boldsymbol{\lambda})$, where $\mathcal{A}(\boldsymbol{\lambda})
=\frac{d}{dx}-\frac{d}{dx}\log\phi_0(x;\boldsymbol{\lambda})$ and
$\mathcal{A}(\boldsymbol{\lambda})^{\dagger}=
-\frac{d}{dx}-\frac{d}{dx}\log\phi_0(x;\boldsymbol{\lambda})$.
These operators act on the eigenstates and virtual states as follows:
\begin{alignat}{2}
  \mathcal{A}(\boldsymbol{\lambda})\phi_n(x;\boldsymbol{\lambda})
  &=f_n(\boldsymbol{\lambda})\phi_{n-1}(x;\boldsymbol{\lambda}+\boldsymbol{\delta}),&\quad
  \mathcal{A}(\boldsymbol{\lambda})^{\dagger}\phi_{n-1}(x;\boldsymbol{\lambda}+\boldsymbol{\delta})
  &=b_{n-1}(\boldsymbol{\lambda})\phi_n(x;\boldsymbol{\lambda}),
  \label{QM:Aphi=}\\
  \mathcal{A}(\boldsymbol{\lambda})\tilde{\phi}_{\text{v}}(x;\boldsymbol{\lambda})
  &=\tilde{f}_{\text{v}}(\boldsymbol{\lambda})
  \tilde{\phi}_{\text{v}}(x;\boldsymbol{\lambda}+\boldsymbol{\delta}),&\quad
  \mathcal{A}(\boldsymbol{\lambda})^{\dagger}
  \tilde{\phi}_{\text{v}}(x;\boldsymbol{\lambda}+\boldsymbol{\delta})
  &=\tilde{b}_{\text{v}}(\boldsymbol{\lambda})
  \tilde{\phi}_{\text{v}}(x;\boldsymbol{\lambda}),
  \label{QM:Aphit=}
\end{alignat}
where
\begin{alignat}{2}
  f_n(\boldsymbol{\lambda})&=\left\{
  \begin{array}{ll}
  -2&:\text{L}\\
  -2(n+g+h)&:\text{J}
  \end{array}\right.,&\quad
  b_{n-1}(\boldsymbol{\lambda})&=-2n:\text{L,\,J}\\
  \tilde{f}_{\text{v}}(\boldsymbol{\lambda})&=\left\{
  \begin{array}{ll}
  2&:\text{L,\,$\text{\Romannumeral{1}}$}\\
  -2(g-\frac12-\text{v})&:\text{L,\,$\text{\Romannumeral{2}}$}\\
  2(h-\frac12-\text{v})&:\text{J,\,$\text{\Romannumeral{1}}$}\\
  -2(g-\frac12-\text{v})&:\text{J,\,$\text{\Romannumeral{2}}$}
  \end{array}\right.,&\quad
  \tilde{b}_{\text{v}}(\boldsymbol{\lambda})&=\left\{
  \begin{array}{ll}
  -2(g+\frac12+\text{v})&:\text{L,\,$\text{\Romannumeral{1}}$}\\
  2&:\text{L,\,$\text{\Romannumeral{2}}$}\\
  -2(g+\frac12+\text{v})&:\text{J,\,$\text{\Romannumeral{1}}$}\\
  2(h+\frac12+\text{v})&:\text{J,\,$\text{\Romannumeral{2}}$}
  \end{array}\right..
\end{alignat}
By taking $\mathcal{A}$ or $\mathcal{A}^{\dagger}$ as $\hat{D}_j$ in
\eqref{QM:Wid} and using shape invariant properties
\eqref{QM:Aphi=}--\eqref{QM:Aphit=}, the Wronskians can be rewritten
in various ways.
Typical two cases are studied in \cite{os37}:
\begin{align}
  \text{A}:\ \ &
  \hat{D}_j=\mathcal{A}\bigl(\boldsymbol{\lambda}+(j-1)\boldsymbol{\delta}\bigr),
  \ \ (j=1,2,\ldots),
  \label{QM:A}\\
  \text{B}:\ \ &
  \hat{D}_{2l-1}=\mathcal{A}(\boldsymbol{\lambda}),
  \ \hat{D}_{2l}=-\mathcal{A}(\boldsymbol{\lambda})^{\dagger},
  \ \ (l=1,2,\ldots).
  \label{QM:B}
\end{align}
Like \S\,\ref{sec:idQM_single}, the equivalence property \eqref{PDn=PD'n}
holds \cite{equiv_miop,t13} and simplifications occur for the index set
with a single type. This calculation is left to readers as an exercise.

\section*{Acknowledgments}

I thank R.\,Sasaki for discussion and reading of the manuscript.
I am supported in part by Grant-in-Aid for Scientific Research
from the Ministry of Education, Culture, Sports, Science and Technology
(MEXT), No.25400395.

redpajama_set_name: RedPajamaArXiv