# Applications of multiple orthogonal polynomials with hypergeometric moment
generating functions

Thomas Wolfs111Department of Mathematics, KU Leuven, Belgium. E-mail adress:
<EMAIL_ADDRESS>
###### Abstract

We investigate several families of multiple orthogonal polynomials associated
with weights for which the moment generating functions are hypergeometric
series with slightly varying parameters. The weights are supported on the unit
interval, the positive real line, or the unit circle and the multiple
orthogonal polynomials are generalizations of the Jacobi, Laguerre or Bessel
orthogonal polynomials. We give explicit formulas for the type I and type II
multiple orthogonal polynomials and study some of their properties. In
particular, we describe the asymptotic distribution of the (scaled) zeros of
the type II multiple orthogonal polynomials via the free convolution.
Essential to our overall approach is the use of the Mellin transform. Finally,
we discuss two applications. First, we show that the multiple orthogonal
polynomials appear naturally in the study of the squared singular values of
(mixed) products of truncated unitary random matrices and Ginibre matrices.
Secondly, we use the multiple orthogonal polynomials to simultaneously
approximate certain hypergeometric series and to provide an explicit proof of
their $\mathbb{Q}$-linear independence.

Keywords: multiple orthogonal polynomials, hypergeometric series, Mellin
transform, random matrix theory, Diophantine approximation

##  1 Introduction

Multiple orthogonal polynomials, compared to (regular) orthogonal polynomials,
satisfy orthogonality conditions with respect to several weights
$(w_{1},\dots,w_{r})$, instead of just one weight (see [12, Chapter 23] and
[31, Chapter 4] for an introduction). They depend on a multi-index
$\vec{n}\in\mathbb{Z}_{\geq 0}^{r}$, of size
$\left|\vec{n}\right|=n_{1}+\dots+n_{r}$, that determines the way in which the
orthogonality conditions are distributed. Suppose that the weights are
supported on $\Lambda\subset\mathbb{C}$ and that their moments are all finite.
The type I multiple orthogonal polynomials are given by vectors of polynomials
$(A_{\vec{n},1},\dots,A_{\vec{n},r})$, with $\deg A_{\vec{n},j}\leq n_{j}-1$,
for which the type I function $F_{\vec{n}}=\sum_{j=1}^{r}A_{\vec{n},j}w_{j}$
satisfies the orthogonality conditions

$\int_{\Lambda}F_{\vec{n}}(x)x^{k}dx=0,\quad
k=0,\dots,\left|\vec{n}\right|-2.$

The type II multiple orthogonal polynomials are polynomials $P_{\vec{n}}$ of
degree at most $\left|\vec{n}\right|$ that satisfy the orthogonality
conditions

$\int_{\Lambda}P_{\vec{n}}(x)x^{k}w_{j}(x)dx=0,\quad k=0,\dots,n_{j}-1,\quad
j=1,\dots,r.$

The existence of non-zero type I functions and type II polynomials is always
guaranteed. The associated multi-index is called normal if the type I and type
II polynomials attain their maximal degree. In that case, the type I functions
and type II polynomials are determined uniquely up to a scalar multiplication.
Systems of weights in which all multi-indices are normal are called perfect
[28]. Many examples of perfect systems (see, e.g., [42]) arise as so-called
Angelesco systems [2] or AT-systems [31, Section 5.7], or more specifically as
Nikishin systems [9, 10].

We will investigate families of multiple orthogonal polynomials for which the
moment generating functions of the weights are (scaled) hypergeometric series
of the form

$f_{0}(z;\vec{a},\vec{b})=\sum_{k\geq
0}\frac{\Gamma(k+\vec{a}+1)}{\Gamma(k+\vec{b}+1)}z^{k},$ (1)

in which we use parameters $\vec{a}\in(-1,\infty)^{p}$ and
$\vec{b}\in(-1,\infty)^{q}$. To simplify notation, we used the convention that
applying the gamma function $\Gamma(s)=\int_{0}^{\infty}e^{-x}x^{s-1}dx$ to a
vector $\vec{s}\in(0,\infty)^{r}$ yields the product
$\prod_{j=1}^{r}\Gamma(s_{j})$ of the values of the gamma function applied to
its components. This convention extends to Pochhammer symbols
$(s)_{k}=\prod_{j=0}^{k-1}(s+j)$ in the natural way:
$(\vec{s})_{k}=\prod_{j=1}^{r}(s_{j})_{k}$. The weight
$w_{0}(x;\vec{a},\vec{b})$ associated with the moment generating function
$f_{0}(z;\vec{a},\vec{b})$ in (1) has moments

$m_{k}(\vec{a},\vec{b})=\int_{\Lambda}x^{k}w_{0}(x;\vec{a},\vec{b})dx=\frac{\Gamma(k+\vec{a}+1)}{\Gamma(k+\vec{b}+1)}$
(2)

and its support $\Lambda$ will be the unit interval $[0,1]$, if $p=q$, the
positive real line $[0,\infty)$, if $p>q$, or the unit circle
$\\{z\in\mathbb{C}:|z|\,=1\\}$, if $p<q$. It will be natural to consider
$r=\max\\{p,q\\}$ weights $w_{j}(x;\vec{a},\vec{b})$, $1\leq j\leq r$, that
arise from $w_{0}(x;\vec{a},\vec{b})$ by slightly varying its parameters.

In the most simple cases, where $r=1$, the weight $w_{0}(x;\vec{a},\vec{b})$
will reduce to one of the three classical orthogonality weights. For an
introduction to orthogonal polynomials, we refer to [41]. Taking
$\vec{a}=(a_{1})$ and $\vec{b}=(b_{1})$, so $(p,q)=(1,1)$, leads to the beta
density $\mathcal{B}^{a_{1},b_{1}}(x)=x^{a_{1}}(1-x)^{b_{1}-a_{1}-1}$ on
$[0,1]$ (assuming $a_{1}<b_{1}$). The associated orthogonal polynomials are
the Jacobi polynomials
$P_{n}^{(a_{1},b_{1}-a_{1}-1)}(x)={}_{2}F_{1}(-n,n+b_{1};a_{1}+1;x)$. The
gamma density $\mathcal{G}^{a_{1}}(x)=x^{a_{1}}e^{-x}$ on $[0,\infty)$
corresponds to the case $(p,q)=(1,0)$ with $\vec{a}=(a_{1})$ and $\vec{b}=()$.
In this case, the orthogonal polynomials are the Laguerre polynomials
$L_{n}^{(a_{1})}(x)={}_{1}F_{1}(-n;a_{1}+1;x)$. Finally, if $(p,q)=(0,1)$, we
need to move to the complex plane and have to consider the weight
$f_{0}(1/z;-,b_{1})/z$ on $\\{z\in\mathbb{C}:|z|\,=1\\}$. The associated
orthogonal polynomials are the Bessel polynomials
$B_{n}^{(b_{1}-1)}(z)={}_{2}F_{0}(-n,n+b_{1};-;z)$. The systems of weights
that we will consider can therefore be seen as generalizations of the
classical settings of Jacobi (Section 2), Laguerre (Section 3) and Bessel
(Section 4).

Some special cases of the multiple orthogonal polynomials have been studied
before; the Jacobi-Piñeiro polynomials, orthogonal with respect to
$x^{\alpha_{j}}(1-x)^{\beta}$ for $1\leq j\leq r$, with $\beta=0$ in [38]
correspond to a subcase of $(p,q)=(r,r)$, the multiple orthogonal polynomials
associated with Meijer-$G$ functions in [22] correspond to $(p,q)=(r,0)$, with
$K$-Bessel functions in [44] to $(p,q)=(2,0)$, with confluent hypergeometric
functions in [25] to $(p,q)=(2,1)$, with the exponential integral in [43] to a
subcase of $(p,q)=(2,1)$ and with Gauss’ hypergeometric functions in [26] to
$(p,q)=(2,2)$. The general $(p,q)$-case was also considered recently in [24]
and [39], where it appeared in connection with the theory of branched
continued fractions. In the former, explicit formulas for the type II
polynomials on the step-line were found and some properties of the zeros were
described.

Essential to our approach, compared to what was done before, will be the use
of the Mellin transform. The Mellin transform of a function $f\in
L^{1}(\mathbb{R}_{\geq 0})$ is given by

$\hat{f}(s)=\int_{0}^{\infty}f(x)x^{s-1}dx.$

Two important properties of the Mellin transform are that there is an inverse
transform

$f(x)=\frac{1}{2\pi i}\int_{c-i\infty}^{c+i\infty}\hat{f}(s)x^{-s}ds$

and that the Mellin convolution of two functions $f,g\in
L^{1}(\mathbb{R}_{\geq 0})$, which is given by

$(f\ast g)(x)=\int_{0}^{\infty}f(t)g(x/t)\frac{dt}{t},$

has Mellin transform

$\int_{0}^{\infty}(f\ast
g)(x)x^{s-1}dx=\left(\int_{0}^{\infty}f(x)x^{s-1}dx\right)\cdot\left(\int_{0}^{\infty}g(x)x^{s-1}dx\right).$

Many results about the type I functions and type II polynomials can be
obtained, very naturally, via the Mellin transform, and for a larger set of
multi-indices than the usual step-line.

The multiple orthogonal polynomials in the Laguerre- and Jacobi-like setting
(i.e. $p\geq q$) will appear in the study of the squared singular values of
products of truncated Haar distributed unitary random matrices and Ginibre
matrices (matrices with independent standard complex Gaussian random variables
as its entries). Such products were for example already investigated in [20,
21]. The connection with multiple orthogonality was only established in the
case where one considers a product of Ginibre matrices, see [22], or the
product of a Ginibre matrix with a truncated unitary random matrix, see [43].
In those cases, the associated weights were respectively Meijer-$G$ functions
generalizing the $K$-Bessel functions and exponential integrals. The role of
multiple orthogonality for other kinds of products has, to our knowledge, not
been investigated before. Other families of multiple orthogonal polynomials
also appear in random matrix theory, see [17, 18].

The multiple orthogonal polynomials in the Bessel- and Jacobi-like setting
(i.e. $p\leq q$) will appear as the solutions of certain Hermite-Padé
approximation problems associated with hypergeometric series. Similar
approximation problems were considered in [32] in connection with the theory
of branched continued fractions, but the overall quality of the approximants
was not studied. We will show that the Diophantine and approximation quality
will be good enough to prove $\mathbb{Q}$-linear independence of the
approximated hypergeometric series and in particular their irrationality. Such
results are also implied by Siegel’s method for $E$-functions (see, e.g., [36]
for an introduction), which was further improved in [37, 7, 45]. However,
since Siegel’s method is based on approximants that are only known implicitly,
they raise the question of whether such results can also be proven by making
use of approximants that can be constructed explicitly. Thus, by giving a
proof using certain explicit Hermite-Padé approximants, we will provide an
answer to this question. Hermite-Padé approximation was also used in Apéry’s
proof of the irrationality of $\zeta(3)$ in [3] as explained in [6].

##  2 Jacobi-like setting

###  2.1 System of weights

Suppose that $\vec{a},\vec{b}\in(-1,\infty)^{r}$ satisfy $a_{j}<b_{j}$ for
$1\leq j\leq r$. We will consider the $r$ weights
$w_{j}(x;\vec{a},\vec{b})=w_{0}(x;\vec{a},\vec{b}+\vec{e}_{j})$, $1\leq j\leq
r$, on $(0,1)$ (here $\vec{e}_{j}$ denotes the unit vector with a $1$ in the
$j$-th component). In this setting, a weight $w_{0}(x;\vec{a},\vec{b})$ with
moments $m_{k}(\vec{a},\vec{b})$ as in (2) can be constructed by taking the
Mellin convolution of the $r$ beta densities $\mathcal{B}^{a_{j},b_{j}}$,
$1\leq j\leq r$. Indeed, a beta density $\mathcal{B}^{a,b}$ has Mellin
transform $\Gamma(s+a)/\Gamma(s+b)$. Since Carleman’s condition
$\sum_{k=0}^{\infty}m_{k}(\vec{a},\vec{b})^{-1/(2k)}=\infty$ is satisfied,
there can be no other (positive) weights on the positive real line with those
moments. The weight $w_{j}(x;\vec{a},\vec{b})$ then has Mellin transform

$\int_{0}^{1}x^{s-1}w_{j}(x;\vec{a},\vec{b})dx=\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+\vec{e}_{j})}$

and arises as the Mellin convolution of $w_{0}(x;\vec{a},\vec{b})$ with the
Jacobi-Piñeiro weight $x^{b_{j}}$. Note that the confluent case
$\vec{a}=\vec{b}$ reduces to the Jacobi-Piñeiro setting in which we use the
weights $x^{b_{j}}$ for $1\leq j\leq r$.

We will study the multiple orthogonal polynomials associated with the weights
$(w_{1},\dots,w_{r})$ for multi-indices near the diagonal. This is a larger
set of multi-indices than the usual step-line. The precise definition of these
sets of multi-indices is stated below.

######  Definition 2.1.

A multi-index $\vec{n}\in\mathbb{N}^{r}$ is said to be near the diagonal if
$\left|n_{i}-n_{j}\right|\leq 1$ for all $i,j$. On the other hand, we say that
a multi-index $\vec{n}\in\mathbb{N}^{r}$ is on the step-line if
$n_{1}\geq\dots\geq n_{r}\geq n_{1}-1$. The sets of all such multi-indices are
denoted by $\mathcal{N}^{r}$ and $\mathcal{S}^{r}$ respectively.

First, we will show that under certain conditions on the parameters the set of
weights $\\{w_{j}(x)\\}_{j=1}^{r}$ is an AT-system for multi-indices near the
diagonal. For this we need to prove that the (extended) set of weights
$\\{x^{k}w_{j}(x)\mid 0\leq k\leq n_{j}-1,\,1\leq j\leq r\\}$ is a T-system,
i.e. that every combination $F(x)=\sum_{j=1}^{r}A_{j}(x)w_{j}(x)$, with $\deg
A_{j}\leq n_{j}-1$, has at most $\left|\vec{n}\right|-1$ zeros in $(0,1)$. We
will prove this by analyzing the Mellin transform of such polynomial
combinations. As a consequence, multi-indices near the diagonal will be normal
and some properties of the zeros of the corresponding multiple orthogonal
polynomials will be known.

######  Lemma 2.2.

Let $\vec{n}\in\mathcal{N}^{r}$. The Mellin transform of
$F(x)=\sum_{j=1}^{r}A_{j}(x)w_{j}(x)$, with $\deg A_{j}\leq n_{j}-1$, is of
the form

$\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+\vec{n})}p_{\left|\vec{n}\right|-1}(s),$

where $p_{\left|\vec{n}\right|-1}(s)$ is a polynomial of degree at most
$\left|\vec{n}\right|-1$.

###### Proof.

The Mellin transform of $F(x)$ is given by

$\int_{0}^{1}F(x)x^{s-1}dx=\sum_{j=1}^{r}\sum_{k=0}^{n_{j}-1}A_{j}[k]\int_{0}^{1}x^{k+s-1}w_{j}(x)dx,$

where

$\int_{0}^{1}x^{k+s-1}w_{j}(x)dx=\frac{\Gamma(s+k+\vec{a})}{\Gamma(s+k+\vec{b}+\vec{e}_{j})}=\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b})}\frac{(s+\vec{a})_{k}}{(s+\vec{b})_{k}(s+b_{j}+k)}.$

It then remains to show that

$\sum_{j=1}^{r}\sum_{k=0}^{n_{j}-1}A_{j}[k]\frac{(s+\vec{a})_{k}}{(s+\vec{b})_{k}(s+b_{j}+k)}=\frac{p_{\left|\vec{n}\right|-1}(s)}{\prod_{i=1}^{r}(s+b_{i})_{n_{i}}}.$

If we multiply the left-hand side by $\prod_{i=1}^{r}(s+b_{i})_{n_{i}}$, we
obtain

$\sum_{j=1}^{r}\sum_{k=0}^{n_{j}-1}A_{j}[k](s+\vec{a})_{k}\prod_{i=1}^{r}\frac{(s+b_{i})_{n_{i}}}{(s+b_{i})_{k+\delta_{i,j}}},$
(3)

which is a polynomial if $n_{i}\geq n_{j}-1+\delta_{i,j}$ for all $i,j$ (here
$\delta_{i,j}$ denotes the Kronecker delta). The latter leads precisely to the
near diagonal condition on $\vec{n}$ that was assumed. The polynomial
$p_{\left|\vec{n}\right|-1}(s)$ is of degree at most $\left|\vec{n}\right|-1$,
because, after comparing with the left-hand side, the ratio
$p_{\left|\vec{n}\right|-1}(s)/\prod_{i=1}^{r}(s+b_{i})_{n_{i}}$ needs to be
$O(s^{-1})$ as $s\to\infty$. ∎

In what follows, we will use the fact the underlying system of Jacobi-Piñeiro
weights $(x^{b_{1}},\dots,x^{b_{r}})$ already forms an AT-system for every
multi-index under the condition that $b_{i}-b_{j}\not\in\mathbb{Z}$ whenever
$i\neq j$ (see [12, Section 23.3.2]). Given an appropriate kernel $K(x,y)$, a
T-system $\\{u_{j}(x)\\}_{j=1}^{n}$ can be lifted to another T-system
$\\{v_{j}(y)\\}_{j=1}^{n}$ defined by
$v_{j}(y)=\int_{\mathbb{R}}u_{j}(x)K(x,y)dx$ by making use of the composition
formula

$\left|v_{j}(y_{k})\right|_{j,k=1}^{n}=\int_{\vec{x}\in
X^{n}(\uparrow)}|u_{j}(x_{l})|_{j,l=1}^{n}|K(x_{l},y_{k})|_{l,k=1}^{n}dx_{1}\dots
dx_{n},$ (4)

where $X^{n}(\uparrow)=\\{\vec{x}\in X^{n}\mid x_{1}<\dots<x_{n}\\}$ (see [14,
Example 8 - Chapter 1]; it is a special case of Andréief’s integration formula
[1]). Indeed, the first determinant on the right-hand side is non-zero and has
a constant sign for every $\vec{x}\in X^{n}(\uparrow)$ because
$\\{u_{j}(x)\\}_{j=1}^{n}$ is a T-system. Hence, if
$|K(x_{l},y_{k})|_{l,k=1}^{n}$ doesn’t change sign for all $\vec{x}\in
X^{n}(\uparrow)$ and $\vec{y}\in Y^{n}(\uparrow)$ and
$|K(x_{k},y_{l})|_{k,l=1}^{n}\neq 0$ on a non-null subset of $X^{n}(\uparrow)$
for fixed $\vec{y}\in Y^{n}(\uparrow)$, the left-hand side shares the same
property and the new system is also a T-system. Using this idea, we can prove
the following extension to more general weights.

######  Proposition 2.3.

Suppose that all $b_{j}-a_{j}\in\mathbb{Z}_{\geq
0}\cup(\left|\vec{n}\right|-1,\infty)$. Then the set of weights
$\\{w_{j}(x)\\}_{j=1}^{r}$ is an AT-system on $(0,1)$ for every
$\vec{n}\in\mathcal{N}^{r}$, if additionally $b_{i}-b_{j}\not\in\mathbb{Z}$
whenever $i\neq j$ or $b_{J}-a_{J}>\left|\vec{n}\right|-n_{J}$ for some $J$.

###### Proof.

Let $F(x)=\sum_{j=1}^{r}A_{j}(x)w_{j}(x)$ be a polynomial combination of the
weights with $\deg A_{j}\leq n_{j}-1$. It follows from Lemma 2.2 that the
Mellin transform of $F(x)$ is of the form

$\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+\vec{n})}p_{\left|\vec{n}\right|-1}(s).$
(5)

We will interpret this result in two ways, each leading to a different set of
conditions on the parameters.

First, since (5) is given by the product of
$p_{\left|\vec{n}\right|-1}(s)/\prod_{i=1}^{r}(s+b_{i})_{n_{i}}$ and
$\Gamma(s+\vec{a})/\Gamma(s+\vec{b})$, it is given by the Mellin convolution
of a polynomial combination $G(x)=\sum_{j=1}^{r}B_{j}(x)x^{b_{j}}$ with $\deg
B_{j}\leq n_{j}-1$ (see Lemma 2.2 with $\vec{a}=\vec{b}$) and $w_{0}(x)$
(which arose as the Mellin convolution of the beta densities
$\mathcal{B}^{a_{j},b_{j}}(x)$ for $1\leq j\leq r$). Hence, $F(x)$ can also be
written as a linear combination of the weights in $\\{x^{k+b_{j}}\ast
w_{0}(x)\mid 0\leq k\leq n_{j}-1,\,1\leq j\leq r\\}$. Since we know that the
set of weights $\\{x^{k+b_{j}}\mid 0\leq k\leq n_{j}-1,\,1\leq j\leq r\\}$ is
a T-system if $b_{i}-b_{j}\not\in\mathbb{Z}$ whenever $i\neq j$, it is
sufficient to show that taking the consecutive Mellin convolution with the
beta densities $\mathcal{B}^{a_{j},b_{j}}(x)$ preserves this property.
Convoluting with a beta density $\mathcal{B}^{a,b}(x)$ corresponds to taking
$K(x,y)=(x/y)^{a}(1-x/y)_{+}^{b-a-1}/y$ on $(0,1)^{2}$ in the composition
formula (4). It is known that $(y-x)_{+}^{b-a-1}$ is totally positive kernel
of order $n=\left|\vec{n}\right|$ if $b-a-1\in\mathbb{Z}_{\geq
0}\cup(n-2,\infty)$, and thus in particular
$\left|(y_{j}-x_{i})_{+}^{b-a-1}\right|_{i,j=1}^{n}\geq 0$ for all $\vec{x}\in
X^{n}(\uparrow)$ and $\vec{y}\in Y^{n}(\uparrow)$ (see [13, Theorem 9.1] and
[14] for an introduction to totally positive kernels and related topics). On
the other hand, it is straightforward to show that
$\left|(y_{j}-x_{i})_{+}^{b-a-1}\right|_{i,j=1}^{n}=(y_{1}-x_{1})^{b-a-1}\dots(y_{n}-x_{n})^{b-a-1}>0$
if $x_{1}<y_{1}\leq x_{2}<y_{2}\leq\dots\leq x_{n}<y_{n}$. Hence the resulting
set of weights will again be a T-system.

Another way to interpret the polynomial combination $F(x)$ is to start from a
function with Mellin transform
$\Gamma(s+a_{J})/\Gamma(s+b_{J}+n_{J})p_{\left|\vec{n}\right|-1}(s)$ and to
take the convolution with the missing beta densities
$\mathcal{B}^{a_{j},b_{j}+n_{j}}(x)$. Such a function exists under the
assumption that $b_{J}+n_{J}-a_{J}>\left|\vec{n}\right|$. Indeed, in that case
we can do a partial fraction decomposition and take the inverse Mellin
transform to show that it is given by a linear combination of the weights in
$\\{x^{k+a_{J}}(1-x)^{b_{J}+n_{J}-a_{J}-\left|\vec{n}\right|-1}\mid 0\leq
k\leq\left|\vec{n}\right|-1\\}$. The latter is a T-system as implied by the
Jacobi(-Piñeiro) case and we can proceed similarly as before. ∎

The bottleneck of the above approach is that in order for $(y-x)_{+}^{b-a-1}$
to be a totally positive kernel of order $\left|\vec{n}\right|$, we have to
assume that $b-a-1\in\mathbb{Z}_{\geq 0}\cup(\left|\vec{n}\right|-2,\infty)$.
It is not possible to relax this condition as it was shown recently in [16,
Theorem 1.7] that these two statements are equivalent (see [5] for an
extensive discussion on the appearance of this kind of condition).

######  Corollary 2.4.

Suppose that all $b_{j}-a_{j}\in\mathbb{Z}_{\geq 0}$ and
$b_{i}-b_{j}\not\in\mathbb{Z}$ whenever $i\neq j$. Then the set of weights
$\\{w_{j}(x)\\}_{j=1}^{r}$ is an AT-system on $(0,1)$ for every multi-index
near the diagonal.

It is known that the classical orthogonality weights satisfy a Pearson
equation of the form

$[\sigma(x)w(x)]^{\prime}=\tau(x)w(x),$

where $\sigma,\tau$ are polynomials with $\deg\sigma\leq 2$ and $\deg\tau=1$
(see, e.g., [12, Chapter 23]). For the beta density $\mathcal{B}^{a,b}(x)$, we
have $\sigma(x)=x(1-x)$ and $\tau(x)=(a+1)-(b+1)x$. Instead of a Pearson
equation, the system of weights is now governed by a matrix Pearson equation.

######  Proposition 2.5.

Suppose that $b_{i}\neq b_{j}$ whenever $i\neq j$. Define
$\vec{W}=(w_{1},\dots,w_{r})^{T}$. Then

$[x(1-x)\vec{W}(x)]^{\prime}=\mathcal{T}(x)\vec{W}(x),$

where $\mathcal{T}(x)$ is a matrix with entries

$\mathcal{T}_{k,j}(x)=-\frac{\prod_{i=1}^{r}(a_{i}-b_{j})}{\prod_{i\neq
j}^{r}(b_{i}-b_{j})}+\begin{dcases}0,\ k\neq j,\\\ x-(b_{j}+1)(1-x),\
k=j.\end{dcases}$

###### Proof.

We aim to find a relation between the Mellin transforms $\hat{w}_{J}(s+1)$ and
the $\hat{w}_{j}(s)$ for $1\leq j\leq r$. Note that

$\sum_{j=1}^{r}[c_{j}+s\delta_{j,J}]\hat{w}_{j}(s)=\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+1)}\sum_{j=1}^{r}[c_{j}+s\delta_{j,J}]\prod_{i\neq
j}(s+b_{i})$

and that the sum on the right-hand side is a monic polynomial of degree $r$
with $r$ parameters $c_{j}$. Hence, we can try to pick the parameters in such
a way that the polynomial becomes $\prod_{i=1}^{r}(s+a_{i})$, which would then
lead to the relation

$\sum_{j=1}^{r}[c_{j}+s\delta_{j,J}]\hat{w}_{j}(s)=(s+b_{J}+1)\hat{w}_{J}(s+1).$
(6)

It is indeed possible to do this. Taking $s=-b_{k}$ with $k\neq J$, we find
$c_{k}\prod_{i\neq k}(b_{i}-b_{k})=\prod_{i=1}^{r}(a_{i}-b_{k})$, while taking
$s=-b_{J}$ gives $(c_{J}-b_{J})\prod_{i\neq
J}(b_{i}-b_{J})=\prod_{i=1}^{r}(a_{i}-b_{J})$. Hence,

$c_{j}=\frac{\prod_{i=1}^{r}(a_{i}-b_{j})}{\prod_{i\neq
j}^{r}(b_{i}-b_{j})}+\begin{dcases}0,\ j\neq J,\\\ b_{J},\ j=J.\end{dcases}$

Formula (6) can be written as

$(s-1)\hat{w}_{J}(s+1)-(s-1)\hat{w}_{J}(s)=\sum_{j=1}^{r}[c_{j}+\delta_{j,J}]\hat{w}_{j}(s)-(b_{J}+2)\hat{w}_{J}(s+1),$

after which it is straightforward to take the inverse Mellin transform and
obtain

$(x^{2}w_{J}(x))^{\prime}-(xw_{J}(x))^{\prime}=\sum_{j=1}^{r}[c_{j}+\delta_{j,J}]w_{j}(x)-(b_{J}+2)xw_{J}(x).$

Combining these relations for $1\leq J\leq r$ then leads to the desired
result. ∎

###  2.2 Type I polynomials

In what follows, we will study the type I multiple orthogonal polynomials
associated with the weights $(w_{1},\dots,w_{r})$ introduced in the previous
section for multi-indices near the diagonal. Essential to our approach will be
the fact that we can find an explicit expression for the Mellin transform of
the associated type I functions.

######  Theorem 2.6.

Let $\vec{n}\in\mathcal{N}^{r}$. Every type I function $F_{\vec{n}}$ has a
Mellin transform of the form

$\mathcal{F}_{\vec{n}}\cdot\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+\vec{n})}(1-s)_{\left|\vec{n}\right|-1},\quad\mathcal{F}_{\vec{n}}\in\mathbb{R}.$

###### Proof.

It follows from Lemma 2.2 that the Mellin transform of $F_{\vec{n}}$ is of the
form

$\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+\vec{n})}p_{\left|\vec{n}\right|-1}(s),$

where $p_{\left|\vec{n}\right|-1}(s)$ is a polynomial of degree at most
$\left|\vec{n}\right|-1$. The orthogonality conditions then imply that this
polynomial must be given by a scalar multiple of
$(1-s)_{\left|\vec{n}\right|-1}$. ∎

It turns out that if $a_{j}=b_{j}$ for some $j$, we can be more generous with
the restrictions on the multi-index, for example in the Jacobi-Piñeiro setting
(where $\vec{a}=\vec{b}$) the above holds for all multi-indices. However, we
won’t go into more detail on this phenomenon here.

From now on, we will assume that the (non-zero) type I functions
$F_{\vec{n}}(x)$ are normalized such that $\mathcal{F}_{\vec{n}}=1$.

Theorem 2.6 allows us to prove normality of multi-indices near the diagonal
under a different set of conditions on the parameters than in Corollary 2.4.
These conditions are less restrictive. Indeed, if
$b_{j}-a_{i}\in\mathbb{Z}_{<0}$ for some $i,j$, then we must also have
$b_{j}-b_{i}\in\mathbb{Z}_{<0}$ because $b_{i}-a_{i}\in\mathbb{Z}_{\geq 0}$.

######  Proposition 2.7.

Suppose that $b_{j}-a_{i}\not\in\mathbb{Z}_{<0}$ and
$b_{i}-b_{j}\not\in\mathbb{Z}$ whenever $i\neq j$. Then every multi-index near
the diagonal is normal for the set of weights $\\{w_{j}(x)\\}_{j=1}^{r}$.

###### Proof.

Under the stated conditions, the leading coefficients of the type I
polynomials can be determined explicitly as the unique residues at
$s=-b_{j}-n_{j}+1$. Indeed, one has,

$A_{\vec{n},j}[n_{j}-1]\frac{(\vec{a}-b_{j}-n_{j}+1)_{n_{j}-1}}{(\vec{b}-b_{j}-n_{j}+1)_{n_{j}-1}}=\frac{(-1)^{n_{j}-1}}{(n_{j}-1)!\prod_{i=1,i\neq
j}^{r}(b_{i}-b_{j}-n_{j}+1)_{n_{i}}}(b_{j}+n_{j})_{\left|\vec{n}\right|-1},$

from which it follows that $A_{\vec{n},j}[n_{j}-1]\neq 0$. ∎

The inverse Mellin transform immediately provides a representation as a
contour integral for the type I functions.

######  Corollary 2.8.

Let $\vec{n}\in\mathcal{N}^{r}$. Then,

$F_{\vec{n}}(x)=\frac{1}{2\pi
i}\int_{c-i\infty}^{c+i\infty}\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+\vec{n})}(1-s)_{\left|\vec{n}\right|-1}x^{-s}ds,$

for some $c>0$ with $-\min\\{a_{j}\mid j\\}<c<1$.

A Rodrigues-type formula also follows from Theorem 2.6.

######  Corollary 2.9.

Let $\vec{n}\in\mathcal{N}^{r}$. Then,

$F_{\vec{n}}(x)=\frac{d^{\left|\vec{n}\right|-1}}{dx^{\left|\vec{n}\right|-1}}\left[x^{\left|\vec{n}\right|-1}w_{0}(x;\vec{a},\vec{b}+\vec{n})\right].$

The type I functions can be interpreted as certain Mellin convolutions as
well.

######  Proposition 2.10.

Let $\vec{n}\in\mathcal{N}^{r}$. Then,

$F_{\vec{n}}=P_{\vec{n}}^{(I\mid\vec{b})}\ast\mathcal{B}^{a_{1},b_{1}}\ast\dots\ast\mathcal{B}^{a_{r},b_{r}},$

in terms of the type I Jacobi-Piñeiro polynomial
$P_{\vec{n}}^{(I\mid\vec{b})}$ associated with $(x^{b_{1}},\dots,x^{b_{r}})$.

###### Proof.

The Mellin transform of the beta density $\mathcal{B}^{a_{j},b_{j}}$ is
$\Gamma(s+a_{j})/\Gamma(s+b_{j})$, while the type I Jacobi-Piñeiro polynomial
$P_{\vec{n}}^{(I\mid\vec{b})}$ has Mellin transform
$(1-s)_{\left|\vec{n}\right|-1}/\prod_{j=1}^{q}(s+b_{j})_{n_{j}}$ (set
$\vec{a}=\vec{b}$ in Theorem 2.6, or see [38]). Hence the result follows from
the multiplicative property of the Mellin convolution and the uniqueness of
the Mellin transform. ∎

We can also express the type I function as a combination of certain
hypergeometric series. Here we used the following notation: given a vector
$\vec{a}\in\mathbb{R}^{r}$, we denote $\vec{a}^{\ast j}$ for the vector
$(a_{1},\dots,a_{j-1},a_{j+1},\dots,a_{r})\in\mathbb{R}^{r-1}$.

######  Proposition 2.11.

Let $\vec{n}\in\mathcal{N}^{r}$ and suppose that
$a_{i}-a_{j}\not\in\mathbb{Z}$ whenever $i\neq j$ and
$b_{i}-a_{j}\not\in\mathbb{Z}$ for all $i,j$. Then,

$F_{\vec{n}}(x)=\sum_{j=1}^{r}\frac{(a_{j}+1)_{\left|\vec{n}\right|-1}\Gamma(\vec{a}^{\ast
j}-a_{j})}{\Gamma(\vec{b}-a_{j}+\vec{n})}{}_{r+1}F_{r}\left(\begin{array}[]{c}\left|\vec{n}\right|+a_{j},-\vec{n}-\vec{b}+a_{j}+1\\\
a_{j}+1,-\vec{a}^{\ast j}+a_{j}+1\end{array};x\right)x^{a_{j}}.$

###### Proof.

We will work out the contour integral in Corollary 2.8 via the residue
theorem. The integrand has simple poles at $s=-a_{j}-k$ for
$k\in\mathbb{Z}_{\geq 0}$ and $1\leq j\leq r$ with residue

$\displaystyle x^{a_{j}+k}\frac{(-1)^{k}}{k!}$
$\displaystyle\frac{\Gamma(\vec{a}^{\ast
j}-a_{j}-k)}{\Gamma(\vec{b}-a_{j}+\vec{n}-k)}(a_{j}+k+1)_{\left|\vec{n}\right|-1}$
$\displaystyle=x^{a_{j}+k}\frac{(-1)^{k}}{k!}\frac{\Gamma(\vec{a}^{\ast
j}-a_{j})}{(\vec{a}^{\ast
j}-a_{j}-k)_{k}}\frac{(\vec{b}-a_{j}+\vec{n}-k)_{k}}{\Gamma(\vec{b}-a_{j}+\vec{n})}\frac{\Gamma(a_{j}+k+\left|\vec{n}\right|)}{\Gamma(a_{j}+k+1)}.$

By making use of the identity $(x)_{n}=(-x-n+1)_{n}$, we can rewrite the
latter as

$x^{a_{j}+k}\frac{1}{k!}\frac{\Gamma(\vec{a}^{\ast
j}-a_{j})}{(a_{j}-\vec{a}^{\ast
j}+1)_{k}}\frac{(a_{j}-\vec{b}-\vec{n}+1)_{k}}{\Gamma(\vec{b}-a_{j}+\vec{n})}\frac{(a_{j}+\left|\vec{n}\right|)_{k}}{(a_{j}+1)_{k}}\frac{\Gamma(a_{j}+\left|\vec{n}\right|)}{\Gamma(a_{j}+1)},$

which then leads to the desired result. ∎

Observe that the terms that appear in this sum (so including the $x^{a_{j}}$)
are all solutions to the standard hypergeometric differential equation, see
[8, Eq. 16.8.6], for

${}_{r+1}F_{r}\left(\begin{array}[]{c}\left|\vec{n}\right|,-\vec{n}-\vec{b}+1\\\
-\vec{a}+1\end{array};x\right).$ (7)

Hence the type I function satisfies the same differential equation.

Finally, we will derive some expressions for the type I polynomials
themselves.

######  Theorem 2.12.

Let $\vec{n}\in\mathcal{N}^{r}$ and suppose that
$b_{j}-a_{i}\not\in\mathbb{Z}_{<0}$ and $b_{i}-b_{j}\not\in\mathbb{Z}_{\geq
0}$ whenever $i\neq j$. The type I polynomials $A_{\vec{n},j}$ are given by

$A_{\vec{n},j}(x)=\frac{(\vec{a}-b_{j})_{1}}{(\vec{b}^{\ast
j}-b_{j})_{1}}\sum_{J=1}^{r}\sum_{K=0}^{n_{J}-1}P_{\vec{n},J}[K]\sum_{k=0}^{K-1+\delta_{j,J}}\frac{(\vec{b}^{\ast
j}-b_{J}-K)_{k+1}(b_{j}-b_{J}-K)_{k}}{(\vec{a}-b_{J}-K)_{k+1}}x^{k},$

where

$P_{\vec{n},J}[K]=\frac{(b_{J}+K+1)_{\left|\vec{n}\right|-1}}{\prod_{i=1,i\neq
J}^{r}(b_{i}-b_{J}-K)_{n_{i}}}\frac{(-1)^{K}}{K!(n_{J}-K-1)!}.$

###### Proof.

We will make use of the connection between multiple orthogonal polynomials and
Hermite-Padé approximation (see, e.g., [12, Chapter 23]). It is known that
there exists a polynomial $B_{\vec{n}}(z)$ such that

$\sum_{j=1}^{r}A_{\vec{n},j}(z)f_{j}(z)-B_{\vec{n}}(z)=\sum_{k\geq
0}\hat{F}_{\vec{n}}(k+1)\frac{1}{z^{k+1}},\quad f_{j}(z)=\sum_{k\geq
0}\frac{\Gamma(k+\vec{a}+1)}{\Gamma(k+\vec{b}+\vec{e}_{j}+1)}\frac{1}{z^{k+1}}.$

We can obtain an expression for the $A_{\vec{n},j}(z)$ by working out the
right-hand side of the above. It follows from Theorem 2.6 that

$\hat{F}_{\vec{n}}(s)=\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b})}\frac{(1-s)_{\left|\vec{n}\right|-1}}{\prod_{j=1}^{q}(s+b_{j})_{n_{j}}}.$

A partial fraction decomposition can be used to expand the second factor as

$\frac{(1-s)_{\left|\vec{n}\right|-1}}{\prod_{j=1}^{r}(s+b_{j})_{n_{j}}}=\sum_{j=1}^{r}\sum_{l=0}^{n_{j}-1}\frac{P_{\vec{n},j}[l]}{s+l+b_{j}},$

in terms of the coefficients $P_{\vec{n},j}[l]$ of the type I Jacobi-Piñeiro
polynomials associated with the weights $(x^{b_{1}},\dots,x^{b_{r}})$. We then
obtain

$\sum_{k\geq
0}\hat{F}_{\vec{n}}(k+1)\frac{1}{z^{k+1}}=\sum_{j=1}^{r}\sum_{l=0}^{n_{j}-1}P_{\vec{n},j}[l]\sum_{k\geq
0}\frac{\Gamma(k+\vec{a}+1)}{\Gamma(k+\vec{b}+1)}\frac{1}{k+l+b_{j}+1}\frac{1}{z^{k+1}},$
(8)

and it remains to find an operation to lower $L$ in

$I_{J,L}(z)=\sum_{k\geq
0}\frac{\Gamma(k+\vec{a}+1)}{\Gamma(k+\vec{b}+1)}\frac{1}{k+L+b_{J}+1}\frac{1}{z^{k+1}},$

because $I_{J,0}(z)=f_{J}(z)$. We can take the first term out of the series to
get

$I_{J,L-1}(z)=\frac{\Gamma(\vec{a}+1)}{\Gamma(\vec{b}+1)}\frac{1}{L+b_{J}}+\sum_{k\geq
0}\frac{\Gamma(k+\vec{a}+1)}{\Gamma(k+\vec{b}+1)}\frac{(k+\vec{a}+1)_{1}}{(k+\vec{b}+1)_{1}}\frac{1}{k+L+b_{J}+1}\frac{1}{z^{k+2}}.$

A partial fraction decomposition in the variable $k$ then gives

$\frac{(k+\vec{a}+1)_{1}}{(k+\vec{b}+1)_{1}}\frac{1}{k+L+b_{J}+1}=\sum_{j=1}^{r}\frac{\alpha_{j}^{L}}{k+b_{j}+1}+\frac{\alpha_{0}^{L}}{k+L+b_{J}+1},$
(9)

where

$\alpha_{j}^{L}=\frac{(\vec{a}-b_{j})_{1}}{(\vec{b}^{\ast
j}-b_{j})_{1}}\frac{1}{L+b_{J}-b_{j}},\quad\alpha_{0}^{L}=\frac{(\vec{a}-b_{J}-L)_{1}}{(\vec{b}-b_{J}-L)_{1}}.$

Therefore,

$I_{J,L-1}(z)=\frac{\Gamma(\vec{a}+1)}{\Gamma(\vec{b}+1)}\frac{1}{L+b_{J}}+\frac{1}{z}\sum_{j=1}^{r}\alpha_{j}^{L}f_{j}(z)+\frac{\alpha_{0}^{L}}{z}I_{J,L}(z),$

which is a recurrence relation of the form
$I_{J,L}(z)=C_{L}(z)+D_{L}(z)I_{J,L-1}(z)$ with

$C_{L}(z)=-\frac{1}{\alpha_{0}^{L}}\left(\sum_{j=1}^{r}\alpha_{j}^{L}f_{j}(z)+\frac{\Gamma(\vec{a}+1)}{\Gamma(\vec{b}+1)}\frac{z}{L+b_{J}}\right),\quad
D_{L}(z)=\frac{z}{\alpha_{0}^{L}}.$

We can obtain a formula for $I_{J,L}(z)$ by consecutive application of this
recurrence relation

$I_{J,L}(z)=\sum_{l=0}^{L-1}C_{L-l}(z)\prod_{j=1}^{l}D_{L-l+j}(z)+\prod_{j=1}^{L}D_{j}(z)f_{J}(z).$

Note that
$\prod_{j=1}^{l}D_{L-l+j}(z)=z^{l}(\vec{b}-b_{J}-L)_{l}/(\vec{a}-b_{J}-L)_{l}$
so that

$\displaystyle I_{J,L}(z)$
$\displaystyle=-\sum_{j=1}^{r}\frac{(\vec{a}-b_{j})_{1}}{(\vec{b}^{\ast
j}-b_{j})_{1}}f_{j}(z)\sum_{l=0}^{L-1}\frac{(\vec{b}-b_{J}-L)_{l+1}}{(\vec{a}-b_{J}-L)_{l+1}}\frac{z^{l}}{L-l+b_{J}-b_{j}}$
$\displaystyle-\frac{\Gamma(\vec{a}+1)}{\Gamma(\vec{b}+1)}\sum_{l=0}^{L-1}\frac{(\vec{b}-b_{J}-L)_{l+1}}{(\vec{a}-b_{J}-L)_{l+1}}\frac{z^{l+1}}{L-l+b_{J}}$
$\displaystyle+z^{L}\frac{(\vec{b}-b_{J}-L)_{L}}{(\vec{a}-b_{J}-L)_{L}}f_{J}(z).$

We can write this as

$\displaystyle I_{J,L}(z)$
$\displaystyle=\sum_{j=1}^{r}\frac{(\vec{a}-b_{j})_{1}}{(\vec{b}^{\ast
j}-b_{j})_{1}}f_{j}(z)\sum_{l=0}^{L-1+\delta_{j,J}}\frac{(\vec{b}^{\ast
j}-b_{J}-L)_{l+1}(b_{j}-b_{J}-L)_{l}}{(\vec{a}-b_{J}-L)_{l+1}}z^{l}$
$\displaystyle-\frac{\Gamma(\vec{a}+1)}{\Gamma(\vec{b}+1)}\sum_{l=0}^{L-1}\frac{(\vec{b}-b_{J}-L)_{l+1}}{(\vec{a}-b_{J}-L)_{l+1}}\frac{1}{L-l+b_{J}}z^{l+1},$
(10)

because the term we added to the sum is

$\frac{(\vec{a}-b_{J})_{1}}{(\vec{b}^{\ast
j}-b_{J})_{1}}f_{J}(z)\frac{(\vec{b}^{\ast
J}-b_{J}-L)_{L+1}(-L)_{L}}{(\vec{a}-b_{J}-L)_{L+1}}z^{L}=f_{J}(z)\frac{(\vec{b}^{\ast
J}-b_{J}-L)_{L}(-L)_{L}}{(\vec{a}-b_{J}-L)_{L}}z^{L}.$

The stated formulas for the type I polynomials then follow. ∎

The polynomial that appears together with the coefficient $P_{\vec{n},J}[K]$
in the above formula can be represented as the hypergeometric polynomial

$\frac{(\vec{b}-b_{J}-K)_{1}^{\ast
j}}{(\vec{a}-b_{J}-K)_{1}}\,{}_{r+1}F_{r}\left(\begin{array}[]{c}1,\vec{b}-b_{J}+1-\vec{e}_{j}-K\\\
\vec{a}-b_{J}+1-K\end{array};x\right).$

The type I polynomials can therefore be constructed with the same coefficients
as in the Jacobi-Piñeiro setting, but using these
${}_{r+1}F_{r}$-hypergeometric polynomials as a basis instead of the
monomials.

###  2.3 Type II polynomials

In this section, we will study the type II multiple orthogonal polynomials
associated with the weights $(w_{1},\dots,w_{r})$ introduced in Section 2.1
for multi-indices near the diagonal. Formulas for these type II polynomials
were already obtained in [24], however, their proof is implicit in the sense
that one simply verifies that the stated polynomials satisfy the desired
orthogonality conditions by making use of a property of the underlying
hypergeometric series. We will provide a constructive method to derive those
formulas. Once an expression is known, a similar approach would also allow us
to verify the orthogonality conditions directly. Additionally, our formulas
are proven for a larger set of multi-indices: near the diagonal
$\mathcal{N}^{r}$ instead of on the step-line $\mathcal{S}^{r}$. Afterwards,
we will describe the asymptotic zero distribution and interpret it through
free probability.

######  Theorem 2.13.

Let $\vec{n}\in\mathcal{N}^{r}$. Every type II polynomial $P_{\vec{n}}(x)$ is
of the form

$\mathcal{P}_{\vec{n}}\cdot{}_{r+1}F_{r}\left(\begin{array}[]{c}-\left|\vec{n}\right|,\vec{b}+\vec{n}+1\\\
\vec{a}+1\end{array};x\right),\quad\mathcal{P}_{\vec{n}}\in\mathbb{R}.$

###### Proof.

Moving into the complex plane allows us to write

$P_{\vec{n}}(x)=\frac{1}{2\pi
i}\sum_{k=0}^{\left|\vec{n}\right|}P_{\vec{n}}[k]\frac{\Gamma(k+\vec{a}+1)}{\Gamma(k+\vec{b}+1)}\int_{\Sigma}\frac{\Gamma(t+\vec{b}+1)}{\Gamma(t+\vec{a}+1)}x^{t}\frac{dt}{t-k},$

in terms of a contour $\Sigma$ encircling the interval
$[0,\left|\vec{n}\right|]$ avoiding the poles of the integrand that are not in
this interval. Hence, we find

$P_{\vec{n}}(x)=\frac{1}{2\pi
i}\sum_{k=0}^{\left|\vec{n}\right|}\int_{\Sigma}\frac{p_{\left|\vec{n}\right|}(t)}{(-t)_{\left|\vec{n}\right|+1}}\frac{\Gamma(t+\vec{b}+1)}{\Gamma(t+\vec{a}+1)}x^{t}dt,$

where $p_{\left|\vec{n}\right|}(t)$ is a polynomial of degree at most
$\left|\vec{n}\right|$ that is determined by the relation

$\frac{p_{\left|\vec{n}\right|}(t)}{(-t)_{\left|\vec{n}\right|+1}}=\sum_{k=0}^{\left|\vec{n}\right|}\frac{\Gamma(k+\vec{a}+1)}{\Gamma(k+\vec{b}+1)}\frac{P_{\vec{n}}[k]}{t-k}.$

Let $s\in\mathbb{N}$. Then,

$\int_{0}^{1}P_{\vec{n}}(x)x^{s}w_{j}(x)dx=\frac{1}{2\pi
i}\int_{\Sigma}\frac{p_{\left|\vec{n}\right|}(t)}{(-t)_{\left|\vec{n}\right|+1}}\frac{(t+\vec{a}+1)_{s}}{(t+\vec{b}+1)_{s}(s+t+b_{j}+1)}dt.$
(11)

Note that the integrand is a rational function of order $\mathcal{O}(t^{-2})$
as $t\to\infty$, so that we can compute the integral by moving the contour
away to infinity and by adding the residues of poles in
$\mathbb{C}\backslash[0,\left|\vec{n}\right|]$. Since the integral needs to
vanish for $0\leq s\leq n_{j}-1$ and $1\leq j\leq r$, we don’t want to have
any poles in that domain. Under the near-diagonal restriction, this implies
that $\prod_{j=1}^{r}(t+b_{j}+1)_{n_{j}}$ needs to divide
$p_{\left|\vec{n}\right|}(t)$. Since the former is a polynomial of degree
$\left|\vec{n}\right|$, we must have
$p_{\left|\vec{n}\right|}(t)=\mathcal{P}_{\vec{n}}\cdot\prod_{j=1}^{r}(t+b_{j}+1)_{n_{j}}$
for some $\mathcal{P}_{\vec{n}}\in\mathbb{R}$. Consequently, after comparing
residues, we find

$P_{\vec{n}}[k]=\mathcal{P}_{\vec{n}}\cdot\frac{\Gamma(k+\vec{b}+1)}{\Gamma(k+\vec{a}+1)}\frac{\prod_{j=1}^{r}(k+b_{j}+1)_{n_{j}}}{(-k)_{k}(1)_{\left|\vec{n}\right|-k}}$

and this leads to the desired result. ∎

Similarly as with the type I functions, if $a_{j}=b_{j}$ for some $j$, we can
be more generous with the restrictions on the multi-index. Also, observe the
similarities with the hypergeometric series in (7).

In what follows, we will assume that the (non-zero) type II polynomials
$P_{\vec{n}}(x)$ are normalized such that
$\mathcal{P}_{\vec{n}}=\Gamma(\vec{b}+\vec{n}+1)/[\Gamma(\vec{a}+1)(\left|\vec{n}\right|)!]$.

From the proof of Theorem 2.13, it is clear that we have the following contour
integral representation for the type II polynomials.

######  Corollary 2.14.

Let $\vec{n}\in\mathcal{N}^{r}$. Then,

$P_{\vec{n}}(x)=\frac{1}{2\pi
i}\int_{\Sigma}\frac{1}{(-t)_{\left|\vec{n}\right|+1}}\frac{\Gamma(t+\vec{b}+\vec{n}+1)}{\Gamma(t+\vec{a}+1)}x^{t}dt,$

where $\Sigma$ is a positively oriented contour that starts and ends at
$\infty$ and encircles the positive real line, without enclosing any poles of
the integrand not in $[0,\infty)$.

Note that there is some kind of duality between the contour integral
representations for the type I functions in Corollary 2.8 and type II
polynomials in Corollary 2.14. Such duality is expected due to the underlying
biorthogonality relations.

We will end this section by discussing some properties of the zeros of the
type II polynomials that were obtained in [24]. By exploiting the connection
with branched continued fractions, one was able to prove the theorem below.
Note that under the, generally more restrictive, condition that all
$b_{j}-a_{j}\in\mathbb{Z}_{\geq 0}\cup(\left|\vec{n}\right|-1,\infty)$ and
$b_{J}-a_{J}>\left|\vec{n}\right|-n_{J}$ for some $J$, it also follows from
the fact that the underlying system of weights is an AT-system (see
Proposition 2.3).

######  Theorem 2.15 (Lima [24], Theorem 4.8).

Let $\vec{n}\in\mathcal{S}^{r}$ and suppose that $b_{j}>a_{i}$ if $i\leq j$
and $b_{j}>a_{i}-1$ if $i>j$. The type II polynomials $P_{\vec{n}}(x)$ have
real, positive, simple zeros.

Using an underlying recurrence relation, it was then shown that the Stieltjes
transform $S(z)$ of the asymptotic zero distribution of $P_{\vec{n}}(x)$ with
$\vec{n}\in\mathcal{S}^{r}$ as $\left|\vec{n}\right|\to\infty$ satisfies the
algebraic equation $(S(z)+r^{r}/(r+1)^{r+1})^{r+1}=zS(z)^{r}$. From this, one
could conclude that $P_{\vec{n}}(x)$ has the same asymptotic zero distribution
as the type II Jacobi-Piñeiro polynomials; it was already proven in [34] that
the Stieltjes transform of the latter satisfies the same algebraic equation.

######  Theorem 2.16 (Lima [24], Theorem 8.6).

Suppose that $b_{j}>a_{i}$ if $i\leq j$ and $b_{j}>a_{i}-1$ if $i>j$. The
asymptotic zero distribution of $P_{\vec{n}}(x)$ with
$\vec{n}\in\mathcal{S}^{r}$ as $\left|\vec{n}\right|\to\infty$ has a density,
with support on $(0,1)$, given by

$v_{r}(x)=\frac{r+1}{\pi
x}\frac{\sin{\varphi}\sin{r\varphi}\sin{(r+1)\varphi}}{(r+1)^{2}(\sin{r\varphi})^{2}-2r(r+1)\sin{(r+1)\varphi\sin{r\varphi}\cos{\varphi}}+r^{2}(\sin{(r+1)\varphi})^{2}},$

after the change of variables

$x=\frac{r^{r}}{(r+1)^{r+1}}\frac{(\sin{(r+1)\varphi})^{r+1}}{\sin\varphi(\sin{r\varphi})^{r}},\quad\varphi\in\left(0,\frac{\pi}{r+1}\right).$

Another way to interpret these results is through the (finite) free
multiplicative convolution. The idea to use this tool to study the zeros of
hypergeometric polynomials was developed recently in [29]. Many known (and
new) results were collected here, some go back to Szegő [40]. Recently, there
has been renewed interest in the finite free convolution of polynomials due to
the connection with free probability, see [27]. Since the type II polynomials
here are hypergeometric polynomials that don’t belong to the class that was
studied in [29], we can’t use their results on the zero behavior and
asymptotic zero distribution directly.

Following the definition in [27], the finite free multiplicative convolution
(called $d$-th symmetric multiplicative convolution here) of two polynomials

$p(x)=\sum_{k=0}^{d}x^{d-k}(-1)^{k}p[k],\quad
q(x)=\sum_{k=0}^{d}x^{d-k}(-1)^{k}q[k],$

of degree at most $d$, is defined as

$p(x)\boxtimes_{d}q(x)=\sum_{k=0}^{d}x^{d-k}(-1)^{k}p[k]q[k].$

An important property is that, given two sequences of polynomials
$(p_{d})_{d\in\mathbb{N}}$ and $(q_{d})_{d\in\mathbb{N}}$ with $\deg
p_{d}=\deg q_{d}=d$, the asymptotic zero distribution of the finite free
multiplicative convolution $p_{d}\boxtimes_{d}q_{d}$ as $d\to\infty$ is given
by the free multiplicative convolution $\mu\boxtimes\nu$ of the asymptotic
zero distribution $\mu$ of $p_{d}$ and $\nu$ of $q_{d}$ as $d\to\infty$ (see
[4, Theorem 1.4]). The multiplicative convolution of two measures is defined
through the $S$-transform as

$S_{\mu\boxtimes\nu}(z)=S_{\mu}(z)\cdot S_{\nu}(z).$

The $S$-transform of a measure $\mu$ is defined by

$S_{\mu}(z)=m_{\mu}^{-1}(z)\frac{z+1}{z},$

where $m_{\mu}^{-1}(z)$ is the (formal) inverse of the formal power series
$m_{\mu}(z)=\sum_{k\geq 1}m_{\mu,k}z^{k}$ with
$m_{\mu,k}=\int_{\mathbb{R}}x^{k}d\mu(x)$. The condition $m_{\mu,1}\neq 0$
ensures that this inverse exists.

Unfortunately, in general, a direct decomposition of the type II polynomials
$P_{\vec{n}}(x;\vec{a},\vec{b})$ isn’t helpful to study its zeros. It would
lead to

$\displaystyle{}_{r+1}F_{r}$
$\displaystyle\left(\begin{array}[]{c}-\left|\vec{n}\right|,\vec{b}+\vec{n}+1\\\
\vec{a}+1\end{array};x\right)$
$\displaystyle={}_{2}F_{1}\left(\begin{array}[]{c}-\left|\vec{n}\right|,b_{1}+n_{1}+1\\\
a_{1}+1\end{array};x\right)\boxtimes_{\left|\vec{n}\right|}\cdots\boxtimes_{\left|\vec{n}\right|}{}_{2}F_{1}\left(\begin{array}[]{c}-\left|\vec{n}\right|,b_{r}+n_{r}+1\\\
a_{r}+1\end{array};x\right),$

but the polynomials that appear here are Jacobi polynomials
$P_{\left|\vec{n}\right|}^{(a_{j},b_{j}-a_{j}+n_{j}-\left|\vec{n}\right|)}(x)$
with a parameter that is typically negative and not an integer, thus having
roots in the complex plane. Under the condition that all
$b_{j}-a_{j}\in\mathbb{Z}_{\geq 0}\cup(\left|\vec{n}\right|-1,\infty)$ and
$b_{J}-a_{J}>\left|\vec{n}\right|-n_{J}$ for some $J$, the Jacobi polynomials
do have real, positive roots, and at least one of them also has simple roots.
It is known that if the roots of a collection of polynomials are real,
positive and at least one polynomial also has simple roots, the same property
holds for their convolution (see [40] and [15, Theorem 6]). Hence under these
conditions, we would be able to obtain the result in Theorem 2.15. Remarkably,
these were exactly the conditions that were required for the system of weights
to be an AT-system for the multi-index $\vec{n}$ (see Proposition 2.3).

By making use of the fact that the type II Jacobi-Piñeiro polynomials have the
same asymptotic zero distribution as the type II polynomials
$P_{\vec{n}}(x;\vec{a},\vec{b})$, it is still possible to give a description
of the asymptotic zero distribution using the free convolution.

######  Theorem 2.17.

Suppose that $b_{j}>a_{i}$ if $i\leq j$ and $b_{j}>a_{i}-1$ if $i>j$. The
asymptotic zero distribution of $P_{\vec{n}}(x)$ for
$\vec{n}\in\mathcal{S}^{r}$ as $\left|\vec{n}\right|\to\infty$ is given by
$((1-1/r)\delta_{1}+(1/r)\mu_{r-1})^{\boxtimes r}$ in terms of the deformed
arcsin distribution $\mu_{r-1}$ from Lemma 2.18.

###### Proof.

It is sufficient to prove this for the type II Jacobi-Piñeiro polynomials.
They are given by

$\displaystyle{}_{r+1}F_{r}$
$\displaystyle\left(\begin{array}[]{c}-\left|\vec{n}\right|,\vec{b}+\vec{n}+1\\\
\vec{b}+1\end{array};x\right)$
$\displaystyle={}_{2}F_{1}\left(\begin{array}[]{c}-\left|\vec{n}\right|,b_{1}+n_{1}+1\\\
b_{1}+1\end{array};x\right)\boxtimes_{\left|\vec{n}\right|}\cdots\boxtimes_{\left|\vec{n}\right|}{}_{2}F_{1}\left(\begin{array}[]{c}-\left|\vec{n}\right|,b_{r}+n_{r}+1\\\
b_{r}+1\end{array};x\right),$

where

${}_{2}F_{1}\left(\begin{array}[]{c}-\left|\vec{n}\right|,b_{j}+n_{j}+1\\\
b_{j}+1\end{array};x\right)=(1-x)^{\left|\vec{n}\right|-n_{j}}P_{n_{j}}^{(b_{j},\left|\vec{n}\right|-n_{j})}(x)$

in terms of a Jacobi polynomial on $[0,1]$. ∎

######  Lemma 2.18.

The asymptotic zero distribution $\mu_{\alpha}$ of the Jacobi polynomial
$P_{n}^{(a,\alpha_{n})}(x)$ with $\alpha_{n}>0$ and $\alpha_{n}/n\to\alpha$ as
$n\to\infty$ has density

$u_{\alpha}(x)=\frac{\sqrt{4(\alpha+1)x-(\alpha+2)^{2}x^{2}}}{2\pi
x(1-x)},\quad 0<x<4(\alpha+1)/(\alpha+2)^{2}.$

###### Proof.

The zeros of $P_{n}^{(a,\alpha_{n})}(x)$ are real, simple and contained in the
interval $[0,1]$, hence we can use the usual method to transform the
differential equation for $P_{n}^{(a,\alpha_{n})}(z)$ to an algebraic equation
for the Stieltjes transform $S_{\alpha}(z)$ of its asymptotic zero
distribution (see, e.g., [23, Section 3.1.4] or [43, Section 3.4.2]). The
differential equation is given by

$z(1-z)y^{\prime\prime}+[a+1-(a+\alpha_{n}+2)z]y^{\prime}+n(n+a+\alpha_{n}+1)y=0$

(see, e.g., [41, Eq. (4.21.1)]). The associated algebraic equation is

$z(1-z)S^{2}-\alpha zS+\alpha+1=0,$

which we may solve explicitly for two solutions

$S_{\pm}(z)=\frac{\alpha
z\pm\sqrt{(\alpha+2)^{2}z^{2}-4(\alpha+1)z}}{2z(1-z)}.$

Since $zS_{-}(z)\to 1$ as $z\to\infty$, we have

$u_{\alpha}(x)=-\lim_{y\to
0^{+}}\frac{1}{\pi}\text{Im}(S_{-}(x+iy))=\frac{\sqrt{4(\alpha+1)x-(\alpha+2)^{2}x^{2}}}{2\pi
x(1-x)},$

for $0<x<4(\alpha+1)/(\alpha+2)^{2}$. ∎

Another interesting decomposition of the type II polynomials
$P_{\vec{n}}(x;\vec{a},\vec{b})$ is in terms of the type II Jacobi-Piñeiro
polynomial $P_{\vec{n}}^{(II\mid\vec{b})}(x)$ associated with
$(x^{b_{1}},\dots,x^{b_{r}})$,

$P_{\vec{n}}(x;\vec{a},\vec{b})=P_{\vec{n}}^{(II\mid\vec{b})}(x)\boxtimes_{\left|\vec{n}\right|}{}_{2}F_{1}\left(\begin{array}[]{c}-\left|\vec{n}\right|,b_{1}+1\\\
a_{1}+1\end{array};x\right)\boxtimes_{\left|\vec{n}\right|}\cdots\boxtimes_{\left|\vec{n}\right|}{}_{2}F_{1}\left(\begin{array}[]{c}-\left|\vec{n}\right|,b_{r}+1\\\
a_{r}+1\end{array};x\right)$

because of its similarities to the Mellin convolution formula for the type I
function in Proposition 2.10.

##  3 Laguerre-like setting

###  3.1 System of weights

Let $p>q$ and suppose that $\vec{a}\in(-1,\infty)^{p}$ and
$\vec{b}\in(-1,\infty)^{q}$ satisfy $a_{j}<b_{j}$ for $1\leq j\leq q$. We will
consider $p$ weights
$\\{w_{j}(x;\vec{a},\vec{b})\\}_{j=1}^{q}\cup\\{v_{j}(x;\vec{a},\vec{b})\\}_{j=1}^{p-q}$
on $(0,\infty)$. The first weights are defined by
$w_{j}(x;\vec{a},\vec{b})=w_{0}(x;\vec{a},\vec{b}+\vec{e}_{j})$, similarly as
before. In this setting, a weight $w_{0}(x;\vec{a},\vec{b})$ with moments
$m_{k}(\vec{a},\vec{b})$ as in (2) can be constructed by taking the Mellin
convolution of $q$ beta densities $\mathcal{B}^{a_{j},b_{j}}$, $1\leq j\leq
q$, and $p-q$ gamma densities $\mathcal{G}^{a_{j}}$, $q+1\leq j\leq p$.
Indeed, a beta density $\mathcal{B}^{a,b}$ has Mellin transform
$\Gamma(s+a)/\Gamma(s+b)$ and a gamma density $\mathcal{G}^{a}$ has Mellin
transform $\Gamma(s+a)$. Since Carleman’s condition
$\sum_{k=0}^{\infty}m_{k}(\vec{a},\vec{b})^{-1/(2k)}=\infty$ is satisfied,
there can be no other (positive) measures on the positive real line with those
moments. The weight $w_{j}(x;\vec{a},\vec{b})$ then has Mellin transform

$\int_{0}^{1}x^{s-1}w_{j}(x;\vec{a},\vec{b})dx=\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+\vec{e}_{j})}.$

The remaining weights arise as $v_{j}(x;\vec{a},\vec{b})=(x\cdot
d/dx)^{j-1}[w_{0}(x;\vec{a},\vec{b})]$ and have Mellin transforms

$\int_{0}^{\infty}x^{s-1}v_{j}(x;\vec{a},\vec{b})dx=\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b})}(s-1)^{j-1}.$

It is natural to set up the weights in this way, because it leads to desirable
limiting relations between the sets of weights in all cases with $p\geq q$. We
will describe those relations in the following lemma. Similar limiting
relations will then also hold for the type I and type II polynomials.

######  Lemma 3.1.

Given a set of weights
$\\{w_{j}(x;\vec{a},\vec{b})\\}_{j=1}^{q}\cup\\{v_{j}(x;\vec{a},\vec{b})\\}_{j=1}^{p-q}$,
we can construct a new set of weights $\\{w_{j}(x;\vec{a},\vec{b}^{\ast
q})\\}_{j=1}^{q-1}\cup\\{v_{j}(x;\vec{a},\vec{b}^{\ast q})\\}_{j=1}^{p-q+1}$
via the limiting relations

  * i)

$\lim_{b_{q}\to\infty}\Gamma(b_{q})w_{j}(x/b_{q};\vec{a},\vec{b})=w_{j}(x;\vec{a},\vec{b}^{\ast
q})$ for $1\leq j\leq q-1$,

  * ii)

$\lim_{b_{q}\to\infty}\Gamma(b_{q}+1)(\lambda_{0}w_{q}(x/b_{q};\vec{a},\vec{b})+\sum_{j=1}^{J}\lambda_{j}v_{j}(x/b_{q};\vec{a},\vec{b}))=v_{J+1}(x;\vec{a},\vec{b}^{\ast
q})$ with $\lambda_{j}=(-1)^{J-j}(b_{q}+1)^{J-j}$ for $0\leq J\leq p-q$.

###### Proof.

Part i) immediately follows by comparing the underlying Mellin transforms. For
Part ii), we note that the coefficients are such that

$\lambda_{0}\hat{w}_{q}(k+1;\vec{a},\vec{b})+\sum_{j=1}^{J}\lambda_{j}\hat{v}_{j}(k+1;\vec{a},\vec{b})=\hat{w}_{q}(k+1;\vec{a},\vec{b})k^{J}.$

Indeed, the left-hand side is of the form

$\frac{\Gamma(k+\vec{a}+1)}{\Gamma(k+\vec{b}+1)}\left(\frac{\lambda_{0}}{k+b_{q}+1}+\sum_{j=1}^{J}\lambda_{j}k^{j-1}\right)=\hat{w}_{q}(k+1;\vec{a},\vec{b})p_{J}(k)$

in terms of a polynomial $p_{J}(k)$ of degree $J$ with $J+1$ unknown
coefficients, which we may solve for $p_{J}(k)=k^{J}$. ∎

In particular, the above procedure describes a way to convert a set of Jacobi-
like weights $\\{w_{j}(x;\vec{a},\vec{b})\\}_{j=1}^{q}$ on $[0,1]$ (so
$\vec{a},\vec{b}\in\mathbb{R}^{q}$) to a set of Laguerre-like weights
$\\{w_{j}(x;\vec{a},\vec{b}^{\ast
q})\\}_{j=1}^{q}\cup\\{v_{1}(x;\vec{a},\vec{b}^{\ast q})\\}$ on $[0,\infty)$;
we have to use the relations

$\lim_{b_{q}\to\infty}\Gamma(b_{q})w_{j}(x/b_{q};\vec{a},\vec{b})=w_{j}(x;\vec{a},\vec{b}^{\ast
q}),\quad\lim_{b_{q}\to\infty}\Gamma(b_{q}+1)w_{q}(x/b_{q};\vec{a},\vec{b})=v_{1}(x;\vec{a},\vec{b}^{\ast
q}).$

The limiting relations indicate that we should be able to obtain formulas for
the multiple orthogonal polynomials for multi-indices of the form
$\vec{n}\sqcup\vec{m}:=(n_{1},\dots,n_{q},m_{1},\dots,m_{p-q})$ that are near
the diagonal, for which the second part, associated with the weights
$\\{v_{j}\\}_{j=1}^{p-q}$, is restricted to the step-line (because we require
weights that correspond to previous components in their construction).

We can show that, under some conditions, the system of weights will be an AT-
system for the described multi-indices. To prove this, we will again analyze
the Mellin transform of polynomial combinations of the weights.

######  Lemma 3.2.

Let $\vec{n}\sqcup\vec{m}\in\mathcal{N}^{p}$ with
$\vec{m}\in\mathcal{S}^{p-q}$. The Mellin transform of a polynomial
combination
$F(x)=\sum_{j=1}^{q}A_{j}(x)w_{j}(x)+\sum_{j=1}^{p-q}A_{q+j}(x)v_{j}(x)$ with
$\deg A_{j}\leq n_{j}-1$ and $\deg A_{q+j}\leq m_{j}-1$ is of the form

$\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+\vec{n})}p_{\left|\vec{n}\right|+\left|\vec{m}\right|-1}(s),$

where $p_{\left|\vec{n}\right|+\vec{m}|-1}(s)$ is a polynomial of degree at
most $\left|\vec{n}\right|+\left|\vec{m}\right|-1$.

###### Proof.

The Mellin transform of $F(x)$ is given by

$\int_{0}^{1}F(x)x^{s-1}dx=\sum_{j=1}^{q}\sum_{k=0}^{n_{j}-1}A_{j}[k]\int_{0}^{1}x^{k+s-1}w_{j}(x)dx+\sum_{j=1}^{p-q}\sum_{k=0}^{m_{j}-1}A_{q+j}[k]\int_{0}^{1}x^{k+s-1}v_{j}(x)dx,$

where

$\int_{0}^{1}x^{k+s-1}w_{j}(x)dx=\frac{\Gamma(s+k+\vec{a})}{\Gamma(s+k+\vec{b}+\vec{e}_{j})}=\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b})}\frac{(s+\vec{a})_{k}}{(s+\vec{b})_{k}(s+b_{j}+k)}$

and

$\int_{0}^{1}x^{k+s-1}v_{j}(x)dx=\frac{\Gamma(s+k+\vec{a})}{\Gamma(s+k+\vec{b})}(k+s-1)^{j-1}=\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b})}\frac{(s+\vec{a})_{k}}{(s+\vec{b})_{k}}(k+s-1)^{j-1}.$

It then remains to show that

$\sum_{j=1}^{q}\sum_{k=0}^{n_{j}-1}A_{j}[k]\frac{(s+\vec{a})_{k}}{(s+\vec{b})_{k}(s+b_{j}+k)}+\sum_{j=1}^{p-q}\sum_{k=0}^{m_{j}-1}A_{q+j}[k]\frac{(s+\vec{a})_{k}}{(s+\vec{b})_{k}}(k+s-1)^{j-1}=\frac{p_{\left|\vec{n}\right|+\left|\vec{m}\right|-1}(s)}{\prod_{i=1}^{q}(s+b_{i})_{n_{i}}}.$

Multiplying the left-hand side by $\prod_{i=1}^{q}(s+b_{i})_{n_{i}}$ yields

$\sum_{j=1}^{q}\sum_{k=0}^{n_{j}-1}A_{j}[k](s+\vec{a})_{k}\prod_{i=1}^{q}\frac{(s+b_{i})_{n_{i}}}{(s+b_{i})_{k+\delta_{i,j}}}+\sum_{j=1}^{p-q}\sum_{k=0}^{m_{j}-1}A_{q+j}[k](s+\vec{a})_{k}\prod_{i=1}^{q}\frac{(s+b_{i})_{n_{i}}}{(s+b_{i})_{k}}(k+s-1)^{j-1},$

which is a polynomial if $n_{i}\geq n_{j}-1+\delta_{i,j}$ and $n_{i}\geq
m_{j}-1$ for all $i,j$. The polynomial produced by the first terms has a
degree of at most $\left|\vec{n}\right|+(p-q)\max\\{n_{j}-1\mid j\\}-1$, which
is $\leq\left|\vec{n}\right|+\left|\vec{m}\right|-1$ if $m_{j}\geq n_{i}-1$
for all $i,j$. The last two sets of conditions lead exactly to the near
diagonal condition on $\vec{n}\sqcup\vec{m}$. The remaining terms produce a
polynomial of degree at most $\left|\vec{n}\right|+\max\\{(p-q)(m_{j}-1)+j\mid
j\\}-1$. The step-line condition on $\vec{m}$ then implies that this value is
$\leq\left|\vec{n}\right|+\left|\vec{m}\right|-1$. Indeed, suppose that
$m_{j}=m$ for $j\leq J$ and $m_{j}=m-1$ for $j>J$, then
$\left|\vec{m}\right|=Jm+(p-q-J)(m-1)=(p-q)(m-1)+J$. If $j\leq J$, then
$(p-q)(m-1)+j\leq\left|\vec{m}\right|$. If $J<j\leq p-q$, then
$(p-q)(m-2)+j\leq(p-q)(m-1)+j-(p-q)\leq\left|\vec{m}\right|$. ∎

Similarly as in the proof of Proposition 2.3, the composition formula (4) and
the fact that certain Jacobi-like systems of weights are AT-systems allow us
to show that the same holds for certain Laguerre-like systems of weights.

######  Proposition 3.3.

Suppose that all $b_{j}-a_{j}\in\mathbb{Z}_{\geq
0}\cup(\left|\vec{n}\right|+\left|\vec{m}\right|-1,\infty)$. Then the set of
weights $\\{w_{j}(x)\\}_{j=1}^{q}\cup\\{v_{j}(x)\\}_{j=1}^{p-q}$ is an AT-
system on $(0,\infty)$ for every $\vec{n}\sqcup\vec{m}\in\mathcal{N}^{p}$ with
$\vec{m}\in\mathcal{S}^{p-q}$.

###### Proof.

Suppose that $F(x)$ is a polynomial combination of the described Laguerre-like
weights that corresponds to the multi-index $\vec{n}\sqcup\vec{m}$. Its Mellin
transform $\hat{F}(s)$ is then of the form

$\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+\vec{n})}p_{\left|\vec{n}\right|+\left|\vec{m}\right|-1}(s),$

see Lemma 3.2. We can decompose $\hat{F}(s)$ as
$\hat{F}(s)/\Gamma(s+\vec{c}+\vec{m})\cdot\Gamma(s+\vec{c}+\vec{m})$ for some
$\vec{c}\in(-1,\infty)^{p-q}$ with
$c_{j}-a_{q+j}>\left|\vec{n}\right|+\left|\vec{m}\right|-1$ for all $1\leq
j\leq p-q$. The first factor corresponds to the Mellin transform of a
polynomial combination of the Jacobi-like weights in
$\\{w_{j}(x;\vec{a},\vec{b}\sqcup\vec{c})\\}_{j=1}^{p}$, see Lemma 2.2. We
know from Proposition 2.3 that the latter is an AT-system for the multi-index
$\vec{n}\sqcup\vec{m}$, hence the extended set of weights is a T-system. The
second factor is the Mellin convolution of the gamma densities
$\mathcal{G}^{c_{j}+m_{j}}(x)$, $1\leq j\leq p-q$. It is therefore sufficient
to show that taking the consecutive convolution with the stated gamma
densities preserves the property of being a T-system. Convoluting with a gamma
density $\mathcal{G}^{a}(x)$ corresponds to taking
$K(x,y)=(x/y)^{a}\exp(x/y)/y$ in the composition formula (4). The desired
result then follows from the fact that the kernel $\exp(xy)$ is strictly
totally positive of any order, see, e.g., [14, Example 1 - Section 1.3]. ∎

######  Corollary 3.4.

Suppose that $b_{j}-a_{j}\in\mathbb{Z}_{\geq 0}$ for all $j$. Then the set of
weights $\\{w_{j}(x)\\}_{j=1}^{q}\cup\\{v_{j}(x)\\}_{j=1}^{p-q}$ is an AT-
system on $(0,\infty)$ for every multi-index
$\vec{n}\sqcup\vec{m}\in\mathcal{N}^{p}$ with $\vec{m}\in\mathcal{S}^{p-q}$.

Note that, compared to Corollary 2.4 (the analogous result in the Jacobi-like
setting), there is no condition of the form $b_{i}-b_{j}\not\in\mathbb{Z}$
whenever $i\neq j$. Such a condition was required in the former because the
weights could reduce to the Jacobi-Piñeiro weights
$(x^{b_{1}},\dots,x^{b_{r}})$ on $[0,1]$, but this is not the case here.

We can again derive a matrix Pearson equation for the system of weights. For
the gamma density $\mathcal{G}^{a}(x)$, we have a Pearson equation
$[\sigma(x)w(x)]^{\prime}=\tau(x)w(x),$ with $\sigma(x)=x$ and
$\tau(x)=(a+1)x$. Due to the limiting nature of the weights, the structure of
the underlying matrix will be more involved than in the Jacobi-like setting.

######  Proposition 3.5.

Suppose that $b_{i}\neq b_{j}$ if $i\neq j$. Define
$\vec{W}=(w_{1},\dots,w_{q},v_{1},\dots,v_{p-q})^{T}$. Then

$[x\vec{W}(x)]^{\prime}=\mathcal{T}(x)\vec{W}(x),\quad\mathcal{T}(x)=\begin{pmatrix}\mathcal{T}_{0}&\vec{1}\,^{T}&0\\\
0&\vec{0}^{T}&I_{p-q-1}\\\ \vec{c}&x+d_{0}&\vec{d}\end{pmatrix},$

where $\mathcal{T}_{0}=-\text{diag}(b_{j}+1\mid 1\leq j\leq q)$,
$c_{k}=-(\vec{a}-b_{k})_{1}/(\vec{b}^{\ast k}-b_{k})_{1}$ for $1\leq k\leq q$
and
$d_{k}=-(d^{k}/ds^{k})[(s+\vec{a})_{1}/(s+\vec{b})_{1}]_{s=1}/k!+(-1)^{k+1}\sum_{j=1}^{q}c_{j}/(b_{j}+1)^{k+1}$
for $0\leq k\leq p-q-1$.

###### Proof.

First, note that $(s+b_{j})\hat{w}_{j}(s)=\hat{v}_{1}(s)$ and thus
$(xw_{j}(x))^{\prime}=-(b_{j}+1)w_{j}(x)+v_{1}(x)$ for all $1\leq j\leq q$. On
the other hand, by definition, we have $(xv_{j}(x))^{\prime}=v_{j+1}(x)$ for
all $1\leq j\leq p-q-1$. It then remains to find the relation between
$(xv_{p-q}(x))^{\prime}$ and the other weights. Observe that

$\hat{v}_{1}(s+1)=\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b})}\frac{(s+\vec{a})_{1}}{(s+\vec{b})_{1}},$

in terms of a rational function that has a numerator of degree $p>q$ and a
denominator of degree $q$. Hence, after a partial fraction decomposition, we
find

$\frac{(s+\vec{a})_{1}}{(s+\vec{b})_{1}}=-\sum_{j=1}^{q}\frac{c_{j}}{s+b_{j}}-\sum_{j=1}^{p-q}d_{j-1}(s-1)^{j-1}+(s-1)^{p-q}.$

Consequently,

$(s-1)\hat{v}_{p-q}(s)=\sum_{j=1}^{q}c_{j}\hat{w}_{j}(s)+\hat{v}_{1}(s+1)+\sum_{j=1}^{p-q}d_{j-1}\hat{v}_{j}(s),$

and thus after taking the inverse Mellin transform

$(xv_{p-q}(x))^{\prime}=\sum_{j=1}^{q}c_{j}w_{j}(x)+(x+d_{0})v_{1}(x)+\sum_{j=2}^{p-q}d_{j-1}v_{j}(x).$

The coefficients in the partial fraction decomposition can be computed as
follows. Comparing residues in $s=-b_{k}$, we obtain
$c_{k}=-(\vec{a}-b_{k})_{1}/(\vec{b}^{\ast k}-b_{k})_{1}$. On the other hand,
by applying $(d^{k}/ds^{k})[\,\cdot\,]_{s=1}/k!$, we get the stated expression
for $d_{k}$. ∎

###  3.2 Type I

In this section, we will investigate the type I multiple orthogonal
polynomials associated with the weights
$(w_{1},\dots,w_{q},v_{1},\dots,v_{p-q})$ introduced in the previous section
for multi-indices $\vec{n}\sqcup\vec{m}$ near the diagonal with $\vec{m}$ on
the step-line. Similarly as in the Jacobi-like setting, we will first derive
an explicit expression for the Mellin transform of the corresponding type I
functions. Essentially, this result also follows from Theorem 2.6 and the
limiting relations between the weights; to remove the last component from
$\vec{b}$ to get $\vec{b}^{\ast q}$, we can do the change of variables
$x\mapsto x/b_{q}$ and let $b_{q}\to\infty$.

######  Theorem 3.6.

Let $\vec{n}\sqcup\vec{m}\in\mathcal{N}^{p}$ with
$\vec{m}\in\mathcal{S}^{p-q}$. Every type I function
$F_{\vec{n}\sqcup\vec{m}}$ has a Mellin transform of the form

$\mathcal{F}_{\vec{n}\sqcup\vec{m}}\cdot\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+\vec{n})}(1-s)_{\left|\vec{n}\right|+\left|\vec{m}\right|-1},\quad\mathcal{F}_{\vec{n}\sqcup\vec{m}}\in\mathbb{R}.$

###### Proof.

We know from Lemma 3.2 that the Mellin transform of $F_{\vec{n}\sqcup\vec{m}}$
is of the form

$\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+\vec{n})}p_{\left|\vec{n}\right|+\left|\vec{m}\right|-1}(s),$

in terms of a polynomial $p_{\left|\vec{n}\right|+\left|\vec{m}\right|-1}(s)$
of degree at most $\left|\vec{n}\right|+\left|\vec{m}\right|-1$. Due to the
orthogonality conditions, we can conclude that this polynomial must be given
by a scalar multiple of $(1-s)_{\left|\vec{n}\right|+\left|\vec{m}\right|-1}$.
∎

It is possible to be more generous with the restrictions on the multi-index if
$a_{j}=b_{j}$ for some $j$, for example, in the exponential integral case
(where $(p,q)=(2,1)$, $\vec{a}=(\alpha,\alpha+\nu)$ and
$\vec{b}=(\alpha+\nu$)) the above holds for multi-indices $(n,m)$ with
$m+1\geq n$ (near and above the diagonal), see [43, Lemma 2.1]. However, we
won’t go into further detail on this here.

In what follows, we will normalize the (non-zero) type I functions
$F_{\vec{n}\sqcup\vec{m}}(x)$ such that
$\mathcal{F}_{\vec{n}\sqcup\vec{m}}=1$.

Taking the inverse Mellin transform immediately yields an integral
representation for the type I functions.

######  Corollary 3.7.

Let $\vec{n}\sqcup\vec{m}\in\mathcal{N}^{p}$ with
$\vec{m}\in\mathcal{S}^{p-q}$. Then,

$F_{\vec{n}\sqcup\vec{m}}(x)=\frac{1}{2\pi
i}\int_{c-i\infty}^{c+i\infty}\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+\vec{n})}(1-s)_{\left|\vec{n}\right|+\left|\vec{m}\right|-1}x^{-s}ds,$

for some $c>0$ with $-\min\\{a_{j}\mid j\\}<c<1$.

We can derive the following Rodrigues-type formula from Theorem 3.6.

######  Corollary 3.8.

Let $\vec{n}\sqcup\vec{m}\in\mathcal{N}^{p}$ with
$\vec{m}\in\mathcal{S}^{p-q}$. Then,

$F_{\vec{n}\sqcup\vec{m}}(x)=\frac{d^{\left|\vec{n}\right|+\left|\vec{m}\right|-1}}{dx^{\left|\vec{n}\right|+\left|\vec{m}\right|-1}}\left[x^{\left|\vec{n}\right|+\left|\vec{m}\right|-1}w_{0}(x;\vec{a},\vec{b}+\vec{n})\right].$

The type I functions can be expressed as certain Mellin convolutions in which
now gamma densities play a role.

######  Proposition 3.9.

Let $\vec{n}\sqcup\vec{m}\in\mathcal{N}^{p}$ with
$\vec{m}\in\mathcal{S}^{p-q}$. Denote $\vec{a_{\ast}}=(a_{q+1}\dots,a_{p})$.
Then,

$F_{\vec{n}\sqcup\vec{m}}=F_{\vec{n}\sqcup\vec{m}}^{(\vec{a},\vec{b}\sqcup\vec{a}_{\ast})}\ast\mathcal{G}^{m_{1}+a_{q+1}}\ast\dots\ast\mathcal{G}^{m_{p-q}+a_{p}},$

in terms the type I function
$F_{\vec{n}\sqcup\vec{m}}(x;\vec{a},\vec{b}\sqcup\vec{a}_{\ast})$ from the
Jacobi-like setting.

###### Proof.

The result follows from the multiplicative property of the Mellin convolution
and the uniqueness of the Mellin transform; the Mellin transform of the type I
functions in the Jacobi-like setting are given in Theorem 2.6, the Mellin
transform of a gamma density $\mathcal{G}^{a}(x)$ is $\Gamma(s+a)$. ∎

By working out the contour integral in Corollary 3.7 via the residue theorem,
we can again find an expansion of the type I functions in terms of certain
hypergeometric series (see Proposition 2.11).

######  Proposition 3.10.

Let $\vec{n}\sqcup\vec{m}\in\mathcal{N}^{p}$ with
$\vec{m}\in\mathcal{S}^{p-q}$ and suppose that $a_{i}-a_{j}\not\in\mathbb{Z}$
whenever $i\neq j$ and $b_{i}-a_{j}\not\in\mathbb{Z}$ for all $i,j$. Set
$N=\left|\vec{n}\right|+\left|\vec{m}\right|$. Then,

$F_{\vec{n}\sqcup\vec{m}}(x)=\sum_{j=1}^{p}\frac{(a_{j}+1)_{N-1}\Gamma(\vec{a}^{\ast
j}-a_{j})}{\Gamma(\vec{b}-a_{j}+\vec{n})}{}_{p+1}F_{q}\left(\begin{array}[]{c}N+a_{j},-\vec{n}-\vec{b}+a_{j}+1\\\
a_{j}+1,-\vec{a}^{\ast j}+a_{j}+1\end{array};x\right)x^{a_{j}}.$

The terms that appear (including the $x^{a_{j}}$) are all solutions to the
differential equation for

${}_{q+1}F_{p}\left(\begin{array}[]{c}\left|\vec{n}\right|+\left|\vec{m}\right|,-\vec{n}-\vec{b}+1\\\
-\vec{a}+1\end{array};x\right),$

see [8, Eq. 16.8.6]. The type I function itself therefore satisfies the same
differential equation.

Unfortunately, we can’t use the same strategy as in the Jacobi-like setting to
obtain formulas for the type I polynomials themselves, because the series for
the error in the underlying Hermite-Padé approximation isn’t convergent
anymore (though formally it would still make sense). In principle, one should
be able to describe the polynomials inductively from the expressions found in
the Jacobi-like setting by making use of the limiting relations that are
induced by the weights. Since doing so would lead to expressions that are
rather involved, we won’t do it here explicitly.

###  3.3 Type II

In what follows, we will study the type II multiple orthogonal polynomials
associated with the weights $(w_{1},\dots,w_{q},v_{1},\dots,v_{p-q})$
introduced in Section 3.1 for multi-indices $\vec{n}\sqcup\vec{m}$ near the
diagonal with $\vec{m}$ on the step-line. Formulas for these type II
polynomials for multi-indices on the step-line were already obtained in [24],
but we will derive them explicitly from the orthogonality conditions and for a
larger set of multi-indices. They also follow immediately from Theorem 2.13
and the limiting relations; to go from $\vec{b}$ to a vector $\vec{b}^{\ast
q}$ in which the last component is missing, we do the change of variables
$x\mapsto x/b_{q}$ and let $b_{q}\to\infty$. Afterwards, we will describe the
asymptotic zero distribution through free probability.

######  Theorem 3.11.

Let $\vec{n}\sqcup\vec{m}\in\mathcal{N}^{p}$ with
$\vec{m}\in\mathcal{S}^{p-q}$. Every type II polynomial
$P_{\vec{n}\sqcup\vec{m}}(x)$ is of the form

$\mathcal{P}_{\vec{n}\sqcup\vec{m}}\cdot{}_{q+1}F_{p}\left(\begin{array}[]{c}-\left|\vec{n}\right|-\left|\vec{m}\right|,\vec{b}+\vec{n}+1\\\
\vec{a}+1\end{array};x\right),\quad\mathcal{P}_{\vec{n}\sqcup\vec{m}}\in\mathbb{R}.$

###### Proof.

Similarly as in the Jacobi-like setting (see Theorem 2.13), it will be more
natural to move into the complex plane and write

$P_{\vec{n}\sqcup\vec{m}}(x)=\frac{1}{2\pi
i}\sum_{k=0}^{\left|\vec{n}\right|+\left|\vec{m}\right|}\int_{\Sigma}\frac{p_{\left|\vec{n}\right|+\left|\vec{m}\right|}(t)}{(-t)_{\left|\vec{n}\right|+\left|\vec{m}\right|+1}}\frac{\Gamma(t+\vec{b}+1)}{\Gamma(t+\vec{a}+1)}x^{t}dt,$

in terms of a contour $\Sigma$ encircling the interval
$[0,\left|\vec{n}\right|+\left|\vec{m}\right|]$ avoiding the poles of the
integrand that are not in this interval and a polynomial
$p_{\left|\vec{n}\right|+\left|\vec{m}\right|}(t)$ of degree at most
$\left|\vec{n}\right|+\left|\vec{m}\right|$ that is determined by the relation

$\frac{p_{\left|\vec{n}\right|+\left|\vec{m}\right|}(t)}{(-t)_{\left|\vec{n}\right|+\left|\vec{m}\right|+1}}=\sum_{k=0}^{\left|\vec{n}\right|+\left|\vec{m}\right|}\frac{\Gamma(k+\vec{a}+1)}{\Gamma(k+\vec{b}+1)}\frac{P_{\vec{n}+\left|\vec{m}\right|}[k]}{t-k}.$

The orthogonality conditions imply that

$\int_{\Sigma}\frac{(t+\vec{a}+1)_{s}}{(t+\vec{b}+1)_{s}(s+t+b_{j}+1)}\frac{p_{\left|\vec{n}\right|+\left|\vec{m}\right|}(t)}{(-t)_{\left|\vec{n}\right|+\left|\vec{m}\right|+1}}dt=0,\quad
s=0,\dots,n_{j}-1,$
$\int_{\Sigma}\frac{(t+\vec{a}+1)_{s}(s+t)^{j-1}}{(t+\vec{b}+1)_{s}}\frac{p_{\left|\vec{n}\right|+\left|\vec{m}\right|}(t)}{(-t)_{\left|\vec{n}\right|+\left|\vec{m}\right|+1}}dt=0,\quad
s=0,\dots,m_{j}-1.$

As a consequence, since $\vec{n}\sqcup\vec{m}\in\mathcal{N}^{p}$ with
$\vec{m}\in\mathcal{S}^{p-q}$, we must have

$\int_{\Sigma}\frac{p_{\left|\vec{n}\right|+\left|\vec{m}\right|}(t)}{\prod_{j=1}^{q}(t+b_{j})_{n_{j}}}\frac{q_{\left|\vec{n}\right|+\left|\vec{m}\right|-1}(t)}{(-t)_{\left|\vec{n}\right|+\left|\vec{m}\right|+1}}x^{t}dt=0,$

for every polynomial $q_{\left|\vec{n}\right|+\left|\vec{m}\right|-1}(t)$ of
degree at most $\left|\vec{n}\right|+\left|\vec{m}\right|-1$. By considering
the polynomials that have roots in all possible
$\left|\vec{n}\right|+\left|\vec{m}\right|-1$ roots of
$(-t)_{\left|\vec{n}\right|+\left|\vec{m}\right|+1}$ and applying the residue
theorem, we can conclude that
$p_{\left|\vec{n}\right|+\left|\vec{m}\right|}(t)/\prod_{j=1}^{q}(t+b_{j})_{n_{j}}$
is constant along the $\left|\vec{n}\right|+\left|\vec{m}\right|+1$ roots of
$(-t)_{\left|\vec{n}\right|+\left|\vec{m}\right|+1}$. We therefore must have
that $p_{\left|\vec{n}\right|+\left|\vec{m}\right|}(t)$ is a scalar multiple
of $\prod_{j=1}^{q}(t+b_{j}+1)_{n_{j}}$, from which the desired result then
follows. ∎

From now on, we will assume that the (non-zero) type II polynomials
$P_{\vec{n}\sqcup\vec{m}}(x)$ are normalized such that
$\mathcal{P}_{\vec{n}\sqcup\vec{m}}=1$.

The proof of Theorem 3.11 provides the following contour integral
representation for the type II polynomials. Observe again the duality with
Corollary 3.7.

######  Corollary 3.12.

Let $\vec{n}\in\mathcal{N}^{p}$ with $\vec{m}\in\mathcal{S}^{p-q}$. Then,

$P_{\vec{n}\sqcup\vec{m}}(x)=\frac{1}{2\pi
i}\int_{\Sigma}\frac{1}{(-t)_{\left|\vec{n}\right|+\left|\vec{m}\right|+1}}\frac{\Gamma(t+\vec{b}+\vec{n}+1)}{\Gamma(t+\vec{a}+1)}x^{t}dt,$

where $\Sigma$ is a positively oriented contour that starts and ends at
$\infty$ and encircles the positive real line, without enclosing any poles of
the integrand not in $[0,\infty)$.

Some properties of the zeros of these type II polynomials for multi-indices on
the step-line were already obtained in [24] by exploiting the connection with
branched continued fractions. We will give a proof by making use of Theorem
2.15 (the analogue in the Jacobi-like setting) and the finite free
multiplicative convolution. Under the condition that all
$b_{j}-a_{j}\in\mathbb{Z}_{\geq
0}\cup(\left|\vec{n}\right|+\left|\vec{m}\right|-1,\infty)$, these properties
also follow from the fact that the underlying system of weights is an AT-
system (see Proposition 3.3).

######  Theorem 3.13 (Lima [24], Theorem 4.8).

Let $\vec{n}\sqcup\vec{m}\in\mathcal{S}^{p}$ and suppose that $b_{j}>a_{i}$ if
$i\leq j$ and $b_{j}>a_{i}-1$ if $i>j$. The type II polynomials
$P_{\vec{n}\sqcup\vec{m}}(x)$ have real, positive, simple zeros.

Alternative proof of Theorem 3.13. Set
$N=\left|\vec{n}\right|+\left|\vec{m}\right|$ and denote
$\vec{a}_{\ast}=(a_{q+1},\dots,a_{p})$, then we have the following
decomposition

$\displaystyle{}_{q+1}F_{p}$
$\displaystyle\left(\begin{array}[]{c}-N,\vec{b}+\vec{n}+1\\\
\vec{a}+1\end{array};x\right)={}_{p+1}F_{p}\left(\begin{array}[]{c}-N,\vec{b}+\vec{n}+1,\vec{a}_{\ast}+\vec{m}+1\\\
\vec{a}+1\end{array};x\right)$ (16)
$\displaystyle\boxtimes_{N}{}_{1}F_{1}\left(\begin{array}[]{c}-N\\\
a_{p-q+1}+m_{1}+1\end{array};x\right)\boxtimes_{N}\cdots\boxtimes_{N}{}_{1}F_{1}\left(\begin{array}[]{c}-N\\\
a_{p}+m_{p-q}+1\end{array};x\right),$ (21)

The first factor is
$P_{\vec{n}\sqcup\vec{m}}(x;\vec{a},\vec{b}\sqcup\vec{a}_{\ast})$, a type II
polynomial in the Jacobi-like setting, and the others are Laguerre polynomials

${}_{1}F_{1}\left(\begin{array}[]{c}-N\\\
a_{q+j}+m_{j}+1\end{array};x\right)=L_{N}^{(a_{q+j}+m_{j})}(x).$

Since each factor has real, positive, simple roots, the same holds for their
finite free convolution (see [40] and [15, Theorem 6]). ∎

We will now describe the asymptotic zero distribution of the (scaled) type II
polynomials as a free convolution of measures. It will be more natural to use
a different decomposition than (16) to achieve this. In principle, we could
also use (16), but then the expression for the asymptotic zero distribution
would contain more factors than necessary and a later interpretation in terms
of random matrices would be lost.

######  Theorem 3.14.

Suppose that $b_{j}>a_{i}$ if $i\leq j$ and $b_{j}>a_{i}-1$ if $i>j$. The
asymptotic zero distribution of $P_{\vec{n}\sqcup\vec{m}}(N^{p-q}x)$ for
$\vec{n}\sqcup\vec{m}\in\mathcal{S}^{p}$ as
$N:=\left|\vec{n}\right|+\left|\vec{m}\right|\to\infty$ is given by
$((1-1/p)\delta_{1}+(1/p)\mu_{p-1})^{\boxtimes
q}\boxtimes\nu^{\boxtimes(p-q)}$ in terms of the deformed arcsin distribution
$\mu_{p-1}$ from Lemma 2.18 and the Marchenko-Pastur distribution $\nu$.

###### Proof.

It is clear from (16) that the described scaling is indeed appropriate: the
zeros of the free multiplicative convolution of polynomials are bounded by the
product of the largest zeros of its factors (see [40]) and each factor $N$ in
the scaling corresponds to an appearance of a Laguerre polynomial.
Additionally, it follows that the asymptotic zero distribution is independent
of the precise value of the parameters. We can therefore assume that
$\vec{b}=\vec{a}^{\ast}$ with $\vec{a}^{\ast}=(a_{1},\dots,a_{p})$. The
decomposition

$\displaystyle{}_{q+1}F_{p}\left(\begin{array}[]{c}-N,\vec{a}+\vec{n}+1\\\
\vec{a}+1\end{array};N^{p-q}x\right)$
$\displaystyle={}_{2}F_{1}\left(\begin{array}[]{c}-N,a_{1}+n_{1}+1\\\
a_{1}+1\end{array};x\right)\boxtimes_{N}\cdots\boxtimes_{N}{}_{2}F_{1}\left(\begin{array}[]{c}-N,a_{q}+n_{q}+1\\\
a_{q}+1\end{array};x\right)$
$\displaystyle\quad\boxtimes_{N}{}_{1}F_{1}\left(\begin{array}[]{c}-N\\\
a_{q+1}+1\end{array};Nx\right)\boxtimes_{N}\cdots\boxtimes_{N}{}_{1}F_{1}\left(\begin{array}[]{c}-N\\\
a_{p}+1\end{array};Nx\right),$

in terms of the Jacobi polynomials
$P_{N}^{(a_{j},n_{j}-N)}(x)=(1-x)^{N-n_{j}}P_{n_{j}}^{(a_{j},N-n_{j})}(x)$ and
scaled Laguerre polynomials $L_{N}^{(a_{q+j})}(Nx)$, together with Lemma 2.18
and the fact that the Marchenko-Pastur is the asymptotic zero distribution of
the scaled Laguerre polynomials, then leads to the desired result. ∎

By making use of the previous result, we can derive an algebraic equation for
the Stieltjes transform of the asymptotic zero distribution. Alternatively, we
can start from the hypergeometric differential equation for the type II
polynomials, see [8, Eq. (16.8.3)] and employ the usual techniques (see, e.g.,
[23, Section 3.1.4] or [43, Section 3.4.2]). Generally, the described
algebraic equation doesn’t seem to be solvable explicitly.

######  Corollary 3.15.

Suppose that $b_{j}>a_{i}$ if $i\leq j$ and $b_{j}>a_{i}-1$ if $i>j$. The
Stieltjes transform of the asymptotic zero distribution of
$P_{\vec{n}\sqcup\vec{m}}(N^{p-q}x)$ as
$N:=\left|\vec{n}\right|+\left|\vec{m}\right|\to\infty$ satisfies the
algebraic equation $(zS(z))^{p+1}=z(zS(z)+1/p)^{q}(zS(z)-1)$.

###### Proof.

We can prove this by working out the $S$-transform of the asymptotic zero
distribution in terms of the $S$-transforms of its factors
$(1-1/p)\delta_{1}+(1/p)\mu_{p-1}$ and $\nu$. Afterwards, we can use the fact
that a Stieltjes transform $S(z)=\sum_{k\geq 0}m_{k}z^{-k-1}=(1+m(1/z))/z$ and
$S$-transform $s(z)=m^{-1}(z)(z+1)/z$ are connected through the moment
generating function $m(z)=\sum_{k\geq 0}m_{k}z^{k}$. The $S$-transform of
$(1-1/p)\delta_{1}+(1/p)\mu_{p-1}$ and $\nu$ is given by
$s_{1}(z)=1+1/(p(z+1))$ and $s_{2}(z)=1/(z+1)$ respectively. This follows from
the algebraic equations for the corresponding Stieltjes transforms, which can
be obtained from the hypergeometric differential equation for
$P_{n}^{(a,(p-1)n)}(x)$ and $L_{n}^{(a)}(nx)$ in the standard way. ∎

##  4 Bessel-like setting

###  4.1 System of weights

Let $p<q$ and suppose that $\vec{a}\in(-1,\infty)^{p}$ and
$\vec{b}\in(-1,\infty)^{q}$ satisfy $a_{j}<b_{j}$ for $1\leq j\leq p$. We will
consider the $q$ weights
$w_{j}(z;\vec{a},\vec{b})=f_{0}(z;\vec{a},\vec{b}+\vec{e}_{j})$, with $f_{0}$
as in (1), on the unit circle $\\{z\in\mathbb{C}:|z|\,=1\\}$. In this setting,
such weights are convergent for every non-zero $z\in\mathbb{C}$. It follows
from properties of the Mellin transform that there is no weight
$w_{j}(z;\vec{a},\vec{b})$ on the positive real line with a Mellin transform
that leads to moments of the form

$\int_{|z|=1}z^{s-1}w_{j}(z;\vec{a},\vec{b})dz=\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+\vec{e}_{j})}.$

The limiting relations between these kinds of weights are straightforward: we
may remove the last component in the vector $\vec{a}$ via

$\lim_{a_{q}\to\infty}\frac{1}{\Gamma(a_{q})}w_{j}(z;\vec{a},\vec{b})=w_{j}(z;\vec{a}^{\ast
p},\vec{b}).$

We can also use it to convert a set of Jacobi-like weights from Section 2.1 to
a set of Bessel-like weights. This is because the Jacobi-like weights can also
be interpreted as weights on the unit circle. Indeed, if
$\vec{a},\vec{b}\in(-1,\infty)^{r}$, we can define
$w_{j}(z;\vec{a},\vec{b})=f_{0}(1/z;\vec{a},\vec{b}+\vec{e}_{j})/z$ (it is
convergent in $\\{z\in\mathbb{C}:0<\left|z\right|\leq 1\\}$) and still have
the desired moments.

Since we are working with weights on the unit circle, the notion of an AT-
system isn’t well-defined anymore. The matrix Pearson equation for this set of
weights is similar to the one in the Jacobi-like setting, see Proposition 2.5.

######  Proposition 4.1.

Define $\vec{W}=(w_{1},\dots,w_{q})^{T}$. Then

$[z(1-z)\vec{W}(x)]^{\prime}=\mathcal{T}(z)\vec{W}(z),$

where

$\mathcal{T}_{k,j}(z)=-\frac{\prod_{i=1}^{p}(a_{i}-b_{j})}{\prod_{i\neq
j}^{q}(b_{i}-b_{j})}+\begin{dcases}0,\ k\neq j,\\\ z-(b_{j}+1)(1-z),\
k=j.\end{dcases}$

###### Proof.

We can use the same strategy as in the Jacobi-like setting (Proposition 2.5).
The results obtained for the Mellin transforms can be interpreted as results
about the coefficients in the expansions of the weights here. ∎

###  4.2 Type I polynomials

Below, we will study the type I multiple orthogonal polynomials associated
with the weights $(w_{1},\dots,w_{q})$ from the previous section for multi-
indices $\vec{n}$ near the diagonal. Since, in this setting, the notion of
Mellin transform is not well-defined anymore, we will only be able to compute
the moments of the type I functions. However, this is sufficient to obtain
formulas for the associated type I polynomials. Other results, like the
differential formula, convolution formula and hypergeometric expansion don’t
necessarily follow anymore.

We will first compute the moments of polynomial combinations of the weights.

######  Lemma 4.2.

Let $\vec{n}\in\mathcal{N}^{q}$. The $(s-1)$-th moment of
$F(z)=\sum_{j=1}^{q}A_{j}(z)w_{j}(z)$ with $\deg A_{j}\leq n_{j}-1$ is of the
form

$\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+\vec{n})}p_{\left|\vec{n}\right|-1}(s),$

where $p_{\left|\vec{n}\right|-1}(s)$ is a polynomial of degree at most
$\left|\vec{n}\right|-1$.

###### Proof.

We can follow the proof of Lemma 2.2 (the analogue in the Jacobi-like
setting). The degree of the polynomial in (3) is at most
$\left|\vec{n}\right|-1$, because the degree of the term that corresponds to
the indices $(j,k)$ is $k(p-q)+\left|\vec{n}\right|-1$ and this is at most
$\left|\vec{n}\right|-1$ for $0\leq k\leq n_{j}-1$ since $p<q$. ∎

We can obtain an explicit expression for the moments of the type I functions
in the same way as in the Jacobi-like setting, see Theorem 2.6. It is also
implied by the limiting relations between the weights: we may remove the last
component from $\vec{a}$ to get $\vec{a}^{\ast p}$ by doing the change of
variables $a_{p}\mapsto a_{p}z$ and letting $a_{p}\to\infty$.

######  Theorem 4.3.

Let $\vec{n}\in\mathcal{N}^{q}$. The $(s-1)$-th moment of every type I
function $F_{\vec{n}}$ is of the form

$\mathcal{F}_{\vec{n}}\cdot\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+\vec{n})}(1-s)_{\left|\vec{n}\right|-1},\quad\mathcal{F}_{\vec{n}}\in\mathbb{R}.$

In what follows, we will assume that the (non-zero) type I functions
$F_{\vec{n}}(x)$ are normalized such that $\mathcal{F}_{\vec{n}}=1$.

Theorem 4.3 allows us to obtain explicit expressions for the type I
polynomials.

######  Theorem 4.4.

Let $\vec{n}\in\mathcal{N}^{q}$ and suppose that
$b_{j}-a_{i}\not\in\mathbb{Z}_{<0}$ and $b_{i}-b_{j}\not\in\mathbb{Z}_{\geq
0}$ whenever $i\neq j$. The type I polynomials $A_{\vec{n},j}$ are given by

$A_{\vec{n},j}(z)=\frac{(\vec{a}-b_{j})_{1}}{(\vec{b}-b_{j})_{1}^{\ast}}\sum_{J=1}^{q}\sum_{K=0}^{n_{J}-1}P_{\vec{n},J}[K]\sum_{k=0}^{K-1+\delta_{j,J}}\frac{(\vec{b}-b_{J}-K)_{k+1}^{\ast
j}(b_{j}-b_{J}-K)_{k}}{(\vec{a}-b_{J}-K)_{k+1}}z^{k},$

where

$P_{\vec{n},J}[K]=\frac{(b_{J}+K+1)_{\left|\vec{n}\right|-1}}{\prod_{i=1,i\neq
J}^{q}(b_{i}-b_{J}-K)_{n_{i}}}\frac{(-1)^{K}}{K!(n_{J}-K-1)!}.$

###### Proof.

The proof is analogous to that in the Jacobi-like setting, see Theorem 2.12.
The essential partial fraction decomposition (9) still goes through. ∎

######  Corollary 4.5.

Suppose that $b_{j}-a_{i}\not\in\mathbb{Z}_{<0}$ and
$b_{i}-b_{j}\not\in\mathbb{Z}_{\geq 0}$ whenever $i\neq j$. Then every multi-
index near the diagonal is normal for the set of weights
$\\{w_{j}(z)\\}_{j=1}^{q}$.

###  4.3 Type II polynomials

In this section, we will determine formulas for the type II multiple
orthogonal polynomials associated with the weights $(w_{1},\dots,w_{q})$
introduced in Section 4.1 for multi-indices near the diagonal. These formulas
were already obtained in [24], but we will provide a constructive proof that
holds for a larger set of multi-indices. They immediately follow from the
limiting relations as well; to go from $\vec{a}$ to $\vec{a}^{\ast p}$ and
remove the last component, we do the change of variables $z\mapsto a_{p}z$ and
let $a_{p}\to\infty$.

######  Theorem 4.6.

Let $\vec{n}\in\mathcal{N}^{q}$. Every type II polynomial $P_{\vec{n}}(z)$ is
of the form

$\mathcal{P}_{\vec{n}}\cdot{}_{q+1}F_{p}\left(\begin{array}[]{c}-\left|\vec{n}\right|,\vec{b}+\vec{n}+1\\\
\vec{a}+1\end{array};z\right),\quad\mathcal{P}_{\vec{n}}\in\mathbb{R}.$

###### Proof.

We can use the same idea as in the Jacobi-like setting, see Theorem 2.13.
Since $p<q$, the integrand of the contour integral (11) is still
$\mathcal{O}(t^{-2})$ as $t\to\infty$. ∎

The zeros of these polynomials have not been studied before in full
generality. In the case of the Bessel polynomials, some properties of the
zeros are known, see [12, Section 4.10] for an overview, but such results are
not straightforward to obtain because the zeros accumulate on a curve in the
complex plane. The approach with the finite free multiplicative convolution
doesn’t seem to lead to new results either, because decompositions always
involve polynomials that can have roots in the complex plane (even in the case
$\vec{b}=\vec{a}\sqcup\vec{c}$). For example, we can write

$\displaystyle{}_{q+1}F_{p}\left(\begin{array}[]{c}-\left|\vec{n}\right|,\vec{b}+\vec{n}+1\\\
\vec{a}+1\end{array};z\right)$
$\displaystyle={}_{2}F_{1}\left(\begin{array}[]{c}-\left|\vec{n}\right|,b_{1}+n_{1}+1\\\
a_{1}+1\end{array};z\right)\boxtimes_{\left|\vec{n}\right|}\cdots\boxtimes_{\left|\vec{n}\right|}{}_{2}F_{1}\left(\begin{array}[]{c}-\left|\vec{n}\right|,b_{p}+n_{p}+1\\\
a_{p}+1\end{array};z\right)$
$\displaystyle\quad\boxtimes_{\left|\vec{n}\right|}{}_{2}F_{0}\left(\begin{array}[]{c}-\left|\vec{n}\right|,b_{p+1}+n_{p+1}+1\\\
-\end{array};z\right)\boxtimes_{\left|\vec{n}\right|}\cdots\boxtimes_{\left|\vec{n}\right|}{}_{2}F_{0}\left(\begin{array}[]{c}-\left|\vec{n}\right|,b_{q}+n_{q}+1\\\
-\end{array};z\right),$

which involves the Jacobi polynomials
$P_{\left|\vec{n}\right|}^{(a_{j},b_{j}-a_{j}+n_{j}-\left|\vec{n}\right|)}(z)$
and the Bessel polynomials
$B_{\left|\vec{n}\right|}^{(b_{j}+n_{j}-\left|\vec{n}\right|)}(z)$, both in a
typically negative non-integer parameter.

##  5 Applications

###  5.1 Products of random matrices

We will establish the connection between the multiple orthogonal polynomials
in the Jacobi/Laguerre-like setting and products of truncated unitary random
matrices and Ginibre matrices. In the most basic cases, where we would only
consider a Ginibre matrix or a truncated unitary random matrix, it is well-
known that there is a connection with the Laguerre and Jacobi orthogonal
polynomials, which are the most basic cases of our multiple orthogonal
polynomials (see, e.g., [19]). The connection arises in the following way. Let
$\Delta_{n}(\vec{x})=\prod_{1\leq i<j\leq n}(x_{j}-x_{i})$ denote the
Vandermonde determinant associated with the vector $\vec{x}\in\mathbb{R}^{n}$.
It is known that the joint probability density of the squared singular values
of these random matrices is given by a polynomial ensemble

$\propto\Delta_{n}(\vec{x})\det\left[w_{k}(x_{j})\right]_{j,k=1}^{n},\quad\vec{x}\in\mathbb{R}^{n}$

(the symbol $\propto$ means proportional to). For a size $(n+a)\times n$ upper
left corner of a unitary matrix of size $(n+b)\times(n+b)$ (assuming that
$b\geq n+a+1$), the weight functions are supported on $(0,1)$ and given by
beta densities $w_{k}(x)=\mathcal{B}^{a+k-1,b-n}(x)$. For a size $(n+a)\times
n$ Ginibre matrix, the weight functions are supported on $(0,\infty)$ and
given by gamma densities $w_{k}(x)=\mathcal{G}^{a+k-1}(x)$. Note that the
latter arises as an appropriate limit of the former by doing a suitable
scaling and letting $b\to\infty$.

Typically one is interested in the correlation kernel of such polynomial
ensembles

$K_{n}(x,y)=\sum_{k=0}^{n-1}P_{k}(x)Q_{k}(y).$

The kernel may be built out of monic polynomials $P_{k}$ of degree $k$ and
functions $Q_{l}$ in the linear span of $w_{1},\dots,w_{l+1}$ that satisfy the
biorthogonality conditions

$\int_{0}^{\infty}P_{k}(x)Q_{l}(x)dx=\delta_{k,l},\quad k,l=0,\dots,n-1.$ (22)

These are satisfied whenever

$\int_{0}^{\infty}P_{k}(x)w_{j}(x)dx=0,\quad j=0,\dots,k-1,$ (23)
$\int_{0}^{\infty}Q_{l}(x)x^{j}dx=\delta_{l,j},\quad j=0,\dots,l.$ (24)

Such a construction may not always be possible, however if it is, the $P_{k}$
and $Q_{l}$ are uniquely determined. Constructions using other types of
biorthogonal functions $P_{k}$ and $Q_{l}$ are also allowed as long as they
span the spaces $\text{span}\\{1,\dots,x^{n-1}\\}$ and
$\text{span}\\{w_{1},\dots,w_{n}\\}$ respectively. Another object of interest
is the monic polynomial $P_{n}(x)$ of degree $n$ that extends the
biorthogonality conditions (22) to $k=n$, because it is the average
characteristic polynomial $\mathbb{E}[\prod_{j=1}^{n}(x-x_{j})]$ over the
particles in the associated polynomial ensemble.

For the described truncated unitary matrix, it can be shown that

$P_{k}(x)=P_{k}^{(a,b-n)}(x),\quad
Q_{k-1}(x)=\mathcal{B}^{a,b-n}(x)P_{k}^{(a,b-n)}(x),\quad k=0,\dots,n.$

On the other hand, for the stated Ginibre matrix, we have

$P_{k}(x)=L_{k}^{(a)}(x),\quad
Q_{k-1}(x)=\mathcal{G}^{a}(x)L_{k}^{(a)}(x),\quad k=0,\dots,n.$

We will show that to study the kernel associated with products of these random
matrices, we essentially have to take finite free convolutions of the
underlying polynomials $P_{k}$ and Mellin convolutions of the underlying
functions $Q_{l}$.

The joint probability density of the squared singular values of products of
truncated unitary random matrices and Ginibre matrices can be described using
the following results.

######  Theorem 5.1 (Kieburg-Kuijlaars-Stivigny [20], Corollary 2.4).

Let $T$ be a size $(n+\nu)\times l$ truncation of a Haar distributed unitary
matrix of size $m\times m$ with $\nu\geq 0$, $m\geq l\geq n\geq 1$ and $m\geq
n+\nu+1$. Let $X$ be a random matrix of size $l\times n$, independent of $T$,
such that the joint probability density of its squared singular values is
given by a polynomial ensemble on $(0,\infty)^{n}$ in terms of functions
$f_{1},\dots,f_{n}$. Then the joint probability density of the squared
singular values of $Y=TX$ is given by a polynomial ensemble on
$(0,\infty)^{n}$ in terms of functions $g_{1},\dots,g_{n}$ with

$g_{k}(y)=\int_{0}^{1}x^{\nu}(1-x)^{m-n-\nu-1}f_{k}\left(\frac{y}{x}\right)\frac{dx}{x},\quad
k=1,\dots,n.$

######  Theorem 5.2 (Kuijlaars-Stivigny [21], Theorem 2.1).

Let $G$ be a Ginibre matrix of size $(n+\nu)\times l$ with $\nu\geq 0$ and
$l\geq n\geq 1$. Let $X$ be a random matrix of size $l\times n$, independent
of $G$, such that the joint probability density of its squared singular values
is given by a polynomial ensemble on $(0,\infty)^{n}$ in terms of functions
$f_{1},\dots,f_{n}$. Then the joint probability density of the squared
singular values of $Y=GX$ is given by a polynomial ensemble on
$(0,\infty)^{n}$ in terms of functions $g_{1},\dots,g_{n}$ with

$g_{k}(y)=\int_{0}^{\infty}x^{\nu}e^{-x}f_{k}\left(\frac{y}{x}\right)\frac{dx}{x},\quad
k=1,\dots,n.$

We will first establish the connection between certain products of truncated
unitary matrices and the multiple orthogonal polynomials in the Jacobi-like
setting (studied in Section 2, see Theorem 2.6 and Theorem 2.13). Consider a
multi-index $\vec{n}\in\mathbb{N}^{r}$ with $n=\left|\vec{n}\right|$. Suppose
that $\vec{a},\vec{b}\in\mathbb{N}^{r}$ with
$b_{j}\geq(n-n_{1})\delta_{j,1}+a_{j}$ for all $j$ and set $a_{0}=0$. Let
$T_{j}$ be a size $(n+a_{j})\times(n+a_{j-1})$ truncation of a size
$(n+n_{j}+b_{j})\times(n+n_{j}+b_{j})$ Haar distributed unitary matrix. Then
the joint probability density of the squared singular values of $T_{r}\cdots
T_{1}$ is given by

$\propto\Delta_{n}(x)\det\left[w_{0}(x_{j};\vec{a}+(k-1)\vec{e}_{1},\vec{b}+\vec{n}+(k-n-1)\vec{e}_{1})\right]_{j,k=1}^{n}.$
(25)

with $w_{0}$ as in Section 2.1.

######  Proposition 5.3.

The biorthogonal functions in (23) and (24) associated with the squared
singular values of the product $T_{r}\cdots T_{1}$ described above are given
by

$P_{k}(x)=P_{\vec{k}}(x;\vec{a},\vec{b}+\vec{n}-\vec{k}),\quad
Q_{k-1}(x)=F_{\vec{k}}(x;\vec{a},\vec{b}+\vec{n}-\vec{k}),\quad k=0,\dots,n,$

in terms of type II multiple orthogonal polynomials and type I functions in
the Jacobi-like setting, for every $\vec{k}\in\mathcal{N}^{r}$ with
$|\vec{k}|=k$.

###### Proof.

It was proven in [20, Proposition 2.7] that

$P_{k}(x)=\frac{1}{2\pi
i}\int_{\Sigma_{k}}\frac{1}{(t-k)_{k+1}}\frac{\Gamma(t+\vec{b}+\vec{n}+1)}{\Gamma(t+\vec{a}+1)}x^{t}dt,\quad
k=0,\dots,n-1$

(it also extends to $k=n$), for a positively oriented contour $\Sigma_{k}$
encircling $[0,k]$ once and not enclosing any pole of the integrand not in
this interval, and

$Q_{l}(x)=\frac{1}{2\pi
i}\int_{C}(s-l)_{l}\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+\vec{n})}x^{-s}ds,\quad
k=0,\dots,n-1,$

for a positively oriented contour $C$ that starts and ends at $-\infty$ and
encircles the negative real line. Hence, the result follows after comparison
with the integral representations for the type II multiple orthogonal
polynomials and type I functions in the Jacobi-like setting, see Corollaries
2.14 and 2.8. ∎

######  Corollary 5.4.

If $\vec{n}\in\mathcal{N}^{r}$, then
$P_{\vec{n}}(x;\vec{a},\vec{b})=\mathbb{E}[\prod_{j=1}^{n}(x-x_{j})]$ with the
expectation taken over the density (25).

It is remarkable that a condition of the form $b_{1}\geq(n-n_{1})+a_{j}$ from
Proposition 2.3 is required here in this construction. Since already
$b_{j}-a_{j}\in\mathbb{Z}_{\geq 0}$, the system of weights associated with
$P_{\vec{n}}(x;\vec{a},\vec{b})$ is therefore necessarily an AT-system.

The multiple orthogonal polynomials in the Laguerre-like setting will become
relevant if we also include Ginibre matrices in the product: the results here
should be compatible with the limiting relations between the multiple
orthogonal polynomials and the random matrices that were mentioned before. Let
$p>q$ and consider multi-indices
$\vec{n}\in\mathbb{N}^{q},\vec{m}\in\mathbb{N}^{p-q}$ with
$n=\left|\vec{n}\right|+\left|\vec{m}\right|$. Suppose that
$\vec{a},\vec{b}\in\mathbb{N}^{q},\vec{c}\in\mathbb{N}^{p-q}$ with $a_{j}\leq
b_{j}$ for all $j$. Set $a_{0}=c_{p-q}$ and $c_{0}=0$. Let $T_{j}$, $1\leq
j\leq q$, be a size $(n+a_{j})\times(n+a_{j-1})$ truncation of a size
$(n+n_{j}+b_{j})\times(n+n_{j}+b_{j})$ Haar distributed unitary matrix. Let
$G_{j}$, $1\leq j\leq p-q$, be a size $(n+c_{j})\times(n+c_{j-1})$ Ginibre
matrix. Then the joint probability density of the squared singular values of
$T_{q}\cdots T_{1}\cdot G_{p-q}\cdots G_{1}$ is given by

$\propto\Delta_{n}(x)\det\left[w_{0}(x_{j};\vec{a}\sqcup\vec{c}+(k-1)\vec{e}_{q+1},\vec{b}+\vec{n})\right]_{j,k=1}^{n}.$
(26)

with $w_{0}$ as in Section 3.1.

######  Proposition 5.5.

The biorthogonal functions in (23) and (24) associated with the squared
singular values of the product $T_{q}\cdots T_{1}\cdot G_{p-q}\cdots G_{1}$
described above are given by

$P_{k}(x)=P_{\vec{k}_{1}\sqcup\vec{k}_{2}}(x;\vec{a}\sqcup\vec{c},\vec{b}+\vec{n}-\vec{k}_{1}),\quad
Q_{k-1}(x)=F_{\vec{k}_{1}\sqcup\vec{k}_{2}}(x;\vec{a}\sqcup\vec{c},\vec{b}+\vec{n}-\vec{k}_{1}),\quad
k=0,\dots,n,$

in terms of type II multiple orthogonal polynomials and type I functions in
the Laguerre-like setting, for every
$\vec{k}_{1}\sqcup\vec{k}_{2}\in\mathcal{N}^{p}$, for which
$\vec{k}_{2}\in\mathcal{S}^{p-q}$, with $|\vec{k}_{1}|\,+\,|\vec{k}_{2}|\,=k$.

###### Proof.

The Mellin transform of $Q_{l}$ is of the form

$\hat{Q}_{l}(s)=\sum_{j=1}^{l+1}c_{l}\hat{w}_{0}(s;\vec{a}\sqcup\vec{c}+(j-1)\vec{e}_{q+1},\vec{b}+\vec{n})=\frac{\Gamma(s+\vec{a}\sqcup\vec{c})}{\Gamma(s+\vec{b}+\vec{n})}q_{l}(s)$

in terms of a polynomial $q_{l}(x)$ of degree at most $l$. The conditions in
(24) then imply that $q_{l}(x)\propto(s-l+1)_{l-1}$ and that
$\hat{Q}_{l}(l)=1$. Comparing with the Mellin transform of the type I
functions in the Laguerre-like setting (Theorem 3.6) then leads to the desired
result for the functions $Q_{l}$. Similarly as in the proofs of [20,
Proposition 2.7] (for a product $T_{r}\dots T_{1}$) and [21, Proposition 4.2
and 4.3] (for a product $G_{M-1}\dots G_{1}T_{1}$), we can then verify that

$P_{k}(x)=\frac{1}{2\pi
i}\int_{\Sigma_{k}}\frac{1}{(t-k)_{k+1}}\frac{\Gamma(t+\vec{b}+\vec{n}+1)}{\Gamma(t+\vec{a}\sqcup\vec{c}+1)}x^{t}dt,\quad
k=0,\dots,n.$

Hence, after comparison with the integral representation for the type II
multiple orthogonal polynomials in the Laguerre-like setting (Corollary 3.12),
we obtain the result about the polynomials $P_{k}$ as well. ∎

######  Corollary 5.6.

If $\vec{n}\sqcup\vec{m}\in\mathcal{N}^{p}$ and $\vec{m}\in\mathcal{S}^{p-q}$,
then
$P_{\vec{n}\sqcup\vec{m}}(x;\vec{a}\sqcup\vec{c},\vec{b})=\mathbb{E}[\prod_{j=1}^{n}(x-x_{j})]$
with the expectation taken over the density (26).

Corollaries 5.4 and 5.6 imply that the zeros of the type II multiple
orthogonal polynomials correspond to the zeros of the average characteristic
polynomial of the squared singular values of a product of truncated unitary
and Ginibre matrices. This connection leads to certain properties of the zeros
of the type II polynomials (like simplicity and their location), at least
under appropriate conditions on the parameters. Moreover, one can interpret
the asymptotic zero distribution of the type II multiple orthogonal
polynomials as the asymptotic distribution of the squared singular values of
the corresponding product of truncated unitary and Ginibre matrices. This was
to be expected because the latter arises as the free convolution of the
asymptotic distribution of the squared singular values of its truncated
unitary and Ginibre matrix factors, which corresponds to the deformed arcsin
distributions and Marchenko-Pastur distributions that appear in the former
(see Theorem 2.17 and Theorem 3.14).

The integral representations in the proofs of Propositions 5.3 and 5.5 allow
us to derive a double integral representation for the associated correlation
kernel.

######  Proposition 5.7.

In the setting of Propositions 5.3 and 5.5, we have

$K_{n}(x,y)=\frac{1}{(2\pi
i)^{2}}\int_{C}\int_{\Sigma}\frac{(s-n)_{n}}{(t-n+1)_{n}}\frac{\Gamma(s+\vec{a}\sqcup\vec{c})}{\Gamma(t+\vec{a}\sqcup\vec{c}+1)}\frac{\Gamma(t+\vec{b}+\vec{n}+1)}{\Gamma(s+\vec{b}+\vec{n})}\frac{x^{t-s}}{s-t-1}dtds,$

where $\Sigma$ is a positively oriented contour that starts and ends at
$\infty$ and encircles the positive real line, without enclosing any poles of
the integrand not in $[0,\infty)$, and $C$ is a positively oriented contour
that starts and ends at $-\infty$ and encircles the negative real line.

###### Proof.

The integral representations in the proofs of Propositions 5.3 and 5.5 show
that

$K_{n}(x,y)=\frac{1}{(2\pi
i)^{2}}\int_{C}\int_{\Sigma}\sum_{k=0}^{n-1}\frac{(s-k)_{k}}{(t-k)_{k+1}}\frac{\Gamma(s+\vec{a}\sqcup\vec{c})}{\Gamma(t+\vec{a}\sqcup\vec{c}+1)}\frac{\Gamma(t+\vec{b}+\vec{n}+1)}{\Gamma(s+\vec{b}+\vec{n})}x^{t-s}dtds.$

Then we can use [21, Eq. (6.5)], which is

$\sum_{k=0}^{n-1}\frac{(s-k)_{k}}{(t-k)_{k+1}}=\frac{1}{s-t-1}\left(\frac{(s-n)_{n}}{(t-n+1)_{n}}-1\right).$

The double integral that corresponds to the second term vanishes, because the
contour $\Sigma$ doesn’t contain any poles of the integrand in the variable
$t$. Hence we immediately obtain the desired result. ∎

In particular, this result allows us to obtain hard edge scaling limits at the
origin as $n=\left|\vec{n}\right|+\left|\vec{m}\right|\to\infty$. It would
lead to the different kernels described in [21, Theorem 2.8] that arose as
scaling limits at the origin of kernels associated with certain products of
truncated unitary matrices.

###  5.2 Diophantine approximation

In this section, we will use the connection between multiple orthogonal
polynomials and Hermite-Padé approximation (see, e.g., [12, Chapter 23]) to
study the Diophantine properties of values of certain hypergeometric series.
Generally, Hermite-Padé approximation is used to simultaneously approximate
several functions $f_{j}(z)$ with an expansion of the form

$f_{j}(z)=\sum_{k\geq 0}\frac{m_{j,k}}{z^{k+1}},\ z\to\infty.$

Consider $p\leq q$ and suppose that $\vec{a}\in(-1,\infty)^{p}$ and
$\vec{b}\in(-1,\infty)^{q}$ satisfy $a_{j}<b_{j}$ for all $1\leq j\leq p$. The
multiple orthogonal polynomials in the Jacobi- and Bessel-like setting will
appear in the Hermite-Padé approximation problems associated with the $q$
functions $f_{j}(z;\vec{a},\vec{b})=f_{0}(1/z;\vec{a},\vec{b}+\vec{e}_{j})/z$,
with $f_{0}$ as in (1). The series $f_{j}(z;\vec{a},\vec{b})$ are convergent
in $\mathbb{C}_{0}$, if $p<q$, and in $\\{z\in\mathbb{C}:0<\left|z\right|\leq
1\\}$, if $p=q$. The associated type I Hermite-Padé approximation problem is
given by

$\sum_{j=1}^{q}A_{\vec{n},j}(z)f_{j}(z;\vec{a},\vec{b})-B_{\vec{n}}(z)=\mathcal{O}(z^{-\left|\vec{n}\right|}),\
z\to\infty,$

in terms of polynomials $A_{\vec{n},j}(z)$ and $B_{\vec{n}}(z)$ of degree at
most $n_{j}-1$ and $\max\\{n_{j}\mid j\\}-2$ respectively. The corresponding
type II Hermite-Padé approximation problem is

$P_{\vec{n}}(z)f_{j}(z;\vec{a},\vec{b})-Q_{\vec{n},j}(z)=\mathcal{O}(z^{-n_{j}-1}),\
z\to\infty,\quad j=1,\dots,q,$

in terms of polynomials $P_{\vec{n}}(z)$ and $Q_{\vec{n},j}(z)$ of degree at
most $\left|\vec{n}\right|$ and $\left|\vec{n}\right|-1$ respectively. The
vector of polynomials $(A_{\vec{n},1},\dots,A_{\vec{n},q})$ and the polynomial
$P_{\vec{n}}$ are precisely the type I and type II multiple orthogonal
polynomials associated with the moment generating functions
$(f_{1}(z;\vec{a},\vec{b}),\dots,f_{q}(z;\vec{a},\vec{b}))$, i.e. with weights
$(w_{1}(z;\vec{a},\vec{b}),\dots,w_{q}(z;\vec{a},\vec{b}))$ as in Sections 2.1
and 4.1. Explicit expressions for the other polynomials and the errors are
then also known.

We will focus on the type I Hermite-Padé approximation problem. In order to
generate interesting examples, it will be more convenient to consider the $q$
functions
$f_{j}^{\ast}(z;\vec{a},\vec{b})=f_{0}(1/z;\vec{a},\vec{b}+\sum_{i=1}^{j}\vec{e}_{i})/z$.
We can move between the two settings using the lemma below. The second
relation in this lemma can be used to obtain a new approximation problem

$\sum_{j=1}^{q}A_{\vec{n},j}^{\ast}(z)f_{j}^{\ast}(z;\vec{a},\vec{b})-B_{\vec{n}}(z)=\mathcal{O}(z^{-\left|\vec{n}\right|}),\
z\to\infty,$ (27)

in terms of polynomials
$A_{\vec{n},k}^{\ast}(z)=\sum_{j=k}^{q}c_{j,k}A_{\vec{n},j}(z)$ of degree at
most $\max\\{n_{k},\dots,n_{q}\\}-1$.

######  Lemma 5.8.

Suppose that the entries of $\vec{b}$ are pairwise distinct. Then,

  * i)

$f_{j}^{\ast}(z;\vec{a},\vec{b})=\sum_{k=1}^{j}\lambda_{j,k}f_{k}(z;\vec{a},\vec{b})$
with $\lambda_{j,k}=1/\prod_{i=1,i\neq k}^{j}(b_{i}-b_{k})$,

  * ii)

$f_{j}(z;\vec{a},\vec{b})=\sum_{k=1}^{j}c_{j,k}f_{k}^{\ast}(z;\vec{a},\vec{b})$
with $c_{j,k}=\prod_{i=1}^{k-1}(b_{i}-b_{j})$.

###### Proof.

The first part follows from a partial fraction decomposition of
$1/\prod_{i=1}^{j}(s+b_{i})$. For the second part, we have to verify that

$\frac{1}{s+b_{j}}=\sum_{k=1}^{j}\frac{\prod_{i=1}^{k-1}(b_{i}-b_{j})}{\prod_{i=1}^{k}(s+b_{i})}.$

We may use the decomposition $(s+b_{j})/(s+b_{k})=1+(b_{j}-b_{k})/(s+b_{k})$
to write

$\sum_{k=1}^{j}\frac{s+b_{j}}{s+b_{k}}\prod_{i=1}^{k-1}\frac{b_{i}-b_{j}}{s+b_{i}}=\sum_{k=1}^{j}\prod_{i=1}^{k-1}\frac{b_{i}-b_{j}}{s+b_{i}}-\sum_{k=1}^{j}\prod_{i=1}^{k}\frac{b_{i}-b_{j}}{s+b_{i}}=1,$

which then immediately leads to the desired result. ∎

It is also possible to approximate other linear combinations of the
$f_{j}(z;\vec{a},\vec{b})$ using combinations of the type I polynomials
$A_{\vec{n},j}(z)$, as long as the resulting $q$ functions are
$\mathbb{C}$-linearly independent.

We can use these approximants to prove the following result.

######  Theorem 5.9.

Consider $p<q$. Let
$\vec{a}\in(\mathbb{Q}_{>-1})^{p},\vec{b}\in(\mathbb{Q}_{>-1})^{q}$ with all
$a_{j}<b_{j}$ and let $z\in\mathbb{Q}_{0}$. Suppose that
$a_{i}-b_{j},b_{i}-b_{j}\not\in\mathbb{Z}_{\geq 0}$ whenever $i\neq j$. Then
$1$ and the constants

$\sum_{k\geq
0}\frac{(\vec{a}+1)_{k}}{(\vec{b}+1)_{k}\prod_{i=1}^{j}(b_{i}+k+1)}z^{k},\quad
j=1,\dots,q,$

are $\mathbb{Q}$-linearly independent. In particular, these constants are
irrational.

In many cases in which some conditions on the parameters are not satisfied, we
can still find a similar result as the above. However, we would have to take
appropriate combinations of the approximated functions such that their
confluent limits are still $\mathbb{C}$-linearly independent (so that, e.g.,
the same function doesn’t appear twice). The Mellin transform of the
corresponding type I function would remain invariant under these operations,
but a priori it is not clear how the underlying type I polynomials would be
affected. Hence, for such cases, a more hands-on approach is required. For
example, in the simplest case, where we take
$\vec{a}=\vec{0}\in\mathbb{R}^{p}$ and $\vec{b}=\vec{0}\in\mathbb{R}^{q}$, we
can show the following result.

######  Theorem 5.10.

Consider $p<q$ and let $z\in\mathbb{Q}_{0}$. Then $1$ and the constants

$\sum_{k\geq 0}\frac{z^{k}}{(k+1)^{j}\cdot k!^{q-p}},\quad j=1,\dots,q,$

are $\mathbb{Q}$-linearly independent. In particular, these constants are
irrational.

Some special cases of Theorems 5.9 and 5.10 have been proven explicitly
before. The case $\vec{a}=(),\vec{b}=(0)$ that corresponds to an approximation
using Bessel polynomials, leads to an irrationality proof of $\exp(z)$ for
$z\in\mathbb{Q}_{0}$ and was studied by Grosswald in [11]. In [35], Rivoal
extended this result to $\vec{a}=(\alpha),\vec{b}=(0,\alpha)$ with
$\alpha\in\mathbb{Q}$ and proved $\mathbb{Q}$-linear independence of $1$,
$\exp(z)$ and $\mathcal{E}_{\alpha}(z)=\sum_{m\geq 0}z^{m}/((m+\alpha+1)\cdot
m!)$. The latter implied that at least one of Euler’s constant $\gamma$ and
Euler-Gompertz’ constant $\delta$ is irrational.

In fact, all the cases covered by Theorems 5.9 and 5.10, follow from Siegel’s
theory on $E$-functions (see, e.g., [36] for an introduction), which was
subsequently improved in [37, 7, 45]; the matrix Pearson equation for the
underlying weights (see Proposition 4.1) can be lifted to a system of
differential equations for the approximated functions. However, the
approximants that are used to prove this are only known implicitly. In our
approach, they are explicit and correspond to Hermite-Padé approximants for
certain hypergeometric moment generating functions.

Note that if we could take $p=q$ in Theorem 5.10, the constants would reduce
to values of consecutive polylogarithms $\text{Li}_{j}(z)$ for $j=1,\dots,q$.
In particular, for $z=-1$, we would recover a logarithm $\ln 2$ and the zeta
values $\zeta(j)$ for $j=2,\dots,q$. In this case, our proof will not be valid
anymore and a more careful analysis of the approximants is required; they will
grow geometrically and will not have a behavior determined by $n!$ anymore.
Typically, this will lead to a condition on the values that $z\in\mathbb{Q}$
can take, see [30] where this case was studied.

The proofs of Theorems 5.9 and 5.10 will consist of two main parts. First, we
will analyze the Diophantine quality of the approximants in (27), i.e. the
growth and Diophantine properties of the denominator polynomials
$A_{\vec{n},j}^{\ast}(z)$ (that arose as linear combinations of the type I
polynomials $A_{\vec{n},j}(z)$ from Theorem 4.4) and the numerator polynomials
$B_{\vec{n}}(z)$. Then, we will determine the approximation quality, i.e. the
growth of the error. Afterwards, since the overall quality of the approximants
will be good enough, we will be able to apply Nesterenko’s criterion for
$\mathbb{Q}$-linear independence of real numbers, stated below, in its
strongest form (with $\tau_{1}=\tau_{2}=N-1$).

######  Theorem 5.11 (Nesterenko [33]).

Let $x_{1},\dots,x_{N}\in\mathbb{R}$ and suppose that there exists sequences
of integers $(p_{j,n})_{n\in\mathbb{N}}$, $j=1,\dots,N$, and an increasing
function $\sigma:\mathbb{R}\to(0,\infty)$, with
$\lim_{t\to\infty}\sigma(t)=\infty$ and
$\limsup_{t\to\infty}\sigma(t+1)/\sigma(t)=1$, such that

  * i)

$\max\\{p_{j,n}\mid j=1,\dots,N\\}\leq e^{\sigma(n)}$,

  * ii)

$c_{1}e^{-\tau_{1}\sigma(n)}\leq\left|\sum_{j=1}^{N}p_{j,n}x_{j}\right|\leq
c_{2}e^{-\tau_{2}\sigma(n)}$ for some $c_{1},c_{2},\tau_{1},\tau_{2}>0$.

Then
$\dim_{\mathbb{Q}}\text{span}_{\mathbb{Q}}\\{x_{1},\dots,x_{N}\\}\geq(\tau_{1}+1)/(1+\tau_{1}-\tau_{2}).$

We will start by analyzing the Diophantine quality of the approximants in the
setting of Theorem 5.9.

######  Lemma 5.12.

(Diophantine quality) Consider $p<q$. Let
$\vec{a}\in(\mathbb{Q}_{>-1})^{p},\vec{b}\in(\mathbb{Q}_{>-1})^{q}$ with all
$a_{j}<b_{j}$ and let $z\in\mathbb{Q}_{0}$. Suppose that
$a_{i}-b_{j},b_{i}-b_{j}\not\in\mathbb{Z}_{\geq 0}$ whenever $i\neq j$. Take
$\vec{n}=(n,\dots,n)\in\mathbb{N}^{q}$. Then $A_{\vec{n},j}(z)\leq
C^{n}n!^{q-p}$ as $n\to\infty$ for some $C=C(z;\vec{a},\vec{b})>0$. There
exists a sequence of integers $(D_{n})_{n\in\mathbb{N}}$ and
$D=D(z;\vec{a},\vec{b})>0$ such that $D_{n}A_{\vec{n},j}(z)\in\mathbb{Z}$ and
$D_{n}\sim D^{n}$ as $n\to\infty$. As a consequence, similar properties hold
for $A^{\ast}_{\vec{n},j}(z)$. It also holds for
$B_{\vec{n}}(z)\Gamma(\vec{b}+1)/\Gamma(\vec{a}+1)$.

###### Proof.

We may use a similar idea as in the proof of [36, Lemma 2] in which one shows
that certain ${}_{p}F_{q}$-series with $p<q$ are $E$-functions. It was shown
that $(a)_{k}/(b)_{k}\sim k^{c}$ as $k\to\infty$ and the least common
denominator of $(a)_{k}/(b)_{k}$ for $k=0,\dots,n$ is of the order $O(d^{n})$
as $n\to\infty$. A suitable representation for $A_{\vec{n},j}(z)$ to analyze
can be found in Theorem 4.4. Since there is always a deficit of $q-p$
Pochhammer symbols in the denominators of the terms of $A_{\vec{n},j}(z)$, the
desired properties immediately follow. The same then holds for
$A_{\vec{n},j}^{\ast}(z)$ due to Lemma 5.8 ii). An appropriate representation
for $B_{\vec{n}}(z)$ follows from the proof of Theorem 4.4. After using
formula (2.2) to expand $I_{J,L}(z)$ in the error (8), we find

$B_{\vec{n}}(z)=\frac{\Gamma(\vec{a}+1)}{\Gamma(\vec{b}+1)}\sum_{J=1}^{r}\sum_{K=0}^{n_{J}-1}P_{\vec{n},J}[K]\sum_{k=0}^{K-1}\frac{(\vec{b}-b_{J}-K)_{k+1}}{(\vec{a}-b_{J}-K)_{k+1}}\frac{1}{K-k+b_{J}}z^{k+1}$

and we can repeat the same argument. ∎

A similar result holds in the simpler setting of Theorem 5.10.

######  Lemma 5.13.

(Diophantine quality) Consider $p<q$ and set
$\vec{a}=\vec{0}\in\mathbb{R}^{p},\vec{b}=\vec{0}\in\mathbb{R}^{q}$. Suppose
that $z\in\mathbb{Q}_{0}$ and take $\vec{n}=(n,\dots,n)$. Then
$A^{\ast}_{\vec{n},j}(z)\leq C^{n}n!^{q-p}$ as $n\to\infty$ for some
$C=C(z;\vec{a},\vec{b})>0$. There exists a sequence of integers
$(D_{n})_{n\in\mathbb{N}}$ and $D=D(z;\vec{a},\vec{b})>0$ such that
$D_{n}A^{\ast}_{\vec{n},j}(z)\in\mathbb{Z}$ and $D_{n}\sim D^{n}$ as
$n\to\infty$. It also holds for $B_{\vec{n}}(z)$.

###### Proof.

We will use the same idea as in the proof of Theorem 4.4 to obtain expressions
for the polynomials $A^{\ast}_{\vec{n},j}(z)$ in (27). The error in the
approximation problem is

$\sum_{j=1}^{q}A^{\ast}_{\vec{n},j}(z)f_{j}(z)-B_{\vec{n}}(z)=\sum_{k\geq
0}\hat{F}_{\vec{n}}(k+1)\frac{1}{z^{k+1}},\quad f_{j}(z)=\sum_{k\geq
0}\frac{1}{(k+1)^{j}k!^{q-p}}\frac{1}{z^{k+1}}.$

We will analyze the polynomials $A^{\ast}_{\vec{n},j}(z)$ by working out the
right-hand side of the above. We know from Theorem 4.3 that

$\hat{F}(s)=\frac{(1-s)_{qn-1}}{\Gamma(s)^{q-p}(s)_{n}^{q}},\quad
s\in\mathbb{Z}_{\geq 1}.$

The ratio $(1-s)_{qn-1}/(s)_{n}^{q}$ can be seen as the Mellin transform of
the (limiting) type I Jacobi-Piñeiro polynomial
$G_{\vec{n}}(x)=\sum_{j=1}^{q}P^{\ast}_{\vec{n},j}(x)(\ln x)^{j-1}$ with $\deg
P^{\ast}_{\vec{n},j}\leq n-1$, see [38] or [30]. It was proven in [30] that
there exists $c=c(z)>0$ such that $P^{\ast}_{\vec{n},j}[k]\leq c^{n}$, for
$0\leq k\leq n_{j}-1$, as $n\to\infty$ and that there exists a sequence of
integers $(d_{n})_{n\in\mathbb{N}}$ and $d=d(z)>0$ such that
$d_{n}P^{\ast}_{\vec{n},j}[k]\in\mathbb{Z}$, for $0\leq k\leq n_{j}-1$ and
$1\leq j\leq q$, and $d_{n}\sim d^{n}$ as $n\to\infty$. The associated partial
fraction decomposition gives

$\frac{(1-s)_{qn-1}}{(s)_{n}^{q}}=\sum_{j=1}^{q}\sum_{l=0}^{n-1}\frac{P^{\ast}_{\vec{n},j}[l]}{(s+l)^{j}}.$

Hence,

$\sum_{k\geq
0}\hat{F}_{\vec{n}}(k+1)\frac{1}{z^{k+1}}=\sum_{j=1}^{q}\sum_{l=0}^{n-1}P^{\ast}_{\vec{n},j}[l]\sum_{k\geq
0}\frac{1}{k!^{q-p}}\frac{1}{(k+l+1)^{j}}\frac{1}{z^{k+1}},$

and it remains to find a way to lower $L$ in

$I_{J,L}^{[0]}(z)=\sum_{k\geq
0}\frac{1}{k!^{q-p}}\frac{1}{(k+L+1)^{J}}\frac{1}{z^{k+1}},$

because $I_{J,0}^{[0]}(z)=f_{J}(z)$. We may do so by making use of the two
operations

$I_{J,L}^{[t]}(z)=I_{J-1,L}^{[t+1]}(z)-L\cdot I_{J,L}^{[t+1]}(z),\quad
I_{J,L}^{[q-p]}(z)=zI_{J,L-1}^{[0]}(z)-\frac{1}{L^{J}},$

where

$I_{J,L}^{[t]}(z)=\sum_{k\geq
0}\frac{1}{k!^{q-p}}\frac{1}{(k+1)^{t}}\frac{1}{(k+L+1)^{J}}\frac{1}{z^{k+1}},\quad
I_{0,L}^{[t]}(z)=f_{t}(z).$

For a given $L$, the first operation needs to be applied at most $q-p$ times
to reach $t=q-p$ which in total leads to an additional factor that is at most
$L^{q-p}$. We can then decrease $L$ by $1$ using the second operation, doing
so adds a factor $z$. Repeating this $L$ times leads to an additional factor
that is at most $z^{L}L!^{q-p}$. Since we have to do this for every $1\leq
L\leq n-1$, the first part of the lemma follows. For the second part, we
observe that for a given $L$, the first operation doesn’t introduce any
denominators, while the second one introduces a denominator $L^{J}$ to the
coefficients of the $A_{\vec{n},j}^{\ast}(z)$. Hence, after reducing $L$ to
$0$, we get a denominator that is at most $\text{lcm}(1,\dots,L)^{J}$. To
handle all $L$, we thus require a factor of at most
$\text{lcm}(1,\dots,n-1)^{q-p}$ to cancel all denominators in the coefficients
of $A_{\vec{n},j}^{\ast}(z)$. As $z\in\mathbb{Q}$, it then remains to multiply
with $\text{denom}(z)^{n-1}$. ∎

Since the error in the approximation problem only depends on the Mellin
transform of the type I functions (and not on its precise decomposition), the
approximation quality can be derived for both cases simultaneously.

######  Lemma 5.14.

(Approximation quality) Consider $p<q$. Let
$\vec{a}\in(-1,\infty)^{p},\vec{b}\in(-1,\infty)^{q}$ with all $a_{j}<b_{j}$.
Take $\vec{n}=(n,\dots,n)$. Then there exists $E=E(z;\vec{a},\vec{b})>0$ such
that $\sum_{k\geq 0}\hat{F}_{\vec{n}}(k+1)/z^{k+1}=E^{n(1+o(1))}n!^{-q(q-p)},\
n\to\infty.$

###### Proof.

It follows from Theorem 4.3 and the orthogonality conditions that

$\sum_{k\geq 0}\hat{F}_{\vec{n}}(k+1)\frac{1}{z^{k+1}}=z^{-qn+1}\sum_{k\geq
0}\hat{F}_{\vec{n}}(qn+k)z^{-k},\quad\hat{F}_{\vec{n}}(s)=\frac{\Gamma(s+\vec{a})}{\Gamma(s+\vec{b}+n)}(1-s)_{qn-1}.$

Hence,

$\sum_{k\geq
0}\hat{F}_{\vec{n}}(k+1)\frac{1}{z^{k+1}}=(-z)^{-qn+1}\frac{(1)_{qn-1}(\vec{a}+1)_{qn-1}}{(\vec{b}+1)_{(q+1)n-1}}\frac{\Gamma(\vec{a}+1)}{\Gamma(\vec{b}+1)}{}_{p+1}F_{q}\left(\begin{array}[]{c}qn,\vec{a}+qn\\\
\vec{b}+(q+1)n\end{array};\frac{1}{z}\right).$

The factor in front of the hypergeometric series is asymptotic to
$E^{n(1+o(1))}n!^{-q(q-p)}$ as $n\to\infty$ and

$\lim_{n\to\infty}{}_{p+1}F_{q}\left(\begin{array}[]{c}qn,\vec{a}+qn\\\
\vec{b}+(q+1)n\end{array};\frac{1}{z}\right)=\begin{dcases}\sum_{k\geq
0}\frac{q^{(p+1)k}}{(q+1)^{qk}}\frac{z^{-k}}{k!}=\exp(q^{q}/(q+1)^{q}z^{-1}),\text{
if }p=q-1,\\\ 1,\text{ if }p<q-1,\end{dcases}$

hence the stated result follows. ∎

We will end this section by noting that the structure of the type II
approximation problem, in which individual approximants for the constants
appear, suggests that it may be used to determine bounds for the irrationality
measure of the approximated constants. At least from trivial estimates using
determinant relations between the type I and type II systems, we would only
obtain the trivial lower bound $\mu(x)\geq 2$ for these irrationality measures
(equivalent to the constants being irrational).

## Acknowledgements

I would like to thank Walter Van Assche for the many discussions on the topics
presented here and for providing feedback on the drafts, Hélder Lima for
discussing the results about the zeros of the multiple orthogonal polynomials,
Andrei Martínez-Finkelshtein for explaining the connection between
hypergeometric polynomials and the finite free multiplicative convolution of
polynomials, Arno Kuijlaars for providing feedback on the section about random
matrices and Wadim Zudilin for suggesting to investigate explicit Diophantine
approximation of the constant $\sum_{k\geq 1}1/(k\cdot k!)$.

## References

  * [1] C. Andréief, Note sur une relation les intégrates définies des produits des fonctions, M´em. de la Soc. Sci. Bordeaux 2 (1883), 1–11.
  * [2] A. Angelesco, Sur deux extensions des fractions continues algébriques, Comptes Rendus Acad. Sci. Paris 168 (1919), 262–265.
  * [3] R. Apéry, Irrationalité de $\zeta(2)$ et $\zeta(3)$, Astérisque 61 (1979), 11–13.
  * [4] O. Arizmendi, Octavio, J. Garza-Vargas, D. Perales, Finite free cumulants: multiplicative convolutions, genus expansion and infinitesimal distributions, Trans. Amer. Math. Soc. 376 (2023), no. 6, 4383–4420.
  * [5] A. Belton, D. Guillot, A. Khare, M. Putinar, Preservers of totally positive kernels and Pólya frequency functions, Mathematics Research Reports 3 (2022), 35–56.
  * [6] F. Beukers, Padé‐approximations in number theory, in Padé approximation and its applications (Amsterdam, 1980), Lecture Notes in Math. 888, Springer, Berlin‐New York, 1981, pp. 90–99.
  * [7] G. V. Chudnovsky, On some applications of diophantine approximations, Proc. Nat. Acad. Sci. USA 81 (1984), 1926–1930.
  * [8] NIST Digital Library of Mathematical Functions. https://dlmf.nist.gov/, Release 1.1.12 of 2023-12-15. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
  * [9] U. Fidalgo Prieto, G. López Lagomasino, Nikishin systems are perfect, Constr. Approx. 34 (2011), no. 3, 297–356.
  * [10] U. Fidalgo Prieto, G. López Lagomasino, Nikishin systems are perfect. The case of unbounded and touching supports, J. Approx. Theory 163 (2011), no. 6, 779–811.
  * [11] E. Grosswald, The Bessel Polynomials, Lecture Notes in Mathematics, vol. 698, Berlin: Springer, 1978.
  * [12] M.E.H. Ismail, Classical and Quantum Orthogonal Polynomials in One Variable, Encyclopedia of Mathematics and its Applications 98, Cambridge University Press, 2008.
  * [13] S. Karlin, Total positivity, absorption probabilities and applications, Trans. Amer. Math. Soc. 111 (1964), 33–107.
  * [14] S. Karlin, W.J. Studden, Tchebycheff Systems: with Applications In Analysis and Statistics, New York (N.Y.): Interscience, 1966.
  * [15] O. Katkova, B. Shapiro, A. Vishnyakova, Multiplier sequences and logarithmic mesh, C. R. Math. Acad. Sci. Paris 349 (2011), no. 1-2, 35–38.
  * [16] A. Khare, Multiply positive functions, critical exponent phenomena, and the Jain–Karlin–Schoenberg kernel. Preprint available at arXiv:2008.05121 (2021).
  * [17] A.B.J. Kuijlaars, Multiple orthogonal polynomials in random matrix theory, Proceedings of the International Congress of Mathematicians. Volume III, 1417–1432, Hindustan Book Agency, New Delhi, 2010.
  * [18] A.B.J. Kuijlaars, Multiple orthogonal polynomial ensembles, in ’Recent trends in orthogonal polynomials and approximation theory’, Contemp. Math. 507, Amer. Math. Soc., Providence, RI, 155–176, 2010.
  * [19] A.B.J. Kuijlaars, Transformations of polynomial ensembles, in ’Modern trends in constructive function theory’, Contemp. Math. 661, Amer. Math. Soc., Providence, RI, 253–268, 2016.
  * [20] M. Kieburg, A.B.J. Kuijlaars, D. Stivigny, Singular value statistics of matrix products with truncated unitary matrices, Int. Math. Res. Not. 11 (2016), 3392–3424.
  * [21] A.B.J. Kuijlaars, D. Stivigny, Singular values of products of random matrices and polynomial ensembles, Random Matrices Theory Appl. 3 (2014), no. 3., 1450011.
  * [22] A.B.J. Kuijlaars, L. Zhang, Singular values of products of Ginibre random matrices, multiple orthogonal polynomials and hard edge scaling limits, Commun. Math. Phys. 332 (2014), 759–781.
  * [23] M. Leurs, W. Van Assche, Laguerre-Angelesco multiple orthogonal polynomials on an $r$-star, J. Approx. Theory 250 (2020), 105324.
  * [24] H. Lima, Multiple orthogonal polynomials associated with branched continued fractions for ratios of hypergeometric series, Adv. Appl. Math. 147 (2023), 102505.
  * [25] H. Lima, A. Loureiro, Multiple orthogonal polynomials associated with confluent hypergeometric functions, J. Approx. Theory 260 (2020), 105484.
  * [26] H. Lima, A. Loureiro, Multiple orthogonal polynomials with respect to Gauss’ hypergeometric function, Stud. in Appl. Math. 148 (2021), no. 1., 154–185.
  * [27] A.W. Marcus, D.A. Spielman, N. Srivastava, Finite free convolutions of polynomials. Probab. Theory Relat. Fields 182 (2022), 807–848.
  * [28] K. Mahler, Perfect systems, Compos. Math. 19 (1968), 95–166.
  * [29] A. Martínez-Finkelshtein, R. Morales, D. Perales, Real roots of hypergeometric polynomials via finite free convolution. Preprint available at arXiv:2309.10970 (2023).
  * [30] E.M. Nikishin, On irrationality of the values of the functions $F(x,s)$, Math. USSR Sb. 37 (1980), 381.
  * [31] E.M. Nikishin, V.N. Sorokin, Rational Approximations and Orthogonality, Translations of Mathematical Monographs 92, Amer. Math. Soc., Providence RI, 1991.
  * [32] Yu. Nesterenko, Hermite-Padé approximants of generalized hypergeometric functions, 1995 Sb. Math. 83 (1995), 189.
  * [33] Yu. Nesterenko, Linear independence of numbers, Vestnik Moskov. Univ. Ser. I Mat. Mekh. 108 (1985), 46–49 (Russian); English translation in Moscow Univ. Math. Bull. 40 (1985), 69–74.
  * [34] T. Neuschel, W. Van Assche, Asymptotic zero distribution of Jacobi-Piñeiro and multiple Laguerre polynomials, J. Approx. Theory 205 (2016), 114–132.
  * [35] T. Rivoal, On the Arithmetic Nature of the Values of the Gamma Function, Euler’s Constant, and Gompertz’s Constant, Michigan Math. J. 61 (2012), 239–254.
  * [36] A.B. Shidlovskii, Transcendental Numbers, De Gruyter Studies in Mathematics, vol. 12, Berlin, New York: De Gruyter, 1989.
  * [37] A.B. Shidlovskii, On estimates of the measure of transcendence of the values of E-functions, Mat. Zametki 2 (1967), 33–44; English translation in Math. Notes 2 (1967).
  * [38] C. Smet, W. Van Assche, Mellin transforms for multiple Jacobi-Piñeiro polynomials and a $q$-analogue, J. Approx. Theory 162 (2010), no. 4, 782–806.
  * [39] A.D. Sokal, Multiple orthogonal polynomials, $d$-orthogonal polynomials, production matrices, and branched continued fractions. Preprint available at arXiv:2204.11528 (2022), to appear in Trans. Amer. Math. Soc.
  * [40] G. Szegő, Bemerkungen zu einem Satz von J. H. Grace über die Wurzeln algebraischer Gleichungen, Math. Z. 13 (1922), 28–55.
  * [41] G. Szegő, Orthogonal Polynomials, Amer. Math. Soc. Colloq. Publ. 23, Amer. Math. Soc., Providence RI, 1939 (fourth edition 1975).
  * [42] W. Van Assche, E. Coussement, Some classical multiple orthogonal polynomials, Numerical Analysis 2000, vol. V, Quadrature and Orthogonal Polynomials, J. Comput. Appl. Math. 127 (2001), 317–347.
  * [43] W. Van Assche, T. Wolfs, Multiple orthogonal polynomials associated with the exponential integral, Stud. Appl. Math. 151 (2023), 411–449.
  * [44] W. Van Assche, S. Yakubovich Multiple orthogonal polynomials associated with Macdonald functions, Integral Transf. Spec. Funct. 9 (2000), no. 3, 229–244.
  * [45] W. Zudilin, On rational approximations of values of a certain class of entire functions, Sb. Math. 186 (1995), 555.