# Matrix-valued Cauchy bi-orthogonal polynomials and a novel noncommutative
integrable lattice

Shi-Hao Li  Department of Mathematics, Sichuan University, Chengdu, 610064, PR
China<EMAIL_ADDRESS>,  Ying Shi  School of Science, Zhejiang
University of Science and Technology, Hangzhou, 310023, PR China
<EMAIL_ADDRESS>,  Guo-Fu Yu  School of Mathematical Sciences, Shanghai
Jiaotong University, Shanghai, PR China<EMAIL_ADDRESS>and  Jun-Xiao Zhao
School of Mathematics and Statistics, University of Glasgow, University Place,
Glasgow G12 8QQ, UK<EMAIL_ADDRESS>
###### Abstract.

Matrix-valued Cauchy bi-orthogonal polynomials were proposed in this paper,
together with its quasideterminant expression. It is shown that the
coefficients in four-term recurrence relation for matrix-valued Cauchy bi-
orthogonal polynomials should satisfy a novel noncommutative integrable
system, whose Lax pair is given by fractional differential operators with non-
abelian variables.

###### Key words and phrases:

matrix-valued Cauchy bi-orthogonal polynomials; noncommutative C-Toda lattice;
quasideterminants; direct method

###### 2020 Mathematics Subject Classification:

39A36, 15A15

##  1\. introduction

Noncommutative generalizations of integrable systems have been widely studied
with noncommutative space-time variables, which arise from quantization of the
phase space. In kinds of literature, there have been a large number of
noncommutative generalizations of integrable systems, including KP and its
related equations [13, 22, 23], Toda lattice and related equations [17, 18,
29, 30], Painlevé equations [7, 35], and so on.

The study of closed expressions for solutions is one of the most important
topics in soliton theory and Hirota’s bilinear method (or so-called direct
method) is one of the powerful ways of finding solutions [25]. It has been
widely shown that solutions of classical integrable systems could be written
in terms of determinantal or Pfaffian $\tau$-functions. The reason owes to
Sato’s observation, demonstrating that the Plücker relation on Grassmannian is
an integrable system. The Plücker relations are actually identities for
determinants if local coordinates are taken [33]. However, as determinants and
Pfaffians are usually defined by commutative elements and they could hardly
act as $\tau$-functions for noncommutative systems, a new algebraic tool is
needed. The achievement was made by quasideterminants, which was proposed by
Gelfand et al [20, 21, 19], and it was shown that quasideterminants could act
as solutions for noncommutative integrable systems. Besides, quasideterminants
take advantages in various fields such as matrix-valued orthogonal
polynomials, noncommutative combinatorics, noncommutative symmetric functions,
continued fractions with noncommutative coefficients, and so on [12].
Therefore, quasideterminants establish ties for all these subjects.

We mainly consider the connection between a new family of matrix-valued
orthogonal polynomials called matrix-valued Cauchy bi-orthogonal polynomials,
and related noncommutative integrable systems in this paper. The study of
matrix-valued orthogonal polynomials could date back to Krein about matrix-
valued moment problem [28]. In recent years, it was found that matrix-valued
orthogonal polynomials could work as wave functions for noncommutative Toda
and Lotka-Volterra lattices [32, 8, 26], as well as Painlevé equations [9].
Such connections were also applied to noncommutative Hermite-Padé
approximation [14] and Wynn recurrence [15]. By introducing a matrix-valued
inner product

$\displaystyle\langle f(x),g(x)\rangle=\int_{\mathbb{R}}f(x)W(x)g(x)dx$

where $f(x),\,g(x)\in\mathbb{R}^{p\times p}[x]$ are matrix-valued polynomials
and $W(x):\,\mathbb{R}\to\mathbb{R}^{p\times p}$ is a symmetric matrix-valued
weight function satisfying $W(x)=W^{\top}(x)$, one could define a family of
matrix-valued orthogonal polynomials $\\{P_{n}(x)\\}_{n\in\mathbb{N}}$ through
the orthogonal relation $\langle
P_{n}(x),P^{\top}_{m}(x)\rangle=H_{n}\delta_{n,m}$ for some nonsingular
normalization factor $H_{n}$. It was shown that matrix-valued orthogonal
polynomials satisfy a three-term recurrence relation

$\displaystyle xP_{n}(x)=P_{n+1}(x)+a_{n}P_{n}(x)+b_{n}P_{n-1}(x),\quad
a_{n},\,b_{n}\in\mathbb{R}^{p\times p},\quad n\in\mathbb{N}.$

Moreover, if time parameters $\\{t_{i},i\in\mathbb{N}\\}$ are introduced into
the weight function $W$ such that
$\partial_{t_{i}}W(x;\mathbf{t})=x^{i}W(x;\mathbf{t})$, then
$\partial_{t_{1}}P_{n}(x;\mathbf{t})=-b_{n}P_{n-1}(x;\mathbf{t}),$ and a
noncommutative version of Toda hierarchy could be obtained from compatibility
condition
$\partial_{t_{1}}(xP_{n}(x;\mathbf{t}))=x\partial_{t_{1}}P_{n}(x;\mathbf{t})$
[32].

We plan to study a matrix-valued generalization of Cauchy bi-orthogonal
polynomials and its related noncommutative integrable systems in this paper.
In recent years, with the development of interdisciplinary research on random
matrices, orthogonal polynomials and integrable systems, connections between
variants of orthogonal polynomials and integrable systems have been made, see
e.g. [1, 2, 11, 31] and references therein. Cauchy bi-orthogonal polynomials
were proposed in the study of the Degasperis-Procesi peakon problem and later
found its applications in random matrix model [34, 5, 6]. With the Cauchy
kernel $(x+y)^{-1}$ and nonnegative weights $\omega_{1},\,\omega_{2}$, an
inner product over $\mathbb{R}[x]\times\mathbb{R}[y]\to\mathbb{R}$ could be
introduced

$\displaystyle\langle
f(x),g(y)\rangle=\int_{\mathbb{R}_{+}\times\mathbb{R}_{+}}\frac{f(x)g(y)}{x+y}\omega_{1}(x)\omega_{2}(y)dxdy.$
(1.1)

From the inner product, one can define Cauchy bi-orthogonal polynomials
$\\{P_{n}(x),Q_{n}(y)\\}_{n\in\mathbb{N}}$, which are dual to each other, by
orthogonality

$\displaystyle\langle P_{n}(x),Q_{m}(y)\rangle=h_{n}\delta_{n,m}$

for some nonzero constant $h_{n}\in\mathbb{R}$. In [31], it was found that
symmetric Cauchy bi-orthogonal polynomials (i.e. when weight functions $w_{1}$
and $w_{2}$ are equivalent) are related to an integrable lattice equation,
which enjoys the same $\tau$-function with CKP equation111CKP equation is a
2-dimensional generalization of Kaup-Kuperschmidt equation. It is called CKP
equation because it is related to the infinite dimensional Lie algebra of type
$C$ according to the classification of classical integrable systems [27]., and
thus we call the corresponding lattice equation as the Toda lattice of CKP
type or C-Toda lattice for brevity. Therefore, it is of interest to consider a
matrix-valued version of Cauchy bi-orthogonal polynomials and a noncommutative
generalization of C-Toda lattice. For this purpose, we introduce a symmetric
matrix-valued inner product

$\displaystyle\langle\cdot,\cdot\rangle:\mathbb{R}^{p\times
p}[x]\times\mathbb{R}^{p\times p}[y]\to\mathbb{R}^{p\times p},\quad\langle
P(x),Q(y)\rangle=\int_{\mathbb{R}_{+}\times\mathbb{R}_{+}}\frac{P(x)W(x)W(y)Q(y)}{x+y}dxdy$

with respect to a matrix-valued weight function $W(x)$.

Therefore, one main purpose of this paper is to study the matrix-valued Cauchy
bi-orthogonal polynomials, including their quasi-determinant expressions and
corresponding recurrence relations. We show that the coefficients of the
recurrence relations can also be expressed in terms of quasideterminants. The
other purpose is to introduce proper time flows, from which a novel
noncommutative integrable system is obtained by compatibility condition. This
paper is organized as follows. In Section 2, we give a formal definition for
matrix-valued Cauchy bi-orthogonal polynomials and their quasi-determinant
expressions. We derive a four-term recurrence relation for the matrix-valued
Cauchy bi-orthogonal polynomials and show the coefficients in the recurrence
relation could also be written in terms of quasideterminants. In Section 3, we
introduce proper time flows and hence derive a novel noncommutative integrable
lattice. We provide a Lax pair to this novel noncommutative integrable
lattice, and verify it through a direct method by quasideterminant identities.
Some concluding remarks are given in Section 4.

##  2\. Matrix-valued Cauchy bi-orthogonal polynomials

The theory of matrix orthogonal polynomials has attracted much attention in
different aspects of mathematics such as approximation problems on Riemann
sphere [4], quasi-birth-and-death processes [24], random tiling problems [16]
and non-commutative integrable systems [3, 32]. Such matrix-valued orthogonal
polynomial is a non-commutative generalization of standard orthogonal
polynomial and related matrix-valued moment problem was proposed by Krein in
late 1940s [28]. In this setting, a bounded matrix-valued Borel measure $d\mu$
on the real line is considered, which maps from the real line $\mathbb{R}$ to
a matrix ring $\mathbb{R}^{p\times p}$. According to the Radon-Nikodym
theorem, such a measure could be expressed in terms of a weight function
$W(x)dx$, and a matrix-valued inner product could be written as

$\displaystyle\langle\cdot,\cdot\rangle:\mathbb{R}^{p\times
p}[x]\times\mathbb{R}^{p\times p}[x]\rightarrow\mathbb{R}^{p\times
p},\quad\langle P(x),Q(x)\rangle=\int_{\mathbb{R}}P(x)W(x)Q(x)dx.$

To connect it with non-abelian integrable lattices, we require that $W(x)$ to
be symmetric and there exists a family of matrix-valued orthogonal polynomials
$\\{P_{n}(x)\\}_{n\in\mathbb{N}}$ satisfying orthogonal relation

$\displaystyle\langle P_{n}(x),P_{m}^{\top}(x)\rangle=H_{n}\delta_{n,m},$

where $H_{n}$ is a nonsingular matrix in $\mathbb{R}^{p\times p}$ and $\top$
represents the transpose of the matrix-valued orthogonal polynomials.

Here we propose a bi-orthogonal generalization of the matrix-valued orthogonal
polynomials. In general, we could consider an inner product over
$\mathbb{R}[x]\times\mathbb{R}[y]$ such that

$\displaystyle\langle
P(x),Q(y)\rangle=\int_{\mathbb{R}\times\mathbb{R}}P(x)W_{1}(x)K(x,y)W_{2}(y)Q(y)dxdy,$

where $K(x,y)$ is a matrix-valued kernel or a scalar kernel. In this paper, we
consider a specific scalar Cauchy kernel $(x+y)^{-1}$ and an inner product
with respect to this kernel is defined by

$\displaystyle\langle\cdot,\cdot\rangle:\,\mathbb{R}^{p\times
p}[x]\times\mathbb{R}^{p\times p}[y]\to\mathbb{R}^{p\times p},\quad\langle
P(x),Q(y)\rangle=\int_{\mathbb{R}_{+}\times\mathbb{R}_{+}}\frac{P(x)W_{1}(x)W_{2}(y)Q(y)}{x+y}dxdy.$
(2.1)

Here we require the integration space is a subset of
$\mathbb{R}_{+}\times\mathbb{R}_{+}$ to avoid the singularity from the kernel.
$W_{1}(x)$ and $W_{2}(y)$ are thus some matrix-valued weight functions from
$\mathbb{R}_{+}$ to $\mathbb{R}^{p\times p}$. For example, $W_{1}$ and $W_{2}$
could be taken as some matrix-valued Laguerre weights. Clearly, we have the
following properties for the inner product

  1. (1)

Linearity. For any $L_{1},\,L_{2},\,R_{1},\,R_{2}\in\mathbb{R}^{p\times p}$,
we have

$\displaystyle\langle L_{1}P_{1}(x)+L_{2}P_{2}(x),Q(y)\rangle=L_{1}\langle
P_{1}(x),Q(y)\rangle+L_{2}\langle P_{2}(x),Q(y)\rangle,$ $\displaystyle\langle
P(x),Q_{1}(y)R_{1}+Q_{2}(y)R_{2}\rangle=\langle P(x),Q_{1}(y)\rangle
R_{1}+\langle P(x),Q_{2}(y)\rangle R_{2}.$

  2. (2)

If $W_{1}$ and $W_{2}$ are symmetric, then $\langle
P(x),Q(y)\rangle^{\top}=\langle Q^{\top}(x),P^{\top}(y)\rangle$.

Besides, we need two requirements on the inner product. One is the existence
and finiteness of the moments

$\displaystyle m_{i,j}=\langle
x^{i}\mathbb{I}_{p},y^{j}\mathbb{I}_{p}\rangle=\int_{\mathbb{R}_{+}\times\mathbb{R}_{+}}\frac{x^{i}y^{j}}{x+y}W_{1}(x)W_{2}(y)dxdy\in\mathbb{R}^{p\times
p},\quad i,\,j\in\mathbb{N}.$

We need to require that $W_{1}(x)$ and $W_{2}(y)$ should be properly taken
such that all $\\{m_{i,j}\\}_{i,j\in\mathbb{N}}$ exist and are finite. The
second is the invertibility of principle minors of moment matrices. For any
$n\in\mathbb{N}$, principle minors of moment matrices
$\left(m_{i,j}\right)_{i,j=0}^{n}$ should be invertible. We call these
requirements as moment conditions.

Now we could give the following definition for matrix-valued Cauchy bi-
orthogonal polynomials.

######  Definition 2.1.

With inner product (2.1) satisfying moment conditions, a family of monic
matrix-valued Cauchy bi-orthogonal polynomials
$\\{P_{n}(x),Q_{n}(x)\\}_{n=0}^{\infty}$ are defined by222A monic matrix-
valued polynomial in the matrix-valued polynomial ring $\mathbb{R}^{p\times
p}[x]$ is set to be like
$\mathbb{I}_{p}x^{n}+a_{n,n-1}x^{n-1}+\cdots+a_{n,0}$, where
$a_{n,i}\in\mathbb{R}^{p\times p}$ for all $i=0,\cdots,n-1$, and
$\mathbb{I}_{p}$ is the unity in the matrix ring.

$\displaystyle\langle P_{n}(x),Q_{m}(y)\rangle=H_{n}\delta_{n,m}$ (2.2)

for some $H_{n}\in\mathbb{R}^{p\times p}$.

In fact, the orthogonal relation (2.2) is equivalent to

$\displaystyle\langle P_{n}(x),y^{i}\mathbb{I}_{p}\rangle=0,\quad 0\leq i\leq
n-1$

according to the linearity of the inner product. Moreover, if we assume that

$\displaystyle
P_{n}(x)=\mathbb{I}_{p}x^{n}+\xi_{n,n-1}x^{n-1}+\cdots+\xi_{n,0}\in\mathbb{R}^{p\times
p}[x],$

then the above orthogonal relation is in fact a combination of linear systems
with matrix elements coefficients

$\displaystyle\xi_{n,0}m_{0,j}+\xi_{n,1}m_{1,j}+\cdots+\xi_{n,n-1}m_{n-1,j}+m_{n,j}=0,\quad
j=0,\cdots,n-1.$

Therefore, the existence and uniqueness of bi-orthogonal polynomials are
equivalent to the existence and uniqueness for solutions in noncommutative
linear systems. According to the quasideterminant theory proposed by Gelfand
et al [19] (a quick introduction is given in the appendix), if moment
conditions are satisfied, then the coefficients of matrix-valued Cauchy bi-
orthogonal polynomials are solved to be

$\displaystyle\xi_{n,i}=-\left(m_{n,0},\cdots,m_{n,n-1}\right)\left(\begin{array}[]{ccc}m_{0,0}&\cdots&m_{0,n-1}\\\
\vdots&&\vdots\\\
m_{n-1,0}&\cdots&m_{n-1,n-1}\end{array}\right)^{-1}e_{i+1}^{\top},$ (2.6)

where $e_{i}=(0,\cdots,\mathbb{I}_{p},\cdots,0)$ is the block unit vector,
whose $i$-th element is the unity and the others are zeros. By substituting
those coefficients into the expression of $P_{n}(x)$, we obtain

$\displaystyle P_{n}(x)$
$\displaystyle=x^{n}\mathbb{I}_{p}-\left(m_{n,0},\cdots,m_{n,n-1}\right)\left(\begin{array}[]{ccc}m_{0,0}&\cdots&m_{0,n-1}\\\
\vdots&&\vdots\\\
m_{n-1,0}&\cdots&m_{n-1,n-1}\end{array}\right)^{-1}\left(\begin{array}[]{c}\mathbb{I}_{p}\\\
\vdots\\\ x^{n-1}\mathbb{I}_{p}\end{array}\right)$
$\displaystyle=\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&\mathbb{I}_{p}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&x^{n-1}\mathbb{I}_{p}\\\
m_{n,0}&\cdots&m_{n,n-1}&\boxed{x^{n}\mathbb{I}_{p}}\end{array}\right|.$

In a similar manner, we could get an explicit formula for $Q_{n}(x)$ and

$\displaystyle
Q_{n}(x)=\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&m_{0,n}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&m_{n-1,n}\\\
\mathbb{I}_{p}&\cdots&x^{n-1}\mathbb{I}_{p}&\boxed{x^{n}\mathbb{I}_{p}}\end{array}\right|.$

To find its connection with integrable systems, we need to further require
that weight functions $W_{1}$ and $W_{2}$ are identical and symmetric, i.e.
$W_{1}=W_{2}=W_{1}^{\top}$. With this assumption, we know that moments are
symmetric, i.e. $m_{i,j}=m_{j,i}^{\top}$, and we have the following lemma.

######  Lemma 2.2.

With weight functions set to be identical and symmetric, matrix-valued Cauchy
bi-orthogonal polynomials defined in (2.2) have the relation

$\displaystyle Q_{n}^{\top}(x)=P_{n}(x).$

###### Proof.

This is a direct verification that

$\displaystyle
Q_{n}^{\top}(x)=x^{n}\mathbb{I}_{p}-(m_{0,n}^{\top},\cdots,m_{n-1,n}^{\top})\left(\begin{array}[]{ccc}m_{0,0}^{\top}&\cdots&m_{n-1,0}^{\top}\\\
\vdots&&\vdots\\\
m_{0,n-1}^{\top}&\cdots&m_{n-1,n-1}^{\top}\end{array}\right)^{-1}\left(\begin{array}[]{c}\mathbb{I}_{p}\\\
\vdots\\\ x^{n-1}\mathbb{I}_{p}\end{array}\right)=P_{n}(x).$

∎

Since $Q_{n}(y)$ could be expressed in terms of $P_{n}(y)$, there is only one
family of symmetric Cauchy bi-orthogonal polynomials. We call this family of
polynomials as matrix-valued Cauchy orthogonal polynomials from now on.

To summarize, we have the following proposition.

######  Proposition 2.3.

With inner product (2.1) in which $W(x)$ is set to be symmetric, a family of
Cauchy orthogonal polynomials $\\{P_{n}(x)\\}_{n\in\mathbb{N}}$ is given by
the orthogonal relation

$\displaystyle\langle P_{n}(x),P_{m}^{\top}(y)\rangle=H_{n}\delta_{n,m},$

where $P_{n}(x)$ and $H_{n}$ have the following quasideterminant expressions

$\displaystyle
P_{n}(x)=\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&\mathbb{I}_{p}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&x^{n-1}\mathbb{I}_{p}\\\
m_{n,0}&\cdots&m_{n,n-1}&\boxed{x^{n}\mathbb{I}_{p}}\end{array}\right|,\quad
H_{n}=\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&m_{0,n}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&m_{n-1,n}\\\
m_{n,0}&\cdots&m_{n,n-1}&\boxed{m_{n,n}}\end{array}\right|$ (2.15)

with $H_{n}=H_{n}^{\top}$.

###  2.1. Recurrence relation for matrix-valued Cauchy orthogonal polynomials

The most important characterization in the inner product (2.1) is the Cauchy
kernel $(x+y)^{-1}$. From this inner product, we have

$\displaystyle\begin{aligned} \langle
xf(x),g(y)\rangle&=\int_{\mathbb{R}_{+}\times\mathbb{R}_{+}}\frac{xf(x)W(x)W(y)g(y)}{x+y}dxdy\\\
&=\int_{\mathbb{R}_{+}}f(x)W(x)dx\int_{\mathbb{R}_{+}}W(y)g(y)dy-\langle
f(x),yg(y)\rangle,\end{aligned}$ (2.16)

which holds for arbitrary $f(x)\in\mathbb{R}^{p\times p}[x]$,
$g(y)\in\mathbb{R}^{p\times p}[y]$. This formula is useful when inducing
corresponding recurrence relations.

######  Proposition 2.4.

The matrix-valued Cauchy orthogonal polynomials
$\\{P_{n}(x)\\}_{n\in\mathbb{N}}$ satisfy the following recurrence relation

$\displaystyle
x(P_{n+1}(x)+a_{n}P_{n}(x))=P_{n+2}(x)+b_{n}P_{n+1}(x)+c_{n}P_{n}(x)+d_{n}P_{n-1}(x)$
(2.17)

for some certain $a_{n}$, $b_{n}$, $c_{n}$, $d_{n}\in\mathbb{R}^{p\times p}$.

###### Proof.

From the formula (2.16), we know that if we take

$\displaystyle
a_{n}=-\left(\int_{\mathbb{R}_{+}}P_{n+1}(x)W(x)dx\right)\left(\int_{\mathbb{R}_{+}}P_{n}(x)W(x)dx\right)^{-1},$
(2.18)

then

$\displaystyle\langle x(P_{n+1}(x)+a_{n}P_{n}(x)),Q(y)\rangle=-\langle
P_{n+1}(x)+a_{n}P_{n}(x),yQ(y)\rangle.$

This is a key ingredient to deriving the recurrence formula. Moreover, since
$\\{P_{n}(x)\\}_{n\in\mathbb{N}}$ form a basis in the Hilbert space, we have
the following expansion formula

$\displaystyle
x(P_{n+1}(x)+a_{n}P_{n}(x))=P_{n+2}(x)+\sum_{i=0}^{n+1}\alpha_{n,i}P_{i}(x),$

where $\alpha_{n,i}\in\mathbb{R}^{p\times p}$ are some coefficients.
Therefore, by taking $a_{n}$ as in (2.18), we have

$\displaystyle\langle
P_{n+2}(x)+\sum_{i=0}^{n+1}\alpha_{n,i}P_{i}(x),P_{m}^{\top}(y)\rangle$
$\displaystyle=\langle x(P_{n+1}(x)+a_{n}P_{n}(x)),P_{m}^{\top}(y)\rangle$
$\displaystyle=-\langle P_{n+1}(x)+a_{n}P_{n}(x),yP_{m}^{\top}(y)\rangle.$

According to the orthogonal relation, we know that when $m<n-1$, the right-
hand side is definitely zero, which means that
$\alpha_{n,0}=\cdots=\alpha_{n,n-2}=0$. Therefore, the four-term recurrence
relation (2.17) is obtained.

Moreover, we could write down the recurrence coefficients with the help of the
inner product. By Taking the inner product with $P_{n-1}(y)$, $P_{n}(y)$ and
$P_{n+1}(y)$ respectively, we have

$\displaystyle d_{n}$ $\displaystyle=-\left\langle
P_{n+1}(x)+a_{n}P_{n}(x),yP_{n-1}^{\top}(y)\right\rangle
H_{n-1}^{-1}=-a_{n}H_{n}H_{n-1}^{-1},$ $\displaystyle c_{n}$
$\displaystyle=-\left\langle
P_{n+1}(x)+a_{n}P_{n}(x),yP_{n}^{\top}(y)\right\rangle
H_{n}^{-1}=-\left(H_{n+1}+a_{n}\left\langle
P_{n}(x),yP_{n}^{\top}(y)\right\rangle\right)H_{n}^{-1},$ $\displaystyle
b_{n}$ $\displaystyle=\left\langle
x(P_{n+1}(x)+a_{n}P_{n}(x)),P_{n+1}^{\top}(y)\right\rangle
H_{n+1}^{-1}=\left\langle xP_{n+1}(x),P_{n+1}^{\top}(y)\right\rangle
H_{n+1}^{-1}+a_{n}.$

∎

It should be emphasized that those coefficients could be written in terms of
quasideterminants. It is obvious that

$\displaystyle
V_{n}:=\int_{\mathbb{R}_{+}}P_{n}(x)W(x)dx=\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&\phi_{0}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&\phi_{n-1}\\\
m_{n,0}&\cdots&m_{n,n-1}&\boxed{\phi_{n}}\end{array}\right|,\quad\phi_{i}=\int_{\mathbb{R}_{+}}x^{i}W(x)dx\in\mathbb{R}^{p\times
p}$ (2.23)

according to the quasi-determinant formula (2.15) and the linearity. To ensure
that $V_{n}$ is well-defined, we assume that $\\{\phi_{i}\\}_{i\in\mathbb{N}}$
exist and are finite. Moreover, according to the exact formula of $a_{n}$ in
(2.18), we should assume that $V_{n}$ is nonsingular to ensure that it has a
proper inverse.

Moreover, if we expand
$P_{n}(x)=x^{n}\mathbb{I}_{p}+\xi_{n,n-1}x^{n-1}+\cdots+\xi_{n,0}$, then we
have

$\displaystyle b_{n}=Z^{\top}_{n+1}H_{n+1}^{-1}+\xi_{n+1,n}+a_{n},\quad
c_{n}=-H_{n+1}H_{n}^{-1}-a_{n}Z_{n}H_{n}^{-1}-a_{n}H_{n}\xi_{n,n-1}^{\top}H_{n}^{-1}$

where

$\displaystyle Z_{n}=\langle
P_{n}(x),y^{n+1}\mathbb{I}_{p}\rangle=\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&m_{0,n+1}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&m_{n-1,n+1}\\\
m_{n,0}&\cdots&m_{n,n-1}&\boxed{m_{n,n+1}}\end{array}\right|.$

It is of interest to notice that $Z_{n}$ could be written in terms of
$\xi_{n+1,n}$.

######  Proposition 2.5.

It holds that

$\displaystyle Z_{n}^{\top}=-\xi_{n+1,n}H_{n}.$

###### Proof.

According to (2.6), it is known that

$\displaystyle\xi_{n+1,n}=\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n}&0\\\
\vdots&&\vdots&\vdots\\\ m_{n,0}&\cdots&m_{n,n}&\mathbb{I}_{p}\\\
m_{n+1,0}&\cdots&m_{n+1,n}&\boxed{0}\end{array}\right|,$

then application of the noncommutative Jacobi identity (A.12) gives the
desired formula. ∎

Therefore, coefficients in the four-term recurrence relation (2.17) could be
expressed in terms of quasideterminants and we have

$\displaystyle\begin{aligned} &a_{n}=-V_{n+1}V_{n}^{-1},\quad
b_{n}=-\xi_{n+2,n+1}+\xi_{n+1,n}+a_{n},\\\
&c_{n}=-H_{n+1}H_{n}^{-1}+a_{n}H_{n}\left(\xi_{n+1,n}^{\top}-\xi_{n,n-1}^{\top}\right)H_{n}^{-1},\quad
d_{n}=-a_{n}H_{n}H_{n-1}^{-1}.\end{aligned}$ (2.24)

##  3\. Time evolutions and non-commutative C-Toda lattice

This part is devoted to the discussion of how to introduce time flows into the
matrix-valued Cauchy bi-orthogonal polynomials. In fact, formula (2.17)
indicates the following spectral problem

$\displaystyle\left(\begin{array}[]{c}P_{n+2}(x)\\\ P_{n+1}(x)\\\
P_{n}(x)\end{array}\right)=\left(\begin{array}[]{ccc}x-b_{n}&a_{n}x-c_{n}&-d_{n}\\\
\mathbb{I}_{p}&0&0\\\
0&\mathbb{I}_{p}&0\end{array}\right)\left(\begin{array}[]{c}P_{n+1}(x)\\\
P_{n}(x)\\\ P_{n-1}(x)\end{array}\right),$

and it is expected to find proper time flows such that an integrable system
could be derived.

Let’s assume that there is a family of odd time variables
$\mathbf{t}=(t_{1},t_{3},\cdots)$ and put them into the weight functions such
that333Since the commutative case belongs to the CKP hierarchy, we constrain
ourselves to the odd flows.

$\displaystyle
W(x;\mathbf{t})=\exp\left(\sum_{i=1}^{\infty}t_{2i+1}x^{2i+1}\right)W(x).$
(3.1)

In literature, those time variables could be possibly infinite many. For
convenience, let’s consider $t_{1}$-flow, which results in the first non-
trivial example in the hierarchy. The role of $t_{1}$-flow is equivalent to a
ladder operator in the orthogonal polynomial theory in some sense since
$\partial_{t_{1}}W(x;\mathbf{t})=xW(x;\mathbf{t})$.

Therefore, the inner product (2.1) becomes a time-deformed inner product

$\displaystyle\langle
f(x),g(y)\rangle=\int_{\mathbb{R}_{+}\times\mathbb{R}_{+}}\frac{f(x)W(x;\mathbf{t})W(y;\mathbf{t})g(y)}{x+y}dxdy,$

and matrix-valued Cauchy bi-orthogonal polynomials are time-dependent, denoted
by $\\{P_{n}(x;\mathbf{t})\\}_{n\in\mathbb{N}}$. Since time variables are
added in the weight function, it doesn’t make any influence on the spectral
part. Moreover, the moments $\\{m_{i,j}\\}_{i,j\in\mathbb{N}}$ satisfy the
time evolutions

$\displaystyle\partial_{t_{1}}m_{i,j}=m_{i+1,j}+m_{i,j+1}=\phi_{i}\phi_{j},$

where $\phi_{i}$ is given in (2.23). With this assumption, $\xi_{n,n-1}$ and
$H_{n}$ have the following relation.

######  Proposition 3.1.

It holds that

$\displaystyle\partial_{t_{1}}H_{n}=\left(\xi_{n,n-1}-\xi_{n+1,n}\right)H_{n}+H_{n}\left(\xi_{n,n-1}^{T}-\xi_{n+1,n}^{T}\right).$
(3.2)

###### Proof.

With $\partial_{t_{1}}m_{i,j}=m_{i+1,j}+m_{i,j+1}$ and applications of
derivative formula (A.16), we notice that

$\displaystyle\partial_{t_{1}}H_{n}$
$\displaystyle=\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&m_{0,n+1}+m_{1,n}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&m_{n-1,n+1}+m_{n,n}\\\
m_{n,0}&\cdots&m_{n,n-1}&\boxed{m_{n,n+1}}\end{array}\right|+\text{(transposed
part)}$ $\displaystyle+mM^{-1}(A_{1}+A_{2})M^{-1}m^{\top}$

where

$\displaystyle m=(m_{n,0},\cdots,m_{n,n-1}),\quad
M=\left(m_{i,j}\right)_{i,j=0}^{n-1},\quad
A_{1}=\left(m_{i+1,j}\right)_{i,j=0}^{n-1},\quad
A_{2}=\left(m_{i,j+1}\right)_{i,j=0}^{n-1}.$

By inserting an identity matrix, we have

$\displaystyle
mM^{-1}\left(\sum_{j=1}^{n}e_{j}e_{j}^{\top}\right)A_{1}M^{-1}m^{\top}=\sum_{j=1}^{n}\left|\begin{array}[]{cc}M&e_{j}\\\
m&\boxed{0}\end{array}\right|\cdot\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&m_{0,n}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&m_{n-1,n}\\\
m_{j,0}&\cdots&m_{j,n-1}&\boxed{0}\end{array}\right|.$

Noting that if two rows are identical, then the corresponding quasideterminant
is zero [19, Prop. 1.4.6], we know that the above formula is equal to

$\displaystyle-\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&m_{1,n}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&m_{n,n}\\\
m_{n,0}&\cdots&m_{n,n-1}&\boxed{0}\end{array}\right|+\left|\begin{array}[]{cc}M&e_{n}\\\
m&\boxed{0}\end{array}\right|\cdot\left|\begin{array}[]{cc}M&m^{\top}\\\
m&\boxed{m_{n,n}}\end{array}\right|.$

Therefore,

$\displaystyle\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&m_{0,n+1}+m_{1,n}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&m_{n-1,n+1}+m_{n,n}\\\
m_{n,0}&\cdots&m_{n,n-1}&\boxed{m_{n,n+1}}\end{array}\right|+mM^{-1}A_{1}M^{-1}m^{\top}=Z_{n}+\xi_{n,n-1}H_{n}.$

Similarly, the transposed part and $mM^{-1}A_{2}M^{-1}m^{\top}$ give the rest
part, and the proof is complete. ∎

Moreover, we have the following derivative formula for the polynomials.

######  Proposition 3.2.

Time-dependent Cauchy orthogonal polynomials
$\\{P_{n}(x;\mathbf{t})\\}_{n\in\mathbb{N}}$ satisfy the evolution equation

$\displaystyle\partial_{t_{1}}\left(P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t})\right)=\partial_{t_{1}}\left(\xi_{n+1,n}+a_{n}\right)P_{n}(x;\mathbf{t}),$
(3.3)

where $a_{n}$ is given in (2.24).

###### Proof.

From the orthogonal relation (2.2), we have the following formula

$\displaystyle\langle
P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t}),P_{m}^{\top}(y;\mathbf{t})\rangle=H_{n+1}\delta_{n+1,m}+a_{n}H_{n}\delta_{n,m}.$
(3.4)

By taking the $t_{1}$-derivative on both sides and noting that

$\displaystyle\langle
x(P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t})),P_{m}^{\top}(y;\mathbf{t})\rangle+\langle
P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t}),yP_{m}^{\top}(y;\mathbf{t})\rangle=0.$
(3.5)

we obtain the formula

$\displaystyle\begin{aligned}
\partial_{t_{1}}H_{n+1}\delta_{n+1,m}+\partial_{t_{1}}(a_{n}H_{n})\delta_{n,m}&=\langle\partial_{t_{1}}(P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t})),P_{m}^{\top}(y;\mathbf{t})\rangle\\\
&+\langle
P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t}),\partial_{t_{1}}P_{m}^{\top}(y;\mathbf{t})\rangle.\end{aligned}$
(3.6)

Therefore, when $m<n$, it is known that the left-hand side of the above
formula is equal to zero. Moreover, we have
$\text{deg}(\partial_{t_{1}}P_{m}^{\top})<\text{deg}P_{m}^{\top}$, and thus

$\displaystyle\langle\partial_{t_{1}}(P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t})),P_{m}^{\top}(y;\mathbf{t})\rangle=0,\quad
m<n.$ (3.7)

On the other hand, since $\\{P_{n}(x;\mathbf{t})\\}_{n\in\mathbb{N}}$ span a
polynomial basis in the Hilbert space, we have

$\displaystyle\partial_{t_{1}}(P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t}))=\sum_{k=0}^{n}\gamma_{n,k}P_{k}(x;\mathbf{t}).$

By taking it into (3.7), we could get $\gamma_{n,0}=\cdots=\gamma_{n,n-1}=0$.
By comparing the coefficients of $x^{n}$ on both sides, we get that
$\gamma_{n,n}=\partial_{t_{1}}(\xi_{n+1,n}+a_{n}).$ ∎

######  Remark 3.3.

The reason why we consider
$\partial_{t_{1}}(P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t}))$ is due to
the equation (3.5).

Below, we derive an integrable lattice from (3.6), which is equivalent to the
compatibility condition of spectral problem (2.17) and time evolution (3.3).
From (3.6), it is known that when $m=n$, we have

$\displaystyle\partial_{t_{1}}(a_{n}H_{n})=\langle\gamma_{n,n}P_{n}(x;\mathbf{t}),P_{n}^{\top}(y;\mathbf{t})\rangle=\partial_{t_{1}}(\xi_{n+1,n}+a_{n})H_{n},$

and when $m=n+1$, we have

$\displaystyle\partial_{t_{1}}H_{n+1}=\langle
P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t}),\partial_{t_{1}}P_{n+1}^{\top}(y;\mathbf{t})\rangle=a_{n}H_{n}\partial_{t_{1}}\xi_{n+1,n}^{\top}.$

To conclude, we have the following proposition.

######  Theorem 3.4.

The non-commutative C-Toda lattice is

$\displaystyle\begin{aligned}
\left\\{\begin{array}[]{l}a_{n}\cdot\partial_{t_{1}}H_{n}=\partial_{t_{1}}\xi_{n+1,n}\cdot
H_{n},\\\ \partial_{t_{1}}H_{n+1}=a_{n}\cdot
H_{n}\cdot\partial_{t_{1}}\xi_{n+1,n}^{\top},\\\
\partial_{t_{1}}H_{n}=\left(\xi_{n,n-1}-\xi_{n+1,n}\right)H_{n}+H_{n}\left(\xi_{n,n-1}^{\top}-\xi_{n+1,n}^{\top}\right).\end{array}\right.\end{aligned}$
(3.8)

with Lax pair

$\displaystyle
x\left(P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t})\right)=P_{n+2}(x;\mathbf{t})+\left(-\xi_{n+2,n+1}+\xi_{n+1,n}+a_{n}\right)P_{n+1}(x;\mathbf{t})$
$\displaystyle\qquad-\left(H_{n+1}H_{n}^{-1}+a_{n}H_{n}\left(-\xi_{n+1,n}^{\top}+\xi_{n,n-1}^{\top}\right)H_{n}^{-1}\right)P_{n}(x;\mathbf{t})-a_{n}H_{n}H_{n-1}^{-1}P_{n-1}(x;\mathbf{t}),$
$\displaystyle\partial_{t_{1}}(P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t}))=\partial_{t_{1}}(\xi_{n+1,n}+a_{n})P_{n}(x;\mathbf{t}).$

It should be noticed that there is also a matrix Lax pair. If we denote

$\displaystyle\Phi=\left(\begin{array}[]{c}P_{0}(x)\\\ P_{1}(x)\\\
\vdots\end{array}\right),\quad
A=\left(\begin{array}[]{cccc}a_{0}&\mathbb{I}_{p}&&\\\
&a_{1}&\mathbb{I}_{p}&\\\ &&\ddots&\ddots\end{array}\right),\quad
B=\left(\begin{array}[]{ccccc}c_{0}&b_{0}&\mathbb{I}_{p}&&\\\
d_{1}&c_{1}&b_{1}&\mathbb{I}_{p}&\\\
&\ddots&\ddots&\ddots&\ddots\end{array}\right),$

and $C=\text{diag}(e_{0},e_{1},\cdots)$ where
$e_{i}=\partial_{t_{1}}\xi_{i+1,i}$, then the four-term recurrence relation
and the derivative formula is equivalent to

$\displaystyle Ax\Phi=B\Phi,\quad A\partial_{t}\Phi=C\Phi,$

and the compatibility condition gives

$\displaystyle\partial_{t_{1}}\mathcal{L}=[\mathcal{N},\mathcal{L}]=\mathcal{N}\mathcal{L}-\mathcal{L}\mathcal{N},$
(3.9)

where $\mathcal{L}=A^{-1}B$ and $\mathcal{N}=A^{-1}C$.

######  Remark 3.5.

In the commutative case, if we denote $\tau$-function $\tau_{n}$ and auxiliary
function $\sigma_{n}$ as

$\displaystyle\tau_{n}=\det\left(\begin{array}[]{ccc}m_{0,0}&\cdots&m_{0,n-1}\\\
\vdots&&\vdots\\\
m_{n-1,0}&\cdots&m_{n-1,n-1}\end{array}\right),\quad\sigma_{n}=\det\left(\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&\phi_{0}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&\phi_{n-1}\\\
m_{n,0}&\cdots&m_{n,n-1}&\phi_{n}\end{array}\right),$

then we have

$\displaystyle H_{n}=\frac{\tau_{n+1}}{\tau_{n}},\quad
a_{n}=\frac{\tau_{n}\sigma_{n+1}}{\sigma_{n}\tau_{n+1}},\quad\xi_{n+1,n}=-\frac{1}{2}\frac{\partial_{t_{1}}\tau_{n+1}}{\tau_{n+1}},$

and the non-commutative equation (3.8) is equivalent to444$D_{t}$ is the
Hirota’s bilinear derivative operator defined by $D_{t}f(t)\cdot
g(t)=\frac{\partial}{\partial s}f(t+s)g(t-s)|_{s=0}$.

$\displaystyle D_{t}\tau_{n+1}\cdot\tau_{n}=\sigma_{n}^{2},\quad
D_{t}^{2}\tau_{n+1}\cdot\tau_{n+1}=4\sigma_{n+1}\sigma_{n},$ (3.10)

which is the so-called C-Toda lattice whose $\tau$-function acts as that of
CKP hierarchy [31].

Below, we verify noncommutative variables in (3.8) admit quasi-determinant
forms. For this purpose, we observe that the first two equations in (3.8)
could be equivalently written as

$\displaystyle\partial_{t_{1}}H_{n+1}=a_{n}\cdot\partial_{t_{1}}H_{n}\cdot
a_{n}^{\top},\quad\partial_{t_{1}}\xi_{n+1,n}=a_{n}\cdot\partial_{t_{1}}H_{n}\cdot
H_{n}^{-1},$

and we need following lemmas.

######  Lemma 3.6.

It holds that

$\displaystyle\partial_{t_{1}}H_{n}=V_{n}V_{n}^{\top}.$ (3.11)

###### Proof.

In light of the condition on the derivative of moments as follows

$\displaystyle\partial_{t_{1}}m_{i,j}=\phi_{i}\phi_{j},\quad i,j=0,1,\dots,n,$

we can do computations on the derivative of quasideterminant following (A.27)
and

$\displaystyle\partial_{t_{1}}H_{n}=$
$\displaystyle\partial_{t_{1}}m_{n,n}+\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&\phi_{0}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&\phi_{n-1}\\\
m_{n,0}&\cdots&m_{n,n-1}&\boxed{0}\end{array}\right|~{}\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&m_{0,n}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&m_{n-1,n}\\\
\phi_{0}&\cdots&\phi_{n-1}&\boxed{0}\end{array}\right|$
$\displaystyle+\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&m_{0,n}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&m_{n-1,n}\\\
\partial_{t_{1}}m_{n,0}&\cdots&\partial_{t_{1}}m_{n,n-1}&\boxed{0}\end{array}\right|+\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&\partial_{t_{1}}m_{0,n}\\\
\vdots&&\vdots&\vdots\\\
m_{n-1,0}&\cdots&m_{n-1,n-1}&\partial_{t_{1}}m_{n-1,n}\\\
m_{n,0}&\cdots&m_{n,n-1}&\boxed{0}\end{array}\right|.$

Then together with the following quasi-determinant relations

$\displaystyle
V_{n}-\phi_{n}=\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&\phi_{0}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&\phi_{n-1}\\\
m_{n,0}&\cdots&m_{n,n-1}&\boxed{0}\end{array}\right|,\quad
V_{n}^{\top}-\phi_{n}=\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&m_{0,n}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&m_{n-1,n}\\\
\phi_{0}&\cdots&\phi_{n-1}&\boxed{0}\end{array}\right|.$

Moreover, by noting that

$\displaystyle
V_{n}\phi_{n}-\partial_{t_{1}}m_{n,n}=\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&\partial_{t_{1}}m_{0,n}\\\
\vdots&&\vdots&\vdots\\\
m_{n-1,0}&\cdots&m_{n-1,n-1}&\partial_{t_{1}}m_{n-1,n}\\\
m_{n,0}&\cdots&m_{n,n-1}&\boxed{0}\end{array}\right|$

and its transposed version, as well as that
$\partial_{t_{1}}m_{n,n}=\phi_{n}\phi_{n}$, we know the result is true. ∎

######  Lemma 3.7.

Regarding with the derivative of $\xi_{n+1,n}$, it holds that

$\displaystyle\partial_{t_{1}}\xi_{n+1,n}=-V_{n+1}V_{n}^{\top}H_{n}^{-1}.$
(3.12)

###### Proof.

Following the derivative formula (A.27), we obtain

$\displaystyle\\!\\!\partial_{t_{1}}\xi_{n+1,n}\\!\\!=\\!\\!\left|\begin{array}[]{cccc}\\!\\!m_{0,0}&\cdots&\\!\\!\\!\\!m_{0,n}&\\!\\!\\!\\!0\\\
\vdots&&\vdots&\vdots\\\
\\!\\!m_{n,0}&\cdots&\\!\\!\\!\\!m_{n,n}&\\!\\!\\!\\!\mathbb{I}_{p}\\\
\\!\\!\partial_{t_{1}}m_{n+1,0}&\\!\\!\\!\\!\cdots&\partial_{t_{1}}m_{n+1,n}&\\!\\!\\!\\!\boxed{0}\end{array}\right|+\left|\begin{array}[]{cccc}\\!\\!m_{0,0}&\cdots&\\!\\!\\!\\!m_{0,n}&\\!\\!\\!\\!\phi_{0}\\\
\vdots&&\vdots&\vdots\\\
\\!\\!m_{n,0}&\cdots&\\!\\!\\!\\!m_{n,n}&\\!\\!\\!\\!\phi_{n}\\\
\\!\\!m_{n+1,0}&\cdots&\\!\\!\\!\\!m_{n+1,n}&\\!\\!\\!\\!\boxed{0}\end{array}\right|\left|\begin{array}[]{cccc}\\!\\!m_{0,0}&\cdots&\\!\\!\\!\\!m_{0,n}&\\!\\!\\!\\!0\\\
\vdots&&\vdots&\vdots\\\
\\!\\!m_{n,0}&\cdots&\\!\\!\\!\\!m_{n,n}&\\!\\!\\!\\!\mathbb{I}_{p}\\\
\\!\\!\phi_{0}&\cdots&\\!\\!\\!\\!\phi_{n}&\\!\\!\\!\\!\boxed{0}\end{array}\right|.$

To simplify the above equation, we need to rewrite the first term as

$\displaystyle\left|\begin{array}[]{cccc}\\!\\!m_{0,0}&\cdots&\\!\\!\\!\\!m_{0,n}&\\!\\!\\!\\!0\\\
\vdots&&\vdots&\vdots\\\
\\!\\!m_{n,0}&\cdots&\\!\\!\\!\\!m_{n,n}&\\!\\!\\!\\!\mathbb{I}_{p}\\\
\\!\\!\partial_{t_{1}}m_{n+1,0}&\\!\\!\\!\\!\cdots&\partial_{t_{1}}m_{n+1,n}&\\!\\!\\!\\!\boxed{0}\end{array}\right|=\phi_{n+1}\left|\begin{array}[]{cccc}\\!\\!m_{0,0}&\cdots&\\!\\!\\!\\!m_{0,n}&\\!\\!\\!\\!0\\\
\vdots&&\vdots&\vdots\\\
\\!\\!m_{n,0}&\cdots&\\!\\!\\!\\!m_{n,n}&\\!\\!\\!\\!\mathbb{I}_{p}\\\
\\!\\!\phi_{0}&\\!\\!\\!\\!\cdots&\phi_{n}&\\!\\!\\!\\!\boxed{0}\end{array}\right|$

and the first part in second term as $V_{n+1}-\phi_{n+1}$. Therefore we have
the formula

$\displaystyle\\!\\!\partial_{t_{1}}\xi_{n+1,n}\\!\\!=V_{n+1}\left|\begin{array}[]{ccccc}\\!\\!m_{0,0}&\cdots&\\!\\!\\!\\!m_{0,n-1}&\\!\\!\\!\\!m_{0,n}&\\!\\!\\!\\!0\\\
\vdots&&\vdots&\vdots&\vdots\\\
\\!\\!m_{n-1,0}&\cdots&\\!\\!\\!\\!m_{n-1,n-1}&\\!\\!\\!\\!m_{n-1,n}&\\!\\!\\!\\!0\\\
\\!\\!m_{n,0}&\cdots&\\!\\!\\!\\!m_{n,n-1}&\\!\\!\\!\\!m_{n,n}&\\!\\!\\!\\!\mathbb{I}_{p}\\\
\\!\\!\phi_{0}&\cdots&\\!\\!\\!\\!\phi_{n-1}&\\!\\!\\!\\!\phi_{n}&\\!\\!\\!\\!\boxed{0}\end{array}\right|.$

According to the homological relation of quasideterminant (LABEL:homo), we
arrive at

$\displaystyle\left|\begin{array}[]{ccccc}\\!\\!m_{0,0}&\cdots&\\!\\!\\!\\!m_{0,n-1}&\\!\\!\\!\\!m_{0,n}&\\!\\!\\!\\!0\\\
\vdots&&\vdots&\vdots&\vdots\\\
\\!\\!m_{n-1,0}&\cdots&\\!\\!\\!\\!m_{n-1,n-1}&\\!\\!\\!\\!m_{n-1,n}&\\!\\!\\!\\!0\\\
\\!\\!m_{n,0}&\cdots&\\!\\!\\!\\!m_{n,n-1}&\\!\\!\\!\\!m_{n,n}&\\!\\!\\!\\!\mathbb{I}_{p}\\\
\\!\\!\phi_{0}&\cdots&\\!\\!\\!\\!\phi_{n-1}&\\!\\!\\!\\!\phi_{n}&\\!\\!\\!\\!\boxed{0}\end{array}\right|$
$\displaystyle\\!\\!=\\!\\!$
$\displaystyle-\left|\begin{array}[]{cccc}\\!\\!m_{0,0}&\cdots&\\!\\!\\!\\!m_{0,n-1}&\\!\\!\\!\\!m_{0,n}\\\
\vdots&&\vdots.&\vdots\\\
\\!\\!m_{n-1,0}&\cdots&\\!\\!\\!\\!m_{n-1,n-1}&\\!\\!\\!\\!m_{n-1,n}\\\
\\!\\!\phi_{0}&\cdots&\\!\\!\\!\\!\phi_{n-1}&\\!\\!\\!\\!\boxed{\phi_{n}}\end{array}\right|\left|\begin{array}[]{cccc}m_{0,0}&\cdots&m_{0,n-1}&m_{0,n}\\\
\vdots&&\vdots&\vdots\\\ m_{n-1,0}&\cdots&m_{n-1,n-1}&m_{n-1,n}\\\
m_{n,0}&\cdots&m_{n,n-1}&\boxed{m_{n,n}}\end{array}\right|^{-1}$
$\displaystyle=$ $\displaystyle-V_{n}^{\top}{H_{n}}^{-1},$

which means the equation (3.12) is satisfied. ∎

Therefore, by substituting equations (3.11) and (3.12) into (3.8), the
noncommutative C-Toda lattice could be directly verified.

At the end of this section, we give a discussion on the $t_{3}$-flow and
higher order flows could be similarly computed. By noting that

$\displaystyle\partial_{t_{3}}\langle f(x;\mathbf{t}),g(y;\mathbf{t})\rangle$
$\displaystyle=\langle\partial_{t_{3}}f(x;\mathbf{t}),g(y;\mathbf{t})\rangle+\langle
f(x;\mathbf{t}),\partial_{t_{3}}g(y;\mathbf{t})\rangle$
$\displaystyle\quad+\langle
x^{3}f(x;\mathbf{t}),g(y;\mathbf{t})\rangle+\langle
f(x;\mathbf{t}),y^{3}g(y;\mathbf{t})\rangle,$

where the last two terms could be alternatively written as

$\displaystyle\int_{\mathbb{R}_{+}}x^{2}f(x;\mathbf{t})W(x;\mathbf{t})dx\int_{\mathbb{R}_{+}}W(y;\mathbf{t})g(y;\mathbf{t})dy$
$\displaystyle-\int_{\mathbb{R}_{+}}xf(x;\mathbf{t})W(x;\mathbf{t})dx\int_{\mathbb{R}_{+}}yW(y;\mathbf{t})g(y;\mathbf{t})dy$
$\displaystyle+\int_{\mathbb{R}_{+}}f(x;\mathbf{t})W(x;\mathbf{t})\int_{\mathbb{R}_{+}}y^{2}W(y;\mathbf{t})g(y;\mathbf{t})dy.$

If we choose $f(x;\mathbf{t})$ as
$P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t})$, then we could take
$g(y;\mathbf{t})$ as a proper linear combination for
$\\{P_{m}^{\top}(y;\mathbf{t})\\}_{m\in\mathbb{N}}$ such that the summand is
zero. We have the following proposition.

######  Proposition 3.8.

Denote
$U_{m}^{\top}=\int_{\mathbb{R}_{+}}yW(y;\mathbf{t})P_{m}^{\top}(y;\mathbf{t})dy$,
and if $\eta_{m}$ and $\zeta_{m}$ are solutions of equations

$\displaystyle V_{m}^{\top}+V_{m-1}^{\top}\eta_{m}+V_{m-2}^{\top}\zeta_{m}=0,$
$\displaystyle U_{m}^{\top}+U_{m-1}^{\top}\eta_{m}+U_{m-2}^{\top}\zeta_{m}=0,$

then
$g(y;\mathbf{t})=P_{m}^{\top}(y;\mathbf{t})+P_{m-1}^{\top}(y;\mathbf{t})\eta_{m}+P_{m-2}^{\top}(y;\mathbf{t})\zeta_{m}$
satisfies

$\displaystyle\int_{\mathbb{R}_{+}}W(y;\mathbf{t})g(y;\mathbf{t})dy=\int_{\mathbb{R}_{+}}yW(y;\mathbf{t})g(y;\mathbf{t})dy=0.$

To solve the linear system, $\eta_{m}$ and $\zeta_{m}$ have quasideterminant
expressions

$\displaystyle\eta_{m}=\left|\begin{array}[]{ccc}V_{m-2}^{\top}&V_{m-1}^{\top}&V_{m}^{\top}\\\
U_{m-2}^{\top}&U_{m-1}^{\top}&U_{m}^{\top}\\\
0&\mathbb{I}_{p}&\boxed{0}\end{array}\right|,\quad\zeta_{m}=\left|\begin{array}[]{ccc}V_{m-2}^{\top}&V_{m-1}^{\top}&V_{m}^{\top}\\\
U_{m-2}^{\top}&U_{m-1}^{\top}&U_{m}^{\top}\\\
\mathbb{I}_{p}&0&\boxed{0}\end{array}\right|.$

Moreover, if we take $t_{3}$-derivative to equation

$\displaystyle\left\langle
P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t}),P^{\top}_{m}(y;\mathbf{t})+P_{m-1}^{\top}(y;\mathbf{t})\eta_{m}+P_{m-2}^{\top}(y;\mathbf{t})\zeta_{m}\right\rangle$
$\displaystyle\quad=H_{n+1}\left(\delta_{n+1,m}+\eta_{m}\delta_{n+1,m-1}+\zeta_{m}\delta_{n+1,m-2}\right)+a_{n}H_{n}\left(\delta_{n,m}+\eta_{m}\delta_{n,m-1}+\zeta_{m}\delta_{n,m-2}\right),$

then when $m<n$,

$\displaystyle\left\langle\partial_{t_{3}}(P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t})),P^{\top}_{m}(y;\mathbf{t})+P_{m-1}^{\top}(y;\mathbf{t})\eta_{m}+P_{m-2}^{\top}(y;\mathbf{t})\zeta_{m}\right\rangle=0,$
(3.13)

and when $m=n$,

$\displaystyle\left\langle\partial_{t_{3}}(P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t})),P^{\top}_{n}(y;\mathbf{t})+P_{n-1}^{\top}(y;\mathbf{t})\eta_{n}+P_{n-2}^{\top}(y;\mathbf{t})\zeta_{n}\right\rangle=\partial_{t_{3}}(a_{n}H_{n}).$
(3.14)

By writing

$\displaystyle\partial_{t_{3}}\left(P_{n+1}(x;\mathbf{t})+a_{n}P_{n}(x;\mathbf{t})\right)=\sum_{i=0}^{n}\gamma_{n,i}^{(3)}P_{i}(x;\mathbf{t}),$

above formulas (3.13) and (3.14) give an inhomogeneous linear system with
$(n+1)$ variables and $(n+1)$ equations. Therefore, an evolution of the wave
function with respect to $t_{3}$ could be obtained, and the compatibility
condition leads to a higher order noncommutative integrable system.

##  4\. Concluding remarks

In this paper, we generalized the Cauchy bi-orthogonal polynomials to the
matrix-valued version, and obtain a novel noncommutative integrable system. We
provide a Lax pair with noncommutative variables and make use of a quasi-
determinant technique to verify the solutions. However, there are some
interesting questions to be studied. One is to study the Pfaffian version. It
is noted that in [31], the solutions of C-Toda lattice (3.10) were verified by
Pfaffian identities. Therefore, it is expected to obtain a Pfaffian version to
verify these integrable systems. Another question is a noncommutative version
of fractional differential operator. In [10], a fractional differential
operator is an operator that can be represented in the form $A^{-1}B$ for two
differential operators $A$ and $B$. Therefore, the Lax equation (3.9) is valid
for a fractional differential operator with noncommutative variables. It is
interesting to develop a Sato-Krichever theory in noncommutative
circumstances.

## Appendix: A quick introduction to quasideterminants

Quasideterminants, introduced by Gelfand et al in the early 1990s [20], are
important and effective tools for noncommutative algebra. Here we give a quick
introduction to quasideterminants, including the basic definitions and some
identities that we used in previous sections. Consider an $n\times n$ matrix
$A$ over a noncommutative ring $\mathcal{R}$. In general, it has $n^{2}$
quasideterminants. One may calculate every quasideterminant $|A|_{i,j}$,
$1\leq i,j\leq n$, as follows

$\displaystyle|A|_{i,j}=a_{i,j}-r_{i}^{j}\left(A^{i,j}\right)^{-1}c_{j}^{i},$

where $A^{i,j}$, denoting the matrix obtained from $A$ by deleting the $i$th
row and $j$th column, is an invertible square matrix over $\mathcal{R}$;
$r_{i}^{j}$ is the row vector over $\mathcal{R}$ obtained form the $i$th row
of $A$ by deleting the $j$th entry; $c_{j}^{i}$ is the column vector over
$\mathcal{R}$ obtained from the $j$the column of $A$ by deleting the $i$th
entry. Here the noncommutative ring is supposed to be a matrix ring, with the
purpose to be in the frame of matrix-valued measure. Through reordering rows
and columns so that the "expansion point" is bottom right, a more graphic
notation for the quasideterminant by boxing the element $a_{i,j}$ can be
adopted as follows

$\displaystyle|A|_{i,j}$
$\displaystyle=\left|\begin{array}[]{ccccccc}a_{1,1}&\cdots&a_{1,j-1}&a_{1,j+1}&\cdots&a_{1,n}&a_{1,j}\\\
\vdots&&\vdots&\vdots&\vdots&&\vdots\\\
a_{i-1,1}&\cdots&a_{i-1,j-1}&a_{i-1,j+1}&\cdots&a_{i-1,n}&a_{i-1,j}\\\
a_{i+1,1}&\cdots&a_{i+1,j-1}&a_{i+1,j+1}&\cdots&a_{i+1,n}&a_{i+1,j}\\\
\vdots&&\vdots&\vdots&\vdots&&\vdots\\\
a_{n,1}&\cdots&a_{n,j-1}&a_{n,j+1}&\cdots&a_{n,n}&a_{n,j}\\\
a_{i,1}&\cdots&a_{i,j-1}&a_{i,j+1}&\cdots&a_{i,n}&\boxed{a_{i,j}}\\\
\end{array}\right|$
$\displaystyle=\left|\begin{array}[]{cc}A^{i,j}&c_{j}^{i}\\\
r_{i}^{j}&\boxed{a_{i,j}}\\\
\end{array}\right|=a_{i,j}-r_{i}^{j}\left(A^{i,j}\right)^{-1}c_{j}^{i}.$

One important application of quasideterminants is to solve the linear system
with noncommutative coefficients. Let $A=(a_{i,j})$ be an $n\times n$ matrix
over a ring $\mathcal{R}$. Assume that all the quasideterminants $|A|_{i,j}$
are well-defined and invertible. Then

$\displaystyle\left\\{\begin{array}[]{c}a_{1,1}x_{1}+\cdots+a_{1,n}x_{n}=\xi_{1}\\\
\vdots\\\ a_{n,1}x_{1}+\cdots+a_{n,n}x_{n}=\xi_{n}\end{array}\right.$

have solution $x_{i}\in\mathcal{R}$ if and only if

$\displaystyle x_{i}=\sum_{j=1}^{n}|A|_{j,i}^{-1}\xi_{j}.$

Below we show some important identities of quasideterminants [20, 21, 19],
which are useful in the verifications of solutions for noncommutative
integrable systems. One is the non-commutative Jacobi identity, which can be
written as

$\displaystyle\left|\begin{array}[]{ccc}A&B&C\\\ D&f&g\\\ E&h&\boxed{i}\\\
\end{array}\right|=\left|\begin{array}[]{cc}A&C\\\ E&\boxed{i}\\\
\end{array}\right|-\left|\begin{array}[]{cc}A&B\\\ E&\boxed{h}\\\
\end{array}\right|\left|\begin{array}[]{cc}A&B\\\ D&\boxed{f}\\\
\end{array}\right|^{-1}\left|\begin{array}[]{cc}A&C\\\ D&\boxed{g}\\\
\end{array}\right|.$ (A.12)

where $A$ should be invertible. By utilizing the noncommutative Jacobi
identity, one could obtain the so-called homological relation

$\displaystyle\begin{aligned} &\left|\begin{array}[]{ccc}A&B&C\\\ D&f&g\\\
E&\boxed{h}&i\end{array}\right|=\left|\begin{array}[]{ccc}A&B&C\\\ D&f&g\\\
E&h&\boxed{i}\end{array}\right|\left|\begin{array}[]{ccc}A&B&C\\\ D&f&g\\\
0&\boxed{0}&1\end{array}\right|,\\\ &\left|\begin{array}[]{ccc}A&B&C\\\
D&f&\boxed{g}\\\ E&{h}&i\end{array}\right|=\left|\begin{array}[]{ccc}A&B&0\\\
D&f&\boxed{0}\\\ E&h&1\end{array}\right|\left|\begin{array}[]{ccc}A&B&C\\\
D&f&g\\\ E&h&\boxed{i}\end{array}\right|.\end{aligned}$ (A.13)

Besides, it is necessary to consider the derivative formulae for
quasideterminants. Considering the expansion of quasideterminant

$\displaystyle\left|\begin{array}[]{cc}A&B\\\
C&\boxed{d}\end{array}\right|=d-CA^{-1}B,$

and assume that $A,\,B,C$ and $d$ are $t$-dependent, then we have the general
derivative formula

$\displaystyle\left|\begin{array}[]{cc}A&B\\\
C&\boxed{d}\end{array}\right|^{\prime}=d^{\prime}-C^{\prime}A^{-1}B-CA^{-1}B^{\prime}+CA^{-1}A^{\prime}A^{-1}B,$
(A.16)

where ′ means the derivative with respect to $t$. Especially, if the
quasideterminant is well structured, then derivative formulae will be
modified. One is the noncommutative version of Grammian determinant. If $A$ is
a Grammian-like matrix satisfying $A^{\prime}=\sum_{i=1}^{k}E_{i}F_{i}$ where
$E_{i}$ (resp. $F_{i}$) is a column (resp. row) vector of appropriate length.
Then we can factorise the RHS of (A.16) to obtain

$\displaystyle\left|\begin{array}[]{cc}A&B\\\
C&\boxed{d}\end{array}\right|^{\prime}=d^{\prime}+\left|\begin{array}[]{cc}A&B\\\
C^{\prime}&\boxed{0}\end{array}\right|+\left|\begin{array}[]{cc}A&B^{\prime}\\\
C&\boxed{0}\end{array}\right|+\sum_{i=1}^{k}\left|\begin{array}[]{cc}A&E_{i}\\\
C&\boxed{0}\end{array}\right|\left|\begin{array}[]{cc}A&B\\\
F_{i}&\boxed{0}\end{array}\right|.$ (A.27)

Another is the noncommutative version of the Wronskian determinant. By
inserting the identity matrix expressed as $\sum_{j=1}^{n}e_{j}^{\top}e_{j}$,
where $e_{j}=(0,\cdots,\mathbb{I}_{p},\cdots,0)$ is the block unit vector,
then the above derivative for quasideterminants can be written as

$\displaystyle\left|\begin{array}[]{cc}A&B\\\
C&\boxed{d}\end{array}\right|^{\prime}=\left|\begin{array}[]{cc}A&B\\\
C^{\prime}&\boxed{d^{\prime}}\end{array}\right|+\sum_{j=1}^{n}\left|\begin{array}[]{cc}A&e_{j}^{\top}\\\
C&\boxed{0}\end{array}\right|\left|\begin{array}[]{cc}A&B\\\
\left(A^{j}\right)^{\prime}&\boxed{\left(B^{j}\right)^{\prime}}\end{array}\right|$
(A.36)

and

$\displaystyle\left|\begin{array}[]{cc}A&B\\\
C&\boxed{d}\end{array}\right|^{\prime}=\left|\begin{array}[]{cc}A&B^{\prime}\\\
C&\boxed{d^{\prime}}\end{array}\right|+\sum_{j=1}^{n}\left|\begin{array}[]{cc}A&\left(A_{j}\right)^{\prime}\\\
C&\boxed{\left(C_{j}\right)^{\prime}}\end{array}\right|\left|\begin{array}[]{cc}A&B\\\
e_{j}&\boxed{0}\end{array}\right|$ (A.45)

where $A^{j}=e_{j}A$ and $A_{j}=Ae_{j}^{\top}$ denote the $j$th row and column
of matrix $A$ respectively.

## Acknowledgement

This work is partially funded by grants (NSFC12101432, NSFC12175155,
NSFC11971322).

## References

  * [1] M. Adler and P. van Moerbeke.  Recursion relations for unitary integrals, combinatorics and the Toeplitz lattice.  Commun. Math. Phys., 237 (2003), 397-440.
  * [2] M. Adler, P. van Moerbeke and P. Vanhaecke.  Moment matrices and multi-component KP, with applications to random matrix theory.  Commun. Math. Phys., 286 (2009), 1-38.
  * [3] C. Álvarez-Fernández, G. Ariznabarreta, J. Garćia-Ardila, M. Mañas and F. Marcellán.  Christoffel transformations for matrix orthogonal polynomials in the real line and the non-Abelian 2D Toda lattice hierarchy.  _Int. Math. Res. Not._ , (2017), 1285-1341.
  * [4] M. Bertola.  Abelianization of matrix orthogonal polynomials.  arXiv: 2107.12998.
  * [5] M. Bertola, M. Gekhtman and J. Szmigielski.  The Cauchy two-matrix model.  Commun. Math. Phys., 287 (2009), 983-1014.
  * [6] M. Bertola, M. Gekhtman and J. Szmigielski.  Cauchy biorthogonal polynomials.  J. Approx. Theory, 162 (2010), 832-867.
  * [7] M. Bertola, M. Cafasso and V. Rubtsov.  Noncommutative Painlevé equations and systems of Calogero type.  Commun. Math. Phys., 363 (2018), 503-530.
  * [8] A. Branquinho, A. Foulquié-Moreno and J. Garciá-Ardila.  Matrix Toda and Volterra lattices.  _Appl. Math. Comp._ , 365 (2020), 124722.
  * [9] M. Cafasso and D. Manuel.  Non-commutative Painlevé equations and Hermite-type matrix orthogonal polynomials.  _Commun. Math. Phys._ , 326 (2014), 559-583.
  * [10] W. Casper, E. Horozov, P. Iliev and M. Yakimov.  A Sato-Krichever theory for fractional differential operators.  arXiv: 2108.12010.
  * [11] X. Chang, Y. He, X. Hu and S. Li. Partial-skew-orthogonal polynomials and related integrable lattices with Pfaffian tau-functions.  Commun. Math. Phys., 364 (2018), 1069-1119.
  * [12] P. Di Francesco and R. Kedem.  Non-commutative integrability, paths and quasi-determinants.  _Adv. Math._ , 10 (2011), 97-152.
  * [13] A. Dimakis and F. Müller-Hoissen.  An algebraic scheme associated with the non-commutative KP hierarchy and some of its extensions.  J. Phys. A, 38 (2005), 5453-5505.
  * [14] A. Doliwa.  Non-commutative Hermite–Padé approximation and integrability.  Lett. Math. Phys., 112 (2022), article no. 68.
  * [15] A. Doliwa and A. Siemaszko.  Integrability and geometry of the Wynn recurrence.  Numer. Algor., 2022, DOI: 10.1007/s11075-022-01344-5.
  * [16] M. Duits and A. Kuijlaars.  The two periodic Aztec diamond and matrix valued orthogonal polynomials.  to appear in _JEMS_ , arXiv: 1712.05636.
  * [17] P. Etingof, I. Gelfand, and V. Retakh.  Factorization of differential operators, quasideterminants, and nonabelian Toda field equations.  Math. Res. Lett., 4 (1997), 413-425.
  * [18] P. Etingof, I. Gelfand, and V. Retakh.  Nonabelian integrable systems, quasideterminants, and Marchenko lemma.  Math. Res. Lett., 5 (1998), 1-12.
  * [19] I. Gelfand, S. Gelfand, V. Retakh and R. Wilson.  quasideterminants.  _Adv. Math._ , 193 (2005), 56-141.
  * [20] I. Gelfand and V. Retakh.  Determinants of matrices over commutative rings.  Funct. Anal. Appl., 25 (1991), 91-102.
  * [21] I. Gelfand and V. Retakh.  A theory of noncommutative determinants and characteristic functions of graphs.  Funct. Anal. Appl., 26 (1992), 1-20.
  * [22] C. Gilson and J. Nimmo.  On a direct approach to quasideterminant solutions of a noncommutative equation.  J. Phys. A, 40 (2007), 3839-3850.
  * [23] C. Gilson, J. Nimmo and Y. Ohta.  Quasideterminant solutions of a non-Abelian Hirota-Miwa equation.  J. Phys. A, 40 (2007), 12607-12617.
  * [24] F. Grünbaum and M. de la Iglesia.  Matrix valued orthogonal polynomials arising from group representation theory and a family of quasi-birth-and-death processes.  _SIAM J. Mat. Anal. Appl._ , 30 (2008), 741-761.
  * [25] R. Hirota.  The direct method in soliton theory.  Edited and translated by A. Nagai, J. Nimmo and C. Gilson, Cambridge Tracts in Mathematics 155, Cambridge University Press, 2004.
  * [26] M. Ismail, E. Keolink and P. Román.  Matrix valued Hermite polynomial, Burchnall formulas and non-abelian Toda lattice.  _Adv. Appl. Math.,_ 110 (2019), 235-269.
  * [27] M. Jimbo and T. Miwa.  Solitons and infinite dimensional Lie algebras.  Publ. RIMS, Kyoto Univ., 19 (1983), 943-1001.
  * [28] M. Krein.  Infinite J-matrices and a matrix-moment problem.  Dokl. Akad. Nauk. SSSR, 69 (1949), 125-128.
  * [29] I. Krichever.  The periodic non-abelian Toda chain and its two-dimensional generalization.  Russ. Math. Surv., 35 (1980), 53-79.
  * [30] C. Li, J. Nimmo and K. Tamizhmani.  On solutions to the non-abelian Hirota-Miwa equation and its continuum limits.  Proc. R. Soc. A, 465 (2009), 1441-1451.
  * [31] C. Li and S. Li.  The Cauchy two-matrix model, C-Toda lattice and CKP hierarchy.  _J. Nonlinear Sci._ , 29 (2019), 3-27.
  * [32] S. Li.  Matrix orthogonal polynomials, non-abelian Toda lattice and Bäcklund transformation.  arXiv: 2109.00671.
  * [33] T. Miwa, M. Jimbo, and E. Date.  Solitons: Differential equations, symmetries and infinite dimensional algebras.  Cambridge Tracts in Mathematics 135, Cambridge University Press, 1999.
  * [34] H. Lundmark and J. Szmigielski.
  * [35] V. Retakh and V. Rubtsov.  Noncommutative Toda chains, Hankel quasideterminants and the Painlevé II equation.  J. Phys. A, 43 (2010), 505204.