\section{Introduction}\label{intro}
The properties of {\em Jacobi polynomials}, $P_n^{(\alpha, \beta)}$ $(n=1,2,\dots; \alpha>-1, \beta>-1)$, may well consult in \cite[Chapter IV]{S75}; they are orthogonal on $[-1,1]$ with the weight function $(1-X)^\alpha(1+X)^\beta$, and
$$
y=P_n^{(\alpha, \beta)}(X)=\pFq{2}{1}{-n, n+\alpha+\beta+1}{1+\alpha}{\frac{1-X}{2}}
$$
satisfy the homogeneous differential equation of second order
\begin{align}\label{edo}
a \,y''+ b \, y'+c\,y=0,
\end{align}
where $a=a(X)=-X^2+1$ and $b=b(X; \alpha, \beta)= -(\alpha+\beta+2)X-\alpha+\beta$. (Having the explicit expression of $P_n^{(\alpha, \beta)}$, the parameter $c$ does not play any interesting role and it can be easily calculated.) Recall that the hypergeometric function ${}_{i}F_{j}$ is formally defined by the series (see \cite[(2.1.2)]{AAR99})
$$
\pFq{i}{j}{\alpha_1,\dots, \alpha_i}{\beta_1\dots, \beta_j}{X}=\sum_{k=0}^\infty \frac{(\alpha_1,\dots, \alpha_i)_k}{(\beta_1\dots, \beta_j)_k}\frac{X^k}{k!},
$$
where
$$
(\alpha_1,\dots, \alpha_i)_k=(\alpha_1)_k\cdots(\alpha_i)_k, \quad (\alpha_1)_k=\prod_{j=1}^k (\alpha_1+j-1);
$$
by convention, the empty product is $1$. From \eqref{edo}, Stieltjes proves an important statement concerning the dependence of the zeros of Jacobi polynomials on the parameters $\alpha$ and $\beta$  (see \cite{S87}). Soon after his work has been accepted for publication, Stieltjes implicitly acknowledges, in a note added at the end of the manuscript itself and also in a letter of February 3, 1887 to Hermite (see \cite[Lettre 106]{SH}), that the monotonicity of the zeros of Jacobi polynomials was previously proved by A. Markov in \cite{M86}. However, in Stieltjes' work, {\em ut in multis aliis rebus}, the ``How” is more important than the ``What” and, as he wrote to Hermite, {\em ``la d$\acute{e}$monstration que j'ai d$\acute{e}$velopp$\acute{e}$e pour les Acta Mathematica est diff$\acute{e}$rente de celle de M. Markoff''}. Historically, Stieltjes proves that given a positive definite real symmetric matrix with non-positive off-diagonal elements\footnote{These matrices are now known as Stieltjes' matrices (see \cite[Definition 3.4]{V62}).}, its inverse is also positive definite.  As a consequence, putting aside a clever manipulation of the differential equation \eqref{edo}, he shows that, since
\begin{align}\label{s1}
\frac ba=\frac{-(\alpha+\beta+2)X-\alpha+\beta}{-X^2+1}
\end{align}
 is a strictly decreasing function of $X\in (-1,1)$, and
\begin{align}\label{s2}
\frac{\partial}{\partial \alpha}\left(\frac ba\right)=\frac{1}{X-1}<0, \quad \frac{\partial}{\partial \beta}\left(\frac ba\right)=\frac{1}{X+1}>0,
\end{align}
for each $X\in (-1,1)$, the zeros of Jacobi polynomials are strictly decreasing functions of $\alpha$ on $(-1, \infty)$ and strictly increasing functions of $\beta$ on $(-1, \infty)$ (see \cite[(6)-(7)]{S87}). The proof of Markov is entirely different from that of Stieltjes and is based on the weight function. While it is true that Stieltjes worked directly with Jacobi polynomials ---and Markov proves his result through a general theorem---, his argument furnishes similar results for a more general differential equation (see \cite[Section $6.22$]{S75}). Stieltjes himself considered the ultraspherical case $\alpha=\beta$ and Szeg\H{o} noted that {\em ``the same method applies to Laguerre polynomials''} (see \cite[p. 123]{S75}). In particular, by definition, the {\em old classical orthogonal polynomials on the real line} (Hermite, Jacobi, and Laguerre)\footnote{We write ``classical orthogonal polynomials on the real line" rather than simply  ``classical orthogonal polynomials" because, for instance, from the algebraic point of view of Maroni (see \cite{M91}), the Jacobi polynomials exist and are ``classical'' even when $-\alpha, -\beta, -\alpha-\beta+1\not\in \mathbb{N}$  (see \cite[Chapters $8$ and $9$]{P18} for a recent survey on the subject). Moreover, the Bessel polynomials are classical in the same sense as the other three systems. As Maroni says, ``comme dans le roman d'Alexandre Dumas, les trois mousquetaires \'etaient quatre en r\'ealité''.} are solutions of a differential equation of the same type of \eqref{edo} (see \cite[Section 4.2]{Si15O}). However, only Jacobi and Laguerre polynomials depend on a real parameter and, therefore, the Stieltjes result is no longer applicable in this framework.

From a truly practical point of view ---for example in Physics, which was traditionally the birthplace of some of the most beautiful families (see \cite{VK91})---, it rarely will require more than the classical orthogonal polynomials. For classical (discrete) orthogonal polynomials on a uniform grid \footnote{A non-empty (totally) ordered set of equidistant (respectively, non-equidistant) points is called uniform (respectively, non-uniform)  grid, which in turn is an elementary example of lattice.} (Charlier, Krawtchouk, Hahn, and Meixner polynomials), Markov's theorem can be used (see \cite[Chapter $7$]{I05}) and, of course, for other discrete families \footnote{For a ``continuous'' case as the Askey-Wilson polynomials, Askey and Wilson used a consequence of Markov's theorem, which goes back to Szeg\H{o} (see \cite[Theorem 6.12.2]{S75}), to study the monotonicity of zeros of these polynomials (see \cite[Section $7$]{AW85}).}. Here we purpose an alternative approach, establishing a bridge between two works separated in time by almost one century, on one hand the work of Stieltjes and, on the other hand, the work \cite{NU83} by Nikiforov and Uvarov. This allows one to carry out a systematic study of the monotonicity of zeros of classical orthogonal polynomials on linear, quadratic, q-linear, and q-quadratic grids. Indeed, the purpose of this note is to prove that under suitable regularity conditions the Stieltjes result for Jacobi polynomials remains valid if we replace the differential equation \eqref{edo}, with $a$ and $b$ arbitrary polynomials of degree at most $2$ and $1$, respectively, and $c$ an arbitrary constant (from now on, whenever we refer to \eqref{edo}, we are assuming these conditions), by the following difference equation introduced in \cite[(5)]{NU83} (see also \cite[p. $127$]{NU84} and \cite[p. $71$]{NSU85}):
\begin{align}\label{edd}
a(X)\frac{\Delta}{\Delta x\left(s-1/2\right)}
\left(\frac{\nabla y(X)}{\nabla X}\right)
+\frac{b(X)}{2}\left(\frac{\Delta y(X)}{\Delta X}+\frac{\nabla y(X)}{\nabla X}\right)+c\,y(X)=0,
\end{align}
or, equivalently,
\begin{align}\label{edd2}
a(s)\frac{\Delta}{\Delta x\left(s-1/2\right)}
\left(\frac{\nabla y(X)}{\nabla X}\right)
+b(X)
\frac{\Delta y(X)}{\Delta X}+c\, y(X)=0,
\end{align}
where
\begin{align*}
a(s)=a(X)-\frac{1}{2}\,b(X)\Delta x\left(s-\frac12\right),
\end{align*}
$X=x(s)$\footnote{In first edition of ``{S}pecial {F}unctions of {M}athematical {P}hysics'' (see \cite{NU78}) the author only consider the case $x(s)=s$. The second edition (see \cite{NU84}) was significantly enriched with the equation \eqref{edd}; although A. A. Samarskii's preface is the same in both editions.}  defines a class of grids with, generally nonuniform, step-size $\Delta X=\Delta x(s)=x(s+1)-x(s)$ and $\nabla X=\nabla x(s)=\Delta x(s-1)$. (By abuse of notation, we use the same letter $a$ for the function $a(s)$ and the polynomial $a(X)$.) In what follows, we assume that $x$ is a real-valued function defined on an interval of the real line. Any solution of \eqref{edd} can be brought in correspondence with the solution of \eqref{edo} by replacing $s$ by $s/h$ and then taking limit $h \to 0$, whenever it exists. It is important to highlight that \eqref{edd} has polynomial solutions, in $X$, whose difference-derivatives satisfy equations of the same kind if and only if, for $q\neq1$ fixed, $x$ is a linear, quadratic, q-linear, or q-quadratic grid of the form
 \begin{align*}
x(s)=
\left\{
\begin{array}{ll}
C_1s^2+C_2s,\\[7pt]
C_3q^{-s}+C_4q^s,
\end{array}
\right.
\end{align*}
where $(C_1,C_2)\neq(0,0)$ and $(C_3,C_4)\neq(0,0)$ (see \cite[(1.68)]{ARS95}). The grids that depend on ``q'' are called q-linear if $C_3$ or $C_4$ is zero; otherwise it is q-quadratic. By using linear transformations (see \cite[(3.4.1)]{NSU}) we can reduce the expressions for the grids to simpler forms. In what follows, we assume that the grid $x$ takes on the following canonical forms (see \cite[p. 74]{NSU}):
 \begin{align}\label{grid}
x(s)=
\left\{
\begin{array}{lll}
s & \qquad \qquad (\mathrm{I})\\[7pt]
s(s+1) & \qquad \qquad (\mathrm{II})\\[7pt]
q^{s} \qquad \quad \quad \ \ \ \ \ \ \ (q>1) & \qquad \qquad (\mathrm{III})\\[7pt]
\displaystyle\frac12(q^s-q^{-s}) \qquad (q>1) & \qquad \qquad (\mathrm{IV}) \\[7pt]
\displaystyle\frac12(q^s+q^{-s}) \qquad (q>1) & \qquad \qquad (\mathrm{V})\\[7pt]
\displaystyle\frac12(q^s+q^{-s}) \qquad \, (q=e^{2i\theta},\, 0<\theta<\pi/2). & \qquad \qquad (\mathrm{VI})
\end{array}
\right.
\end{align}
Under the above notation, we adopt the following definition of classical discrete orthogonal polynomials on the real line (COPRL), which is enough for our purpose:

\begin{definition}\label{def1}
Fix $\mathrm{a}\in\mathbb{R}\cup\{-\infty\}$ and $N\in\mathbb{N}\cup \{\infty\}$ and define $\mathrm{b}=\mathrm{a}+N$. Fix $q$ and let $X=x(s)$ be a real-valued function given by \eqref{grid}, where the variable $s$ ranges over the finite interval $[\mathrm{a}, \mathrm{b}]$ or the infinity interval $[\mathrm{a}, \infty)$. A sequence of polynomials, $(P_n(X))_{n=0}^{N-1}$,  is said to be sequence of classical discrete orthogonal polynomials on the set $\{x(\mathrm{a}), x(\mathrm{a+1}), \dots, x(\mathrm{b-1})\}$ or, simply, COPRL if $:$\\
\begin{enumerate}
\item  $P_n$ satisfy \eqref{edd}, $x$ being a strictly  monotone function on $[\mathrm{a}, \mathrm{b}]$ or $[\mathrm{a}, \infty)$ given, up to a linear transformation, by \eqref{grid};\\[1pt]
\item   a positive weight function $\omega$ satisfying the boundary conditions \footnote{If $\mathrm{a}$ and $\mathrm{b}$ are finite, \eqref{bc} can be written in the form $\omega(\mathrm{a})a(\mathrm{a})=0$ and $\omega(\mathrm{b})a(\mathrm{b})=0$. If $\mathrm{a}=-\infty$ and/or $\mathrm{b}=\infty$, then \eqref{bc} must hold for each $k$ in the limiting sense.}
\begin{align}\label{bc}
\omega(s) a(s) x^k\left.\left(s-\frac12\right)\right|_{\mathrm{a}, \mathrm{b}}=0 \quad (k=0,1,\dots)
\end{align}
exists;
\item the difference equation
\begin{align}\label{pearson}
\frac{\Delta}{\Delta x \left(s-\displaystyle\frac12 \right)}\big(\omega(s)a(s)\big)=\omega(s)b(X)
\end{align}
holds.
\end{enumerate}
\end{definition}

Subsequently, when we say that a certain polynomial is a COPRL, we are assuming the definition and notation given in Definition \ref{def1}. From \eqref{edd2}, \eqref{bc}, and \eqref{pearson}, we conclude that the COPRL satisfy the orthogonality condition  (see \cite[(3.3.4)]{NSU})
\begin{align}\label{orth}
\sum_{s=\mathrm{a}}^{\mathrm{b}-1} P_m(X)P_n(X)\omega(s)\Delta x\left(s-\frac12\right)=0 \quad (m\not=n).
\end{align}
The polynomials $(P_n(X))_{n=0}^{N-1}$ given in Definition \ref{def1} are called simply discrete orthogonal polynomials, in $X$, on the set $\{x(\mathrm{a}), x(\mathrm{a+1}), \dots, x(\mathrm{b-1})\}$ with respect to a positive weight function $\omega$ if they satisfy the relation \eqref{orth} instead of the requirements $\mathrm{i})-\mathrm{iii})$. Of course, COPRL are a special case of discrete orthogonal polynomials.
From \eqref{orth}, we can see that the zeros of discrete orthogonal polynomials on $\{x(\mathrm{a}), x(\mathrm{a+1}), \dots, x(\mathrm{b-1})\}$ are real and distinct and are located in $(\min\{x(\mathrm{a}), x(\mathrm{b-1})\}, \max\{x(\mathrm{a}), x(\mathrm{b-1})\})$ (see \cite[Theorem 3.3.1]{S75}). In concluding this section we remark that it is possible to obtain a series representation of COPRL (see \cite[(4.19)]{ARS95} and \cite[Section 3]{NU93}):
\begin{align}\label{explicit}
&P_n(X)=(-1)^n \gamma_n \sum_{j=0}^n \displaystyle\binom{n}{j}_q \big(-q^{-n/2}\big)^j\frac{\Delta x\big(s-(n-1)/2+j\big)}{\displaystyle \prod_{k=0}^n \Delta x\big(s+(j-k+1)/2\big)}\\[7pt]
\nonumber&\times  \prod_{k=0}^n a(s-n+j+k)\prod_{l=0}^j \frac{a(x(s+l-1))+1/2\, b(x(s+l-1))\Delta x(s+l-1/2)}{a(x(s+l))-1/2\, b(x(s+l)) \Delta x(s+l+1/2)},
\end{align}
where $\gamma_n$ is a constant and
$$
\binom{n}{j}_q=\frac{(q; q)_n}{(q; q)_j (q; q)_{n-j}}.
$$
Furthermore (see \cite{NU93}), COPRL represent special cases of hypergeometric series or q-hypergeometric series, the latter formally defined by (see \cite[(10.9.4)]{AAR99})
$$
\pPq{i}{j}{\alpha_1,\dots, \alpha_i}{\beta_1\dots, \beta_j}{q, X}=\sum_{k=0}^\infty \frac{(\alpha_1,\dots, \alpha_i; q)_k}{(\beta_1\dots, \beta_j; q)_k}(-1)^{(1-i+j)k\displaystyle\binom{k}{2}}\frac{X^k}{(q; q)_k}.
$$

The outline of this note is as follows. In Section \ref{ext} we present an extension of the Stieltjes work in the framework of COPRL. In Section \ref{app} we study the variation of zeros of some families of COPRL according to the type of underlying grid.
\section{Main results}\label{ext}

Unless otherwise stated we assume that $a$, $b$, and $c$ appearing in \eqref{edd} depend on a parameter $t$ varying in a non-degenerate open interval of the real line. Rewrite  \eqref{edd} in the more suggestive form
\begin{align}\label{DE2}
A\, y(s-1)+B\, y(s+1)+C\, y(s)=0,
\end{align}
where $y(s)=y(X; t)$ and
\begin{align}
\label{AB} A=A(s; t)&=\frac{a(s; t)}{\nabla X\Delta x(s-1/2)}, \quad B=B(s; t)=\frac{a(s; t)+b(X; t)\Delta x(s-1/2)}{\Delta X\Delta x(s-1/2)},\\[7pt]
\nonumber C=C(s; t)&=c(t)-A(s; t)-B(s; t).
\end{align}
The following example will help motivate our main result.

\begin{eje}\label{eje1}
In 1960, Karlin and McGregor proved  (see \cite[(1.3)]{KM60} and \cite[(5.1)]{L67}) that the COPRL known as Hahn polynomials (see \cite[Section 9.5]{KLS10})
$$
y(s)=H^{(\alpha, \beta)}_n(X)=\pFq{3}{2}{-n, -X, \alpha+\beta+n+1}{\beta+1, 1-N}{1}
$$
$(n=1,\dots, N-1; \mathrm{a}=0, \mathrm{b}=N; \alpha>-1, \beta>-1)$, ---which constitute the finite discrete analogue of Jacobi polynomials considered by Stieltjes in \cite{S87}--- satisfy \eqref{DE2} with $X=x(s)=s$, and $A$ and $B$ given by
\begin{align*}
A=A(X; \alpha)=X(-X+\alpha+N),\quad B=B(X; \beta)=(X+\beta+1)(-X+N-1).
\end{align*}
The function
\begin{align}\label{S1}
\frac BA=\frac{(X+\beta+1)(-X+N-1)}{X(-X+\alpha+N)}
\end{align}
 is a positive and strictly decreasing function of $X \in(\mathrm{a}, \mathrm{b}-1)$, and
\begin{align}
\label{S21}\frac{\partial}{\partial \alpha}\left(\frac BA\right)&=\frac{(X+\beta+1)(X-N+1)}{X(-X+\alpha+N)^2}<0,\\[7pt]
\label{S22} \frac{\partial}{\partial \beta}\left(\frac BA\right)&=\frac{-X+N-1}{X(-X+\alpha+N)}>0.
\end{align}
for each $X \in(\mathrm{a}, \mathrm{b}-1)$. Since this example corresponds to the  grid $(\mathrm{I})$, from Markov's theorem, Ismail proves that the zeros of $H^{(\alpha, \beta)}_n$ are decreasing functions of $\alpha$ on $(-1, \infty)$ and increasing functions of $\beta$ on $(-1, \infty)$ (see \cite[Theorem $7.1.2$]{I05}).  Comparison of \eqref{s1} and \eqref{S1}, and \eqref{s2} and \eqref{S21}-\eqref{S22} suggests that, as for Jacobi polynomials, the information on the monotonicity of the zeros of Hahn polynomials is stored in the rational function $B/A$. Therefore, it is not unreasonable to conjecture that the same happens with any COPRL.
\end{eje}

The following two mathematical objects play a central role in our exposition.

\begin{definition}
Let $A$ and $B$ be given by \eqref{AB}. The function $f$, from now on called monotonicity function, is defined by
\begin{align}\label{f}
f(s; t)=\frac{B(s; t)}{A(s; t)}=\frac{a(s;t)+b(X; t)\Delta x(s-1/2)}{a(s; t)}\frac{\nabla X}{\Delta X}.
\end{align}
\end{definition}

\begin{definition}
Let $X=x(s)$ be given by \eqref{grid} and let $x(y_j(t))$ $(j=1,\dots,n)$ be the zeros of a COPRL of degree $n$, say $P_n(X; t)$, depending on a parameter $t$ taking values on $\mathrm{J}\subseteq \mathbb{R}$. For any set $\mathrm{I}\subseteq \mathbb{R}$, the (nonempty) subset $\mathrm{S}_\mathrm{I}^{(t)}(P)$ of $(\mathrm{a}, \mathrm{b}-1)$ is defined by
$$
\mathrm{S}_\mathrm{I}^{(t)}(P_n)=\Big\{y\in (\mathrm{a}, \mathrm{b}-1)\,\big|\, t\in \mathrm{I}\cap \mathrm{J}\,\wedge\, \big(\forall j\in \{1,\dots,n \}\big)\big[y=y_j(t)\big]\Big\}.
$$
We will write it simply $\mathrm{S}_\mathrm{I}(P_n)$, $\mathrm{S}_\mathrm{I}^{(t)}$ or $\mathrm{S}_\mathrm{I}$ when no confusion can arise\footnote{Note that $\mathrm{S}_\mathrm{J}=\mathrm{S}_\mathbb{R}$ and  $\mathrm{S}_\mathrm{I}\subseteqq\mathrm{S}_\mathrm{J}$ whenever $\mathrm{I}\subseteqq \mathrm{J}$.}.
\end{definition}

To prove our main result we need two lemmas. The first one was proved for the Hahn polynomials by Levit  (see \cite[Theorem $3$]{L67}). Here we reproduce in a more general framework, {\em mutatis mutandis}, his arguments.

\begin{lemma}\label{lemma}
Let $X=x(s)$ be given by \eqref{grid}. Let $x(y(t))$ and $x(z(t))$ be consecutive zeros of a COPRL depending on a real parameter $t$. Let $f$ be the monotonicity function given by \eqref{f}. Then $|z(t)-y(t)|>1$ for those values of $t$ such that $f(\cdot; t)>0$ on $\mathrm{S}_{\mathbb{R}}$.
\end{lemma}
\begin{proof}
We give the proof only for the case in which $x$ is a strictly increasing function. The same arguments apply to the case in which $x$ is a strictly decreasing function. There is no loss of generality in assuming that $x(y(t))<x(z(t))$ are the greatest pair of consecutive zeros of a COPRL, say $P(X; t)$. (Recall that the zeros of $P$ are distinct and are located in $(x(\mathrm{a}), x(\mathrm{b-1}))$. Moreover, under our assumptions, $y(t)<z(t)$.) Suppose that there exists $t_0$ such that $z(t_0)=y(t_0)+1$. Replacing $s$ by $z(t_0)$ and $t$ by $t_0$ in (\ref{DE2}) we have $P(x(z(t_0)+1);t_0)=0$, because $f(z(t_0); t_0)\neq0$, which contradicts the assumption that $x(z(t_0))$ is the greatest zero of $P(\cdot; t_0)$. (Here we have also used the fact that $x$ is, in particular, a strictly increasing function on $[\mathrm{a}, \mathrm{b})$, and not only on $[\mathrm{a}, \mathrm{b}-1]$, because the possibility that $z(t_0)+1\in [\mathrm{b}-1,\mathrm{b})$ is not excluded.) Suppose now that $z(t_0)<y(t_0)+1$. Replacing again $s$ by $z(t_0)$ and $t$ by $t_0$ in (\ref{DE2}) we have
$$
-\frac{P(x(z(t_0)-1); t_0)}{P(x(z(t_0)+1); t_0)}=f(z(t_0); t_0)>0.
$$
Since $z(t_0)\not=y(t_0)+1$, there is an odd number of zero of $P(\cdot; t_0)$ in $[x(z(t_0)-1),x(z(t_0)+1)]$. We next claim that there exists at least one integer $m$ $(1\leq m\leq N-1)$ such that $y(t_0)\leq \mathrm{a}+m \leq z(t_0)$ or, what is the same, $x(y(t_0)) \leq x(\mathrm{a}+m) \leq x(z(t_0))$. Suppose the assertion were false. Hence
$$
\sum_{s=\mathrm{a}}^{\mathrm{b}-1} \frac{\displaystyle P^2(X; t_0)}{\big(X-x(y(t_0))\big)\big(X-x(z(t_0))\big)}\omega(s; t_0)\Delta x\left(s-\frac12\right)>0,
$$
a contradiction with \eqref{orth}. Since $z(t_0)<y(t_0)+1$, $m$ is unique. Thus the only zeros of $P(\cdot; t_0)$ in $[x(z(t_0)-1),x(z(t_0)+1)]$ are $x(y(t_0))$ and $x(z(t_0))$. This contradicts the fact that there is an odd number of zeros of $P(\cdot; t_0)$ in $[x(z(t_0)-1),x(z(t_0)+1)]$, and the lemma follows.
\end{proof}

\begin{lemma}\label{Lemma}
Let $X=x(s)$ be given by \eqref{grid}. Let $P$ be a non-constant COPRL, in $X$, satisfying \eqref{DE2} and depending on a parameter $t$ varying in a non-degenerate open interval of the real line containing $t_0$. Suppose that $A(\cdot; t)$ and $B(\cdot; t)$ given by \eqref{AB} admit partial derivatives with respect to $t$ on a neighbourhood of $t_0$. Assume that $P(Y_0; t_0)=0$. Then there exist $\epsilon>0$ and $\delta>0$ such that  $(Y_0-\delta, Y_0+\delta) \times (t_0-\epsilon, t_0+\epsilon)$  is on the neighbourhood where $P$ is defined, and there exists $Y: (t_0-\epsilon, t_0+\epsilon) \to (Y_0-\delta, Y_0+\delta)$, such that
\begin{align}\label{zero}
P(Y(t); t)=0
\end{align}
and, for each $t \in (t_0-\epsilon, t_0+\epsilon)$, $Y$ is the unique solution of \eqref{zero} with $Y(t)\in (Y_0-\delta, Y_0+\delta)$. Moreover, $Y$  possess a continuous derivative on $(t_0-\epsilon, t_0+\epsilon)$.
\end{lemma}
\begin{proof}
From \eqref{explicit} we see that the coefficients of $P(\cdot; t)$ are differentiable functions of $t$. Moreover, $P(Y_0; t_0)=0$; from this it follows that
$$
\left. \frac{\partial P}{\partial X}(X; t)\right|_{X=Y_0, t=t_0}\not=0,
$$
because the zeros of $P(\cdot; t_0)$ are distinct. Thus, the result is a direct consequence of the implicit function theorem (see \cite[Theorem 3.4.2]{Si15}).
\end{proof}

We shall refer to Theorem \ref{main} below as discrete Stieltjes theorem.

\begin{theorem}\label{main}
Assume the hypotheses and notation of Lemma \ref{Lemma}. Let $f$ be the monotonicity function given by \eqref{f}. Denote $f_1(s; \cdot)=(\partial f/\partial s) (s;\cdot)$ and $f_2(\cdot;t)=(\partial f/\partial t)(\cdot;t)$. Suppose that $f(s;t)>0$ and $f_1(s;t)<0$ on $\mathrm{S}_{(t_0-\epsilon, t_0+\epsilon)}$. In the case of the grid $(\mathrm{IV})$, assume also that $n f(s; \cdot)+f_1(s; \cdot)\leq 0$. Suppose furthermore that $f_2(s;t)>0$ (respectively, $f_2(s;t)<0$) on $\mathrm{S}_{(t_0-\epsilon, t_0+\epsilon)}$. Then $Y$ is a strictly increasing (respectively, decreasing) function on $(t_0-\epsilon, t_0+\epsilon)$ if $x$ is a strictly increasing function, or else $Y$ is a strictly decreasing (respectively, increasing) function on $(t_0-\epsilon, t_0+\epsilon)$.
\end{theorem}
\begin{proof}
We give the proof only for the case in which $x$ is a strictly increasing function. The same arguments apply to the case in which $x$ is a strictly decreasing function. Assume that $P$ is monic and has fixed degree $n$. Let $Y_j$ $(j=1,\dots,n)$ denote the zeros of $P$. Since the indeterminate of $P$ takes values on the real (open) interval $x((\mathrm{a}, \mathrm{b}-1))$, there exist functions $y_j$ defined on a neighbourhood of $t_0$ and taking values on $(\mathrm{a}, \mathrm{b}-1)$ such that $Y_j(t)=x(y_j(t))$. Moreover, since $x$ is strictly  increasing on $(\mathrm{a}, \mathrm{b}-1)$, $x^{-1}$ is differentiable on $x((\mathrm{a}, \mathrm{b}-1))$. Hence, by Lemma \ref{Lemma}, there exists a neighbourhood of $t_0$ where $y_j$ is differentiable.

Let
\begin{align*}
P(X; t)=\prod_{j=1}^{n}(X-Y_j(t))
\end{align*}
defined in a neighbourhood of $t_0$. Replacing $s$ by $y_j(t)$ in (\ref{DE2}) we have
\begin{align}\label{eq1}
f(y_j(t); t)=-\frac{P(x(y_j(t)-1); t)}{P(x(y_j(t)+1); t)}.
\end{align}
Taking the partial derivative of \eqref{eq1} with respect to $t$ on the neighbourhood of $t_0$ leads to
\begin{align*}
&y'_j(t)f_1(y_j(t); t)+f_2(y_j(t); t)\\[7pt]
\nonumber &\quad=\frac{P(x(y_j(t)-1);t)\displaystyle \frac{\partial P}{\partial t}(x(y_j(t)+1); t)-P(x(y_j(t)+1);t)\displaystyle \frac{\partial P}{\partial t}(x(y_j(t)-1); t) }{P^2(x(y_j(t)+1);t)},
\end{align*}
where
\begin{align*}
&\frac{\partial P}{\partial t}(x(y_j(t)\pm1); t)\\[7pt]
\nonumber &\quad=P(x(y_j(t)\pm1); t)\displaystyle \sum_{k=1}^n \frac{\left.\displaystyle\frac{\mathrm{d}x}{\mathrm{d}s}(s\pm1)\right|_{s=y_j(t)} y'_j(t)-  \left.\displaystyle\frac{\mathrm{d}x}{\mathrm{d}s}(s)\right|_{s=y_k(t)} y'_k(t)}{x(y_j(t)\pm1)-Y_k(t)}.
\end{align*}
Hence
\begin{align}\label{f2}
f_{2}(y_j(t); t)=\sum\limits_{k=1}^{n} a_{j k}(t)y'_k(t),
\end{align}
where
\begin{align}
\label{aii}a_{jj}(t)=-f_{1}(y_j(t); t)+f(y_j(t); t)\sum_{k=1}^{n}b_{j k}(t)-\sum\limits_{\substack{k=1 \\ k\neq j}}^{n}a_{jk},
\end{align}
with
\begin{align*}
b_{jk}(t)= \frac{\left.\displaystyle\frac{\mathrm{d}x}{\mathrm{d}s}(s-1)\right|_{s=y_j(t)}-\left.\displaystyle\frac{\mathrm{d}x}{\mathrm{d}s}(s)\right|_{s=y_k(t)}}{x(y_j(t)-1)-Y_k(t)}-\frac{\left.\displaystyle\frac{\mathrm{d}x}{\mathrm{d}s}(s+1)\right|_{s=y_j(t)}-  \left.\displaystyle\frac{\mathrm{d}x}{\mathrm{d}s}(s)\right|_{s=y_k(t)}}{x(y_j(t)+1)-Y_k(t)},
\end{align*}
and $a_{jk}(t)=f(y_j(t); t)\,c_{jk}(t)$ $(j\not=k)$, where
\begin{align*}
c_{jk}(t)=\left(\frac{1}{x(y_j(t)+1)-Y_k(t)}-\frac{1}{x(y_j(t)-1)-Y_k(t)}  \right)\left.\displaystyle\frac{\mathrm{d}x}{\mathrm{d}s}(s)\right|_{s=y_k(t)}.
\end{align*}
We claim that $a_{jk}(t)<0$ $(j\not=k)$. From now on, without loss of generality, we assume $Y_j(t)<Y_k(t)$. Then, by Lemma \ref{lemma}, $x(y_j(t)+1)<Y_k(t)$. Since $\left.(\mathrm{d}x/\mathrm{d}s)(s)\right|_{s=y_k(t)}$ $>0$, $a_{jk}(t)<0$ as required. Consider now the canonical forms \eqref{grid} to conclude that $a_{jj}(t)>0$. We next claim that, for the grids $(\mathrm{I})-(\mathrm{III})$ and $(\mathrm{V})-(\mathrm{VI})$, $b_{jk}(t)>0$. Indeed, for the grids $(\mathrm{I})$ and $(\mathrm{III})$, the proof is straightforward. For the grid $(\mathrm{II})$, we have
$$
b_{jk}(t)=\frac{4}{\big(y_j(t)+y_k(t)\big)\big(y_j(t)+y_k(t)+2\big)}
$$
and $y_j(t)\in \mathrm{S}_{(t_0-\epsilon, t_0+\epsilon)}\subset (-1/2, \infty)$, the latter because $x$ is strictly  increasing on $(-1/2, \infty)$. Hence, by Lemma \ref{lemma}, $b_{jk}(t)>0$ as claimed. Since $f(y_j(t); t)>0$, $f_{1}(y_j(t); t)<0$, $a_{jk}(t)<0$, and $b_{jk}(t)\geq 0$, \eqref{aii} implies $a_{jj}(t)>0$. The corresponding result for the grid $(\mathrm{V})$ follows in the same way after noting that
$$
b_{jk}(t)=2\theta \sinh (2 \theta) \csch\big((y_k(t)+y_j(t)-1)\theta\big)\csch\big((y_k(t)+y_j(t)+1)\theta \big),
$$
for $q=e^{2\theta}$ ($\theta>0$), and $y_j(t)\in \mathrm{S}_{(t_0-\epsilon, t_0+\epsilon)}\subset (0, \infty)$. And the same goes for  the grid $(\mathrm{VI})$ using hyperbolic identities and noting that $y_j(t)\in \mathrm{S}_{(t_0-\epsilon, t_0+\epsilon)}\subset (-\pi/(2\theta), 0)$. Define
\begin{align*}
&\quad\mathbf{f}(t)=(f_2(y_1(t); t), \dots, f_2(y_n(t); t))^\mathsf{T},\\[7pt]
\mathbf{A}&(t)=(a_{jk}(t)),\quad \mathbf{y}(t)=\big(y'_1(t),\dots, y'_n(t) \big)^\mathsf{T},
\end{align*}
and rewrite \eqref{f2} as $\mathbf{f}(t)=\mathbf{A}(t)\mathbf{y}(t)$. Observe that $\mathbf{A}(t)$ has positive diagonal entries and negative of off-diagonal entries. Moreover, from \eqref{aii} we get
$$
 |a_{jj}(t)|>\sum\limits_{\substack{k=1 \\ k\neq j}}^{n}|a_{jk}(t)|.
$$
Hence $\mathbf{A}(t)$ is a real irreducibly diagonally dominant matrix, and so, by \cite[Corollary 1, p. 85]{V62} all the entries of $\mathbf{A}^{-1}(t)$ are positive \footnote{Indeed, $\mathbf{A}(t)$ is a Stieltjes matrix.}. Thus all the entries of $\mathbf{y}(t)=\mathbf{A}^{-1}(t) \mathbf{f}(t)$ are positive, and the theorem is proved for the grids $(\mathrm{I})-(\mathrm{III})$ and $(\mathrm{V})-(\mathrm{VI})$. The above argument does not work for the grid $(\mathrm{IV})$, because, for $q=e^{2\theta}$ ($\theta>0$),
$$
b_{jk}(t)=-2\theta \sinh(2\theta)\sech\big((y_j(t)+y_k(t)+1)\theta\big)\sech\big((y_j(t)+y_k(t)-1)\theta\big)<0.
$$
Indeed, $-1<b_{jk}(t)<0$ and, therefore, under the additional hypothesis for this grid, the theorem follows from \eqref{aii}.
\end{proof}

\begin{obs}
In some cases, the hypotheses of the discrete Stieltjes theorem are fulfilled in $(\mathrm{a}, \mathrm{b}-1)\supset  \mathrm{S}_\mathbb{R}$. However, there are several, and important, examples where this is only true on a subset of $(\mathrm{a}, \mathrm{b}-1)$ containing all the elements of $\mathrm{S}_\mathbb{R}$; see, for instance, the Racah polynomials in Section \ref{latticeII} below. Of course, since the precise location of the zeros of $P$ is not known, these cases require a more careful analysis.
\end{obs}

In Example \ref{eje1}, we can write $H^{(\alpha, \beta, N)}_n$ instead of $H^{(\alpha, \beta)}_n$. The interlacing property between the zeros of $H^{(\alpha, \beta, N)}_n$ and $H^{(\alpha, \beta, N+1)}_n$ was observed by Levit (see \cite[Theorem 6]{L67}). It was proved independently, and almost simultaneously, by Mesztenyi that  \cite[Theorem 6]{L67} is always true for discrete orthogonal polynomials on the grid $(\mathrm{I})$ when the variable $s$ ranges over a finite interval (see \cite[Lemma 3]{Ma66}). In this way, the following general result might be of interest to the reader.

\begin{theorem}\label{N}
Fix $\mathrm{a}\in\mathbb{R}$ and $N\in\mathbb{N}$ and define $\mathrm{b}=\mathrm{a}+N$. Let $X=x(s)$ be a real-valued function, where the variable $s$ ranges over the finite interval $[\mathrm{a}, \mathrm{b}+1]$.  Let $(P_n(X; N))_{n=1}^{N-1}$ be the sequence of discrete orthogonal polynomials, in $X$, on the set $\{x(\mathrm{a}), x(\mathrm{a}+1), \dots,$ $ x(\mathrm{b}-1)\}$. Let $(P_n(X; N+1))_{n=1}^{N}$ be the sequence of discrete orthogonal polynomials, in $X$, on the set $\{x(\mathrm{a}), x(\mathrm{a}+1), \dots, x(\mathrm{b})\}$. Assume that both sequence are orthogonal with respect to the same positive weight function. Let $Y_1<Y_2<\cdots <Y_n$ be the zeros of $P_n(\cdot; N)$. Then one of the following situations holds:
\begin{enumerate}
\item The zeros of $P_n(\cdot; N+1)$ are those of $P_n(\cdot; N)$, if $P_n(x(\mathrm{b}); N)=0$;
\item  $P_n(\cdot; N+1)$ has a zero on $(Y_k, Y_{k+1})$, if $P_n(x(\mathrm{b}); N)\not=0$  and $x(\mathrm{b})\not\in (Y_k, Y_{k+1})$ for fixed $k\in \{1,\dots, n-1\}$;
\item  $P_n(\cdot; N+1)$ has a zero on each of the intervals $(Y_k, x(\mathrm{b}))$ and $(x(\mathrm{b}), Y_{k+1})$, if $P_n(x(\mathrm{b}); N)\not=0$ and $x(\mathrm{b})\in (Y_k, Y_{k+1})$ for fixed $k\in \{1,\dots, n-1\}$.
\end{enumerate}
\end{theorem}
\begin{proof}
To shorten notation, write $P_n^{(N)}$ and $P_n^{(N+1)}$ instead of $P_n(\cdot; N)$ and $P_n(\cdot; N+1)$, respectively. The reader may check for himself that by expressing $P_n^{(N+1)}$ as a linear combination of the elements of the set $\big\{1, P_1^{(N)}, \dots, P_{n}^{(N)}\big\}$ and using the orthogonality property, we obtain $P_n^{(N+1)}(x(\mathrm{b}))=\zeta_n\, P_n^{(N)}(x(\mathrm{b}))$ and
\begin{align}\label{linearcomb}
P_n^{(N+1)}(X)=& P_n^{(N)}(X)\\[7pt]
\nonumber&-\eta_n P_n^{(N+1)}(x(\mathrm{b}))\frac{P_{n-1}^{(N)}(x(\mathrm{b})) P_n^{(N)}(X)-P_n^{(N)}(x(\mathrm{b}))P_{n-1}^{(N)}(X)}{X-x(\mathrm{b})},
\end{align}
where $\eta_n=\omega(\mathrm{b})\Delta x(\mathrm{b}-1/2)\big\|P_{n-1}^{(N)}\big\|^{-2}$ and
\begin{align*}
\zeta_n&=1+\eta_n\left( P_{n-1}^{(N)}(x(\mathrm{b}))  \displaystyle \left.\displaystyle\frac{\mathrm{d}P_n^{(N)}}{\mathrm{d}X}(X)\right|_{X=x(\mathrm{b})}-\displaystyle \left.\displaystyle\frac{\mathrm{d}P_{n-1}^{(N)}}{\mathrm{d}X}(X)\right|_{X=x(\mathrm{b})} P_n^{(N)}(x(\mathrm{b}))\right).
\end{align*}
Following a standard procedure, the rest of the proof follows as \cite[Lemma 3]{Ma66}.
\end{proof}

We end this section with the following consequence of Theorem \ref{N}, which is valid, in particular, for COPRL.

\begin{coro}\label{CN}
Assume the hypotheses and notation of Theorem \ref{N}. Suppose furthermore that $x$ is a strictly monotone function on $[\mathrm{a}, \mathrm{b}+1]$. Set $Y_0=Y_{n+1}=x(\mathrm{b})$. Then $P_n(\cdot; N+1)$ has a zero on each interval $(Y_k, Y_{k+1})$ for all $k\in \{1,\dots, n\}$  if $x$ is a strictly increasing function, or else $P_n(\cdot; N+1)$ has a zero on each interval $(Y_k, Y_{k+1})$ for all $k\in \{0,\dots, n-1\}$.
\end{coro}
\begin{proof}
We give the proof only for the case in which $x$ is a strictly increasing function. The same arguments apply to the case in which $x$ is a strictly decreasing function. From \eqref{linearcomb}, it follows that
$$
\mathrm{sgn}\left(P_n(Y_n; N+1)\right)=-\mathrm{sgn}\left(P_n(Y_{n+1}; N)\right)=-1.
$$
Hence $P_n(\cdot; N+1)$ has a zero on $(Y_n, Y_{n+1})$. The rest of the proof is a direct consequence of Theorem \ref{N} $\mathrm{ii})$.
\end{proof}

\section{Applications}\label{app}
In this section we apply the discrete Stieltjes theorem to specific families of COPRL. (Of course, Corollary \ref{CN} is applicable to any family of COPRL on a finite grid, see for instance the Hahn, Krawtchouk, Racah, dual Hahn, q-Hahn, q-Krawtchouk, affine q-Krawtchouk, quantum q-Krawtchouk, q-Racah, and dual q-Hanh polynomials below.) We prove, or sketch the proof in similar cases, only of some illustrative examples; the other cases are stated and the proofs are left as exercises for the reader. The reader also should satisfy himself that the hypotheses of Lemma \ref{Lemma} are fulfilled. As far as we know, only the monotonicity of zeros of COPRL on the grid $(\mathrm{I})$ has been studied (see \cite[Chapter $7$]{I05}) \footnote{Recall that  for a continuous case on the grid $(\mathrm{VI})$  as the Askey-Wilson polynomials, the monotonicity of their zeros was studied in \cite[Section 7]{AW85}.}. We include this case for the sake of completeness.

\subsection{The grid $(\mathrm{I})$}\label{latticeI}
\subsubsection*{Examples of COPRL on $X=x(s)=s$ (Hahn, Charlier, Krawtchouk, and Meixner polynomials).}

The {\em Hahn polynomials}, $H^{(\alpha, \beta)}_n$, are defined in Example \ref{eje1}.
\begin{proposition}
The zeros of $H^{(\alpha, \beta)}_n$ are strictly decreasing functions of $\alpha$ on $(-1,\infty)$ and strictly increasing functions of $\beta$ on $(-1,\infty)$.
\end{proposition}
\begin{proof}
It suffices to use \eqref{S1} and \eqref{S21}-\eqref{S22} together with the discrete Stieltjes theorem.
\end{proof}

The {\em Charlier polynomials} (see \cite[Section 9.14]{KLS10}),
$$
y(s)=C_n^{(\alpha)}(X)=\pFq{2}{0}{-n, -X}{-}{-\frac{1}{\alpha}}
$$
$(n=1,2,\dots; \mathrm{a}=0, \mathrm{b}=\infty; \alpha>0)$, satisfy \eqref{DE2} with $A$ and $B$ given by
\begin{align*}
A=A(X)=X, \quad B=B(\alpha)=\alpha.
\end{align*}

\begin{proposition}
The zeros of $C_n^{(\alpha)}$ are strictly increasing functions of $\alpha$ on $(0, \infty)$.
\end{proposition}

The {\em Krawtchouk polynomials} (see \cite[Section 9.11]{KLS10}),
$$
y(s)=K_n^{(\alpha)}(X)=\pFq{2}{1}{-n, -X}{1-N}{\frac{1}{\alpha}}
$$
$(n=1,\dots, N-1; \mathrm{a}=0, \mathrm{b}=N; 0<\alpha<1)$, satisfy \eqref{DE2} with $A$ and $B$ given by
\begin{align*}
A=A(X; \alpha)=(1-\alpha)X,\quad B=B(X; \alpha)=\alpha(-X+N-1).
\end{align*}

\begin{proposition}
The zeros of $K_n^{(\alpha)}$ are strictly increasing functions of $\alpha$ on $(0, 1)$.
\end{proposition}

The {\em Meixner polynomials} (see \cite[Section 9.10]{KLS10}),
$$
y(s)=M_n^{(\alpha, \beta)}(X)=\pFq{2}{1}{-n, -X}{\beta}{1-\frac{1}{\alpha}}
$$
$(n=1,2,\dots; \mathrm{a}=0, \mathrm{b}=\infty; 0<\alpha<1, \beta>0)$, satisfy \eqref{DE2} with $A$ and $B$ given by
\begin{align*}
A=A(X)=X, \quad B=B(X; \alpha, \beta)=\alpha(X+\beta).
\end{align*}

\begin{proposition}
The zeros of $M_n^{(\alpha, \beta)}$ are strictly increasing functions of $\alpha$ on $(0, 1)$ and strictly increasing functions of $\beta$ on $(0,\infty)$.
\end{proposition}
\subsection{The grid $(\mathrm{II})$}\label{latticeII}
\subsubsection*{Examples of COPRL on $X=x(s)=s(s+1)$ (Racah and dual Hahn polynomials).}

The {\em Racah polynomials} (see \cite[p. 236]{NU93}),
\begin{align}\label{FRacah}
y(s)=R_n^{(\alpha, \beta)}(X)&=\pFq{4}{3}{-n, \alpha+\beta+n+1, \mathrm{a}-s, s+\mathrm{a}+1}{2\mathrm{a}+\alpha+N+1, \beta+1, 1-N}{1}
\end{align}
$(n=1,\dots, N-1; \mathrm{a}>-1/2, \mathrm{b}=\mathrm{a}+N; \alpha>-1, -1<\beta<2\mathrm{a}+1)$, satisfy \eqref{DE2} with $A$ and $B$ given by
\begin{align*}
A&=A(s; \alpha, \beta)=\frac{(s-\mathrm{a})(s+\mathrm{a}+N)(s-\mathrm{a}-\alpha-N)(s+\mathrm{a}-\beta)}{2s(2s+1)},\\[7pt]
B&=B(s; \alpha, \beta)=\frac{(s+\mathrm{a}+1)(s-\mathrm{a}-N+1)(s+\mathrm{a}+\alpha+N+1)(s-\mathrm{a}+\beta+1)}{2(s+1)(2s+1)}.
\end{align*}

\begin{proposition}\label{Racah}
The zeros of $R_n^{(\alpha,\beta)}$ are strictly decreasing functions of $\alpha$ on $(-1, \infty)$ and strictly increasing function of $\beta$ on $(-1, 2\mathrm{a}+1)$ if $\mathrm{a}\geq 0$, or else the zeros of $R_n^{(\alpha,\beta)}$ are strictly decreasing functions of $\alpha$ on $(-1, \infty)$ for each $\beta\in (\mathrm{a}, 2\mathrm{a}+1)$ and strictly increasing function of $\beta$ on $(\mathrm{a}, 2\mathrm{a}+1)$.
\end{proposition}
\begin{proof}
We give the proof only for the case in which $\mathrm{a}\geq 0$. The proof for $-1/2<\mathrm{a}<0$ is similar. Define the interval $\mathrm{K}=\big(\max\big\{\mathrm{a}, \beta-\mathrm{a}\big\}, \mathrm{a}+N-1\big)$.
The monotonicity function
\begin{align}\label{fR}
f=\frac{B}{A}=\frac{s(s+\mathrm{a}+1)(s-\mathrm{a}-N+1)(s+\mathrm{a}+\alpha+N+1)(s-\mathrm{a}+\beta+1)}{(s+1)(s-\mathrm{a})(s+\mathrm{a}+N)(s-\mathrm{a}-\alpha-N)(s+\mathrm{a}-\beta)}
\end{align}
is a positive and strictly decreasing function of $s \in\mathrm{K}$, and
\begin{align}
\label{DfR1}\frac{\partial f}{\partial \alpha}&=\frac{s(2s+1)(s+\mathrm{a}+1)(s-\mathrm{a}-N+1)(s-\mathrm{a}+\beta+1)}{(s+1)(s-\mathrm{a})(s+\mathrm{a}+N)(s-\mathrm{a}-\alpha-N)^2(s+\mathrm{a}-\beta)}<0,\\[7pt]
\label{DfR2}\frac{\partial f}{\partial \beta}&=\frac{s(2s+1)(s+\mathrm{a}+1)(s-\mathrm{a}-N+1)(s+\mathrm{a}+\alpha+N+1)}{(s+1)(s-\mathrm{a})(s+\mathrm{a}+N)(s-\mathrm{a}-\alpha-N)(s+\mathrm{a}-\beta)^2}>0,
\end{align}
for each $s \in\mathrm{K}$. It is immediate that $\mathrm{S}_{(-1,2\mathrm{a}]}^{(\beta)}\subset \mathrm{K}$. We next claim that $\mathrm{S}_{(2\mathrm{a}, 2\mathrm{a}+1)}^{(\beta)}\subset \mathrm{K}$. Indeed, this is equivalent to prove that $R_n^{(\alpha,\beta)}$ has no zeros on $\big(\mathrm{a}(\mathrm{a}+1), (\beta-\mathrm{a})(\beta-\mathrm{a}+1)\big)$ for $\beta \in (2\mathrm{a}, 2\mathrm{a}+1)$. For $\beta \in (2\mathrm{a}, 2\mathrm{a}+1)$, $R_n^{(\alpha, \beta)}(\mathrm{a}(\mathrm{a}+1))=1$ and
\begin{align*}
 R_n^{(\alpha, \beta)}\big((\beta-\mathrm{a})(\beta-\mathrm{a}+1)\big)&=\pFq{3}{2}{-n, \alpha+\beta+n+1, 2\mathrm{a}-\beta}{2\mathrm{a}+\alpha+N+1, 1-N}{1}\\[7pt]
 &=\frac{(2\mathrm{a}-\beta+N-n)_n(\alpha+\beta+N+1)_n}{(2\mathrm{a}+\alpha+N+1)_n(N-n)_n}>0,
\end{align*}
the last equality being a consequence of Sheppard's identity (see \cite[Corollary 3.3.4]{AAR99}). Hence $R_n^{(\alpha, \beta)}$ has no zeros on $\big(\mathrm{a}(\mathrm{a}+1), (\beta-\mathrm{a})(\beta-\mathrm{a}+1)\big)$ or has at least two zeros there. Suppose the second of these possibilities is true. From the proof of Lemma \ref{lemma}, we have $x(\mathrm{a}+1)<x(\beta-\mathrm{a})$, which is impossible. Thus $\mathrm{S}_{(2\mathrm{a}, 2\mathrm{a}+1)}^{(\beta)}\subset \mathrm{K}$ as claimed. The same proof actually shows that $\mathrm{S}_{(-1,\infty)}^{(\alpha)}\subset \mathrm{K}$. The result follows from the discrete Stieltjes theorem.
\end{proof}

\begin{obs
In \cite[Section 9.2]{KLS10}, the Racah polynomials (interchanging $\alpha$ and $\beta$) are defined by
$$
\widetilde{R}_n^{(\alpha, \beta, \delta)}(\widetilde{x}(s))=\pFq{4}{3}{-n, \alpha+\beta+n+1, -s, s+\delta-N+1}{\alpha+\delta+1, \beta+1, 1-N}{1}
$$
$(n=1,\dots, N-1; \mathrm{a}=0, \mathrm{b}=N; \alpha>-1, -1<\beta<\delta-N+1, \delta>N-1)$, where $\widetilde{x}(s)=s(s+\delta-N+1)$.  In \eqref{FRacah}, we can write $R^{(\alpha, \beta, \mathrm{a})}_n$ instead of $R^{(\alpha, \beta)}_n$; although, implicitly, we have the agreement to omit those parameters that are assumed fixed. Clearly, $\widetilde{x}$ is not one of the canonical forms \eqref{grid}. However, $\widetilde{x}$ and $x$ are related by a linear transformation. Therefore, for  $\mathrm{a}=(\delta-N)/2$ fixed,
$$
R^{(\alpha, \beta, (\delta-N)/2)}_n(x(s))=\widetilde{R}_n^{(\alpha, \beta, \delta)}\left(\widetilde{x}\left(s-\frac{\delta-N}{2}\right)\right),
$$
where
$$
\widetilde{x}\left(s-\frac{\delta-N}{2}\right)=x(s)-\frac{(\delta-N+1)^2-1}{4}.
$$
Consequently, Proposition \ref{Racah} remains valid if we replace $R^{(\alpha, \beta)}_n$ by $\widetilde{R}_n^{(\alpha, \beta, \delta)}$.
\end{obs}

The {\em dual Hahn polynomials} (see \cite[p. 236]{NU93}),
\begin{align}\label{DualHahneq}
y(s)=W_n^{(\alpha)}(X)&=\pFq{3}{2}{-n, \mathrm{a}-s, s+\mathrm{a}+1}{\alpha+1, 1-N}{1}
\end{align}
$(n=1,\dots, N-1; \mathrm{a}>-1/2, \mathrm{b}=\mathrm{a}+N; -1<\alpha<2\mathrm{a}+1)$, satisfy \eqref{DE2} with $A$ and $B$ given by
\begin{align*}
A&=A(s; \alpha)=\frac{(s-\mathrm{a})(s+\mathrm{a}+N)(s+\mathrm{a}-\alpha)}{2s(2s+1)},\\[7pt]
B&=B(s; \alpha)=\frac{(s+\mathrm{a}+1)(-s+\mathrm{a}+N-1)(s-\mathrm{a}+\alpha+1)}{2(s+1)(2s+1)}.
\end{align*}

\begin{proposition}\label{DualHahn}
The zeros of $W_n^{(\alpha)}$ are strictly increasing functions of $\alpha$ on $(-1, 2\mathrm{a}+1)$ if $\mathrm{a}\geq 0$, or else the zeros of $W_n^{(\alpha)}$ are strictly increasing functions of $\alpha$ on $(\mathrm{a}, 2\mathrm{a}+1)$.
\end{proposition}
\begin{proof}
We sketch the proof only for the case in which $\mathrm{a}\geq 0$. The proof for $-1/2<\mathrm{a}<0$ is similar. Define the interval $\mathrm{K}=\big(\max\big\{\mathrm{a}, \alpha-\mathrm{a}\big\}, \mathrm{a}+N-1\big)$. Note that the hypotheses of the discrete Stieltjes theorem are fulfilled in $\mathrm{K}$. Thus we only need to prove that $\mathrm{S}_{(2\mathrm{a}, 2\mathrm{a}+1)}\subset \mathrm{K}$\footnote{It is immediate that $\mathrm{S}_{(-1, 2\mathrm{a}]} \subset \mathrm{K}$.}.  Note that $W_n^{(\alpha)}(\mathrm{a}(\mathrm{a}+1))=1$ and
\begin{align*}
 W_n^{(\alpha)}\big((\alpha-\mathrm{a})(\alpha-\mathrm{a}+1)\big)=\pFq{2}{1}{-n, 2\mathrm{a}-\alpha}{1-N}{1}=\frac{(1-N+\alpha-2\mathrm{a})_n}{(1-N)_n}>0,
 \end{align*}
the last equality being a consequence of Chu-Vandermonde's identity (\cite[Corollary 2.2.3]{AAR99}). The rest of the proof runs as in Proposition \ref{Racah}.
\end{proof}

\begin{obs}
In \cite[Section 9.6]{KLS10}, the dual Hahn polynomials are defined by
$$
\widetilde{W}_n^{(\alpha, \beta)}(\widetilde{x}(s))=\pFq{3}{2}{-n, -s, s+\alpha+\beta+1}{\alpha+1, 1-N}{1}
$$
$(n=1,\dots, N-1; \mathrm{a}=0, \mathrm{b}=N; \alpha>-1, \beta>-1\; \text{or}\; \alpha<1-N, \beta<1-N)$, where $\widetilde{x}(s)=s(s+\alpha+\beta+1)$.  In \eqref{DualHahneq}, we can write $W^{(\alpha, \mathrm{a})}_n$ instead of $W^{(\alpha)}_n$. Hence, for $\mathrm{a}=(\alpha+\beta)/2$ fixed,
$$
W^{(\alpha, (\alpha+\beta)/2)}_n(x(s))=\widetilde{W}_n^{(\alpha, \beta)}\left(\widetilde{x}\left(s-\frac{\alpha+\beta}{2}\right)\right),
$$
where
$$
\widetilde{x}\left(s-\frac{\alpha+\beta}{2}\right)=x(s)-\frac{(\alpha+\beta+1)^2-1}{4}.
$$
Consequently, Proposition \ref{DualHahn} remains valid if we replace $W^{(\alpha)}_n$ by $\widetilde{W}_n^{(\alpha, \beta)}$ and assume that $\alpha+\beta$ is constant.
\end{obs}

\subsection{The grid $(\mathrm{III})$}\label{qlinear}
\subsubsection{Examples of COPRL on $X=x(s)=q^{-s}$ $(0<q<1)$ (q-Meixner, Al-Salam-Carlitz, q-Charlier, q-Hahn, q-Krawtchouk, affine q-Krawtchouk, and quantum q-Krawtchouk) and some related cases (q-Charlier and big q-Laguerre).}
\vspace{.2cm}
The {\em q-Meixner polynomials} (see \cite[Section 14.13]{KLS10}),
$$
y(s)=M_n^{(\alpha, \beta)}(X; q)=\pPq{2}{1}{q^{-n}, X}{\beta q}{q, -\frac{q^{n+1}}{\alpha}}
$$
$(n=1,2,\dots; \mathrm{a}=0, \mathrm{b}=\infty; \alpha>0, 0\leq\beta<q^{-1})$, satisfy \eqref{DE2} with $A$ and $B$ given by
\begin{align*}
A=A(s; \alpha, \beta)=(1-q^s)(1+\alpha \beta q^s), \quad B=B(s; \alpha, \beta)=\alpha q^s(1-\beta q^{s+1}).
\end{align*}

\begin{proposition}
The zeros of $M_n^{(\alpha, \beta)}(\cdot; q)$ are strictly increasing functions of $\alpha$ on $(0, \infty)$ and strictly decreasing functions of $\beta$ on $[0,q^{-1})$.
\end{proposition}
\begin{proof}
The monotonicity function
$$
f=\frac{B}{A}=\frac{\alpha q^s(1-\beta q^{s+1})}{(1-q^s)(1+\alpha \beta q^s)}
$$
is a positive and strictly decreasing function of $s \in (0, \infty)$, and
\begin{align*}
\frac{\partial f}{\partial \alpha}&=\frac{q^s(1-\beta q^{s+1})}{(1-q^s)(1+\alpha \beta q^s)^2}>0,\\[7pt]
\frac{\partial f}{\partial \beta}&=-\frac{\alpha q^{2s} (\alpha+q)}{(1-q^s)(1+\alpha \beta q^s)^2}<0,
\end{align*}
for each $s \in(0, \infty)$. The result follows from the discrete Stieltjes theorem.
\end{proof}
\begin{obs}
The q-Charlier polynomials (see \cite[Section 14.23]{KLS10}) are given by $C_n^{(\alpha)}(\cdot; q)=M_n^{(\alpha, 0)}(\cdot; q)$.
\end{obs}

The {\em second family of Al-Salam-Carlitz polynomials} (see \cite[Section 14.25]{KLS10}),
$$
y(s)=V_n^{(\alpha)}(X; q)=(-\alpha)^n q^{-\binom{n}{2}} \pPq{2}{0}{q^{-n}, X}{-}{q, \frac{q^n}{\alpha}}
$$
$(n=1,2,\dots; \mathrm{a}=0, \mathrm{b}=\infty; 0<\alpha<q^{-1})$, satisfy \eqref{DE2} with $A$ and $B$ given by
\begin{align*}
A=A(s; \alpha)=(1-q^{-s})(\alpha-q^{-s}),\quad B=B(s; \alpha)=\alpha q.
\end{align*}

\begin{proposition}\label{aSC}
The zeros of $V_n^{(\alpha)}(\cdot; q)$ are strictly increasing functions of $\alpha$ on $(0, q^{-1})$.
\end{proposition}
\begin{proof}
Define the interval $\mathrm{K}=\big(\max\{-\log_q \alpha, 0\}, \infty\big)$. Note that the hypotheses of the discrete Stieltjes theorem are fulfilled in $\mathrm{K}$. Thus we only need to prove that $\mathrm{S}_{(1, q^{-1})}\subset \mathrm{K}$\footnote{It is immediate that $\mathrm{S}_{(0, 1]} \subset \mathrm{K}$.}.  The rest of the proof runs as in Proposition \ref{Racah}; although given the simplicity of this case, we do not need to use q-hypergeometric identities.
\end{proof}

\begin{obs}
The first family of Al-Salam-Carlitz polynomials is given by $U_n^{(\alpha)}$ $(\cdot; q^{-1})=V_n^{(\alpha)}(\cdot; q)$.
\end{obs}

The {\em q-Hahn polynomials} (see \cite[Section 14.6]{KLS10}),
$$
y(s)=H_n^{(\alpha, \beta)}(X; q)=\pPq{3}{2}{q^{-n}, \alpha \beta q^{n+1}, X}{\alpha q, q^{1-N}}{q, q}
$$
$(n=1,\dots N-1; \mathrm{a}=0, \mathrm{b}=N; 0<\alpha<q^{-1}, 0<\beta<q^{-1})$, satisfy \eqref{DE2} with $A$ and $B$ given by
\begin{align*}
A&=A(s; \alpha, \beta)=\alpha q (1-q^s)(\beta-q^{s-N}),\\[7pt]
B&=B(s; \alpha)=(1-q^{s-N+1})(1-\alpha q^{s+1}).
\end{align*}
\begin{proposition}
The zeros of $H_n^{(\alpha, \beta)}(\cdot; q)$ are strictly decreasing functions of $\alpha$ on $(0, q^{-1})$ and strictly increasing functions of $\beta$ on $(0,q^{-1})$.
\end{proposition}

The {\em q-Krawtchouk polynomials} (see \cite[Section 14.15]{KLS10}),
$$
y(s)=K_n^{(\alpha)}(X; q)=\pPq{3}{2}{q^{-n}, -\alpha q^{n}, X}{0, q^{1-N}}{q, q}
$$
$(n=1,\dots, N-1; \mathrm{a}=0, \mathrm{b}=N; \alpha>0)$, satisfy \eqref{DE2} with $A$ and $B$ given by
\begin{align*}
A=A(s; \alpha)=\alpha (q^s-1),\quad B=B(s)=1-q^{s-N+1}.
\end{align*}
\begin{proposition}
The zeros of $K_n^{(\alpha)}(\cdot; q)$ are strictly decreasing functions of $\alpha$ on $(0, \infty)$.
\end{proposition}

The {\em affine q-Krawtchouk polynomials} (see \cite[Section 14.16]{KLS10}),
$$
y(s)=\widehat{K}_n^{(\alpha)}(X; q)=\pPq{3}{2}{q^{-n},0, X}{\alpha q, q^{1-N}}{q, q}
$$
$(n=1,\dots, N-1; \mathrm{a}=0, \mathrm{b}=N; 0<\alpha<q^{-1})$, satisfy \eqref{DE2} with $A$ and $B$ given by
\begin{align*}
A=A(s; \alpha)=\alpha q^{s-N+1} (q^s-1),\quad B=B(s; \alpha)=(1-q^{s-N+1})(1-\alpha q^{s+1}).
\end{align*}
\begin{proposition}
The zeros of $\widehat{K}_n^{(\alpha)}(\cdot; q)$ are strictly decreasing functions of $\alpha$ on $(0, q^{-1})$.
\end{proposition}
\begin{obs}\label{bqL1}
The big q-Laguerre polynomials (see \cite[Section 14.16]{KLS10}) are given by
$$
L_n^{(\alpha, \beta)}(X; q)=\pPq{3}{2}{q^{-n}, 0, X}{\alpha q, \beta q}{q, q}
$$
$(n=1,2,\dots ; 0<\alpha<q^{-1}, \beta<0)$. In particular, $\widehat{K}_n^{(\alpha)}(\cdot; q)=L_n^{(\alpha, q^{-N})}(X; q)$ (see also Remark \ref{bqL2}).
\end{obs}

The {\em quantum q-Krawtchouk polynomials} (see \cite[Section 14.14]{KLS10}),
$$
y(s)=\widetilde{K}_n^{(\alpha)}(X; q)=\frac{(q^{-N},q)_n}{\alpha^n q^{n^2}} \pPq{2}{1}{q^{-n}, X}{q^{1-N}}{q, \alpha q^{n+1}}
$$
$(n=1,\dots, N-1; \mathrm{a}=0, \mathrm{b}=N; \alpha>q^{1-N})$, satisfy \eqref{DE2} with $A$ and $B$ given by
\begin{align*}
A=A(s; \alpha)=(1-q^s)(\alpha-q^{s-N}),\quad B=B(s)=-q^s(1-q^{s-N+1}).
\end{align*}
\begin{proposition}
The zeros of $\widetilde{K}_n^{(\alpha)}(\cdot; q)$ are strictly decreasing functions of $\alpha$ on $(q^{1-N}, \infty)$.
\end{proposition}
\begin{proof}
Define the interval $\mathrm{K}=\big(\max\{0, \log_q \alpha+N\}, N-1\big)$. Note that the hypotheses of the discrete Stieltjes theorem are fulfilled in $\mathrm{K}$. Thus we only need to prove that $\mathrm{S}_{(q^{1-N}, q^{-N})}\subset \mathrm{K}$\footnote{It is immediate that $\mathrm{S}_{[q^{-N}, \infty)}\subset \mathrm{K}$.}.  The rest of the proof runs as in Proposition \ref{Racah}; although given the simplicity of this case, we do not need to use q-hypergeometric identities.
\end{proof}

\subsubsection{Examples of COPRL on $X=x(s)=q^{s}$ $(0<q<1)$ (q-Bessel, little q-Jacobi, and little q-Laguerre/Wall) and some related cases (big q-Jacobi, big q-Laguerre, and q-Laguerre).}
\vspace{.2cm}
The {\em q-Bessel} (see \cite[Section 14.22]{KLS10}),
$$
y(s)=B_n^{(\alpha)}(X; q)=\pPq{2}{1}{q^{-n}, -\alpha q^{n}}{0}{q, q X}
$$
$(n=1, 2,\dots; \mathrm{a}=0, \mathrm{b}=\infty; \alpha>0)$, satisfy \eqref{DE2} with $A$ and $B$ given by
\begin{align*}
A=A(s)=q^{s}-1,\quad B=B(s; \alpha)=\alpha.
\end{align*}
\begin{proposition}
The zeros of $B_n^{(\alpha)}(\cdot; q)$ are strictly decreasing functions of $\alpha$ on $(0, \infty)$.
\end{proposition}

The {\em little q-Jacobi polynomials} (see \cite[Section 14.12]{KLS10}),
$$
y(s)=P_n^{(\alpha, \beta)}(X; q)=\pPq{2}{1}{q^{-n}, \alpha \beta q^{n+1}}{\alpha q}{q, q X}
$$
$(n=1, 2,\dots; \mathrm{a}=0, \mathrm{b}=\infty; 0<\alpha<q^{-1}, \beta<q^{-1})$, satisfy \eqref{DE2} with $A$ and $B$ given by
\begin{align*}
A=A(s)=q^{-s} (q^{s}-1),\quad B=B(s; \alpha, \beta)=\alpha q^{-s} (\beta q^{s+1}-1).
\end{align*}
\begin{proposition}
The zeros of $P_n^{(\alpha, \beta)}(\cdot; q)$ are strictly decreasing functions of $\alpha$ on $(0, q^{-1})$ and strictly increasing functions of $\beta$ on $(-\infty,q^{-1})$.
\end{proposition}

\begin{coro}
The zeros of the special case of big q-Jacobi polynomials (see \cite[Section 14.5]{KLS10}),
$$
P_n^{(\alpha, \beta, 0)}(X; q)=\frac{(\beta q; q)_n}{(\alpha q; q)_n} (-1)^n \alpha^n q^{n+\binom{n}{2}} P_n^{(\beta, \alpha)}(\alpha^{-1} q^{-1} X; q)
$$
 $(n=1, 2,\dots; 0<\alpha<q^{-1}, 0<\beta<q^{-1})$, are strictly increasing functions of $\alpha$ on $(0, q^{-1})$ and strictly decreasing functions of $\beta$ on $(0,q^{-1})$.
\end{coro}

\begin{obs}\label{bqL2}
$P_n^{(\alpha, 0, 0)}(\cdot; q)=L_n^{(\alpha, 0)}(\cdot; q)$ (see also Remark \ref{bqL1}).
\end{obs}

The {\em little q-Laguerre/Wall polynomials} (see \cite[Section 14.20]{KLS10}),
$$
y(s)=L_n^{(\alpha)}(X; q)=\pPq{2}{1}{q^{-n}, 0}{\alpha q}{q, q X}
$$
$(n=1, 2,\dots; \mathrm{a}=0, \mathrm{b}=\infty; 0<\alpha<q^{-1})$, satisfy \eqref{DE2} with $A$ and $B$ given by
\begin{align*}
A=A(s)=q^{s}-1,\quad B=B(s; \alpha)=\alpha q^{-s}.
\end{align*}
\begin{proposition}
The zeros of $L_n^{(\alpha)}(\cdot; q)$ are strictly decreasing functions of $\alpha$ on $(0, q^{-1})$.
\end{proposition}

\begin{coro}
The zeros of the q-Laguerre polynomials (see \cite[Section 14.21]{KLS10}),
$$
\widehat{L}_n^{(\alpha)}(\cdot; q)=q^{-\alpha n}\frac{(q^{\alpha+1}; q)_n}{(q; q)_n}L_n^{(q^\alpha)}(\cdot; q)
$$
 $(n=1, 2,\dots; \alpha>-1)$, are strictly decreasing functions of $\alpha$ on $(-1, \infty)$.
\end{coro}

\subsection{The grid $(\mathrm{VI})$}\label{latticeII}
\subsubsection*{Examples of COPRL on $X=x(s)=(q^{s}+q^{-s})/2$ $(0<q<1)$ (q-Racah and dual q-Hahn polynomials).}

The {\em q-Racah polynomials} (see \cite[p. 239]{NU93}),
\begin{align}\label{qFRacah}
y(s)=R_n^{(\alpha, \beta)}(X; q)&=\pPq{4}{3}{q^{-n}, q^{\alpha+\beta+n+1}, q^{\mathrm{a}-s},q^{s+\mathrm{a}}}{q^{2\mathrm{a}+\alpha+N}, q^{\beta+1}, q^{1-N}}{q,q}
\end{align}
$(n=1,\dots, N-1; \mathrm{a}>0, \mathrm{b}=\mathrm{a}+N; \alpha>-1, -1<\beta<2\mathrm{a})$, satisfy \eqref{DE2} with $A$ and $B$ given by
\begin{align*}
A&=A(s; \alpha, \beta)\\[7pt]
&=-\frac{4\,q^{\alpha+\beta+5/2}(q^{s-\mathrm{a}}-1)(q^{s+\mathrm{a}+N-1}-1)(q^{s-\mathrm{a}-\alpha-N}-1)(q^{s+\mathrm{a}-\beta-1}-1)}{(q-1)^2(q^{2s}-1)(q^{2s-1}-1)},\\[7pt]
B&=B(s; \alpha, \beta)\\[7pt]
&=-\frac{4\, q^{3/2}(q^{s+\mathrm{a}}-1)(q^{s-\mathrm{a}-N+1}-1)(q^{s+\mathrm{a}+\alpha+N}-1)(q^{s-\mathrm{a}+\beta+1}-1)}{(q-1)^2(q^{2s}-1)(q^{2s+1}-1)}.
\end{align*}

\begin{proposition}\label{qRacah}
The zeros of $R_n^{(\alpha,\beta)}(\cdot; q)$ are strictly decreasing functions of $\alpha$ on $(-1, \infty)$ and strictly increasing function of $\beta$ on $(-1, 2\mathrm{a})$ if $\mathrm{a}\geq 1/2$, or else the zeros of $R_n^{(\alpha,\beta)}(\cdot; q)$ are strictly decreasing functions of $\alpha$ on $(-1, \infty)$ for each $\beta\in (\mathrm{a}-1/2, 2\mathrm{a})$ and strictly increasing function of $\beta$ on $(\mathrm{a}-1/2, 2\mathrm{a})$.
\end{proposition}
\begin{proof}
We give the proof only for the case in which $\mathrm{a}\geq 1/2$. The proof for $0<\mathrm{a}<1/2$ is similar. Define the interval $\mathrm{K}=\big(\max\big\{\mathrm{a}, \beta-\mathrm{a}+1\big\}, \mathrm{a}+N-1\big)$.
The monotonicity function
\begin{align*}
f&=\frac{B}{A}\\[7pt]
&=\frac{(q^{s+\mathrm{a}}-1)(q^{s-\mathrm{a}-N+1}-1)(q^{s+\mathrm{a}+\alpha+N}-1)(q^{s-\mathrm{a}+\beta+1}-1)}{q^{\alpha+\beta+1}(q^{s-\mathrm{a}}-1)(q^{s+\mathrm{a}+N-1}-1)(q^{s-\mathrm{a}-\alpha-N}-1)(q^{s+\mathrm{a}-\beta-1}-1)}\frac{q^{2s-1}-1}{q^{2s+1}-1}
\end{align*}
is a positive and strictly decreasing function of $s \in\mathrm{K}$, and
\begin{align*}
&0>\frac{\partial f}{\partial \alpha}=\\[7pt]
&\frac{\log q\,(q^{2s}-1)(q^{2s-1}-1)(q^{s+\mathrm{a}}-1)(q^{s-\mathrm{a}-N+1}-1)(q^{s-\mathrm{a}+\beta+1})}{q^{\alpha+\beta+1/}(q^{2s+1}-1)(q^{s-\mathrm{a}}-1)(q^{s+\mathrm{a}+N-1}-1)(q^{s-\mathrm{a}-\alpha-N}-1)^2(q^{s+\mathrm{a}-\beta-1}-1)},\\[7pt]
&0<\frac{\partial f}{\partial \beta}=\\[7pt]
&\frac{\log q\,(q^{2s}-1)(q^{2s-1}-1)(q^{s+\mathrm{a}}-1)(q^{s-\mathrm{a}-N+1}-1)(q^{s+\mathrm{a}+\alpha+N})}{q^{\alpha+\beta+1}(q^{2s+1}-1)(q^{s-\mathrm{a}}-1)(q^{s+\mathrm{a}+N-1}-1)(q^{s-\mathrm{a}-\alpha-N}-1)(q^{s+\mathrm{a}-\beta-1}-1)^2},
\end{align*}
for each $s \in\mathrm{K}$. Thus we only need to prove that $\mathrm{S}^{(\beta)}_{(2\mathrm{a}-1, 2\mathrm{a})}\subset \mathrm{K}$\footnote{It is immediate that $\mathrm{S}^{(\beta)}_{(-1, 2\mathrm{a}-1]}\subset \mathrm{K}$.}. Indeed, this is equivalent to prove that $R_n^{(\alpha,\beta)}(\cdot; q)$ has no zeros on $\big((q^{\mathrm{a}}+q^{-\mathrm{a}})/2, (q^{\beta-\mathrm{a}+1}+q^{\mathrm{a}-\beta-1})/2\big)$ for $\beta \in (2\mathrm{a}-1, 2\mathrm{a})$. For $\beta \in (2\mathrm{a}-1, 2\mathrm{a})$, $R_n^{(\alpha, \beta)}\big((q^{\mathrm{a}}+q^{-\mathrm{a}})/2; q\big)=1$ and
\begin{align*}
 R_n^{(\alpha, \beta)}\big((q^{\beta-\mathrm{a}+1}+q^{\mathrm{a}-\beta-1})/2; q\big)&=\pPq{3}{2}{q^{-n}, q^{\alpha+\beta+n+1}, q^{2\mathrm{a}-\beta-1}}{q^{2\mathrm{a}+\alpha+N}, q^{1-N}}{q,q}\\[7pt]
 &=\frac{(q^{2\mathrm{a}-\beta+N-n-1}; q)_n (q^{\alpha+\beta+N+1}; q)_n}{(q^{2\mathrm{a}+\alpha+N} ;q)_n (q^{N-n}; q)_n}>0,
\end{align*}
the last equality being a consequence of q-Pfaff-Saalsch\"utz's identity (see \cite[(10.10.3)]{AAR99}).  We thus get, as in Proposition \ref{Racah}, $\mathrm{S}^{(\beta)}_{(-1, 2\mathrm{a})}\subset \mathrm{K}$. The same proof actually shows that $\mathrm{S}_{(-1,\infty)}^{(\alpha)}\subset \mathrm{K}$. The result follows from the discrete Stieltjes theorem.
\end{proof}

\begin{obs}
In \cite[Section 14.2]{KLS10}, the q-Racah polynomials (replacing $\alpha$ by $q^{\beta}$ and $\beta$ by $q^{\alpha}$) are defined by
$$
\widetilde{R}_n^{(\alpha, \beta, \delta)}(\widetilde{x}(s); q)=\pPq{4}{3}{q^{-n}, q^{\alpha+\beta+n+1}, q^{-s},q^{s-N+1}}{\delta q^{\alpha+1}, q^{\beta+1}, q^{1-N}}{q,q}
$$
$(n=1,\dots, N-1; \mathrm{a}=0, \mathrm{b}=N; \alpha>-1, -1<\beta<1-N+\log_q \delta, 0<\delta<q^{N-1})$, where $\widetilde{x}(s)=\delta q^{s-N+1}+q^{-s}$.  In \eqref{qFRacah}, we can write $R^{(\alpha, \beta, \mathrm{a})}_n(\cdot; q)$ instead of $R^{(\alpha, \beta)}_n(\cdot; q)$. Hence, for  $\mathrm{a}=(1-N+\log_q \delta)/2$ fixed,
$$
R^{(\alpha, \beta, (1-N+\log_q \delta)/2)}_n(x(s); q)=\widetilde{R}_n^{(\alpha, \beta, \delta)}\left(\widetilde{x}\left(s-(1-N+\log_q \delta)/2\right); q\right),
$$
where
$$
\widetilde{x}\left(s-(1-N+\log_q \delta)/2\right)=2 q^{(1-N+\log_q \delta)/2} x(s).
$$
Consequently, Proposition \ref{qRacah} remains valid if we replace $R^{(\alpha, \beta)}_n(\cdot; q)$ by $\widetilde{R}_n^{(\alpha, \beta, \delta)}$ $(\cdot; q)$.\end{obs}

The {\em dual q-Hahn polynomials} (see \cite[p. 239]{NU93}),
\begin{align}\label{DqH}
y(s)=W_n^{(\alpha)}(X; q)&=\pPq{3}{2}{q^{-n}, q^{\mathrm{a}-s},q^{s+\mathrm{a}}}{q^{\alpha+1}, q^{1-N}}{q,q}
\end{align}
$(n=1,\dots, N-1; \mathrm{a}>0, \mathrm{b}=\mathrm{a}+N; -1<\alpha<2\mathrm{a})$, satisfy \eqref{DE2} with $A$ and $B$ given by
\begin{align*}
A&=A(s; \alpha)=-\frac{4 q^{-\mathrm{a}+\alpha-N+5/2}(q^{s-\mathrm{a}}-1)(q^{s+\mathrm{a}+N-1}-1)(q^{s+\mathrm{a}-\alpha-1}-1)}{(q-1)^2(q^{2s}-1)(q^{2s-1}-1)},\\[7pt]
B&=B(s; \alpha)=-\frac{4 q^{s+3/2}(q^{s+\mathrm{a}}-1)(q^{s-\mathrm{a}-N+1}-1)(q^{s-\mathrm{a}+\alpha+1}-1)}{(q-1)^2(q^{2s}-1)(q^{2s+1}-1)}.
\end{align*}

\begin{proposition}\label{DqHahn}
The zeros of $W_n^{(\alpha)}(\cdot; q)$ are strictly increasing functions of $\alpha$ on $(-1, 2\mathrm{a})$ if $\mathrm{a}\geq 1/2$, or else the zeros of $W_n^{(\alpha)}(\cdot; q)$ are strictly increasing functions of $\alpha$ on $(\mathrm{a}-1/2, 2\mathrm{a})$.
\end{proposition}
\begin{proof}
We sketch the proof only for the case in which $\mathrm{a}\geq 1/2$. The proof for $0<\mathrm{a}<1/2$ is similar. Define the interval $\mathrm{K}=\big(\max\big\{\mathrm{a}, \alpha-\mathrm{a}+1\big\}, \mathrm{a}+N-1\big)$. Note that the hypotheses of the discrete Stieltjes theorem are fulfilled in $\mathrm{K}$. Thus we only need to prove that $\mathrm{S}_{(2\mathrm{a}-1, 2\mathrm{a})}\subset \mathrm{K}$\footnote{It is immediate that $\mathrm{S}_{(-1, 2\mathrm{a}-1]} \subset \mathrm{K}$.}.  Note that $W_n^{(\alpha)}\big((q^{\mathrm{a}}+q^{-\mathrm{a}})/2; q\big)=1$ and
\begin{align*}
 W_n^{(\alpha)}\big((q^{\alpha-\mathrm{a}+1}+q^{\mathrm{a}-\alpha-1})/2; q\big)&=\pPq{2}{1}{q^{-n}, q^{2\mathrm{a}-\alpha-1}}{q^{1-N}}{q, q}\\[7pt]
 &=q^{n(2\mathrm{a}-\alpha-1)}\frac{(q^{-2\mathrm{a}+\alpha-N-2}; q)_n}{(q^{1-N}; q)_n}>0,
 \end{align*}
the last equality being a consequence of q-Chu-Vandermonde's identity (\cite[(1.11.5)]{KLS10}). The rest of the proof runs as in Proposition \ref{qRacah}.
\end{proof}

\begin{obs}
In \cite[Section 14.7]{KLS10}, the dual q-Hahn polynomials are defined by
$$
\widetilde{W}_n^{(\alpha, \beta)}(\widetilde{x}(s); q)=\pPq{3}{2}{q^{-n}, q^{\mathrm{a}-s},q^{s+\mathrm{a}}}{q^{\alpha+1}, q^{1-N}}{q,q}
$$
$(n=1,\dots, N-1; \mathrm{a}=0, \mathrm{b}=N; \alpha>-1, \beta>-1\; \text{or}\; \alpha<-N, \beta<-N)$, where $\widetilde{x}(s)=\delta q^{s+\alpha+\beta+1}+q^{-s}$.  In \eqref{DqH}, we can write $W^{(\alpha, \mathrm{a})}_n(\cdot; q)$ instead of $W^{(\alpha)}_n(\cdot; q)$. Hence, for  $\mathrm{a}=(\alpha+\beta+1)/2$ fixed,
$$
W^{(\alpha, (\alpha+\beta+1)/2)}_n(x(s); q)=\widetilde{W}_n^{(\alpha, \beta)}\left(\widetilde{x}\left(s-(\alpha+\beta+1)/2\right); q\right),
$$
where
$$
\widetilde{x}\left(s-(\alpha+\beta+1)/2\right)=2 q^{(\alpha+\beta+1)/2} x(s).
$$
Consequently, Proposition \ref{DqHahn} remains valid if we replace $W^{(\alpha)}_n(\cdot; q)$ by $\widetilde{W}_n^{(\alpha, \beta)}$ $(\cdot; q)$ and assume that $\alpha+\beta$ is constant.
\end{obs}

\section*{Acknowledgements}
The authors thank the Warsaw University Library for kindly sending them the book \cite{NU78} and the Keldysh Institute of Applied Mathematics for making \cite{NU83} available at \url{https://keldysh.ru/papers/1983/prep1983_17.pdf} after they request. KC is supported by the Centre for Mathematics of the University of Coimbra - UIDB/00324/2020, funded by the Portuguese Government through FCT/ MCTES. FRR and AS are supported by the Funda\c{c}\~ao de Amparo \`a Pesquisa do Estado de Minas Gerais (FAPEMIG) Demanda Universal under the grant APQ-03782-18, Conselho Nacional de Desenvolvimento Cient\'\i fico e Tecnol\'ogico (CNPq), and Coordena\c{c}\~ao de Aperfei\c{c}oamento de Pessoal de N\'\i vel Superior (CAPES).
\bibliographystyle{plain}


\section{Introduction}

Far-field automatic speech recognition (ASR) and speaker diarization are important areas of research and have many real-world applications such as transcribing meetings~\cite{Hain2012TranscribingMW,Hori2012LowLatencyRM,Renals2017DistantSR,Yoshioka2018RecognizingOS} and in-home conversations~\cite{barker2018fifth}. Although deep learning methods (including end-to-end approaches) have achieved promising results for several tasks such as Switchboard~\cite{Xiong2016AchievingHP,wang2019espresso} and LibriSpeech~\cite{Lscher2019RWTHAS, Synnaeve2019EndtoendAF}, their performance remains unsatisfactory for far-field conditions in real environments, such as the CHiME-5 dataset~\cite{barker2018fifth}. This can be attributed to: (i) noise and reverberation in the acoustic conditions, (ii) conversational speech and speaker overlaps, and (iii) challenge-specific restrictions such as insufficient training data.

Several advances have been made in the last decade to tackle the challenges offered by real, far-field speech. For ASR, this improvement can be attributed to improved neural network architectures~\cite{zorilatoshiba,kanda2018hitachi}, effective data augmentation techniques~\cite{zorila2019investigation}, and advances in speech enhancement~\cite{boeddeker2018front}. Previous work has tried tackling reverberation and noise present in the far-field recording by multi-style training with data augmentation via room impulse responses and background noises~\cite{ko2017study,wang2019jhu}. Recently, spectral augmentation has been successfully used for both end-to-end~\cite{park2019specaugment} and hybrid ASR systems~\cite{Zhou2020TheRA}. Adapting the acoustic model to the environment~\cite{seltzer2013investigation} and speaker~\cite{saon2013speaker} has also been studied. Another popular direction is front-end based approaches such as dereverberation~\cite{drude2018nara} and denoising through beamforming~\cite{anguera2007acoustic,Nakatani2019AUC}, which utilize multi-microphone data.
Far-field speaker diarization~\cite{Garca2019SpeakerDI} has also benefited from enhancement methods~\cite{Sun2020ProgressiveMN,Kataria2020FeatureEW} and approaches to handle overlapping speech~\cite{bullock2019overlap}. Recently, guided source separation (GSS)~\cite{boeddeker2018front} was proposed, which makes use of additional information such as time and speaker annotations for mask estimation. However, it requires a strong diarization system to perform good separation.

In this paper, we describe a multi-microphone multi-speaker ASR system developed using many of these methods for the CHiME-6 challenge \cite{watanabe2020chime}. The challenge aims to improve speech recognition and speaker diarization for far-field conversational speech in challenging environments in a multi-microphone setting. The CHiME-6 data \cite{barker2018fifth} contains a total of 20 4-speaker dinner party recordings. Each dinner party is two to three hours long and is recorded simultaneously on the participants' ear-worn microphone and six microphone arrays placed in the kitchen, dining room, and the living room. The challenge consists of two tracks. Track 1 allows the use of oracle start and end times of each utterance, and speaker labels for each segment. This track focuses on core ASR techniques, and measures system performance in terms of transcription accuracy. Track 2 is a ``diarization+ASR'' track. It additionally requires end-pointing speech segments in the recording, and assigning them speaker labels, i.e diarization.  To this end, VoxCeleb2 data~\cite{nagrani2017voxceleb} is permitted for training a diarization system, and concatenated minimum-permutation word error rate (cpWER) is used to measure speaker-attributed transcription accuracy.

Our system for Track 2, as shown in Figure \ref{fig:pipeline}, consists of three main modules: \emph{enhancement}, \emph{diarization} and \emph{recognition}, described in Sections \ref{sec:enhancement}, \ref{sec:diarization} and \ref{sec:recognition} below. The \emph{enhancement module} performs (i) dereverberation using online multichannel weighted prediction error (WPE), followed by (ii) denoising with a weighted delay-and-sum beamformer, and (iii) multi-array guided source separation (GSS), described below. In the \emph{diarization module}, the beamformer outputs, one per array, are used for (i) speech activity detection and (ii) speaker diarization, both of which fuse information across arrays to improve accuracy, and (iii) overlap-aware variational Bayes hidden Markov model (VB-HMM) resegmentation to assign multiple speakers to overlapped speech regions.  Speaker-marks from this diarization module are used in the multi-array GSS (part of the enhancement module) to produce enhanced, speaker-separated waveforms.  The \emph{recognition module} processes the GSS output using (i) an acoustic and n-gram language model for ASR decoding, (ii) an RNN language model for lattice rescoring, and (iii) sMBR lattice combination.  We augment the clean acoustic training data with dereverberated, beamformed and GSS-enhanced far-field data to match the test conditions.

The diarization module is replaced with oracle speech segments and speaker labels in our system for Track 1.

\begin{figure}[t]
\centering
\includegraphics[width=0.7\linewidth]{figs/pipeline.pdf}
\caption{Overview of the decoding pipeline for track 2. For track 1, we use a similar system, with the exception that the diarization module (shown in the dotted box in the figure) is replaced with oracle speech segments and speaker labels.}
\label{fig:pipeline}
\end{figure}

\section{Speech Enhancement}
\label{sec:enhancement}

\subsection{Dereverberation and Denoising}
\label{sec:wpe+beamformit}

We used an online version of the publicly available NARA-WPE~\cite{drude2018nara} implementation of weighted prediction error (WPE) based dereverberation for multi-channel signals~\cite{nakatani2010speech} for all the channels in each array. This was followed by array-level weighted delay-and-sum beamforming using the BeamformIt tool~\cite{anguera2007acoustic}. All further processing was done on the dereverberated and beamformed signals.

\subsection{Guided Source Separation (GSS)}
\label{sec:gss}

Multi-array GSS~\cite{boeddeker2018front,kanda2019guided} was applied to enhance target speaker speech signals. For track 1, we used oracle speech segmentations and speaker labels, while for track 2, we used the segmentation estimated by the speaker diarization module described in Section \ref{sec:diarization}. In GSS, the
source activity pattern for each speaker derived from the segmentation aids in resolving the permutation ambiguity. A context window is used to obtain an extended segment which can bring sufficient sparsity in the activity pattern of the target speaker to further reduce speaker permutation issues. We found a 20-second context to be ideal from our experiments. GSS gave a relative WER improvement of 9.3\% and 6.5\% on Dev and Eval respectively for track 2 over the baseline delay-and-sum beamformer, as seen in Section \ref{sec:chime-results}, Table \ref{tab:asr_results2}.

\section{Speaker Diarization}
\label{sec:diarization}

\begin{figure*}[t]
\centering
\includegraphics[width=0.8\linewidth]{figs/diarization_overlap.pdf}
\caption{The two-pass speaker diarization module.  Synchronous SAD marks across microphone arrays enables PLDA score fusion before AHC in first-pass diarization (dotted rectangle). Overlap detection (bottom) enables the second-pass resegmentation, initialized by the first-pass output (upper RTTM), to assign more than one speaker to overlapped-speech regions in the final output (right RTTM).}
\label{fig:diarization}
\end{figure*}

\subsection{Speech Activity Detection}
\label{sec:SAD}

For speech activity detection (SAD), we first trained a neural network classifier to assign each frame in an utterance a label from $\cal C$ = \{\textit{silence}, \textit{speech}, \textit{garbage}\}. We used the architecture shown in Table~\ref{tab:sad_dnn}, consisting of time-delay neural network (TDNN) layers to capture long temporal contexts \cite{Peddinti2015ATD}, interleaved with stats-pooling layers to aggregate utterance-level statistics. To obtain classifier training targets, we used a speaker-independent GMM-HMM ASR system to align whole recordings with the training transcriptions.
To leverage multiple channels, we carried out posterior fusion of the classifier outputs across all the arrays at test time. i.e. if $p^k_t(i)$ denotes the classifier's probability for class $i\in C$ at frame $t$ based on array $k$, then $p_t(i) = f(p^k_t(i))$, where $f$ is the fusion criterion.

We then post-processed the per-frame classifier outputs to enforce minimum and maximum speech/silence durations, by constructing a simple HMM whose state transition diagram encodes these constraints, treating the per-frame SAD posteriors $p_t(\cdot)$ like emission probabilities, and performing Viterbi decoding to obtain the most likely SAD label-sequence.

\begin{table}[tb]
\centering
\caption{Neural network architecture for speech activity detection (SAD).  $T$ is the input length, $\cal C$ the set of output classes.}
\label{tab:sad_dnn}
\begin{adjustbox}{max width=\linewidth}
\begin{tabular}{cccc}
\toprule
\thead{\textbf{Layer}} & \thead{\textbf{Layer context}} & \thead{\textbf{Total}\\\textbf{context}} & \thead{\textbf{Input$\times$Output}} \\
\midrule
tdnn1 & [$t-2$, $t+2$] & 5 & $150\times 256$ \\
tdnn2 & [$t-1$, $t+2$] & 8 & $1024\times 256$ \\
tdnn3 & \{$t-3$, $t$, $t+3$, $t+6$\} & 17 & $1024\times 256$ \\
stats1 & [0, $T$) & $T$ & \\
tdnn4 & \{$t-6$, $t$, $t+6$, $t+12$\} & $T$ & $1024 \times 256$\\
stats2 & [0, $T$) & $T$ & \\
tdnn5 & \{$t-12$, $t$, $t+12$, $t+24$\} & $T$ & 256 $\times |\cal C|$ \\
\bottomrule
\end{tabular}
\end{adjustbox}
\end{table}

\begin{table}[t]
\centering
\caption{Reducing SAD errors by fusing posterior probabilities from multiple arrays. SAD errors are comprised of missed speech (MS) and false alarms (FA).}
\label{tab:sad_results}
\begin{adjustbox}{max width=\linewidth}
\begin{tabular}{lcccccc}
\toprule
\multicolumn{1}{c}{\multirow{2}{*}{\textbf{System}}} & \multicolumn{3}{c}{\textbf{Dev}}           & \multicolumn{3}{c}{\textbf{Eval}}          \\
                                 & \textbf{MS} & \textbf{FA} & \textbf{Total} & \textbf{MS} & \textbf{FA} & \textbf{Total} \\
\midrule
Baseline (U06) & 2.7 & \textbf{0.6} & 3.3 & 4.4 & \textbf{1.5} & 5.9 \\
Posterior Mean & 1.7 & 0.7 & 2.4 & 3.2 & 1.9 & \textbf{5.1} \\
Posterior Max & \textbf{1.1} & 0.8 & \textbf{1.9} & \textbf{2.4} & 2.8 & 5.2 \\
\bottomrule
\end{tabular}
\end{adjustbox}
\end{table}

For the fusion criterion, we experimented with simple element-wise mean and max, and found them to be effective. From Table~\ref{tab:sad_results}, we can see that applying posterior fusion across arrays improved the SAD error by approximately 34\% relative. Since most of the gain comes from reduction in missed speech, this also positively impacts downstream recognition rate.

We also tried microphone-level posterior fusion jointly across all the arrays, but it did not yield any improvements over using the beamformed signals as described above.

\subsection{First Pass Speaker Diarization}
\label{sec:xvec-diarization}

Our first-pass diarization followed the method described in \cite{sell2018diarization}.  The test recordings were cut into into overlapping 1.5 second \emph{segments} with a 0.25 second\footnote{The CHiME baseline system used a 0.75 second stride.} stride \cite{landini2019but}, an x-vector was extracted from each segment, and agglomerative hierarchical clustering (AHC) was performed on the x-vectors using the probabilistic linear discriminant analysis (PLDA) score of each pair of x-vectors as their pairwise similarity.

The x-vector extractor we used is similar to that of \cite{snyder2018x}: it is comprised of several TDNN layers and stats-pooling, and was trained, per Challenge stipulations, on only the VoxCeleb data \cite{nagrani2017voxceleb}.  Data augmentation was performed by convolving the audio with simulated room impulse responses \cite{ko2017study} and by adding background noises from the CHiME-6 training data.  PLDA parameters were trained on (x-vectors of) 100k speech segments from the ca 40 speakers in the CHiME-6 training data.

Similar to the SAD posterior probability fusion described in Section \ref{sec:SAD}, we investigated improving diarization by leveraging multiple microphone( array)s at test time.  To compute the pairwise similarity of two 1.5 second segments during clustering, we fused the PLDA scores for their x-vectors extracted from different arrays.  We found that multi-array PLDA score fusion, specifically the element-wise maximum across arrays, provided noticeable gains.

\subsection{Overlap-Aware Resegmentation}
\label{sec:vb-reseg}

Since AHC is not designed to handle overlapping speakers, we resegmented the audio using an \emph{overlap-aware} version of the VB-HMM of \cite{Dez2018SpeakerDB}.  Speaker labels from the first-stage diarization of Section \ref{sec:xvec-diarization} were used to initialize the per-frame speaker posterior matrix, also known as the $\boldsymbol{Q}$-matrix, and one iteration of VB-HMM inference was performed to convert this (binary) $\boldsymbol{Q}$-matrix into per-frame speaker-probabilities.  Separately, we trained an \emph{overlap detector}|a 2-layer bidirectional LSTM with SincNet input features \cite{Ravanelli2018SpeakerRF} and binary (non-overlapped/overlapped speech) output labels|using the CHiME-6 training data.  The per-frame decisions of the overlap detector on the test data were then used to assign each frame to either one or two most likely speakers according to $\boldsymbol{Q}$, as described in \cite{bullock2019overlap}.  An unintended attribute of VB resegmentation was a significant number of very short segments.  The computational complexity of the GSS module (cf Section \ref{sec:gss}) is severely impacted by this growth in the number of segments.  We therefore removed all segments shorter than 200ms.

The overall diarization process is shown in Figure \ref{fig:diarization}.

\subsection{Diarization Performance}

The CHiME Challenge provided two ``ground truths'' for diarization, i.e. two NIST-style rich transcription time marks (RTTM): one based on utterance-level time marks by human annotators (named Annotation RTTM), and another on forced-alignment of the acoustics to the transcripts (resp Alignment RTTM).  The diarization output was scored against each RTTM using the DiHARD \texttt{dscore} toolkit\footnote{https://github.com/nryant/dscore}, and diarization error rate (DER) as well as Jaccard error rate (JER) were computed.

\begin{table}
\centering
\caption{Diarization performance on track 2, showing the impact of the modifications of Sections \ref{sec:xvec-diarization} and \ref{sec:vb-reseg} to the baseline x-vector/AHC system. Some of the improvement derives from the improved SAD of Section \ref{sec:SAD}.}
\label{tab:diar_results}
\begin{adjustbox}{max width=\linewidth}
\begin{tabular}{lcccc}
\toprule
\multicolumn{1}{c}{\multirow{2}{*}{\textbf{System}}} & \multicolumn{2}{c}{\textbf{Dev}} & \multicolumn{2}{c}{\textbf{Eval}} \\
\multicolumn{1}{c}{}                                 & \textbf{DER}    & \textbf{JER}   & \textbf{DER}    & \textbf{JER}    \\
\midrule
\multicolumn{5}{c}{\textit{Alignment RTTM}} \\
\midrule
Baseline (U06)  & 63.42  & 70.83 & \textbf{68.20} & 72.54  \\
PLDA Fusion & 63.97 & 71.65 & 71.56 & 71.32 \\
+ 0.25s stride   & 61.00 & 66.23 & 69.64 & 69.81 \\
+ overlap assign. & \textbf{58.18} & \textbf{59.92} & 69.92 & \textbf{65.64}  \\
\midrule
\multicolumn{5}{c}{\textit{Annotation RTTM}} \\
\midrule
Baseline (U06)  & 61.62  & 69.84 & 62.01 & 71.43  \\
PLDA Fusion & 60.09 & 70.31 & 62.97 & 70.09 \\
+ 0.25s stride   & 57.85 & 65.36 & 61.60 & 69.35 \\
+ overlap assign. & \textbf{50.43} & \textbf{57.81} & \textbf{58.26} & \textbf{64.38}  \\
\bottomrule
\end{tabular}
\end{adjustbox}
\end{table}

Table~\ref{tab:diar_results} shows improvements in diarization performance, relative to the AHC baseline, due to PLDA score fusion, using a 0.25s stride (v/s 0.75s), and overlap-aware resegmentation.  The gains from PLDA fusion across arrays appears modest and somewhat inconsistent relative to the single array (U06) DER, but was consistently better than the \emph{average} single-array DER across the six arrays.  The shorter (0.25s) x-vector stride yielded a robust improvement, and the most significant improvement came from overlap detection and multiple-speaker assignment (denoted overlap assign.)

Finally, note that diarization performance seems to degrade (particularly for the Eval data) when scored against the Alignment RTTM, but shows significant improvement with the Annotation RTTM.  The former stipulates tighter speech boundaries by, for instance, marking short pauses between words as non-speech.  This increases the (measured) false alarm errors of our diarization module.  However, retaining such pauses is beneficial for the downstream ASR task, by inducing more appropriate utterance segmentation.

\section{Automatic Speech Recognition}
\label{sec:recognition}
We built a hybrid DNN-HMM system using algorithms and tools available in the \href{https://github.com/kaldi-asr/kaldi}{Kaldi} ASR toolkit.

\subsection{Acoustic Modeling}

Our baseline acoustic model (AM) was a 15-layer factorized time-delay neural network or TDNN-F \cite{PoveyEtAlIS2018}. Training data for this model comprised of 80h clean CHiME-6 audio from the 2 worn-microphones of the speaker of each transcribed utterance, 320h from a 4x distortion of this clean audio using synthetic room impulse responses and CHiME-6 background noises, and 200h raw far-field audio derived by randomly sampling utterances from the many arrays.  These 600h were subject to 0.9x and 1.1x speed-perturbation to yield 1800h of AM training data.  The baseline TDNN-F was trained with the lattice-free maximum mutual information (LF-MMI) objective.  A context-dependent triphone HMM-GMM system was first trained using standard procedures in Kaldi, and frame-level forced-alignments were generated between the speech and reference transcripts to \emph{guide} LF-MMI training \cite{PoveyEtAlIS2016}.

Our CHiME-6 Challenge baseline system used this acoustic model and achieved a word error rate (WER) of 51.8\% and 51.3\% respectively on the CHiME-6 Dev and Eval test sets.  We will designate this AM by (a) in this paper.

We then experimented with several other model architectures and data selection and augmentation methods.

\begin{figure}[t]
\centering
\includegraphics[width=0.7\linewidth]{figs/am.pdf}
\caption{Illustration of the acoustic model architecture, inputs and outputs.  N-targets is the number of leaves (context dependent HMM states) in the bi-phone clustering trees.}
\label{fig:acoustic}
\end{figure}

\subsubsection{Neural Network Architecture}
\label{sec:arch}

We first trained an AM with 6 CNN layers comprising 3$\times$3 convolutional kernels, and 16 TDNN-F layers using the 1800h training set, designated model (b) in this paper, and lowered the WERs to 49.6\% and 49.3\% on the Dev and Eval sets, respectively.  But we discovered that these models take an inordinate amount of time to train, limiting our ability for exploration.

To expedite the turnaround time of AM training experiments, and informed by experiments in data selection and augmentation described in Section \ref{data_augmentation}, we created a 500h training set comprised of the 80h of clean ear-worn microphone audio and 320h of 4x distortions of the clean audio using synthetic room impulse responses and CHiME-6 background noises, as above, supplemented with 60h raw far-field audio derived by randomly sampling utterances from the many arrays and 40h from multi-array GSS enhancement, as described in Section \ref{sec:gss}.

With this 500h training set, we experimented with the following AM architectures.

\begin{enumerate}[wide, labelwidth=!, labelindent=0pt, label=(\alph*), start=3]

\item
\textbf{6 CNN + 12 TDNN-F}, comprised of 6 CNN layers with 3x3 convolutional kernels, followed by 12 TDNN-F layers. Frequency subsampling by a factor of two was applied at the third and fifth convolutional layers.  The TDNN-F layers were 1536-dimensional with a factorized (lower) rank of 160.

\item
\textbf{6 CNN + 10 TDNN-F + 2 Self-attention}, in which the last two TDNN-F layers were replaced by self-attention layers. Each self-attention layer had 15 attention heads with key and value dimensions of 40 and 80, respectively.

\item
\textbf{CNN + TDNN + LSTM}, comprised of 6 CNN layers, followed by 9 TDNN and 3 LSTM layers, interleaved (3+1)$\times$3.  The CNN layers the are same as above, TDNN layers have 1024-dimensional hidden units, and the LSTM has a 1024-dim layer and 256-dim output and recurrent projections.

\end{enumerate}

The input to these AMs was 40-dim log-Mel filterbank coefficients, and the outputs were context-dependent HMM states (Senones) derived from a left-biphone tree.  We determined empirically that a tree with 4500 leaves, an L$_2$-regularization value of 0.03 for the CNN and TDNN-F layers, and 0.04 for the dense layers performed well.  Using  6-epochs instead of 4-epochs (in the baseline) further improved the AMs.

As shown in Table~\ref{tab:architecture}, the CNN + TDNN-F architecture performed better than the others, so we selected it for subsequent experiments.

\begin{table}
\centering
\caption{Track 1 ASR WERs for AM architecture described in Section \ref{sec:arch}.  The CNN + TDNN-F configuration works best.}
\label{tab:architecture}
\begin{adjustbox}{max width=\linewidth}
\begin{tabular}{lcc}
\toprule
\textbf{(Model) Architecture} & \textbf{Dev} & \textbf{Eval} \\
\midrule
(a) 15 TDNN-F & 51.8\% & 51.3\% \\
(b) 6 CNN + 16 TDNN-F & 49.6\% & 49.3\% \\
(c) 6 CNN + 12 TDNN-F & \textbf{48.3}\% & \textbf{48.5}\%  \\
(d) 6 CNN + 10 TDNN-F + 2 SA & 49.9\% & 49.4\%  \\
(e) 6 CNN + 3 $\times$ (3 TDNN + 1 LSTM) & 50.1\%& 49.8\%  \\
\bottomrule
\end{tabular}
\end{adjustbox}
\end{table}

\subsubsection{Training Data Selection and Augmentation}
\label{data_augmentation}

\begin{table}[t]
\centering
\caption{Track 1 WERs when training the 6 CNN + 16 TDNN-F AM on only enhanced (and not reverberated, noisy) speech.}
\label{tab:data}
\begin{adjustbox}{max width=\linewidth}
\begin{tabular}{lcc}
\toprule
\textbf{(Model) Training Data}& \textbf{Dev} & \textbf{Eval} \\
\midrule
(b) 1800h incl. reverb'ed \& raw far-field& 49.6\% & 49.3\%  \\
(f)  675h incl. only clean \& enhanced& \textbf{44.5\%} & \textbf{44.9\%}  \\
(g)  675h of model (f) + sp/sil probs& 45.0\% & 45.3\%  \\
(h)  3 $\times$ 500h from models (c)-(e)& 44.6\% & 45.4\%  \\
\bottomrule
\end{tabular}
\end{adjustbox}
\end{table}

Inspired by the findings of \cite{zorila2019investigation}, where they reported good performance on CHiME data not by using the matched far-field speech, or synthetically reverberated speech, for AM training, but instead using speech enhancement during both AM training and test time, we created a new AM training data set as follows.  We applied beam-forming to the 4 microphone arrays to obtain 4 $\times$ 40h of speech, obtained another 40h from multi-array GSS described in Section \ref{sec:gss}, and combined them with the 80h of ear-worn microphone data described above. Data clean-up was applied to this 280h data set, to remove utterances that failed forced-alignment under a narrow beam, followed by speed perturbation, resulting in a 675h AM training data set.

Note from Table \ref{tab:data} that careful data selection indeed confirmed the findings of \cite{zorila2019investigation}, reducing the WER to 44.5\% and 44.9\% respectively on the Dev and Eval sets.  The model that attained this is designated (f) in this paper.

We also trained two additional AMs, one with the same 675h of data as the model (f), but with improved estimates of the inter-word silence and pronunciation probabilities (see \cite{ChenEtAlIS2015}), and another with speed perturbation of the 500h data set of models (c)-(e).  These models are designated (g) and (h) respectively in Table \ref{tab:data}.  While they do not perform better than the model (f), they were used for system combination in track 1, as described later in Table \ref{tab:asr_results} below.

\subsubsection{Overlap-Aware Training}
\label{overlap_info}

Since the CHiME data have a significant proportion of overlapped speech, we looked into providing the AM a 1-bit input indicating the presence/absence of overlap.  This \emph{overlap bit} was determined during training from the time alignments, and used alongside the 40-dim filter-bank features and 100-dim i-vectors.  We first projected the overlap bit to 40-dim, and the i-vector to 200-dim, and applied batch normalization and L$_2$ regularization.  We then combined the two with the filter-bank features to create a single-channel input to the first CNN layer of model (f) described above.

The resulting model is designated (i) in Table~\ref{tab:overlap}, which illustrates that knowledge of the presence of overlapped speech yields a modest WER improvement in track 1 conditions.
While we could have used model (i) in track 2 by using the output of the overlap detector of Section \ref{sec:vb-reseg} as the overlap bit, we did not have sufficient time to carefully conduct these experiments.

\begin{table}[t]
\centering
\caption{Track 1 WERs when the presence of overlapped speech is known to the acoustic model.}
\label{tab:overlap}
\begin{adjustbox}{max width=\linewidth}
\begin{tabular}{lcc}
\toprule
\textbf{(Model) Input features}& \textbf{Dev} & \textbf{Eval} \\
\midrule
(f) log-Mel filter-bank and i-vector & 44.5\% & 44.9\%  \\
(i) + overlapped speech indicator bit& \textbf{44.4}\% & \textbf{44.5}\%  \\
\bottomrule
\end{tabular}
\end{adjustbox}
\end{table}

\subsection{Language Modeling and Rescoring}

We used the training transcriptions to build our language models (LMs). We used a 3-gram LM trained with the SRILM toolkit~\cite{stolcke2002srilm} in the first pass decoding. For neural LMs, we used Kaldi to train recurrent neural network LMs (RNNLMs)~\cite{xu2018neural}. We performed a pruned lattice rescoring~\cite{xu2018pruned} with a forward and a backward (reversing the text at the sentence level) LSTM. We first rescored with the forward LSTM and then performed another rescoring on top of the rescored lattices using the backward LSTM.
Both LSTMs are 2-layer projected LSTMs with hidden and projection layer dimensions 512 and 128, respectively. We also used L$_2$-regularization on the embedding, hidden, and output layers.

\subsection{Lattice Combination}

For track 1, we used the lattice combination method to combine four CNN + TDNN-F acoustic models, named models (f), (g), (h) and (i) in Section \ref{data_augmentation}, all with WER's in the 44\%-45\% range, as seen in Tables \ref{tab:data} and \ref{tab:overlap}.  While models (f), (g) and (i) were trained with only worn-microphone and enhanced far-field data, model (h) included raw far-field data and synthetically reverberated/noisy data.  For track 2, we used only one acoustic model (f), we performed GSS twice with two individual arrays and once with all arrays together, followed by lattice combination. Diarization output was shared by all input array signals, and minimum Bayes risk decoding \cite{xu2011minimum} was applied on top of the combined lattice.

\section{CHiME Challenge Results}
\label{sec:chime-results}

We show improvement in WER for tracks 1 and 2, obtained using the different modifications described in the previous sections, in Table~\ref{tab:asr_results} and Table~\ref{tab:asr_results2}, respectively.

From Table~\ref{tab:asr_results}, we see that using a larger CNN + TDNN-F acoustic models improved performance over the TDNN-F baseline system by approx. 2\% absolute. For further improvement, we trained the AM on data closer to test conditions by augmenting the clean worn-microphone data with beamformed and GSS-enhanced data \cite{zorila2019investigation}. This provided almost 5\% absolute WER reduction. Adding the overlap bit provided a modest improvement. Lattice combination and LM rescoring were individually effective, and their combination provided a significant WER reduction of 4\% absolute.  Cumulatively, we obtained more than 10\% absolute WER improvement over the Challenge baseline system for track 1.

Table~\ref{tab:asr_results2} shows a similar step-wise WER improvement for track 2. Again, we obtained 2\% improvement using the larger CNN + TDNN-F architecture. Using multi-array fusion techniques in SAD and first-pass diarization reduced missed speech and speaker confusion, resulting in additional 4\% and 2\% improvement on the Dev and Eval sets, respectively. Using GSS and overlap-aware VB-HMM diarization provided significant improvements of 5\%--8\%, since they permit separating the overlapped speech for ASR. Finally, LM rescoring and lattice combination of ASR output from multiple input streams (but one ASR system) provided additional gains, similar to those observed in track 1 from combining multiple ASR systems.  Cumulatively, we obtained an absolute WER improvement of 17\% and 10\% respectively on the Dev and Eval sets.

\begin{table}[t]
\centering
\caption{Stepwise improvement in WER on Track 1.}
\label{tab:asr_results}
\begin{adjustbox}{max width=\linewidth}
\begin{tabular}{lcc}
\toprule
\textbf{System components (acoustic model)}  &\textbf{Dev}   & \textbf{Eval}\\
\midrule
Baseline TDNN-F (a)                            & 51.8\%        & 51.3\%\\
CNN-TDNN-F (b)                          & 49.6\%        & 49.4\%\\
+ Data selection and augmentation (f)    & 44.5\%        & 44.9\%\\
+ Overlap feature (i)                    & 44.4\%        & 44.5\%\\
+ RNN LM rescoring (i)                   & 42.8\%        & 42.9\%\\
+ Lattice combination (f)+(g)+(h)+(i)    & 41.8\%        & 42.1\%\\
+ Lattice combination \& LM rescoring         & 40.3\%        & 40.5\%\\
\bottomrule
\end{tabular}
\end{adjustbox}
\end{table}

\begin{table}
\centering
\caption{Stepwise improvement in WER on Track 2.}
\label{tab:asr_results2}
\begin{adjustbox}{max width=\linewidth}
\begin{tabular}{lcc}
\toprule
\textbf{System components (acoustic model)}   &\textbf{Dev}   &\textbf{Eval}\\
\midrule
Baseline TDNN-F (a)                            & 84.3\%        & 77.9\%\\
CNN-TDNN-F (f)                          & 82.5\%        & 75.8\%\\
+ Multi-array SAD \& PLDA fusion (f)     & 78.3\%        & 73.6\%\\
+ Multi-array GSS (f)                    & 71.0\%        & 68.8\%\\
+ VB Overlap Assignment (f)              & 69.3\%        & 68.8\%\\
+ RNN LM Rescoring (f)                   & 68.7\%        & 67.9\%\\
+ Lattice combination (f)                & 68.3\%        & 68.3\% \\
+ Lattice combination \& LM rescoring (f)   & 67.8\%        & 67.5\%\\
\bottomrule
\end{tabular}
\end{adjustbox}
\end{table}

\section{Conclusion}

We described our system for the sixth CHiME challenge for distant multi-microphone conversational speaker diarization and speech recognition in everyday home environments. We explored several methods to incorporate multi-microphone and multi-array information for speech enhancement, diarization, and ASR. For track 1, most of the improvements in WER were obtained from data selection and augmentation, and language model rescoring. Through careful training data selection, we reduced the training time of the system 3-fold while also improving its performance. In track 2, array fusion and overlap handling in the diarization module provided more accurate speaker segments than the Challenge baseline, resulting in improved speech enhancement via multi-array GSS. The gains from acoustic modeling and RNNLM rescoring developed in track 1 also largely carried over to track 2.

\section{Acknowledgement}
We thank Daniel Povey for the tremendous support throughout the system development, Jing Shi for help with data simulation in our trials with neural beamforming, and Yiming Wang for help with the lattice combination setup. This work was partially supported by grants from the JHU Applied Physics Laboratory, Nanyang Technological University, and the Government of Israel, and an unrestricted gift from \href{https://www.apptek.com/}{Applications Technology (AppTek)} Inc.

\bibliographystyle{IEEEtran}