\section{Introduction}

The  geometry and topology of  mechanical linkages play an important and increasing role in applied problems.
Most of the previous studies were concerned with the workspace and the topology of the configuration space,
which is only known in a few cases summarized in \cite{KM}, \cite{Oh}, \cite{Mo}. The aim of this paper is
to extend the known results to new classes of linkages and enrich them by considering potential functions on
the workspace, with aplications to control of the linkages considered. Our approach is based on Morse theory
which yields an explicit connection between the topology of configuration space and the critical points of potential.

We will be basically concerned with the so-called arachnoid mechanisms, the topology of which is a largely unexplored topic.
Nowadays the same type of objects is often called "spidery linkage". Detailed information on the topology of arachnoid
linkages is important for the design and control of certain types of spider robots. More concretely, we study the simplest
arachnoid mechanism, the $3$-leg spider also known as tripod spider \cite{Oh}, \cite{Mo}.

In Section \ref{s:hooke} we  use Hooke energy
as a potential function and describe in detail the workspace and the critical point theory. This implies, in particular, that
a weighted version of Hooke energy can be used to control such a linkage.

In section \ref{s:coulomb} we study the Coulomb potential of point charges placed at the foots of a $3$-leg spider linkage. We deal with the so-called
stationary charges of the spider's center (the common point of legs)  and the so-called trapping domain of stationary charges
determined in \cite{GK1},  \cite{GK2}. This enables us to determine the domain in the workspace of spider, where the position of
spider's center can be robustly controlled by the values of stationary charges using the Coulomb control scenario developed in \cite{KPS},
\cite{GK1}. For a symmetric spider based on a regular triangle and having a contractible workspace, we show that the domain of
robust Coulomb control is a non-void open subset of the workspace containing the center of the reference triangle.

In the sequel we use Morse theory for manifolds with boundaries and corners. Most of it is "folklore" hidden in the literature.
We mention here \cite{JJT} and \cite{GM} as the most general reference. For critical points of functions on manifolds with
boundary and corners we refer to \cite{Si}. For brevity, we refer to criteria of such critical points and corresponding
topological changes in level surfaces as "standard rules".

In conclusion we mention several related problems and perspectives suggested by our results. In general, this paper may
be considered as a first step in applying our approach to spider linkages and creating a paradigm for further research
in this direction.

The results of this paper were obtained and written up in the framework of a "Research in Residence" project at
the "Centre des Rencontres Mathematiques" (CIRM, Luminy, France) realized in November of 2022. It is our pleasure
to acknowledge the support and excellent working conditions at CIRM which largely facilitated our research.

\section{Hooke Energy as a potential function of tripod spider}\label{s:hooke}

We will consider the Hooke  potential for three points  in  several situations.  First without any constraint, next with constraints on maximal an minimal distances, and finally  for a 3-leg spider.

\subsection{No constraints}
We start with the following simplified situation. Given $3$ points $A,B,C$ (in vector notation $a,b,c$), the foot points  of the  spider,
and its center (joint) the  point $X$. The legs $AX$, $BX$, $CX$ are completely flexible around the foot, their length is allowed to change.

\begin{center}
\includegraphics[width=6cm]{vrijdag2a}
\end{center}

\smallskip
\noindent
The Hooke Energy  is defined by
$$  H(x) =  ||x- a||^2  +||x-b||^2 +  ||x- c||^2$$
The stationary points are determined by:
$\nabla H = 6 x - 2(a+b+c) = 0$, so $x = \frac{1}{3} (a+b+c)$, the center of gravity.
There are no other stationary points.
Note that $H(X) = ||x-z||^2 +  K,$ where $z=\frac{1}{3} (a+b+c)$, the center of gravity $Z$ and $K= ||a||^2 + ||b^||^2+||c||^2 - ||z||^2$.
All level curves are circles.

\subsection{Maximal length constraints}
We require
$$  |XA| \le R_A \; , \; |XB| \le R_B \;  , \; |XC|  \le  R_C.$$
The configuration space is the intersection of 3 discs $D_A, D_B ,D_C$ with centers resp $A,B,C$  and radii $R_A,R_B,R_C$. We study the boundary extrema of $H$.

\begin{lemma}
The map $H$ restricted to the boundary of $D_A$ has an extrema in $Y$ iff $AY$  has the direction of $AB + AC$,  i.e the  two intersection points of the line trough $A$ and $Z$ with the boundary of $D_A$, a minimum and a maximum.

\end{lemma}
\begin{proof}
Suppose $A$ is the origin.
Let $R_A e_{\phi}$ a point on $\delta D_A$. Use Lagrange multiplyers:
$$ \nabla H = 6 R_A e_{\phi} - 2 (b+c) = \lambda e_{\phi}$$
Therefore
$$ (- \lambda + 6 r_A) e_{\phi} = 2(b+c)$$
\end{proof}

\begin{proposition}
 Let  the workspace $W = D_A \cap D_B \cap D_C$ contain an open neighborhood of $Z$ then $Z$ is the only stationary point of $H$, an absolute minimum. There are no boundary singularities.
\end{proposition}
\begin{proof}
The center of gravity $Z$ is clearly a  minimum.
Potential other critical points are the critical points of the restiction of $H$ to the boundary circles (intersections of $ZA$ with the circle around $A$, etc) and also the 3 `corner points `, where 2 boundary circles intersect. The statement follows now from the standaard rules .
For boundary points, they are here  as follows:

\smallskip
\noindent
\begin{tabular}{llll}
type of restriction & direction normal & type & cell attaching \\
\hline
minimum    & inward  &  no change & no\\
minimum  & outward  & minimum & 0-cell\\
maximum    & inward  &  no change & no\\
maximum & outward  &  saddle & 1-cell \\
\end{tabular}

\smallskip
\noindent
With no change we mean, that the topological type of the lower level sets don't change.\\
For corner points there are similar rules. They give in our case: no change. Due to the special (circular) form of level curves and boundaries one can obtain the same result  by `inspection of pictures'.
\end{proof}

\subsection{Two sided length constraints}
We require
$$  0 < R_A^{-} \le |XA| \le R_A \; , \; 0 < R_B^{-} \le |XB| \le R_B \;  , \; 0 < R_C^{-} \le |XC|  < R_C .$$
The configuration space is the intersection of 3 discs with centers $A,B,C$  and radii $R_A,R_B,R_C$, where some smaller (open)  discs have been taken out. From the  many different posibilties, we consider here the case that these 3 discs have a relatively small radii and the workspace $W$ is a disc with 3 holes, containing an open neighborhood of the gravity center  $Z$.

\begin{proposition}\label{p:3holes}
 In this situation the point  $Z$ is an absolute minumum of $H$ on $W$  , moreover there are 3 saddle points on the (outer) intersection points of  the small discs with $ZA$, $ZB$, $ZC$.
\end{proposition}
\begin{proof}
The center of gravity $Z$ is an absolute minimum.
Other potential  critical points are the critical points of the restiction of $H$ to the boundary circle,  the 3 small circles and also the 3 `corner points `, where 2 boundary circles intersect.
The statement follows now from standard rules for boundary and corner singularities as explained above.

\end{proof}
\begin{remark}
The phase portrait  of the gradient of $H$  consists of  straight lines to $Z$, except for those  lines that intersect one of the small discs.
Their trajectories follow from the moment that they intersect these discs  a part of the boundary circle until they become `visible' from $Z$ and then continue as a straight line. There are 3 intervals between the outer and the inner circles, which are conflict strata for the gradient flow.
\end{remark}

\begin{remark}
The Morse theory in the 3 cases above is as follows:\\
In 2.1:  $b_0 =1, b_i = 0 \; (i \ge 1)$  and  $\mu_1=1 ,  \mu_i = 0  \;  (i \ge 1)$,\\
In 2.2:  $b_0 =1, b_i = 0 \; (i \ge 1)$  and  $\mu_1=1 , \mu_i = 0  \;  (i \ge 1)$,\\
In 2.3:  $b_0 =1, b_2= 3  ,  b_i = 0  \; (i \ge 2)$  and  $\mu_1=1, \mu_1 = 3 , \mu_i = 0  \; (i \ge 2)$.\\
In these three cases, $H$ is a perfect Morse function.
\end{remark}

\subsection{Robust control}
A similar study can be made for the weighted Hooke Energy:
$$ H_{\alpha,\beta,\gamma} =  \alpha ||x- a||^2  + \beta ||x-b||^2 + \gamma  ||x- c||^2$$
Assume  $\alpha > 0,  \beta >0,  \gamma > 0$. This potential function has the point $ Z = (\alpha : \beta : \gamma)  $ (barycentric coordinates) as an absolute minimum.  The level curves are circles with the center  at this point $Z$. The critical point theory is similar to the case of $H$, which corresponds to  $ (\alpha: \beta :\gamma ) = (1:1:1)$.

\smallskip
\noindent
A proper choice of the controls $(\alpha,\beta,\gamma)$ can be used to move the 3-leg spider  to any point in the triangle $ABC$ via minimum points of Hooke energy. This procedure yields a \textit{robust control} of the spider.

\subsection{The 3-leg spider}
In this case the telescopic connections are replaced by 2-arms with fixed arm lengths and flexible turning point.
We assume that the two parts of the arm have different lengths.
So $AX$ is replaced by the arm $AP\cdot PX$,  $BX$ by $BQ \cdot QX$ and $CX$ by $CR \cdot RX$. The configuration space $\mathcal C$ is an 8-fold cover of the workspace $W$ with certain identifications at the boundaries.
The topology has been studied in full generality by P. Mounod \cite{Mo}. He showed that an $n$-leg spider with generic arm lengths has a smooth two-dimensional configuration space, and gave a formula for its Euler characteristic.
In his paper he also solved some questions of J. O'hara \cite{Oh}.

The configuration space of 3-leg spider projects now clearly to the workspace $W$ from 2.3 (with the 3 small holes).
$R_A^-$ is equal to the difference of arm-lengths in the  $AP\cdot PX$ and $R_A$ is equal to their sum. Similar for the points $B$ and $C$.
We assume here that the arm lengths are  such, that they give the workspace mentioned above.

\smallskip
\noindent
 According to the formula of \cite{Mo} the Euler characteristic of $\mathcal{C}$ is -22. So we have a smooth  surface with genus 12.

We consider the quadratic distance function $H$ (see above).
\begin{proposition}
The Hooke energy $H$ on the configuration space $\mathcal{C}$ of the 3-leg spider has:
\begin{itemize}
\item 8  minima
\item 36 saddles
\item 6 maxima
\end{itemize}

\end{proposition}
\begin{proof}
The potential $H$ is defined on the configuration space via projection to the work space. We will use the results from Propostion  \ref{p:3holes} and
use (branched) covering arguments. The covering is 8-fold above the open part of the workspace. The branching takes place at the boundaries and corners. The potential critical points are just the  preimages of the special points, which have been studied in Proposition \ref{p:3holes}.

We find in this way $8$ critical pre-images of the center of gravity $Z$, these are all absolute minima (with the same value).
The $36$ saddles come from the $9$ special points on boundary circles, where we have covering degree $4$. Their types can be deduced
by topological study of the gluing of the two boundary pieces. Finally, we have $6$ maxima, which corespond to the three corner points,
where the covering is $2$-fold and at each such point such four pieces are glued together.

There is is a dictionary between the Morse theory on the workspace treated above and the configuration space. Compare the level curves
in both cases via gluing and this results in the Morse indices  mentioned in the  proposition.
\end{proof}

\section{Coulomb Energy as a potential function  of tripod spider} \label{s:coulomb}
From the viewpoint of control theory, it is also interesting and practically important to consider the Coulomb potential of point charges which are placed at the fixed foot points of a spider and can be varied in order to change the position of its body endowed with a fixed charge. A similar scenario was studied in big detail in the papers \cite{GK1}, \cite{GK2} and our exposition in this section relies on the constructions and results given in those two papers.

\smallskip
\noindent
We begin with recalling several concepts and constructions  in the form needed in the sequel. Recall that the Coulomb potential $E = E(Q@A)$ of $n$ point charges of non-zero magnitudes $q_i$ placed at $n$ points $A_i$ is a function of point $X$ in the complement of points $A_i$ by the formula
$$E = \sum \frac{q_i}{d_i}, $$
where $d_i = d(X, A_i)$ is the Euclidean distance between the points $P$ and $A_i$. As usual a point $X$ is called a stationary point (or equilibrium point) of potential $E$ if its gradient $\nabla E$ vanishes at $X$. There exists a huge number of results and problems concerned with the Coulomb potential and equilibrium points of point charges. For us it is important to mention that, as was proven by M. Morse himself, for a generic configuration and generic values of point charges Coulomb potential is a Morse function \cite{CM}. This suggests that Coulomb potential is a reasonable candidate for developing Morse theory for a spider linkage along the same lines as in the case of Hooke’s energy considered above. The aim of this section is to present a number of results in this direction and present an application in the spirit of robust Coulomb control discussed in \cite{KPS}, \cite{GK1}, \cite{GK2}.
Another relevant topic in our context is the so-called Maxwell’s conjecture on point charges stating that the number of isolated equilibrium points of $n$ point charges in $2$-dimensional Euclidean space does not exceed $(n-1)^2$ (see, e.g., a recent review in  \cite{GaNoSh}).

\smallskip
\noindent
We next focus on the case of 3 point charges. Given a triangle $\triangle$ with vertices $A, B, C$, referred to as a reference triangle (or base triangle),
and a point $X$ in the same plane, the triple of normalized stationary charges $Q(X) =  \{q_i(X; \triangle )\}$ is defined  as a triple of non-zero real
numbers $q_i$ such that $\nabla E(X,Q) = 0$. As was proven in \cite{GK1} these charges are uniquely determined by the ratio's
$$q_1:q_2:q_3 = d_1^3 \mathcal{A}_1 : d_2^3 \mathcal{A}_2 : d_3^2 \mathcal{A}_3,$$
where $\mathcal{A}_{1}$ is the area of $\triangle BCX$, $\mathcal{A}_2$ of $\triangle CAX$, $\mathcal{A}_3$ of $\triangle ABX$. We denote the angles at X of these triangles by $\alpha_1,\alpha_2,\alpha_3$ respectively.

\smallskip
\noindent
 In \cite{GK1}  , \cite{GK2}  \textit{ the Coulomb trapping domain} $T(\triangle)$  is defined as the set
of  minimum points of the potential induced by the charges $Q(X)$. The trapping domain is given by the inequality $h(X) > 0$, where $h$ is the Hessian of $E(X,Q(X))$, which we call the \textit{ trapping Hessian}. We have  in geometric terms the following formula:
$$ h(X) =  -2 \mathcal{A} +9 (\prod_{i}^3  \sin \alpha_i  )(\sum_{i=1}^3 d_i^2 \mathcal{A}_i).$$

It was conjectured in \cite{GK2} that
   $T(\triangle)$ is a non-empty subset containing the incenter of  $\triangle$. This conjecture has been verified in a number of cases, including
the regular triangle.

\smallskip
\noindent
As was mentioned, a complete description of equilibria (stationary points) of $E$ and their Morse indices is not known even in the case of three charges.
However rather complete results for special configurations of three point charges have been obtained in \cite{Ts}. In particular, the Maxwell’s conjecture
was proven for arbitrary three non-zero point charges at the vertices of a regular triangle and the Morse indices are computed.

\smallskip
\noindent
For this reason, in the rest of the text we always assume that the reference triangle is regular since it is the most important case for applications
and, at the same time, in this case all necessary background results are rigorously proven in \cite{GK2}.

\smallskip
\noindent
Further, given a $3$-leg spider $S$ with the center at (moving) point $X$, (fixed) foot-points $A, B, C$ and identical legs with the links of lengths $a$
(thigh length) and $b$ (foot length), we denote by $W(S)$ its workspace considered above.  To simplify the discussion we assume that $R >a > b > 0$,
in which case the whole workspace of the spider is the intersection of three  circular annuli with centers $A, B, C$.

\begin{proposition}
The domain of robust Coulomb control of spider $S$ in the above scenario is equal to the intersection $D(S) = W(S) \cap T(\triangle)$.
\end{proposition}

\noindent
As said before the topology of the workspace depends on the triangle and the arm lengths.
This influences the intersection which could even be void. For the regular triangle, we make the choice that the workspace becomes a contractible
region bounded by three circular arcs and containing the incenter.  Also we assume that all charges are positive.
The trapping domain $T(\triangle)$ is explicitly known as the interior of a compact convex region bounded by $h=0$ and containing
the incenter of the triangle. Under these conditions this yields the following conclusion.
\begin{corollary}
The domain of robust Coulomb control of spider $S$ based on a regular triangle is non-void.

\end{corollary}

\medskip
\noindent
For a spider based on the regular triangle one can also carry out a complete Morse analysis in a similar way as for the Hooke potential.
One knows that for given positive charges there is a single minimum of E in the trapping area $T(\triangle)$ and $3$ saddle points
in $\triangle$  ouside  $T(\triangle)$.  The poles of $E$ are not in the workspace and therefore not in the configuration space, which
is an $8$-fold cover of the workspace branched (glued) over boundary components. After obtaining information about the boundary singularities
one gets a complete Morse analysis.

\section{Concluding remarks}

There are several natural problems and research perspectives suggested by our results. We mention
some of them which seem most interesting and feasible.

First of all, note that all the concepts used in our paper make sense for $n$-leg spiders with arbitrary base $n$-gon
and arbitrary lengths of the links. So a natural next step is to generalize the above results to Hooke and Coulomb
potential of an $n$-leg spider.

Next, one can also consider other geometrically or physically meaningful potentials of $n$-leg spider like Riesz energies or
the oriented area of the polygon formed by the moving joints of spider.

Further, one can formulate several extremal problems related to the design of spider linkages with desirable properties
like the area or shape of the workspace.

Moreover, one can try to extend our results to the case of $n$-leg spiders with a moving $n$-gon platform instead of the center
point.

Finally, it is natural to extend our discussion to spatial spiders where the legs can move in a three-dimensional Euclidean space.

It is also possible to consider "legs" with more than two links or with flexible non-extendible tether links. Each of the mentioned
possibilities deserves a closer look and careful exploration which we intend to undertake in the future research.

Returning to $3$-leg spiders we add that an urgent topic is to clarify what happens with the domain of robust Coulomb control for
arbitrary configuration of foots. In this case, there are several feasible problems concerned with the optimization of workspace with
given restrictions on the sum of the lengths of links and circumradius of the base triangle. For example, there is good evidence
that the maximal area of $W(S)$ arises for the most symmetric spider with equal links and regular configuration of foots.


\section{Introduction}
\label{sec:1}
One of the most famous integrable systems is the \yadd{finite} Toda equation \cite{Flaschka_1974_1, Flaschka_1974_2} expressed using the so-called Flaschka's variables:
\begin{align}
\label{eqn:Toda}
\left\{\begin{aligned}
& \dfrac{dx_k(t)}{dt}=y_k(t)-y_{k-1}(t),\quad k=1,2,\dots,m,\\
& \dfrac{dy_k(t)}{dt}=y_k(t)\left(x_{k+1}(t)-x_k(t)\right),\quad k=1,2,\dots,m-1,\\
& y_0(t)\coloneqq 0,\quad y_m(t)\coloneqq 0,
\end{aligned}\right.
\end{align}
where $t$ and $k$ respectively denote continuous time and the spatial index.
The Toda equation \eqref{eqn:Toda} was first considered in a study of nonlinear spring dynamics
and the Toda variables $x_k(t)$ and $y_k(t)$ respectively correspond to the $k$th point-mass variation from its momentum and the equilibrium point at continuous-time $t$ \cite{Toda_1967}.
The Toda equation \eqref{eqn:Toda} was also studied from the perspective of
orthogonal polynomials \cite{Toda_1989} and rational approximations \cite{Nikishin_1991}.
Symes~\cite{Symes_1982} found an interesting relationship where the continuous-time evolution from $t$ to $t+1$ coincides with one step of the well-known $QR$ algorithm \cite{Golub_1996} in the case where target matrices are exponentials of symmetric tridiagonal matrices.
A time-discretization \cite{Hirota_1981} of the Toda equation \eqref{eqn:Toda} leads to
\begin{align}
\label{eqn:qd}
\left\{\begin{aligned}
& Q_k^{(n+1)}+E_{k-1}^{(n+1)}=Q_k^{(n)}+E_k^{(n)},\quad k=1,2,\dots,m,\\
& Q_k^{(n+1)}E_k^{(n+1)}=Q_{k+1}^{(n)}E_k^{(n)},\quad k=1,2,\dots,m-1,\\
& E_0^{(n)}\coloneqq 0,\quad E_m^{(n)}\coloneqq 0,
\end{aligned}\right.
\end{align}
where $Q_k^{(n)}$ and $E_k^{(n)}$ respectively denote the values of $Q_k$ and $E_k$ at discrete-time $n$.
Sogo also showed that the discrete Toda equation \eqref{eqn:qd} is simply the recursion formula of the quotient-difference (qd) algorithm \cite{Rutishauser_1990} for computing symmetric tridiagonal matrices.
Similarity transformations employed in the qd algorithm are $LR$ transformations that decompose symmetric tridiagonal matrices into the products of lower bidiagonal matrices $L^{(n)}$ and upper bidiagonal matrices $R^{(n)}$ and then generate new symmetric tridiagonal matrices as the products $R^{(n)}L^{(n)}$.
In other words, discrete-time evolutions in the discrete Toda equation \eqref{eqn:qd} give a sequence of $LR$ transformations for symmetric tridiagonal matrices.
See also \cite{Gautschi_2002} for the symmetric $LR$ transformations from the perspective of the orthogonal polynomials.
\par
An extension of the Toda equation \eqref{eqn:Toda} is the \yadd{Kostant-Toda} equation \cite{Kostant_1979}, and the hungry Toda equation is a special case of the \yadd{Kostant-Toda} equation \cite{Shinjo_2020}.
The $q$-analogue of the Toda equation \cite{Area_2018}, which involves the so-called $q$-parameter satisfying $0<q<1$, can also be regarded as an extension of the Toda equation \eqref{eqn:Toda}.
This is because taking the limit $q\to 1$ in the $q$-Toda equation leads to the Toda equation \eqref{eqn:Toda}.
The \yadd{Kostant-Toda} equation including the hungry Toda equation is associated with matrix eigenvalue problems,
and a time-discretization of the hungry Toda equation is shown to be applicable to computing eigenvalues of totally nonnegative matrices
whose minors are all positive \cite{Fukuda_2013}.
However, to the best of our knowledge, the relationships of the $q$-Toda equation and its time-discretization
to eigenvalue problems have not yet been clarified.
The main purpose of this paper is thus to associate the $q$-Toda equation -- or more precisely, its further extension --
with eigenvalue problems and then present applications of its time-discretization to computing eigenvalues.
Constructing matrix representations called Lax representations of integrable systems is useful in finding the relationships to eigenvalue problems.
In particular, in the case of discrete integrable systems, the Lax representations are usually related to $LR$ transformations directly.
The key point of this paper is not to be too particular about $LR$ transformations
in designing Lax representations of the extended $q$-Toda equation and its time-discretization.
Based on the resulting Lax representations, we examine distinctive similarity transformations
generated from the extended $q$-discrete Toda equation.
\par
An integrable system closely related to the Toda equation is the Lotka--Volterra (LV) system,
which describes predator--prey interactions
\begin{align}
\label{eqn:LV}
\left\{\begin{aligned}
& \dfrac{du_k(t)}{dt}=u_k(t)\left(u_{k+1}(t)-u_{k-1}(t)\right),\quad k=1,2,\dots,2m-1,\\
& u_0(t)\coloneqq 0,\quad u_{2m}(t)\coloneqq 0,
\end{aligned} \right.
\end{align}
where $k$ denotes species index and $u_{k}(t)$ corresponds to the number of the $k$th species at continuous-time $t$.
The hungry Lotka--Volterra (hLV) system \cite{Bogoyavlenskii_1991,Itoh_1987} is an extension of the LV system \eqref{eqn:LV}.
The Darboux transformations on band matrices were studied in \cite{Adler_2000,Rolania_2015}, and then enable us to relate the \yadd{Kostant-Toda} equation to the discrete Korteweg de Vries equation which is equivalent to the hLV system \cite{Rolania_2019}.
See also \cite{Rolania_2009} for other extensions of the Toda equation \eqref{eqn:Toda} and LV system \eqref{eqn:LV}.
A time-discretization of the LV (dLV) system
\begin{align}
\label{eqn:dLV}
\left\{\begin{aligned}
& u_k^{(n+1)}=u_k^{(n)}\dfrac{1+u_{k+1}^{(n)}}{1+u_{k-1}^{(n+1)}},\quad k=1,2,\ldots,2m-1, \\
& u_0^{(n)}\coloneqq 0,\quad u_{2m}^{(n)}\coloneqq 0.
\end{aligned}\right.
\end{align}
also has a close relationship to the discrete Toda equation.
The dLV system \eqref{eqn:dLV} can be obtained from a variable transformation called the B\"acklund transformation or the Miura transformation of the discrete Toda equation \eqref{eqn:qd}.
We designed an algorithm for computing singular values of bidiagonal matrices,
which is equivalent to eigenvalues of symmetric positive-definite tridiagonal matrices based on the dLV system~\cite{Iwasaki_2002}.
Similarly to the Toda case, the LV case has not been studied from the viewpoint of eigenvalue problems.
In this paper, we thus focus on applications of the $q$-LV case to computing eigenvalues.
\par
The remainder of this paper is organized as follows.
In Section \ref{sec:2}, we briefly explain the $q$-Toda equation and then find Lax representations of its extension in order to associate it with eigenvalue problems.
Next, in Section \ref{sec:3}, we consider a time-discretization of the extended $q$-Toda equation and relate it to a sequence of similarity transformations of matrices.
In Section \ref{sec:4}, we investigate the case of the extended $q$-LV system.
In Section \ref{sec:5}, we clarify a determinantal solution to the extended $q$-discrete Toda equation and its asymptotic behavior as discrete-time goes to infinity.
We also present numerical examples for confirming the convergence to matrix eigenvalues.
Finally, in Section \ref{sec:6}, we give concluding remarks.
\section{Extended $q$-analogue of the Toda equation}
\label{sec:2}
In this section, we first consider an extension of the $q$-analogue of the Toda equation \eqref{eqn:Toda}.
Next, by considering two kinds of Lax representations,
we associate the extended $q$-Toda equation with matrix eigenvalue problems.
\par
The $q$-derivative of the function $f(t)$
is a $q$-analogue of the ordinary derivative $df(t)/dt$ defined by:
\begin{align}
\label{eqn:q_derivative}
D_qf(t)\coloneqq \left\{
\begin{aligned}
& \dfrac{f(t)-f(qt)}{(1-q)t},\quad t\neq 0,\\
& \lim_{n\to\infty}\dfrac{f(q^nt)-f(0)}{q^nt},\quad t=0,
\end{aligned}\right.
\end{align}
where $0<q<1$. Area et al. \cite{Area_2018} presented a $q$-analogue of the Toda equation:
\begin{align}
\label{eqn:qToda}
\left\{\begin{aligned}
& D_qx_k(t)=g_k(qt)-g_{k-1}(qt),\quad k=1,2,\dots,m,\\
& D_qy_k(t)=g_k(qt)\left(x_{k+1}(qt)-x_k(t)\right),\quad k=1,2,\dots,m-1,\\
& y_0(t)\coloneqq 0,\quad y_m(t)\coloneqq 0,
\end{aligned}\right.
\end{align}
where $g_k(qt)$ are auxiliary variables given by:
\begin{align}
\label{eqn:qToda_g}
\left\{\begin{aligned}
& g_1(qt)=\dfrac{y_1(qt)}{1+(1-q)tx_1(qt)},\\
& g_k(qt)=\dfrac{y_k(qt)}{y_{k-1}(t)}g_{k-1}(qt),\quad k=2,3,\dots,m-1,\\
& g_0(qt)\coloneqq 0,\quad g_m(qt)\coloneqq 0.
\end{aligned}\right.
\end{align}
Since $D_qx_k(t)\to dx_k(t)/dt$, $D_qy_k(t)\to dy_k(t)/dt$ and $g_k(qt)\to y_k(t) $ as $q \to 1$,
the $q$-Toda equation \eqref{eqn:qToda} with $q\to 1$ leads to the Toda equation \eqref{eqn:Toda}.
\par
Now we introduce an arbitrary integer $M$ into the $q$-Toda equation \eqref{eqn:qToda} as
\begin{align}
\label{eqn:extended_qToda}
\left\{\begin{aligned}
& D_qx_k(t)=g_k(qt)-g_{k-M-1}(qt),\quad k=1,2,\dots,M_m+M,\\
& D_qy_k(t)=g_k(qt)\left(x_{k+M+1}(qt)-x_k(t)\right),\quad k=1,2,\dots,M_m-1,\\
& y_{-M}(t)\coloneqq 0,\quad y_{-M+1}(t)\coloneqq 0,\quad\dots,\quad y_0(t)\coloneqq 0,\\
& y_{M_m}(t)\coloneqq 0,\quad y_{M_m+1}(t)\coloneqq 0,\quad\dots,\quad y_{M_m+M}(t)\coloneqq 0,
\end{aligned}\right.
\end{align}
where $M_i\coloneqq (M+1)i-M$ and $g_k(qt)$ satisfy
\begin{align}
\label{eqn:extended_qToda_g}
\left\{\begin{aligned}
& g_k(qt)=\dfrac{y_k(qt)}{1+(1-q)tx_k(qt)},\quad k=1,2,\dots,M+1,\\
& g_k(qt)=\dfrac{y_k(qt)}{y_{k-M-1}(t)}g_{k-M-1}(qt),\quad k=M+2,M+3,\dots,M_m-1,\\
& g_{-M}(qt)\coloneqq 0,\quad g_{-M+1}(qt)\coloneqq 0,\quad\dots,\quad g_0(qt)\coloneqq 0,\\
& g_{M_m}(qt)\coloneqq 0,\quad g_{M_m+1}(qt)\coloneqq 0,\quad\dots,\quad g_{M_m+M}(qt)\coloneqq 0.
\end{aligned}\right.
\end{align}
Equation \eqref{eqn:extended_qToda} with $M=0$ is simply the $q$-Toda equation \eqref{eqn:qToda}.
Thus, we hereinafter refer to \eqref{eqn:extended_qToda} as the extended $q$-Toda equation.
Similarly to the simple Toda equation \eqref{eqn:Toda}, we can regard that the extended $q$-Toda equation \eqref{eqn:extended_qToda} describes a system of nonlinear springs. The extended $q$-Toda case however differs from the Toda equation \eqref{eqn:Toda} in that the $k$th point-mass links with the $(k-M-1)$st and $(k+M+1)$st point-masses rather than the $(k-1)$st and $(k+1)$st point-masses.
\par
We prepare $(M_m+M)\times(M_m+M)$ matrices involving the extended $q$-Toda variables
$x_k(t)$, $y_k(t)$, and $g_k(t)$ as:
\[\begin{aligned}
& A(t)\coloneqq\left(\begin{array}{ccccc}
x_1(t) & \overbrace{0\ \cdots\ 0}^M & 1 & & \\
\begin{array}{c} 0 \\ \vdots \\ 0 \end{array}
& \ddots & \ddots & \ddots & \\
y_1(t) & \ddots & & \ddots & 1 \\
 & \ddots & \ddots & \ddots &
\begin{array}{c} 0 \\ \vdots \\ 0 \end{array} \\
 & & y_{M_m-1}(t) & \underbrace{0\ \cdots\ 0}_M & x_{M_m+M}(t)
\end{array}\right),\\
& G(t)\coloneqq\left(\begin{array}{ccccc}
0 & & \\\begin{array}{c} 0 \\ \vdots \\ 0 \end{array}
 & \ddots & & \\
g_1(t) & \ddots & \ddots & \\
 & \ddots & \ddots & \ddots \\
 & & g_{M_m-1}(t) & \underbrace{0\ \cdots\ 0}_M & 0
\end{array}\right).
\end{aligned}\]
We hereinafter call $A(t)$ and $G(t)$ the $(M+1)$-tridiagonal and $(M+1)$-subdiagonal matrices, respectively.
The following proposition then gives the Lax representation
for the extended $q$-Toda equation \eqref{eqn:extended_qToda}.
\begin{proposition}
\label{prop:Lax_extended_qToda}
For the $(M+1)$-tridiagonal matrices $A(t)$ and the $(M+1)$-subdiagonal matrices $G(t)$, it holds that
\begin{align}
\label{eqn:Lax_extended_qToda}
D_qA(t)=A(qt)G(qt)-G(qt)A(t).
\end{align}
\end{proposition}
\begin{proof}
Observing the $(k,k)$, $(k+M+1,k)$ and $(k+2M+2,k)$ entries of $A(qt)G(qt)-G(qt)A(t)$, we easily derive
$g_k(qt)-g_{k-M-1}(qt)$, $g_k(qt)x_{k+M+1}(qt)-g_k(qt)x_k(t)$, and
$g_k(qt)y_{k+M+1}(qt)-g_{k+M+1}(qt)y_k(t)$, respectively.
The other entries of $A(qt)G(qt)$ and $G(qt)A(t)$ are all $0$.
From the first and second equations of the extended $q$-Toda equation \eqref{eqn:extended_qToda},
it is obvious that $g_k(qt)-g_{k-M-1}(qt)=D_qx_k(t)$ and $g_k(qt)x_{k+M+1}(qt)-g_k(qt)x_k(t)=D_qy_k(t)$.
Moreover, from the second equation of \eqref{eqn:extended_qToda_g},
we see that $g_k(qt)y_{k+M+1}(qt)-g_{k+M+1}(qt)y_k(t)=0$.
Taking into account that the $(k+2M+2,k)$ entries of $D_qA(t)$ are all 0,
we thus have \eqref{eqn:Lax_extended_qToda}.
\end{proof}
For convenience, we here replace $x_{M_k+s}(t)$, $y_{M_k+s}(t)$, and $g_{M_k+s}(qt)$
with $x_{k,s}(t)$, $y_{k,s}(t)$, and $g_{k,s}(qt)$, respectively.
We also obtain another Lax representation for the extended $q$-Toda equation \eqref{eqn:extended_qToda} in terms of $m\times m$ tridiagonal and $1$-subdiagonal matrices:
\[\begin{aligned}
& {\cal A}_s(t)\coloneqq\left(\begin{array}{cccc}
x_{1,s}(t) & 1 & & \\
y_{1,s}(t) & x_{2,s}(t) & \ddots & \\
 & \ddots & \ddots & 1 \\
 & & y_{m-1,s}(t) & x_{m,s}(t)
 \end{array} \right), \\
& {\cal G}_{s}(t)\coloneqq\left(\begin{array}{ccccc}
0 & & \\
g_{1,s}(t) & \ddots & \\
 & \ddots & \ddots \\
 & & g_{m-1,s}(t) & 0\
\end{array} \right).
\end{aligned}\]
\begin{proposition}
\label{prop:Lax_extended_qToda_1}
For the block diagonal matrices
$ {\cal A}(t)\coloneqq {\rm diag}({\cal A}_0(t),{\cal A}_1(t),\dots$, ${\cal A}_M(t))$ and
$ {\cal G}(t)\coloneqq {\rm diag}({\cal G}_0(t),{\cal G}_1(t),\dots,{\cal G}_M(t))$, it holds that
\begin{align}
\label{eqn:Lax_extended_qToda_1_all}
D_q{\cal A}(t)={\cal A}(qt){\cal G}(qt)-{\cal G}(qt){\cal A}(t).
\end{align}
\end{proposition}
\begin{proof}
Obviously, $D_q{\cal A}(t)$ and ${\cal A}(qt){\cal G}(qt)-{\cal G}(qt){\cal A}(t)$ are both block diagonal matrices.
Observing the $s$th diagonal blocks of $D_q{\cal A}(t)$ and ${\cal A}(qt){\cal G}(qt)-{\cal G}(qt){\cal A}(t)$, we derive
\begin{align}
\label{eqn:Lax_extended_qToda_1}
D_q{\cal A}_s(t)={\cal A}_s(qt){\cal G}_s(qt)-{\cal G}_s(qt){\cal A}_s(t),\quad s=0,1,\dots,M.
\end{align}
It is worth noting here that ${\cal A}_s(t)$ and ${\cal G}_s(t)$ have the same form as $A(t)$ and $G(t)$, respectively.
Thus, from Proposition \ref{prop:Lax_extended_qToda},
we immediately obtain equalities of nonzero entries in \eqref{eqn:Lax_extended_qToda_1}:
\[\left\{\begin{aligned}
& D_qx_{k,s}(t)=g_{k,s}(qt)-g_{k-1,s}(qt),\quad k=1,2,\dots,m,\quad s=0,1,\dots,M,\\
& D_q y_{k,s}(t)= g_{k,s}(qt)\left(x_{k+1,s}(qt)-x_{k,s}(t)\right),\quad k=1,2,\dots,m-1,\quad s=0,1,\dots,M,
\end{aligned}\right.\]
which is equivalent to the extended $q$-Toda equation \eqref{eqn:extended_qToda} under the replacements
$x_{k,s}(t)=x_{M_k+s}(t)$, $y_{k,s}(t)=y_{M_k+s}(t)$, and $g_{k,s}(t)=g_{M_k+s}(t)$.
Therefore we conclude that \eqref{eqn:Lax_extended_qToda_1_all} is also
the Lax representation of the extended $q$-Toda equation \eqref{eqn:extended_qToda}.
\end{proof}
We here denote an eigenvalue of the $(M+1)$-tridiagonal matrix $A(t)$ by $\lambda_k$
and the corresponding eigenvector by $\Phi_k(t)$.
From the $q$-derivative of the equality $A(t)\Phi_k(t)=\lambda_k\Phi_k(t)$, we then obtain:
\[D_qA(t)\cdot\Phi_k(t)+A(qt)\cdot D_q\Phi_k(t)=\lambda_kD_q\Phi_k(t).\]
Assuming that $D_q\Phi_k(t)=-G(qt)\Phi_k(t)$,
we derive the Lax representation \eqref{eqn:Lax_extended_qToda} again.
As the so-called Lax pair for the extended $q$-Toda equation \eqref{eqn:extended_qToda}, we thus have:
\begin{align}
\label{eqn:Lax_pair_qToda}
\begin{cases}
A(t)\Phi_k(t)=\lambda_k\Phi_k(t),\quad k=1,2,\dots,M_m+M,\\
D_q\Phi_k(t)=-G(qt)\Phi_k(t),\quad k=1,2,\dots,M_m+M.
\end{cases}
\end{align}
Similarly, we obtain the Lax pair associated with the Lax representation \eqref{eqn:Lax_extended_qToda_1}:
\begin{align}
\label{eqn:Lax_pair_qToda_1}
\begin{cases}
{\cal A}_s(t)\Psi_{k,s}(t)=\sigma_{k,s}\Psi_{k,s}(t),\quad k =1,2,\dots,m,\\
D_q\Psi_{k,s}(t)=-{\cal G}_s(qt)\Psi_{k,s}(t),\quad k=1,2,\dots,m,
\end{cases}
\end{align}
where $(\sigma_{k,s},\Psi_{k,s}(t))$ denote eigenpairs of the block diagonal matrices ${\cal A}_s(t)$.
According to Nagata et al.~\cite{Nagata_2016},
for the $(M+1)$-tridiagonal matrices $A(t)$ and the block diagonal matrices $\Psi(t)$,
there exists a permutation matrix $P$ such that $P^{\top}A(t)P=\Psi(t)$.
This implies that eigenvalues of $A(t)$ and $\Psi(t)$ are equal to each other.
Namely, every $\sigma_{k,s}$ coincides with one of $\lambda_1,\lambda_2,\dots,\lambda_{M_m+M}$.
\section{$q$-discretization and similarity transformations}
\label{sec:3}
In this section, we present a time-discretization of the extended $q$-Toda equation \eqref{eqn:extended_qToda},
and then relate it to similarity transformations of $M$-tridiagonal matrices and block diagonal matrices.
\par
We introduce new variables $t^{(0)},t^{(1)},\dots$ sequentially given by $t^{(n)}\coloneqq q^{-1} t^{(n-1)}$, where $t^{(0)}>0$.
The extended $q$-Toda equation \eqref{eqn:extended_qToda} with $t=t^{(n)}$ then leads to:
\begin{align}
\label{eqn:extended_qdToda}
\left\{\begin{aligned}
& \dfrac{x_{k,s}^{(n+1)}-x_{k,s}^{(n)}}{(1-q)t^{(n)}}=g_{k,s}^{(n)}-g_{k-1,s}^{(n)},\quad k=1,2,\dots,m,\quad s=0,1,\dots,M,\\
& \dfrac{y_{k,s}^{(n+1)}-y_{k,s}^{(n)}}{(1-q)t^{(n)}}=g_{k,s}^{(n)}\left(x_{k+1,s}^{(n)}-x_{k,s}^{(n+1)}\right),\\
& \qquad k=1,2,\dots,m-1,\quad s=0,1,\dots,M,\\
& y_{0,s}^{(n)}\coloneqq 0,\quad y_{m,s}^{(n)}\coloneqq 0,
\end{aligned}\right.
\end{align}
where $x_{k,s}^{(n+1)}\coloneqq x_{k,s}(t^{(n)})$, $y_{k,s}^{(n+1)}\coloneqq y_{k,s}(t^{(n)})$,
$g_{k,s}^{(n+1)}\coloneqq g_{k,s}(t^{(n)})$, and $g_{k,s}^{(n)}$ satisfy
\begin{align}
\label{eqn:extended_qdToda_g}
\left\{\begin{aligned}
& g_{1,s}^{(n)}=\dfrac{y_{1,s}^{(n)}}{1+(1-q)t^{(n)}x_{1,s}^{(n)}},\quad s=0,1,\dots,M,\\
& g_{k,s}^{(n)}=\dfrac{y_{k,s}^{(n)}}{y_{k-1,s}^{(n+1)}}g_{k-1,s}^{(n)},
\quad k=2,3,\dots,m-1,\quad s=0,1,\dots,M,\\
& g_{0,s}^{(n)}\coloneqq 0,\quad g_{m,s}^{(n)}\coloneqq 0.
\end{aligned}\right.
\end{align}
We hereinafter refer to \eqref{eqn:extended_qdToda} as the extended $q$-discrete Toda equation.
It is emphasized here that the values of $x_{k,s}^{(n+1)}$, $y_{k,s}^{(n+1)}$ and $g_{k,s}^{(n+1)}$
are uniquely computed from those of $x_{k,s}^{(n)}$, $y_{k,s}^{(n)}$ and $g_{k,s}^{(n)}$
in the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda}.
Figure \ref{fig:diagram_time_evolution} shows the discrete-time evolution from $n$ to $n+1$
in the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda}.
\begin{figure}[tb]
\centering
\includegraphics[width=0.7\textwidth]{q-dToda_diagram}
\caption{Discrete-time evolution in the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda}. \label{fig:diagram_time_evolution}}
\end{figure}
\par
By setting $A^{(n+1)}\coloneqq A(t^{(n)}) $ and $G^{(n+1)}\coloneqq G(t^{(n)})$
in the Lax representation \eqref{eqn:Lax_extended_qToda},
we easily derive a Lax representation for the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda}.
\begin{theorem}
\label{thm:similar_transformation_qdToda}
For the $(M+1)$-tridiagonal matrices $A^{(n)}$ and the $(M+1)$-subdiagonal matrices $G^{(n)}$, it holds that
\begin{align*}
(I_{M_m+M}+(1-q)t^{(n)}G^{(n)})A^{(n+1)}=A^{(n)}(I_{M_m+M}+(1-q)t^{(n)}G^{(n)}),
\end{align*}
where $I_{M_m+M}$ denotes the $(M_m+M)\times(M_m+M)$ identity matrix.
\end{theorem}
\noindent
Since $\det(I_{M_m+M}+(1-q)t^{(n)}G^{(n)})=1$, the inverse of $I_{M_m+M}+(1-q)t^{(n)}G^{(n)}$ always exist.
Theorem \ref{thm:similar_transformation_qdToda} thus implies that
the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda}
generates a sequence of similarity transformations of the initial $A^{(0)}$.
\par
Now we prepare auxiliary variables given using $x_{k,s}^{(n)}$, $y_{k,s}^{(n)}$, $g_{k,s}^{(n)}$, $t^{(n)}$, and $q$ as
\begin{align}
\label{eqn:auxiliary_variable_qdToda}
\left\{\begin{aligned}
& Q_{k,s}^{(n)}\coloneqq x_{k,s}^{(n)}+\dfrac{1}{(1-q)t^{(n)}}-(1-q)t^{(n)}g_{k-1,s}^{(n)},\\
& \qquad k=1,2,\dots,m,\quad s=0,1,\dots,M,\\
& E_{k,s}^{(n)}\coloneqq (1-q)t^{(n)}g_{k,s}^{(n)},\quad k=1,2,\dots,m,\quad s=0,1,\dots,M,
\end{aligned}\right.
\end{align}
where $E_{0,s}^{(n)}\coloneqq 0$ and $E_{m,s}^{(n)}\coloneqq 0 $.
The following proposition then describes discrete-time evolution from $n$ to $n+1$ of $Q_{k,s}^{(n)}$ and $E_{k,s}^{(n)}$.
\begin{proposition}
\label{prop:auxiliary_recursion_formula_qdToda}
For $s=0,1,\dots,M$, it holds that
\begin{align}
\label{eqn:auxiliary_recursion_formula_qdToda}
\left\{\begin{aligned}
& Q_{k,s}^{(n+1)}+E_{k-1,s}^{(n+1)}-\dfrac{1}{(1-q)t^{(n+1)}}=Q_{k,s}^{(n)}+E_{k,s}^{(n)}-\dfrac{1}{(1-q)t^{(n)}},
 \quad k=1,2,\dots,m,\\
& Q_{k,s}^{(n+1)}E_{k,s}^{(n+1)}=Q_{k+1,s}^{(n)}E_{k,s}^{(n)},\quad k=1,2,\dots,m-1.
\end{aligned}\right.
\end{align}
\end{proposition}
\begin{proof}
Observing the first equation in the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda}, we obtain
\begin{align}
\label{eqn:auxiliary_recursion_formula_qdToda_Q}
\left\{\begin{aligned}
& Q_{k,s}^{(n)}=x_{k,s}^{(n)}-E_{k-1,s}^{(n)}+\dfrac{1}{(1-q)t^{(n)}}, \\
& Q_{k,s}^{(n)}=x_{k,s}^{(n+1)}-E_{k,s}^{(n)}+\dfrac{1}{(1-q)t^{(n)}}.
\end{aligned}\right.
\end{align}
Equation \eqref{eqn:auxiliary_recursion_formula_qdToda_Q} immediately
leads to the first equation of \eqref{eqn:auxiliary_recursion_formula_qdToda}.
\par
By dividing both sides of the second equation of the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda}
by $g_{k,s}^{(n)}$, recalling that $g_{k+1,s}^{(n)}y_{k,s}^{(n+1)}=g_{k,s}^{(n)}y_{k+1,s}^{(n)}$,
and using \eqref{eqn:auxiliary_variable_qdToda}, we derive:
\begin{align}
\label{proof:E_q_dToda_qd}
Q_{k+1,s}^{(n)}-\dfrac{y_{k+1,s}^{(n)}}{E_{k+1}^{(n,s)}}=Q_{k,s}^{(n)}-\dfrac{y_{k,s}^{(n)}}{E_{k,s}^{(n)}}.
\end{align}
From \eqref{proof:E_q_dToda_qd}, it follows that:
\[Q_{k+1,s}^{(n)}-\dfrac{y_{k+1,s}^{(n)}}{E_{k+1}^{(n,s)}}=Q_{1,s}^{(n)}-\dfrac{y_{1,s}^{(n)}}{E_1^{(n,s)}}.\]
Since it is obvious from the first equation in \eqref{eqn:extended_qdToda_g}
and \eqref{eqn:auxiliary_variable_qdToda} with $k=1$ that $Q_{1,s}^{(n)}-y_{1,s}^{(n)}/E_1^{(n,s)}=0$, we obtain:
\[Q_{k,s}^{(n)}-\dfrac{y_{k,s}^{(n)}}{E_{k,s}^{(n)}}=0.\]
The second equation in \eqref{eqn:extended_qdToda_g}
and the second equation in \eqref{eqn:auxiliary_variable_qdToda} also yield
$y_{k,s}^{(n+1)}/E_{k+1,s}^{(n)}=y_{k+1,s}^{(n)}/E_{k,s}^{(n)}$.
Thus, we have:
\begin{align}
\label{eqn:auxiliary_recursion_formula_qdToda_E}
\left\{\begin{aligned}
& Q_{k,s}^{(n)}-\dfrac{y_{k,s}^{(n)}}{E_{k,s}^{(n)}}=0,\\
& Q_{k,s}^{(n)}-\dfrac{y_{k-1,s}^{(n+1)}}{E_{k-1,s}^{(n)}}=0.
\end{aligned}\right.
\end{align}
Equations \eqref{eqn:auxiliary_recursion_formula_qdToda_E} therefore
lead to the second equation of \eqref{eqn:auxiliary_recursion_formula_qdToda}.
\end{proof}
\par
The following theorem also gives another matrix representation expressed using $(M_m+M)\times(M_m+M)$ lower and upper block bidiagonal matrices:
\[L^{(n)}\coloneqq\left(\begin{array}{cccc}
I_{M+1} & & & \\
L_1^{(n)} & I_{M+1} & & \\
 & \ddots & \ddots & \\
 & & L_{m-1}^{(n)} & I_{M+1}
\end{array}\right),\quad
R^{(n)}\coloneqq\left(\begin{array}{cccc}
R_1^{(n)} & I_{M+1} & & \\
 & R_2^{(n)} & \ddots & \\
 & & \ddots & I_{M+1} \\
 & & & R_m^{(n)}
\end{array}\right),\]
where $L_k^{(n)}\coloneqq {\rm diag}(E_{k,0}^{(n)},E_{k,1}^{(n)},\dots,E_{k,M}^{(n)})$
and $R_k^{(n)}\coloneqq {\rm diag}(Q_{k,0}^{(n)},Q_{k,1}^{(n)},\dots,$ $Q_{k,M}^{(n)})$.
\begin{theorem}
\label{thm:LR_transformation}
Under discrete-time evolution from $n$ to $n+1$ in \eqref{eqn:auxiliary_recursion_formula_qdToda}, it holds that
\begin{align}
\label{eqn:LU_extended_q_dToda}
L^{(n+1)}R^{(n+1)}-\dfrac{1}{(1-q)t^{(n+1)}}I_{M_m+M}=R^{(n)}L^{(n)}-\dfrac{1}{(1-q)t^{(n)}}I_{M_m+M}.
\end{align}
\end{theorem}
\begin{proof}
All blocks in $L^{(n)}$ and $R^{(n)}$ are diagonal matrices of size $(M+1)\times(M+1)$,
and their sums and products are also diagonal matrices.
The $k$th diagonal block and lower off-diagonal block in $L^{(n+1)}R^{(n+1)}-\{1/[(1-q)t^{(n+1)}]\}I_{M_m+M}$
are, respectively, $R_k^{(n+1)}+L_{k-1}^{(n+1)}-\{1/[(1-q)t^{(n+1)}]\}I_m$ and $L_k^{(n+1)}R_k^{(n+1)}$.
Moreover, the $(s+1,s+1)$ entries of the blocks $R_k^{(n+1)}+L_{k-1}^{(n+1)}-\{1/[(1-q)t^{(n+1)}]\}I_m$
and $L_k^{(n+1)}R_k^{(n+1)}$ are, respectively, $Q_{k,s}^{(n+1)}+E_{k-1,s}^{(n+1)}-1/[(1-q)t^{(n+1)}]$
and $Q_{k,s}^{(n+1)}E_{k,s}^{(n+1)}$, which are respectively equal to
$Q_{k,s}^{(n)}+E_{k,s}^{(n)}-1/[(1-q)t^{(n)}]$ and $Q_{k+1,s}^{(n)}E_{k,s}^{(n)}$.
Namely, the $(s+1,s+1)$ entries of the blocks $R_k^{(n+1)}+L_{k-1}^{(n+1)}-\{1/[(1-q)t^{(n+1)}]\}I_m$
and $L_k^{(n+1)}R_k^{(n+1)}$ coincide with the $(s+1,s+1)$ entries of the blocks
$R_k^{(n)}+L_k^{(n)}-\{1/[(1-q)t^{(n)}]\}I_m$ and $R_{k+1}^{(n)}L_k^{(n)}$,
which are, respectively, the $k$th diagonal block and lower off-diagonal block
in $R^{(n)} L^{(n)} - \{ 1/[(1-q)t^{(n)}] \}I_{M_m+M}$.
Taking into account that the upper off-diagonal blocks of $L^{(n+1)}R^{(n+1)}-\{1/[(1-q)t^{(n+1)}]\}I_{M_m+M}$ and $R^{(n)}L^{(n)}-\{1/[(1-q)t^{(n)}]\}I_{M_m+M}$ are all $I_{M+1}$,
we therefore have \eqref{eqn:LU_extended_q_dToda}.
\end{proof}
\noindent
Obviously, discrete-time evolution from $n$ to $n+1$
in the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda}
completes that in \eqref{eqn:LU_extended_q_dToda}.
Theorem \ref{thm:LR_transformation} thus suggests that
the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda} implicitly generates
similarity transformations, or more precisely, shifted $LR$ transformations of block tridiagonal matrices.
\section{Extended $q$-discrete Lotka--Volterra system and similarity transformations}
\label{sec:4}
In this section, along the same lines as the case of the Toda equation \eqref{eqn:Toda},
we extend $q$-analogues of the LV system \eqref{eqn:LV}, and clarify the relationship to matrix similarity transformations.
\par
Area et al. \cite{Area_2018} presented a $q$-analogue of the LV ($q$-LV) system \eqref{eqn:LV}:
\begin{align}
\label{eqn:qLV}
\left\{\begin{aligned}
& D_qu_k(t)=w_k(qt)\left(u_{k+1}(qt)-u_{k-1}(t)\right),\quad k=1,2,\dots,2m-1,\\
& u_0(t)\coloneqq 0,\quad u_{2m}(t)\coloneqq 0,
\end{aligned} \right.
\end{align}
where the auxiliary variables $w_k(qt)$ satisfy:
\begin{align}
\label{eqn:qLV_w}
\left\{\begin{aligned}
& w_1(qt)=\dfrac{u_1(qt)}{1+(1-q)tu_1(qt)},\\
& w_k(qt)=\dfrac{u_k(qt)}{u_{k-1}(t)}w_{k-1}(qt),\quad k=2,3,\dots,2m-1.
\end{aligned}\right.
\end{align}
Of course, the $q$-LV system \eqref{eqn:qLV} with $q\to 1$ immediately leads to the LV system \eqref{eqn:LV}.
An extension of the $q$-LV system \eqref{eqn:qLV} is given involving arbitrary integer $M$ as:
\begin{align}
\label{eqn:extended_qLV}
\left\{\begin{aligned}
& D_qu_k(t)=w_k(qt)\left(u_{k+M+1}(qt)-u_{k-M-1}(t)\right),\quad k=1,2,\dots,M_{2m}-1,\\
& u_{-M}(t)\coloneqq 0,\quad u_{-M+1}(t)\coloneqq 0,\quad\ldots,\quad u_{0}(t)\coloneqq 0,\\
&u_{M_{2m}}(t)\coloneqq 0,\quad u_{M_{2m}+1}(t)\coloneqq 0,\quad\dots,\quad u_{M_{2m}+M}(t)\coloneqq 0,
\end{aligned}\right.
\end{align}
where $w_k(qt)$ satisfy:
\begin{align}
\label{eqn:extended_qLV_w}
\left\{\begin{aligned}
& w_k(qt)=\dfrac{u_k(qt)}{1+(1-q)tu_k(qt)},\quad k=1,2,\dots,M+1,\\
& w_k(qt)=\dfrac{u_k(qt)}{u_{k-M-1}(t)}w_{k-M-1}(qt),\quad k=M+2,M+3,\dots,M_{2m}-1.
\end{aligned}\right.
\end{align}
Obviously, \eqref{eqn:extended_qLV} with $M=0$ is the $q$-LV system \eqref{eqn:qLV}.
We can regard that the extended $q$-LV system \eqref{eqn:extended_qLV} describes predator-prey interactions, assuming that the $k$th species preys on the $(k+M+1)$st species and is food for the $(k-M-1)$st species. In the simple LV case, the $k$th species preys on the $(k+1)$st species and is food for the $(k-1)$st species.
\par
Given $(M_{2m}+M)\times(M_{2m}+M)$ matrices involving the extended $q$-LV variables $u_{k}(t)$ and $w_{k}(t)$:
\[\begin{aligned}
& U(t)\coloneqq\left(\begin{array}{ccccc}
0 & \overbrace{0\ \cdots\ 0}^M & 1 & & \\
\begin{array}{c} 0 \\ \vdots \\ 0 \end{array}
 & \ddots & \ddots & \ddots & \\
u_{M_1}(t) & \ddots & \ddots & \ddots & 1 \\
 & \ddots & \ddots & \ddots &
 \begin{array}{c} 0 \\ \vdots \\ 0 \end{array} \\
 & & u_{M_{2m}-1}(t) & \underbrace{0\ \cdots\ 0}_{M} & 0
\end{array}\right),\\
& W(t)\coloneqq\left(\begin{array}{ccccc}
0 & & & & \\
\begin{array}{c} 0 \\ \vdots \\ 0 \end{array} & \ddots & & & \\
w_{M_1}(t)u_{M_2}(t) & \ddots & \ddots & & \\
 & \ddots & \ddots & \ddots & \\
 & & w_{M_{2m}-M-2}(t)u_{M_{2m}-1}(t) & \underbrace{0\ \cdots\ 0}_{2M+1} & 0
\end{array}\right),
\end{aligned}\]
we derive a Lax representation for the extended $q$-LV system \eqref{eqn:extended_qLV}.
\begin{proposition}
\label{prop:Lax_extended_qLV}
For $U(t)$ and $W(t)$, it holds that
\begin{align}
\label{eqn:Lax_extended_qLV}
D_qU(t)=U(qt)W(qt)-W(qt)U(t).
\end{align}
\end{proposition}
\begin{proof}
The $(k+M+1,k)$ and $(k+2M+2,k)$ entries of $U(qt)W(qt)$ are, respectively, $w_{k}(qt)u_{k+M+1}(qt)$ and $w_{k}(qt) u_{k+M+1}(qt) u_{k+2M+2}(qt)$.
Those of $W(qt) U(t)$ are, respectively, $w_{k-M-1}(qt)u_{k}(qt)$ and $w_{k+M+1}(qt)u_{k+2M+2}(qt) u_{k}(t)$.
Using the third equation in \eqref{eqn:extended_qLV}, we see that the $(k+M+1,k)$ and $(k+2M+2,k)$ entries of $U(qt)W(qt)-W(qt)U(t)$ are $w_k(qt)(u_{k+M+1}(qt)-u_{k-M-1}(t))$ and $0$, respectively.
The other entries of $U(qt)W(qt)$ and $W(qt)U(t)$ are all $0$.
Thus we have \eqref{eqn:Lax_extended_qLV}.
\end{proof}
Moreover, by preparing $2m\times2m$ matrices involving the extended $q$-LV variable $u_{k,s}(t)\coloneqq u_{M_k+s}(t)$ and $w_{k,s}(t)\coloneqq w_{M_k+s}(t)$:
\[\begin{aligned}
& {\cal U}_s(t)\coloneqq\left(\begin{array}{cccc}
0 & 1 & & \\
u_{1,s}(t) & 0 & \ddots & \\
 & \ddots & \ddots & 1 \\
 & & u_{2m-1,s}(t) & 0
\end{array}\right), \\
& {\cal W}_s(qt)\coloneqq\left( \begin{array}{ccccc}
0 & & & & \\
0 & \ddots & & & \\
w_{1,s}(t) u_{2,s}(t) & \ddots & \ddots & & \\
 & \ddots & \ddots & \ddots & \\
 & & w_{2m-2,s}(t) u_{2m-1,s}(t) & 0 & 0
\end{array}\right),
\end{aligned}\]
we obtain another Lax representation for the extended $q$-LV system \eqref{eqn:extended_qLV}.
\begin{proposition}
\label{prop:Lax_extended_qLV_1}
For the block diagonal matrices ${\cal U}(t)\coloneqq {\rm diag}({\cal U}_0(t),{\cal U}_1(t),\dots,$ ${\cal U}_m(t))$ and ${\cal W}(t)\coloneqq {\rm diag}({\cal W}_{0}(t),{\cal W}_1(t),\dots,{\cal W}_m(t))$, it holds that
\begin{align}
\label{eqn:Lax_extended_qLV_1_all}
D_q{\cal U}(t)={\cal U}(qt){\cal W}(qt)-{\cal W}(qt){\cal U}(t).
\end{align}
\end{proposition}
\begin{proof}
The proof is similar to that of Proposition \ref{prop:Lax_extended_qToda_1}.
Taking into account that$D_q{\cal U}(t)$ and ${\cal U}(qt){\cal W}(qt)-{\cal W}(qt){\cal U}(t)$ are both block diagonal matrices and focusing on the $s$th diagonal blocks of the block diagonal matrices $D_q{\cal U}(t)$ and ${\cal U}(qt){\cal W}(qt)-{\cal W}(qt){\cal U}(t)$, we find directly that:
\begin{align}
\label{eqn:Lax_extended_qLV_1_sub}
D_q{\cal U}_s(t)={\cal U}_s(qt){\cal W}_s(qt)-{\cal W}_s(qt){\cal U}_s(t),\quad s=0,1,\dots,M.
\end{align}
Noting that ${\cal U}_s(t)$ and ${\cal W}_s(t)$ have the same form as $U(t)$ and $W(t)$
in Proposition \ref{prop:Lax_extended_qLV}, respectively,
we obtain equalities of nonzero entries in \eqref{eqn:Lax_extended_qLV_1_sub}:
\[D_qu_{k,s}(t)=w_{k,s}(qt)\left(u_{k+1,s}(qt)-u_{k-1,s}(t)\right),\]
which is equivalent to the extended $q$-LV system \eqref{eqn:extended_qLV} under the replacements
$u_{k,s}(t)=u_{M_k+s}(t)$ and $w_{k,s}(t)=w_{M_k+s}(t)$.
Thus, we have \eqref{eqn:Lax_extended_qLV_1_all}.
\end{proof}
Comparing Propositions \ref{prop:Lax_extended_qToda} and \ref{prop:Lax_extended_qToda_1}
with Propositions \ref{prop:Lax_extended_qLV} and \ref{prop:Lax_extended_qLV_1},
we also find that the extended $q$-Toda equation \eqref{eqn:extended_qToda} and the $q$-LV system \eqref{eqn:extended_qLV}
both generate similarity transformations of $(M+1)$-tridiagonal matrices and block diagonal matrices.
This implies that the extended $q$-LV system \eqref{eqn:extended_qLV} has a close relationship
to the extended $q$-Toda equation \eqref{eqn:extended_qToda}.
The following proposition actually gives the B\"acklund transformation between the extended $q$-Toda equation \eqref{eqn:extended_qToda}
and the extended $q$-LV system \eqref{eqn:extended_qLV}.
\begin{proposition}
\label{prop:Backlund}
The extended $q$-Toda variables $x_{k,s}(t)$, $y_{k,s}(t)$, and $g_{k,s}(t)$
and the extended $q$-LV variables $u_{2k-1,s}(t)$, $u_{2k,s}(t)$, and $w_{2k-1,s}(qt)$ satisfy:
\begin{align}
\label{eqn:Backlund}
\left\{\begin{aligned}
& x_{k,s}(t)=u_{2k-1,s}(t)+u_{2k-2,s}(t),\quad k=1,2,\dots,m,\quad s=0,1,\dots,M,\\
& y_{k,s}(t)=u_{2k,s}(t)u_{2k-1,s}(t),\quad k=1,2,\dots,m, \quad s=0,1,\dots,M,
\end{aligned}\right.
\end{align}
and:
\begin{align}
\label{eqn:Backlund_auxiliary}
g_{k,s}(qt)=w_{2k-1,s}(qt)u_{2k,s}(qt),\quad k=1,2,\dots,m,\quad s=0,1,\dots,M.
\end{align}
\end{proposition}
\begin{proof}
Using \eqref{eqn:extended_qLV_w} and \eqref{eqn:Backlund}, we can easily rewrite the variables $g_{k,s}(qt)$
in \eqref{eqn:extended_qToda_g} as \eqref{eqn:Backlund_auxiliary}.
Considering \eqref{eqn:Backlund} and \eqref{eqn:Backlund_auxiliary}
in the extended $q$-Toda equation \eqref{eqn:extended_qToda}, we obtain:
\[\left\{\begin{aligned}
& D_qu_{2k-1,s}(t)+D_qu_{2k-2,s}(t)=w_{2k-1,s}(qt)u_{2k,s}(qt)-w_{2k-3,s}(qt)u_{2k-2,s}(qt),\\
& \qquad k=1,2,\dots,m, \quad s=0,1,\dots,M,\\
& (D_qu_{2k,s}(t))u_{2k-1,s}(t)+u_{2k,s}(qt)D_qu_{2k-1,s}(t)\\
& \quad =w_{2k-1,s}(qt)u_{2k,s}(qt)(u_{2k+1,s}(qt)+u_{2k,s}(qt)-u_{2k-1,s}(t)-u_{2k-2,s}(t)),\\
& \qquad k=1,2,\dots,m-1,\quad s=0,1,\dots,M.
\end{aligned}\right.\]
Recalling here that $w_{k+1,s}(qt)u_{k,s}(t)=w_{k,s}(qt)u_{k+1,s}(qt)$,
we thus derive the extended $q$-LV system \eqref{eqn:extended_qLV}.
\end{proof}
\begin{figure}[tb]
\centering
\includegraphics[width=0.7\textwidth]{q-dLV_diagram}
\caption{Discrete-time evolution in the extended $q$-discrete Lotka-Volterra system \eqref{eqn:extended_q_dLV}. \label{fig:diagram_time_evolution_qdLV}}
\end{figure}
\par
Now we consider a time-discretization of the extended $q$-LV system \eqref{eqn:extended_qLV}.
Again, by considering a discrete-time sequence $t^{(0)},t^{(1)},\dots$ such that $t^{(n)}\coloneqq q^{-1} t^{(n-1)}$\yadd{,} where $t^{(0)}>0$,
we can rewrite the extended $q$-LV system \eqref{eqn:extended_qLV} with $t=t^{(n)}$ as:
\begin{align}
\label{eqn:extended_q_dLV}
\left\{\begin{aligned}
& \dfrac{u_{k,s}^{(n+1)}-u_{k,s}^{(n)}}{(1-q)t^{(n)}}=w_{k,s}^{(n)}\left(u_{k+1,s}^{(n)}-u_{k-1,s}^{(n+1)}\right),\quad k=1,2,\dots,2m-1,\\
& u_{0,s}^{(n)}\coloneqq 0,\quad u_{2m,s}^{(n)}\coloneqq 0,
\end{aligned}\right.
\end{align}
where $u_{k,s}^{(n+1)}\coloneqq u_{k,s}(t^{(n)})$ and $w_{k,s}^{(n+1)}\coloneqq w_{k,s}(t^{(n)})$ and $w_{k,s}^{(n)}$ satisfy:
\begin{align}
\label{eqn:extended_q_dLV_w}
\left\{\begin{aligned}
& w_{1,s}^{(n)}=\dfrac{u_{1,s}^{(n)}}{1+(1-q)t^{(n)}u_{1,s}^{(n)}},\quad
w_{k,s}^{(n)}=\dfrac{u_{k,s}^{(n)}}{u_{k-1,s}^{(n+1)}}w_{k-1,s}^{(n)},\quad k=2,3,\dots,2m-1,\\
& w_{0,s}^{(n)}\coloneqq 0,\quad w_{2m,s}^{(n)}\coloneqq 0.
\end{aligned}\right.
\end{align}
We can easily check that the sequences $\{u_{k,s}^{(n)}\}_{n=0,1,\dots}$ and $\{w_{k,s}^{(n)}\}_{n=0,1,\dots}$ are uniquely determined under the initial settings of $u_{k,s}^{(n)}$ and $w_{k,s}^{(n)}$.
Figure \ref{fig:diagram_time_evolution_qdLV} shows the discrete-time evolution from $n$ to $n+1$ in the $q$-discrete LV system \eqref{eqn:extended_q_dLV}.
Using $U^{(n)}\coloneqq U(t^{(n-1)})$ and $W^{(n)}\coloneqq W(t^{(n-1)})$, we obtain a matrix representation for the extended $q$-discrete LV system \eqref{eqn:extended_q_dLV}.
\begin{theorem}
\label{thm:similar_transformation_dqLV}
For the $(M+1)$-tridiagonal matrix $U^{(n)}$ and $(2M+2)$-subdiagonal matrix $W^{(n)}$, it holds that
\[(I_{M_{2m}+M}+(1-q)t^{(n)}W^{(n)})U^{(n+1)}=U^{(n)}(I_{M_{2m}+M}+(1-q)t^{(n)}W^{(n)}).\]
\end{theorem}
\noindent
Since $\det(I_{M_{2m}+M}+(1-q)t^{(n)}W^{(n)})=1$, it follows from Theorem \ref{thm:similar_transformation_dqLV} that
$U^{(n+1)}=(I_{M_{2m}+M}+(1-q)t^{(n)}W^{(n)})^{-1}U^{(n)}(I_{M_{2m}+M}+(1-q)t^{(n)}W^{(n)})$.
Theorem \ref{thm:similar_transformation_dqLV} thus suggests that the extended $q$-discrete LV system \eqref{eqn:extended_q_dLV}
generates a similarity transformation from $U^{(n)}$ to $U^{(n+1)}$.
The replacements $x_{k,s}^{(n)}=x_{k,s}(t^{(n-1)})$, $y_{k,s}^{(n)}=y_{k,s}(t^{(n-1)})$, $g_{k,s}^{(n)}=g_{k,s}(t^{(n-1)})$,
$u_{2k-1,s}^{(n)}=u_{2k-1,s}(t^{(n-1)})$, $u_{2k,s}^{(n)}=u_{2k,s}(t^{(n-1)})$, and $w_{2k-1,s}^{(n)}=w_{2k-1,s}(t^{(n-1)})$
in Proposition \ref{prop:Backlund} yields the B\"acklund transformation
between the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda}
and the extended $q$-discrete LV system \eqref{eqn:extended_q_dLV}.
\begin{proposition}
\label{prop:discrete_Backlund}
The extended $q$-discrete Toda variables $x_{k,s}^{(n)}$, $y_{k,s}^{(n)}$, and $g_{k,s}^{(n)}$
and the extended $q$-discrete LV variables $u_{2k-1,s}^{(n)}$, $u_{2k,s}^{(n)}$, and $w_{2k-1,s}^{(n)}$ satisfy:
\begin{align}
\label{eqn:discrete_Backlund}
\left\{\begin{aligned}
& x_{k,s}^{(n)}=u_{2k-1,s}^{(n)}+u_{2k-2,s}^{(n)},\quad k=1,2,\dots,m,\quad s=0,1,\dots,M, \\
& y_{k,s}^{(n)}=u_{2k,s}^{(n)}u_{2k-1,s}^{(n)},\quad k=1,2,\dots,m,\quad s=0,1,\dots,M,
\end{aligned}\right.
\end{align}
and:
\begin{align}
\label{eqn:discrete_Backlund_auxiliary}
g_{k,s}^{(n)}=w_{2k-1,s}^{(n)}u_{2k,s}^{(n)},\quad k=1,2,\dots,m,\quad s=0,1,\dots,M.
\end{align}
\end{proposition}
\noindent
With the help of Proposition \ref{prop:discrete_Backlund},
we can express $Q_{k,s}^{(n)}$ and $E_{k,s}^{(n)}$ appearing in Proposition \ref{prop:auxiliary_recursion_formula_qdToda}
using the extended $q$-discrete LV variables $u_{2k-1,s}^{(n)}$, $u_{2k,s}^{(n)}$, and $w_{2k-1,s}^{(n)}$ as:
\begin{align}
\label{eqn:auxiliary_variable_extended_qLV}
\left\{\begin{aligned}
& Q_{k,s}^{(n)}\coloneqq u_{2k-1,s}^{(n)}+u_{2k-2,s}^{(n)}+\dfrac{1}{(1-q)t^{(n)}}-(1-q)t^{(n)}w_{2k-3,s}^{(n)}u_{2k-2,s}^{(n)},\\
& \qquad k = 1,2,\dots,m,\quad s=0,1,\dots,M,\\
& E_{k,s}^{(n)}\coloneqq (1-q)t^{(n)}w_{2k-1,s}^{(n)}u_{2k,s}^{(n)},\quad k=1,2,\dots,m,\quad s=0,1,\dots,M.
\end{aligned}\right.
\end{align}
This implies that the extended $q$-discrete LV system \eqref{eqn:extended_q_dLV} implicitly gives a sequence
of shifted $LR$ transformations which are the same as shown in Theorem \ref{thm:LR_transformation}.
\section{Asymptotic convergence}
\label{sec:5}
In this section, we describe the properties of the Hankel determinants related to an infinite moment sequence,
and then clarify the determinantal solution expressed using the Hankel determinants
to the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda}.
Moreover, by considering asymptotic expansions of the Hankel determinants as $n\to\infty$,
we clarify asymptotic convergence as $n\to\infty$ in the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda}.
\par
We consider monic polynomials with respect to $z$ having distinct roots $\lambda_{1,s},\lambda_{2,s},$ $\dots,\lambda_{m,s}$:
\begin{align}
\label{eqn:polynomial}
p_s(z)\coloneqq (z-\lambda_{1,s})(z-\lambda_{2,s})\cdots (z-\lambda_{m,s}),\quad s=0,1,\dots,M.
\end{align}
Equation \eqref{eqn:polynomial} can be expanded as:
\begin{equation}
\label{eqn:polynomial_extended}
p_s(z)=z^m+a_{1,s}z^{m-1}+\cdots +a_{m-1,s}z+a_{m,s},\quad s=0,1,\dots,M,
\end{equation}
where $a_{1,s},a_{2,s},\dots,a_{m,s}$ are given in terms of $\lambda_{1,s},\lambda_{2,s},\dots,\lambda_{m,s}$.
Now, we introduce an infinite sequence $\{ f_{k,s}^{(n)} \}_{n=0}^{\infty}$
associated with the constants $\lambda_{1,s},\lambda_{2,s},\dots,\lambda_{m,s}$.
The sequence $\{f_{k,s}^{(n)}\}_{n=0}^{\infty}$ is sometimes called a moment sequence.
Let us assume that $\{ f_{k,s}^{(n)} \}_{n=0}^{\infty}$ satisfies the linear equations involving $a_{1,s},a_{2,s},\dots,a_{m,s}$:
\begin{align}
\label{eqn:moment_linear_equation}
f_{k+m,s}^{(n)}+\sum_{i=1}^ma_{i,s}f_{k+m-i,s}^{(n)}=0,\quad k=0,1,\dots,m,\quad s=0,1,\dots,M,\quad n=0,1,\dots,
\end{align}
and discrete-time evolutions from $n$ to $n+1$:
\begin{align}
\label{eqn:moment_time_evolution}
f_{k,s}^{(n+1)}=f_{k+1,s}^{(n)}-\mu^{(n)}f_{k,s}^{(n)},\quad k=0,1,\dots,m-1,\quad s=0,1,\dots,M,
\end{align}
where $f_{2m+1,s}^{(n)}\coloneqq 0$ and $\mu^{(n)}$ are arbitrary constants.
It is obvious that the moment sequence $\{f_{k,s}^{(n)}\}_{n=0}^{\infty}$ is uniquely determined
from the values of $f_{0,s}^{(0)},f_{1,s}^{(0)},\dots,f_{m-1,s}^{(0)}$ for each $s$.
It is worth noting that, for each $s$, the moment sequence $\{f_{k,s}^{(n)}\}_{n=0}^{\infty}$ is equivalent to
the moment sequence $\{f_{k}^{(n)}\}_{k,n=0}^{\infty}$ given by \cite{Shinjo_2018}.
Thus, by extending their discussion, we can easily derive a proposition for the moments $f_{k,s}^{(n)}$.
\begin{proposition}(cf.~\cite{Shinjo_2018})
For each $s$, the moments $f_{k,s}^{(n)}$ can be expressed as:
\begin{align}
\label{eqn:moment_expressed}
f_{k,s}^{(n)}=\sum_{\ell=1}^mc_{\ell,s}\lambda_{\ell,s}^k\rho_{\ell,s}^{(n)}\quad k=0,1,\dots,2m,\quad n=0,1,\dots,
\end{align}
where:
\begin{align}
\label{eqn:rho}
\rho_{\ell,s}^{(0)}\coloneqq 1,\quad\rho_{\ell,s}^{(n)}\coloneqq\prod_{j=0}^{n-1}(\lambda_{\ell,s}-\mu^{(j)}),\quad\ell =1,2,\dots,m,
\end{align}
and constants $c_{1,s},c_{2,s},\dots,c_{m,s}$ satisfy:
\begin{align}
\label{eqn:c_condition}
\left(\begin{array}{c}
c_{1,s} \\ c_{2,s} \\ \vdots \\ c_{m,s}
\end{array}\right)
=\left(\begin{array}{cccc}
1 & 1 & \cdots & 1 \\
\lambda_{1,s} & \lambda_{2,s} & \cdots & \lambda_{m,s}\\
\vdots & \vdots & & \vdots\\
\lambda_{1,s}^{m-1} & \lambda_{2,s}^{m-1} & \cdots & \lambda_{m,s}^{m-1}
\end{array}\right)^{-1}
\left(\begin{array}{c}
f_{0,s}^{(0)} \\ f_{1,s}^{(0)} \\ \vdots \\ f_{m-1,s}^{(0)}
\end{array}\right).
\end{align}
\end{proposition}
Similarly to Ref.\cite{Shinjo_2018}, by considering Hankel determinants of degree $k$ given using the moment $f_{k,s}^{(n)}$ as:
\begin{align}
\label{eqn:determinant}
\begin{aligned}
& H_{-1,s}^{(n)}\coloneqq 0,\quad H_{0,s}^{(n)}\coloneqq 1,\\
& H_{k,s}^{(n)}\coloneqq\left\vert \begin{array}{cccc}
 f_{0,s}^{(n)} & f_{1,s}^{(n)} & \cdots & f_{k-1,s}^{(n)}\\
 f_{1,s}^{(n)} & f_{2,s}^{(n)} & \cdots & f_{k,s}^{(n)}\\
 \vdots & \vdots & \ddots & \vdots\\
 f_{k-1,s}^{(n)} & f_{k,s}^{(n)} & \cdots & f_{2k-2,s}^{(n)}
\end{array}\right\vert,\quad k=1,2,\dots,m,\quad n=0,1,\dots,
\end{aligned}
\end{align}
we also obtain two propositions for the Hankel determinants $H_{k,s}^{(n)}$.
\begin{proposition}(cf.~\cite{Shinjo_2018})
For each $s$, the Hankel determinants $H_{k,s}^{(n)}$ satisfy:
\begin{align}
\label{prop:determinant_extended}
H_{k,s}^{(n)}=
\left\{
\begin{aligned}
&\sum_{1\leq i_1<i_2<\cdots <i_k\leq m} { \cal C}_{i_1,i_2,\dots,i_k;s}\rho_{i_1,s}^{(n)}\rho_{i_2,s}^{(n)}\cdots\rho_{i_k,s}^{(n)},\quad k=1,2,\dots,m,\\
&0,\quad k=m+1,\\
\end{aligned}
\right.
\end{align}
where:
\begin{align}
\label{eqn:determinant_extended_C}
{\cal C}_{i_1,i_2,\dots,i_k;s}\coloneqq\left\vert\begin{array}{cccc}
1 & 1 & \cdots & 1 \\
\lambda_{i_1,s} & \lambda_{i_2,s} & \cdots & \lambda_{i_k,s} \\
\vdots & \vdots & \ddots & \vdots \\
\lambda_{i_1,s}^{k-1} & \lambda_{i_1,s}^{k-1} & \cdots & \lambda_{i_k,s}^{k-1}
\end{array}\right\vert\left\vert\begin{array}{cccc}
c_{i_1,s} & c_{i_1,s} \lambda_{i_1,s} & \cdots & c_{i_1,s} \lambda_{i_1,s}^{k-1} \\
c_{i_2,s} & c_{i_2,s} \lambda_{i_2,s} & \cdots & c_{i_2,s} \lambda_{i_2,s}^{k-1} \\
\vdots & \vdots & \ddots & \vdots \\
c_{i_k,s} & c_{i_k,s} \lambda_{i_k,s} & \cdots & c_{i_k,s} \lambda_{i_k,s}^{k-1}
\end{array}\right\vert.
\end{align}
\end{proposition}
Let us assume that $H_{k,s}^{(n)}\neq 0$ for $k=1,2,\dots,m$ and $s=0,1,\dots,M$.
Let us define the $k$th degree polynomials with respect to $z$ as:
\begin{align}
\label{eqn:orthogonal_polynomial}
{\cal H}_{k,s}^{(n)}(z)=\dfrac{H_{k,s}^{(n)}(z)}{H_{k,s}^{(n)}},\quad k=1,2,\dots,m,
\end{align}
where
\begin{align}
\label{eqn:determinant_func}
\begin{aligned}
& H_{-1,s}^{(n)}(z)\coloneqq 0,\quad H_{0,s}^{(n)}(z)\coloneqq 1,\\
& H_{k,s}^{(n)}(z)\coloneqq\left\vert\begin{array}{cccc}
f_{0,s}^{(n)} & f_{1,s}^{(n)} & \cdots & f_{k,s}^{(n)} \\
f_{1,s}^{(n)} & f_{2,s}^{(n)} & \cdots & f_{k+1,s}^{(n)} \\
\vdots & \vdots & \ddots & \vdots \\
f_{k-1,s}^{(n)} & f_{k,s}^{(n)} & \cdots & f_{2k-1,s}^{(n)} \\
1 & z & \cdots & z^k
\end{array}\right\vert,\quad k=1,2,\dots,m,\quad n=0,1,\dots.
\end{aligned}
\end{align}
\noindent
With the help of \yadd{\cite{Shinjo_2018}}, we immediately derive two lemmas concerning the polynomials ${\cal H}_{k,s}^{(n)}(z)$.
\begin{lemma}(cf.~\cite{Shinjo_2018})
Let us assume that $H_{k,s}^{(n)}\neq 0$ and $\mu^{(n)}=-1/ [(1-q)t^{(n)}]$.
For each $s$, the Hadamard polynomials ${\cal H}_{k,s}^{(n)}$ satisfy:
\begin{align}
\left\{\begin{aligned}
& (z-\mu^{(n)}){\cal H}_{k-1,s}^{(n+1)}(z)={\cal H}_{k,s}^{(n)}(z)+Q_{k,s}^{(n)}{\cal H}_{k-1,s}^{(n)}(z),\quad k=1,2,\dots,m,\\
& {\cal H}_{k,s}^{(n)}(z)={\cal H}_{k,s}^{(n+1)}(z)+E_{k,s}^{(n)}{\cal H}_{k-1,s}^{(n+1)}(z),\quad k=0,1,\dots,m,
\end{aligned}\right.
\end{align}
where
\begin{align}
\left\{\begin{aligned}
\label{eqn:auxiliary_variable_determinantal_solution}
& Q_{k,s}^{(n)}=\dfrac{H_{k,s}^{(n+1)}H_{k-1,s}^{(n)}}{H_{k,s}^{(n)}H_{k-1,s}^{(n+1)}},\quad k=1,2,\dots,m,\quad s=0,1,\dots,M,\\
& E_{k,s}^{(n)}=\dfrac{H_{k+1,s}^{(n)}H_{k-1,s}^{(n+1)}}{H_{k,s}^{(n)}H_{k,s}^{(n+1)}},\quad k=0,1,\dots,m,\quad s=0,1,\dots,M.
\end{aligned}\right.
\end{align}
Moreover, ${\cal H}_{m,s}^{(n)}(z)$ are just $p_s(z)$, which are characteristic polynomials of matrices with eigenvalues $\lambda_{1,s},\lambda_{2,s},\ldots,\lambda_{m,s}$.
\end{lemma}
Combining \eqref{eqn:auxiliary_variable_determinantal_solution} with \eqref{eqn:auxiliary_recursion_formula_qdToda_Q} and \eqref{eqn:auxiliary_recursion_formula_qdToda_E}, we thus obtain the determinantal solution to the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda}.
\begin{theorem}
\label{thm:solutions_extended_qdToda}
The extended $q$-discrete Toda variables $x_{k,s}^{(n)}$ and $y_{k,s}^{(n)}$ can be expressed using the Hankel determinants $H_{k,s}^{(n)}$ as
\begin{align}
\label{eqn:qdToda_solution}
\begin{aligned}
& x_{k,s}^{(n)}=\dfrac{H_{k,s}^{(n+1)}H_{k-1,s}^{(n)}}{H_{k,s}^{(n)}H_{k-1,s}^{(n+1)}}
+\dfrac{H_{k,s}^{(n)}H_{k-2,s}^{(n+1)}}{H_{k-1,s}^{(n)}H_{k-1,s}^{(n+1)}}-\dfrac{1}{(1-q)t^{(n)}},\quad k=1,2,\dots,m,\\
& y_{k,s}^{(n)}=\dfrac{H_{k-1,s}^{(n)}H_{k+1,s}^{(n)}}{(H_{k,s}^{(n)})^2},\quad k=0,1,\dots,m.
\end{aligned}
\end{align}
\end{theorem}
\noindent
According to \yadd{\cite{Shinjo_2018}}, asymptotic expansions as $n\to\infty$ of the Hankel determinants are given as the following lemma.
\begin{lemma}(cf.~\cite{Shinjo_2018})
\label{lemma:asymptotic_behavior}
Let us assume that $\vert\lambda_{1,s}-\mu^{(n)}\vert >\vert\lambda_{2,s}-\mu^{(n)}\vert >\cdots >\vert\lambda_{m,s}-\mu^{(n)}\vert$,
and $\varrho_{k,s}$ are constants that satisfy $0\leq\varrho_{k,s}<1$ and
$\varrho_{k,s}> \vert\lambda_{k+1,s}-\mu^{(n)}\vert /\vert\lambda_{k,s}-\mu^{(n)}\vert$.
For sufficiently large $n$, it holds that
\[H_{k,s}^{(n)}={\cal C}_{1,2,\dots,k;s}\rho_{1,s}^{(n)}\rho_{2,s}^{(n)}\cdots\rho_{k,s}^{(n)}\left(1+O(\varrho_{k,s}^{n})\right).\]
\end{lemma}
From Theorem \ref{thm:solutions_extended_qdToda} and Lemma \ref{lemma:asymptotic_behavior},
we therefore have the convergence theorem in the extended $q$-Toda equation \eqref{eqn:extended_qdToda}.
\begin{theorem}
For each $s$, let us assume that the moment sequence $\{f_{k,s}^{(n)}\}_{n}^{\infty}$ satisfies $H_{k,s}^{(n)}\neq0$,
and $\vert\lambda_{1,s}-\mu^{(n)}\vert >\vert\lambda_{2,s}-\mu^{(n)}\vert >\cdots >\vert\lambda_{m,s}-\mu^{(n)}\vert$.
Then, it holds that
\begin{align}
\label{eqn:qdToda_solution_x_converge}
& \lim_{n\to\infty}x_{k,s}^{(n)}=\lambda_{k,s},\quad k=1,2,\dots,m,\\
\label{eqn:qdToda_solution_y_converge}
& \lim_{n\to\infty}y_{k,s}^{(n)}= 0,\quad k=1,2,\dots,m-1.
\end{align}
\end{theorem}
\begin{proof}
Using Lemma \ref{lemma:asymptotic_behavior} and
considering $1/[(1-q)t^{(n)}]\to 0$ as $n\to\infty$ in Theorem \ref{thm:solutions_extended_qdToda}, we obtain
\begin{align*}
x_{k,s}^{(n)} & =\dfrac{\rho_{k,s}^{(n+1)}}{\rho_{k,s}^{(n)}}
\dfrac{\left(1+O(\varrho_{k,s}^{n+1})\right)\left(1+O(\varrho_{k-1,s}^{n})\right)}
{\left(1+O(\varrho_{k,s}^{n})\right)\left(1+O(\varrho_{k-1,s}^{n+1})\right)}\\
& \quad +\dfrac{{\cal C}_{1,2,\dots,k;s}{\cal C}_{1,2,\dots,k-2;s}}{{\cal C}_{1,2,\dots,k-1;s}^2}\dfrac{\rho_{k,s}^{(n)}}{\rho_{k-1,s}^{(n+1)}}
\dfrac{\left(1+O(\varrho_{k,s}^{n})\right)\left(1+O(\varrho_{k-2,s}^{n+1})\right)}
{\left(1+O(\varrho_{k-1,s}^{n})\right)\left(1+O(\varrho_{k-1,s}^{n+1})\right)}, \\
y_{k,s}^{(n)} & =\dfrac{{\cal C}_{1,2,\dots,k+1;s}{\cal C}_{1,2,\dots,k-1;s}}{{\cal C}_{1,2,\dots,k;s}^2}
\dfrac{\rho_{k+1,s}^{(n)}}{\rho_{k,s}^{(n)}}
\dfrac{\left(1+O(\varrho_{k+1,s}^{n})\right)\left(1+O(\varrho_{k-1,s}^{n})\right)}{\left(1+O(\varrho_{k,s}^{n})\right)^2}.
\end{align*}
Thus, we derive \eqref{eqn:qdToda_solution_x_converge} and \eqref{eqn:qdToda_solution_y_converge}.
\end{proof}
Considering $x_{k,s}^{(n)}\to\lambda_{k,s}$ and $y_{k,s}^{(n)}\to0$ as $n\to\infty$ in Proposition \ref{prop:discrete_Backlund}, we also have a convergence theorem in the extended $q$-discrete LV system \eqref{eqn:extended_q_dLV}.
\begin{theorem}
For each $s$, let us assume that the moment sequence $\{f_{k,s}^{(n)}\}_{n}^{\infty}$ satisfies $H_{k,s}^{(n)}\neq0$, and $\vert\lambda_{1,s}-\mu^{(n)}\vert >\vert\lambda_{2,s}-\mu^{(n)}\vert >\cdots >\vert\lambda_{m,s}-\mu^{(n)}\vert$.
Then, it holds that
\begin{align*}
& \lim_{n\to\infty}u_{2k-1,s}^{(n)}=\lambda_{k,s},\quad k=1,2,\dots,m,\\
\label{eqn:qdToda_solution_y_converge}
& \lim_{n\to\infty}u_{2k,s}^{(n)}= 0,\quad k=1,2,\dots,m-1.
\end{align*}
\end{theorem}
\par
In the remainder of this section, we present numerical examples to show the convergence to matrix eigenvalues in the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda}.
We used floating point arithmetic on a computer with a Windows 10 Professional
operating system with an Intel(R) Core(TM) i5-7200U CPU @ 2.50 GHz 2.71 GHz, and employed the software Python 3.8.2.
\begin{figure}[tb]
\centering
\includegraphics[width=0.7\textwidth]{example_fig1}
\caption{Discrete-time $n$ ($x$-axis) versus values of $x_{1,0}^{(n)}$, $x_{1,1}^{(n)}$, $x_{2,0}^{(n)}$, and $x_{2,1}^{(n)}$ ($y$-axis).
Pluses: $\{x_{1,0}^{(n)}\}_{n=0,1,\dots,20}$; crosses: $\{x_{1,1}^{(n)}\}_{n=0,1,\dots,20}$; circles: $\{x_{2,0}^{(n)}\}_{n=0,1,\dots,20}$; and triangles: $\{x_{2,1}^{(n)}\}_{n=0,1,\dots,20}$.\label{fig:convergence_eigenvalue}}
\end{figure}
\par
We first study the computation of eigenvalues of a $4\times4$ $2$-tridiagonal matrix
\[A^{(0)}=\left(\begin{array}{cccc}
7 & & 1 & \\
 & -3 & & 1 \\
6 & & 2 & \\
 & 5 & & 1
\end{array}\right).\]
It is easy to check that the eigenvalues of $A^{(0)}$ are $8,2,1$ and $-4$ because $A^{(0)}$ is a similar matrix to the block diagonal matrix with tridiagonal blocks
\[{\cal A}^{(0)}=\left(\begin{array}{cccc}
7 & 1 & & \\
6 & 2 & & \\
 & & -3 & 1 \\
 & & 5 & 1
\end{array}\right).\]
We set the parameters in the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda} as $m=2$, $M=1$, $t^{(0)}=1$ and $q=1/2$.
The initial values of the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda} are directly
given from entries of $A^{(0)}$ as $x_{1}^{(0)}=x_{1,0}^{(0)}=7$, $x_{2}^{(0)}=x_{1,1}^{(0)}=-3$, $x_{3}^{(0)}=x_{2,0}^{(0)}=2$,
$x_{4}^{(0)}=x_{2,1}^{(0)}=1$, $y_{1}^{(0)}=y_{1,0}^{(0)}=6$, and $y_{2}^{(0)}=y_{1,1}^{(0)}=5$.
This differs from the case of the qd recursion formula, namely, the original discrete Toda equation that
requires the decomposition of the target matrix into a product of lower and upper diagonal matrices.
Figure \ref{fig:convergence_eigenvalue} draws the approach of the extended $q$-discrete Toda variables
$x_{1,0}^{(n)}, x_{2,0}^{(n)}$, $x_{1,1}^{(n)}$ and $x_{2,1}^{(n)}$
to the eigenvalues $8$, $1$, $-4$ and $2$, respectively, as $n$ grows larger.
\begin{figure}[tb]
\centering
\includegraphics[width=0.7\textwidth]{example_fig2}
\caption{Discrete-time $n$ ($x$-axis) versus values of $x_{1,0}^{(n)}$, $x_{1,1}^{(n)}$, $x_{2,0}^{(n)}$, and $x_{2,1}^{(n)}$ ($y$-axis).
Pluses: $\{x_{1,0}^{(n)}\}_{n=0,1,\dots,20}$; crosses: $\{x_{1,1}^{(n)}\}_{n=0,1,\dots,20}$;
circles: $\{x_{2,0}^{(n)}\}_{n=0,1,\dots,20}$; and triangles: $\{x_{2,1}^{(n)}\}_{n=0,1,\dots,20}$.\label{fig:not_convergence_eigenvalue}}
\end{figure}
\par
Next, we prepare the $4\times4$ matrices
\[A^{(0)}=\left(\begin{array}{cccc}
3 & & 1 & \\
 & -5 & & 1 \\
6 & & -2 & \\
 & -5 & & -1
\end{array}\right)
\quad and \quad
{\cal A}^{(0)}=\left(\begin{array}{cccc}
3 & 1 & & \\
6 & -2 & & \\
 & & -5 & 1 \\
 & & -5 & -1
\end{array}\right)\]
with real eigenvalues $4$ and $-3$ and complex eigenvalues $-3+i$ and $-3-i$.
In the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda}, we set the initial values and parameters as
$x_{1}^{(0)}=x_{1,0}^{(0)}=3$, $x_{2}^{(0)}=x_{1,1}^{(0)}=-5$, $x_{3}^{(0)}=x_{2,0}^{(0)}=-2$, $x_{4}^{(0)}=x_{2,1}^{(0)}=-1$,
$y_{1}^{(0)}=y_{1,0}^{(0)}=6$, and $y_{2}^{(0)}=y_{1,1}^{(0)}=-5$, as well as $m=2$, $M=1$, $t^{(0)}=1$ and $q=1/2$.
Figure \ref{fig:not_convergence_eigenvalue} implies that $x_{1,0}^{(n)}$ and $x_{2,0}^{(n)}$ respectively converge to real eigenvalues $4$ and $-3$,but $x_{2,0}^{(n)}$ and $x_{2,1}^{(n)}$ do not converge as $n\to\infty$.
In other words, we numerically verified that the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda}
can be applied to computing only real eigenvalues even in the case where target matrices have complex eigenvalues.
Though complex eigenvalues are computed using the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda},
those of the block that is not diagonalized are almost unchanged from $-3+i$ and $-3-i$.
\begin{figure}[tb]
\centering
\includegraphics[width=0.7\textwidth]{example_fig3}
\caption{Discrete-time $n$ ($x$-axis) versus values of $x_{1,0}^{(n)},x_{1,1}^{(n)},x_{1,2}^{(n)},x_{1,3}^{(n)},x_{2,0}^{(n)},x_{2,1}^{(n)},x_{2,2}^{(n)}$ and $x_{2,3}^{(n)}$ ($y$-axis).
Circle: $\{x_{1,0}^{(n)}\}_{n=0,1,\dots,20}$; squares: $\{x_{1,1}^{(n)}\}_{n=0,1,\dots,20}$;
triangles: $\{x_{1,2}^{(n)}\}_{n=0,1,\dots,20}$; diamonds: $\{x_{1,3}^{(n)}\}_{n=0,1,\dots,20}$;
filled-circle: $\{x_{2,0}^{(n)}\}_{n=0,1,\dots,20}$; filled-squares: $\{x_{2,1}^{(n)}\}_{n=0,1,\dots,20}$;
filled-triangles: $\{x_{2,2}^{(n)}\}_{n=0,1,\dots,20}$; filled-diamonds: $\{x_{2,3}^{(n)}\}_{n=0,1,\dots,20}$.
\label{fig:M_parameter_eigenvalue}}
\end{figure}
\par
The last example concerns the case of $8\times8$ matrices
\[A^{(0)}=\left(\begin{array}{cccccccc}
4 & & & & 1 \\
& 3 & & & & 1 \\
& & 6 & & & & 1 \\
& & & 1 & & & & 1 \\
1& & & & 4 \\
& 1 & & & & 5 \\
& & 1 & & & & 2 \\
& & & 1 & & & & 7
\end{array}\right)
\quad and \quad
{\cal A}^{(0)}=\left(\begin{array}{cccccccc}
4 & 1 & & \\
1 & 4 & & \\
 & & 3 & 1 \\
 & & 1 & 5 \\
 & & & & 6 & 1 \\
 & & & & 1 & 2 \\
 & & & & & & 1 & 1 \\
 & & & & & & 1 & 7 \\
\end{array}\right)\]
with eigenvalues $4\pm1,4\pm\sqrt{2},4\pm\sqrt{5}$ and $4\pm\sqrt{10}$.
We set the initial values and parameters as
$x_{1}^{(0)}=x_{1,0}^{(0)}=4$, $x_{2}^{(0)}=x_{1,1}^{(0)}=3$, $x_{3}^{(0)}=x_{1,2}^{(0)}=6$, $x_{4}^{(0)}=x_{2,1}^{(0)}=1$,
$x_{5}^{(0)}=x_{2,0}^{(0)}=4$, $x_{6}^{(0)}=x_{2,1}^{(0)}=5$, $x_{7}^{(0)}=x_{2,2}^{(0)}=2$, $x_{8}^{(0)}=x_{2,1}^{(0)}=7$,
$y_{1}^{(0)}=y_{1,0}^{(0)}=1$, $y_{2}^{(0)}=y_{1,1}^{(0)}=1$, $y_{3}^{(0)}=y_{1,2}^{(0)}=1$, $y_{4}^{(0)}=y_{1,3}^{(0)}=1$
as well as $m=2$, $M=3$, $t^{(0)}=1$ and $q=1/2$.
in the extended $q$-discrete Toda equation \eqref{eqn:extended_qdToda}.
Figure \ref{fig:M_parameter_eigenvalue} shows that $x_{1,0}^{(n)},x_{1,1}^{(n)},x_{1,2}^{(n)},x_{1,3}^{(n)},x_{2,0}^{(n)},x_{2,1}^{(n)},x_{2,2}^{(n)}$ and $x_{2,3}^{(n)}$, respectively, \yadd{converge} to eigenvalues $5,4+\sqrt{2},4+\sqrt{5},4+\sqrt{10},3,4-\sqrt{2},4-\sqrt{5}$ and $4-\sqrt{10}$ as $n\to\infty$. Comparing Figure \ref{fig:M_parameter_eigenvalue} with Figure \ref{fig:convergence_eigenvalue}, we can see that $M=3$ case has the asymptotic behavior as $n$ grows larger similarly to the $M=1$ case.
\section{Concluding remarks}
\label{sec:6}
In this paper, we first considered an extension of the $q$-Toda equation which is a $q$-analogue of the famous Toda equation, and related it to eigenvalue problems of $(M+1)$-tridiagonal matrices and block diagonal matrices whose blocks are tridiagonal.
We next showed that time-discretization of the extended $q$-Toda equation can generate similarity transformations of the $(M+1)$-tridiagonal matrices.
We also found the relationship between the Toda case and LV case and related the LV case to the same eigenvalue problem as the Toda case.
Finally, we proved convergence to eigenvalues in the extended $q$-discrete Toda equation and presented numerical examples for verifying it.
\par
Our future works \yadd{from the numerical analysis perspective} are to examine the numerical stability of the similarity transformations
generated by the extended $q$-discrete Toda equation
and accelerating the convergence speed by introducing explicit shifts, rather than implicit shifts.
From the viewpoint of the study of integrable systems,
we plan to relate a $q$-analogue of hungry integrable systems to eigenvalue problems.
\section*{Acknowledgements}
The authors thank the reviewer for his/her careful reading and constructive suggestions.
This work was partially supported by the joint project of Kyoto University and Toyota Motor Corporation, titled ``Advanced Mathematical Science for Mobility Society''.
\section*{Disclosure statement}
No potential conflict of interest was reported by the authors.
\section*{ORCID}
R. Watanabe: 0000-0002-5758-9587
\bibliographystyle{tfs}


\section{Introduction}

Long before the discovery of the Higgs boson in 2012, vacuum decay had been used as a theoretical tool to constrain the masses of the Standard Model (SM) particles~\cite{isidori2001metastability, casas1995improved, espinosa1995improved, altarelli1994lower, sher1993precise}. The critical reason for vacuum instability in the SM is that with loop corrections resummed via the renormalization group, the quartic coefficient of the effective potential becomes negative as the renormalization scale approaches about $10^{10}$ \ensuremath{\,\text{GeV}}\xspace \cite{buttazzo2013investigating}. This means that the effective potential must be ultraviolet completed such that it is bounded from below and a second, deeper minimum arises~\cite{sher1989electroweak, krive1976vacuum}. Fortunately, detailed calculation shows that the tunneling time of the desired electroweak symmetry breaking~(EW-breaking) vacuum into this deeper vacuum is longer than the age of the universe~\cite{andreassen2018scale}.

In beyond the SM (BSM) scenarios, like supersymmetric models, extra scalar fields are usually introduced, which can get non-zero vacuum expectation values and destabilise the EW-breaking vacuum. If another lower energy minimum exists, the EW-breaking vacuum could transition into such a vacuum via quantum tunneling or thermal excitation. Since the EW-breaking vacuum is compatible with experimental results, we naturally demand that the transition time is longer than the age of the universe, which constrains the model to describe Nature. In simple models such constraints for stability at zero temperature can be derived analytically
~\cite{ferreira2004stability, barroso2007neutral, barroso2013metastability}. To deal with multi-fields or high temperature effects, several numerical tools were developed ~\cite{camargo2013vevacious, masoumi2017efficient, sato2021simplebounce}.
Both analytic and numerical calculations assume that the set of considered parameters ensure the realization of the EW-breaking vacuum during the evolution of the universe.

This assumption, however, is not always correct from the perspective of finite temperature field theory. Previously, Cline \emph{et al.}~\cite{cline1999electroweak} investigated whether the transition process from a color-breaking vacuum to the EW-breaking vacuum can happen and found that the nucleation rate is many orders of magnitude too small even if the parameters are optimized. Subsequently, Sebastian \emph{et al.}~\cite{baum2021nucleation} stressed the importance of nucleation temperature in determining whether the tunneling process can happen and also pointed out that the two-step phase transition will cease due to lack of nucleation temperature and hence the vacuum will be trapped at the trivial minimum. Meanwhile, Thomas \emph{et al.} also found that in the 2HDM~\cite{2hdmtrapped} or N2HDM~\cite{biekotter2021fate} the transition from the electroweak symmetry restored (EW-restored) minimum to the desired EW-breaking vacuum may not be achieved. These studies strongly indicate that when considering the vacuum stability problem, the thermal history may play an important role in addition to traditional analysis. Besides, there is a case that the non-desired vacuum is always the global minimum during the evolution of the universe, which is also need to be excluded.

In this letter, we study the vacuum stability problem incorporating cosmological evolution, and show that it can provide a new limit on new physics models. By comparing results in the singlet extension of the SM (xSM), we reveal that checking the evolution of vacuum from high temperature to low temperature is more reasonable than that from bottom up. Additionally, our method could be computationally economic for some models.

\section{Vacuum stability analysis with cosmological history}

\begin{figure*}[!ht]
\centering
\includegraphics[width=0.49\textwidth]{figs/No_Tc1.pdf}
\includegraphics[width=0.49\textwidth]{figs/No_Tc2.pdf}
\caption{Situations in which the universe always stays in the non-desired vacuum during the thermal history. The solid green ball indicates the vacuum where the universe stays at the corresponding temperature, while the hollow gray ball stands for the SM EW-breaking vacuum at zero-temperature. The arrow with cross means that the transition can not take place. They can be realized in the xSM using parameters shown by the labels. The scaled effective potential ($\Delta V = [V(h,s)-V(0,0)]$) is along the direct path between the two minimums $(0,v_s)$ and $(v_h,0)$.}
\label{fig:bk1}
\end{figure*}

Traditionally, focusing on the effective potential of today's universe, the EW-braking vacuum of new physics models can be divided into two categories, as illustrated by the blue lines in Fig.\ref{fig:bk1}. In the left frame, the desired EW-breaking vacuum of zero-temperature is indicated by gray circle, which is a global minimum. Thus, without considering thermal effect, this desired vacuum is completely safe. A meta-sable desired vacuum is shown in the right frame. If the tunneling time from it to the global minimum, indicated by a green dot, is longer than the age of the universe, we can still retain this model.

However, we find that this analysis is not sufficient. Fig.~\ref{fig:bk1} presents two situations in which today's desired EW-breaking vacuum is stable or long-lived but the universe is trapped in a non-desired vacuum. In the early hot universe, the global minimum (green dots) of effective potential (red lines) is EW-restored. With decreasing temperature a new minimum appears in the potential (orange lines). In the left frame, the new minimum turns into the desired EW-breaking vacuum and global minimum at zero-temperature, but the transition from the origin to it may not happen. In the right frame, the new minimum also evolves into the desired EW-breaking vacuum but always has higher free energy than the origin, so there is no transition at all. Therefore, the correct way to checking the vacuum of today's universe to follow the thermal history, starting from high temperature, calculating every possible transition
and tracing the minimums of potential as they evolve to zero-temperature.

In Ref.~\cite{camargo2013vevacious} thermal transitions are taken into account in the study of vacuum stability, but in the traditional view. For the right frame of Fig.~\ref{fig:bk1}, with increasing temperature, the transition probability of the desired EW-breaking vacuum increases, which makes it a short-lived vacuum. Thus, this method can also correctly exclude this situation. However, if there are more than one minima that the desired EW-breaking vacuum can transition into at finite temperature, this method becomes invalid. Meanwhile, as shown later, the calculation of transition probability is fairly time consuming, while our new method only needs to trace the free energies of the minima for this situation, which precedes calculating transition probability.

An important and time-consuming part to obtain our new limit is determining the occurrence of the thermal transition. Previous studies~\cite{baum2021nucleation, 2hdmtrapped, biekotter2021fate} of the situation shown in the left frame of Fig.~\ref{fig:bk1} adopted the existence of nucleation temperature as a criterion.
The transition process is identified with the nucleation process, in which the surface tension of the vacuum bubble wall and the energy of the bubble expansion compete against each other. As the universe cools, the thermal transition rate increases dramatically. The transition can be regarded as happening once the probability to nucleate one critical bubble in the Hubble volume is of order one, at the so-called nucleation temperature $T_n$, which can be roughly estimated from $S(T_{n})/T_{n} \sim 140$ \ensuremath{\,\text{GeV}}\xspace. However, the nucleation temperature may be ill-defined for slow or supercooled transitions~\cite{cai2017gravitational, Athron:2022mmm}.
Here we calculate the transition probability following \texttt{Vevacious}~\cite{camargo2014constraining}.

The transition rate, per unit time per unit volume, of the meta-stable vacuum at finite temperature can be expressed as~\cite{linde1983decay}
\begin{equation}
    \frac{\Gamma}{V} = B(T)e^{-S_{3}(T)/T} ,
\end{equation}
where $S_{3}$ denotes the Euclidean bounce action over three dimensions and $B(T)$ is a dimensionful parameter, which is often estimated as $T^4$.
The probability $P(T_i, T_f)$ of non-transitioning (surviving) during the period from an initial temperature $T_i$ to a final temperature $T_f$ can be obtained by trading time $t$ with temperature $T$ in the differential equation
\begin{equation}
    P(t+{\rm d}t) = P(t)(1-\Gamma {\rm d}t).
\end{equation}
As a result, we obtain
\small
\begin{equation}
    P(T_i, T_f) = \mathrm{exp}\left[-\int^{T_{f}}_{T_{i}}\frac{{\rm d}t}{{\rm d}T}V(T)B(T)e^{-S_{3}(T)/T}{\rm d}T\right],
\end{equation}
\normalsize
where $V(T)$ is the observable volume.
Assuming the universe is radiation dominated during the evolution from $T_i$ to $T_f$, entropy is approximately conserved and $S_{3}(T)$ increase monotonically with $T$ increasing. The upper bound of $P(T_i, T_f)$ is then expressed as
\small
\begin{eqnarray} \label{eq:vevacious_P}
    P(T_i = T, T_f = 0) \leq  \mathrm{exp}(-1.581 \times 10^{106} e^{-S_{3}/T}/S_{3}).
\end{eqnarray}
\normalsize
A temperature $T_{opt}$ to maximize the right-hand side of Eq.~(\ref{eq:vevacious_P}) can be found by fitting the thermal action to $b/(T-T_c)^2$, where $T_{c}$ is the critical temperature at which the two minima have degenerate free energy. \texttt{Cosmotranstions}~\cite{wainwright2012cosmotransitions} is used to calculate the action. Then we can set a threshold, $20\%$ in this work, on the maximal surviving probability to judge whether the tunneling process occurs or not. The results are not particularly sensitive to the value of the threshold.

\section{Example Toy Model}

The situations displayed in Fig.\ref{fig:bk1} can be realized in a very simple model, the $Z_{2}$ symmetric real scalar extension of the SM (xSM).  In this section we use this model to demonstrate the above generic discussion.

The tree-level effective potential of the xSM can be written in the form of
\begin{equation}\label{eq:xsm_tree}
     V_{0}(h,s) = -\frac{\mu_{h}^{2}}{2}h^2 + \frac{\lambda_{h}}{4}h^4 - \frac{\mu_{s}^2}{2}s^2 + \frac{\lambda_{s}}{4}s^4 + \frac{\lambda_{hs}}{4}h^2s^2,
\end{equation}
where $h$ and $s$ are background field configurations for the Higgs and extra scalar fields, respectively. We build the one-loop effective potential $V(h,s;T)$ in Landau gauge using the on-shell-like renormalization scheme and the Parwani method for daisy resummation~\cite{parwani1992resummation}. See Ref.~\cite{athron2022arbitrary} for detailed expressions and uncertainties induced by these choices. We set $m_h=125~\ensuremath{\,\text{GeV}}\xspace$ and $v_{\rm EW}=246~\ensuremath{\,\text{GeV}}\xspace$, and fix the Lagrangian parameters by tadpole conditions. The remaining free parameters are chosen to be $\lambda_{s}$, $\lambda_{hs}$ and the singlet mass $m_s$.

Fig.\ref{fig:bk1} was, in fact, generated using two benchmark points of the xSM, which have $m_{s} = 234$ and $220~\ensuremath{\,\text{GeV}}\xspace$ respectively, together with fixed $\lambda_{hs}=2.0$ and $\lambda_{s} = 0.2$. To visualize the two-dimensional potential in one-dimension, we only show the potential along the straight line between the two minima in the $(h,s)$ plane.

For the benchmark point in the left frame of Fig.\ref{fig:bk1}, when $T>267.9~\ensuremath{\,\text{GeV}}\xspace$ there is only one minimum in the potential located at the origin. Below that temperature
the singlet field develops a non-zero expectation value smoothly increasing from $0$ to $168.3~\ensuremath{\,\text{GeV}}\xspace$ until zero-temperature.
A new minimum of $(v_h=54.7~\ensuremath{\,\text{GeV}}\xspace,v_s=0)$ appears at $T=147.6~\ensuremath{\,\text{GeV}}\xspace$ and evolves into the desired EW-breaking vacuum $(v_h=246~\ensuremath{\,\text{GeV}}\xspace,v_s=0)$ as the temperature goes to zero. The two minima are degenerated at $T_c=120.9~\ensuremath{\,\text{GeV}}\xspace$, so after that the universe may transition from $(0,v_s)$ to $(v_h,0)$. We find, however, that the maximal probability of non-transitioning $P(T_{opt},0)$ is approximately 100\%. Therefore, it is unlikely that the universe will evolve into the desired EW-breaking vacuum for this benchmark point. Just looking at the zero-temperature potential, however, we see a stable EW-breaking vacuum which has no stability problem. This demonstrates the importance of calculating transitions in the cosmological history.

The benchmark point of the right frame of Fig.\ref{fig:bk1} has similar minimum structure as above. The difference is that there is no critical temperature between the two minima $(0,v_s)$ and $(v_h,0)$, i.e. the free energy of $(0,v_s)$ is larger than that of $(v_h, 0)$ throughout.
As a result, the universe always stays in the EW-restored minimum. In the zero-temperature potential, the desired EW-breaking vacuum is meta-stable with a lifetime much longer than the universe's age, so it is acceptable in the traditional method. For this kind of situation, checking the history to exclude it is time-saving, because it only need to trace the position and free energy of minima, no transition calculation is needed.

Thus, it is necessary to modify the procedure for accessing vacuum stability.  Below, we assume that this is the case for even more complicated models, and the universe is located at the EW-restored vacuum at some high temperature.  Then, we can assess vacuum stability in following steps staring from the EW-restored vacuum.
\begin{itemize}
\item[(1)] Map out the position and free energy of the minima of the effective potential as the function of temperature (e.g. using numerical tools~\cite{wainwright2012cosmotransitions,Athron:2020sbe}).
\item[(2)] As the temperature is lowered,
whenever degenerate minima occur, i.e. if there is a critical temperature $T_{c}$, calculate the transition probability and transition temperature below $T_c$.
\item[(3)] If the transition probability is larger than a certain threshold, set the vacuum of the universe to be the new minimum, and go to Step (2) until $T=0$.
\item[(4)] Check whether the vacuum of the universe at $T=0$ is the desired EW-breaking vacuum. If not, this point is excluded. If yes, it is allowed when the vacuum is stable or long-lived, and is excluded when the vacuum is short-lived.
\end{itemize}

\begin{figure}[t]
\centering
\includegraphics[width=0.5\textwidth]{figs/scan.pdf}
\caption{Classifications of vacuum structure in parameter space of the xSM. In the red region, the desired EW-breaking vacuum is meta-stable with a life-time shorter than the age of the universe. In the orange region, the desired vacuum is long-lived but always has higher free energy than the EW-restored vacuum in the cosmological history. The desired vacuum in the yellow region is stable but the early universe can not transition to it. Only in the green region the universe can safely reach to the stable desired SM vacuum.
Therefore, the dotted and solid curves represent respectively limits from checking minima in the zero-temperature potential and the evolution history of the universe. }
\label{fig:bk3}
\end{figure}

Following this procedure, we explore the parameter space of the xSM. For instance, the result of fixing  $\lambda_{s} = 0.2$ is displayed in Fig.\ref{fig:bk3}, where we compare the limit from only checking zero-temperature potential and the limit from considering  evolution history of the universe.  The parameter space is split by the dashed curve, at which the two minima $(v_h,0)$ and $(0,v_s)$ are degenerate at zero-temperature, and the parameters satisfy
\begin{equation}
    \lambda_{hs} = \frac{v_{h}^2}{2} \left(m_{s}^2 + \mu_{h} \sqrt{\frac{\lambda_{s}}{\lambda_{h}}} \right)
\end{equation}
in the on-shell-like renormalization scheme. The yellow and orange regions correspond the situations described in the left and right frames of Fig.~\ref{fig:bk1}, respectively.

Above the dashed curve, the EW-restored vacuum is the  global minimum at zero-temperature, as well as at finite temperature. The energy gap between the two vacua increases with $\lambda_{hs}$ increasing. Without considering the cosmological history, the desired EW-breaking vacuum is meta-stable, but can be long-lived, such as in the orange region of Fig.~\ref{fig:bk3}. Therefore, the transitional limit is the upper bound of the orange region. Thermal history tells us that the universe starts from the EW-restored vacuum and remains there as it is global minimum throughout. So the orange region is actually ruled out.

Below the dashed curve, the desired EW-breaking vacuum at zero-temperature is the global minimum and thus stable. In the yellow region, the energy gap between the vacua is not large enough to achieve the transition that breaks electroweak symmetry.  Only in the green region, the early universe transits from $(0,v_s)$ or $(0,0)$ to $(v_h,0)$ and then evolves to $(246~\ensuremath{\,\text{GeV}}\xspace,0)$ at $T=0$. The $\lambda_{hs}$ of the new limit~(pink solid curve) is about 0.2 and 0.9 lower than that of the dashed curve and the old limit~(blue dotted curve).

We see that the traditional method accepts some regions of meta-stable EW-breaking vacuum, but checking the cosmological history proves that incorrect. The exclusion of such regions is computationally fast, because it only needs to trace the minima of the potential with temperature. On the other hard, it is more time-consuming to exclude the yellow region, as it involves calculation of transition probability at finite temperature.

The improvement of the new limit on vacuum stability is significant in the xSM, let alone in complicated models. With more than one first-order phase transitions in the thermal history, such as in a supersymmetric model~\cite{Athron:2019teq}, every transition has to be examined carefully. The relationship between the old and new limits may be tangled, and deserve case by case investigation.

\section{Conclusion}

We argue that the study of vacuum stability should involve full consideration of the cosmological history, which is typically ignored. With the experimental discovery of gravitational waves, the impact of BSM on the cosmological history draws attention because a strong first order electroweak phase transition may generate detectable gravitational wave signals. Our study shows that there are parameter space regions where the transition to the desired EW-breaking vacuum is of low probability or does not even exist, which misleads the traditional vacuum stability analysis.

In this letter we studied vacuum stability taking into account the thermal history of the universe. After demonstrating the lack of transition to the desired EW-breaking vacuum in the xSM, we gave a general procedure which incorporates cosmological history in the vacuum stability study, and compared the traditional and new limits. We found that checking the cosmological history can provide a much stringent and sometimes computationally cheaper limit of vacuum stability. In turn, this motivates further studies of phase transitions in new physics models.

{\em Acknowledgments:}
This work was supported by the National Natural Science Foundation of China (2105248, 11821505, 12075300), the Australian Research Council (DP210101636),  Peng-Huan-Wu Theoretical Physics Innovation Center (12047503),
Key R\&D Program of Ministry of Science and Technology (2017YFA0402204), and by the Key Research Program of the Chinese Academy of Sciences (XDPB15).