# Multiple orthogonal polynomials associated with the exponential integral

Walter Van Assche111Department of Mathematics, KU Leuven, Celestijnenlaan 200B
box 2400, BE-3001 Leuven, Belgium. Supported by FWO project G.0C9819N.
<EMAIL_ADDRESS>Thomas
Wolfs11footnotemark: 1

(August 27, 2024)

###### Abstract

We introduce a new family of multiple orthogonal polynomials satisfying
orthogonality conditions with respect to two weights $(w_{1},w_{2})$ on the
positive real line, with $w_{1}(x)=x^{\alpha}e^{-x}$ the gamma density and
$w_{2}(x)=x^{\alpha}E_{\nu+1}(x)$ a density related to the exponential
integral $E_{\nu+1}$. We give explicit formulas for the type I functions and
type II polynomials, their Mellin transform, Rodrigues formulas,
hypergeometric series and recurrence relations. We determine the asymptotic
distribution of the (scaled) zeros of the type II multiple orthogonal
polynomials and make a connection to random matrix theory. Finally, we also
consider a related family of mixed type multiple orthogonal polynomials.

Keywords: multiple orthogonal polynomials, exponential integral, Mellin
transform, random matrices

##  1 Introduction

Multiple orthogonal polynomials extend the notion of (regular) orthogonal
polynomials by distributing orthogonality conditions over several measures
instead of just one measure. Each multiple orthogonal polynomial corresponds
to a multi-index $\vec{n}=(n_{1},\dots,n_{r})$ in $\mathbb{Z}_{\geq 0}^{r}$
that determines the way in which the associated orthogonality conditions are
distributed. There are two main types of multiple orthogonal polynomials, see
for example [17, Chapter 23] and [27, Chapter 4] for an introduction.

Let $\mu_{1},\dots,\mu_{r}$ be $r$ positive Borel measures, with finite
moments, supported on the real line. We start by introducing the type I
multiple orthogonal polynomials. These are vectors of polynomials
$(A_{\vec{n},1},\dots,A_{\vec{n},r})$, where the degree of $A_{\vec{n},j}$ is
at most $n_{j}-1$, that satisfy the orthogonality conditions

$\sum_{j=1}^{r}\int_{\mathbb{R}}A_{\vec{n},j}(x)x^{k}d\mu_{j}(x)=0,\quad
k=0,\dots,\left|\vec{n}\right|-2.$

The orthogonality conditions result in a system of $\left|\vec{n}\right|-1$
homogeneous linear equations for the $\left|\vec{n}\right|$ coefficients of
the $A_{\vec{n},j}$. The existence of non-zero vectors of polynomials is
therefore guaranteed. Another way to set up the orthogonality conditions is
via the type II multiple orthogonal polynomials. These are polynomials
$P_{\vec{n}}$ of degree at most $\left|\vec{n}\right|=n_{1}+\dots n_{r}$ that
satisfy the orthogonality conditions

$\int_{\mathbb{R}}P_{\vec{n}}(x)x^{k}d\mu_{j}(x)=0,\quad k=0,\dots,n_{j}-1,\
j=1,\dots,r.$

The orthogonality conditions result in a system of $\left|\vec{n}\right|$
homogeneous linear equations for the $\left|\vec{n}\right|+1$ coefficients of
$P_{\vec{n}}$. The existence of non-zero polynomials is therefore again
guaranteed. It is also possible to combine type I and type II orthogonality
conditions in order to obtain so-called mixed type multiple orthogonal
polynomials. They appear in the work of Sorokin [31] and in the study of non-
intersecting Brownian motions with several starting points and endpoints [7].
Later, they were also studied in [1, 11, 12]. Other applications of multiple
orthogonal polynomials can be found in analytic number theory (Diophantine
approximation), for example in Apéry’s proof [3] of the irrationality of
$\zeta(3)$ as explained by Beukers in [5], or in random matrix theory [18,
19].

In [26], Mahler introduced the notion of a perfect system. A system of
measures is said to be perfect if, for all multi-indices, every type I and
type II polynomial attains its maximal degree. In that case, the multiple
orthogonal polynomials are uniquely determined by the orthogonality conditions
up to a non-zero scalar multiplication. The system of equations then has a
unique solution if we add a normalization condition. For the type I
polynomials, we can impose that additionally

$\sum_{j=1}^{r}\int_{\mathbb{R}}A_{\vec{n},j}(x)x^{\left|\vec{n}\right|-1}d\mu_{j}(x)=1.$

The (unique) type I polynomial $(A_{\vec{n},1},\dots,A_{\vec{n},r})$ for which
this holds will be denoted by
$(A_{\vec{n},1}^{(I)},\dots,A_{\vec{n},r}^{(I)})$. We can normalize the type
II polynomials by demanding that they are monic. The (unique) monic type II
polynomial $P_{\vec{n}}$ will be denoted by $P_{\vec{n}}^{(II)}$.

Most of the well-known examples of perfect systems arise via the so-called
Angelesco systems and AT-systems. Angelesco systems were introduced in [2] and
can be produced by considering measures supported on intervals that may touch
at the endpoints but are otherwise disjoint. For an AT-system, one considers
multiple measures supported on the same interval that are sufficiently
different. Particular interest goes to the so-called Nikishin systems, see
[27, Chapter 5, §7], where this is called an MT-system, and [9, 10] where it
is proven that Nikishin systems are perfect.

We will investigate a family of multiple orthogonal polynomials that will be
an extension of the Laguerre polynomials, a family of classical orthogonal
polynomials. Laguerre polynomials arise by considering orthogonality
conditions with respect to the gamma density $x^{\alpha}e^{-x}$, where
$\alpha>-1$, on $[0,\infty)$, see, e.g., [33, Chapter 5] for some properties.
Several extensions of these polynomials have been studied before. For example,
one can form an Angelesco system by considering the weights
$x^{\alpha}e^{-\left|x\right|}$ or $x^{\alpha}e^{-x^{r}}$, where $\alpha>-1$,
on $r$ intervals that form an $r$-star, see [32] and [23] respectively. Known
extensions to AT-systems are $(x^{\alpha}e^{-x},x^{\beta}e^{-x})$, where
$\alpha,\beta>-1$ and $\alpha-\beta\not\in\mathbb{Z}$, and
$(x^{\alpha}e^{-c_{1}x},x^{\alpha}e^{-c_{2}x})$, where $\alpha>-1$ and
$c_{1},c_{2}>0$ are distinct, both on $[0,\infty)$. The resulting multiple
orthogonal polynomials are then called multiple Laguerre polynomials of the
first and second kind and were introduced in [32] and [27]. Their properties
were investigated in [36] and [4]. The type I and type II multiple orthogonal
polynomials in this article (Section 2 and 3 respectively) can be seen as
multiple Laguerre polynomials of the third kind.

Observe that for the multiple Laguerre polynomials of the first and second
kind the additional weight function is still a gamma density, but with
different parameters. We propose to add a different kind of orthogonality
weight as the second weight: one that is related to the exponential integral.
The (generalized) exponential integral is defined on $\mathbb{R}_{>0}$ as

$E_{\nu}(x)=\int_{1}^{\infty}\frac{e^{-xt}}{t^{\nu}}dt,$

for every $\nu\in\mathbb{R}$. This function is connected to the incomplete
gamma function, given by

$\Gamma_{\nu}(x)=\int_{x}^{\infty}t^{\nu-1}e^{-t}dt,$

via the relation $E_{\nu}(x)=x^{\nu-1}\Gamma_{-\nu+1}(x)$. Orthogonal
polynomials with respect to the exponential integral have been investigated
before via a numerical approach in [13]. We believe that it is more natural to
consider multiple orthogonality and include a related gamma density as a
weight in the system as well.

We will consider a system of two weights $(w_{1},w_{2})$ on $[0,\infty)$ with

$w_{1}(x)=x^{\alpha}e^{-x},\quad w_{2}(x)=x^{\alpha}E_{\nu+1}(x).$

In order for all the moments of $w_{1}$ and $w_{2}$ to exist, we will assume
that $\alpha>-1$ and $\alpha+\nu>-1$. The moments are then given by

$\int_{0}^{\infty}x^{k}w_{1}(x)dx=\Gamma(k+\alpha+1),\quad\int_{0}^{\infty}x^{k}w_{2}(x)dx=\frac{\Gamma(k+\alpha+1)}{k+\alpha+\nu+1},$

where

$\Gamma(s)=\int_{0}^{\infty}x^{s-1}e^{-x}dx$

denotes the gamma function. We will also exclude $\nu\in\mathbb{Z}_{\leq-1}$
because for these values the weights will not be linearly independent over the
polynomials and the system won’t be perfect.

The described system is perfect if $\nu>-1$, because it can be written as a
Nikishin system; the ratio $w_{2}/w_{1}$ is given by the Stieltjes transform
of a gamma density on the negative real line. Indeed, if $\nu>-1$, we can use
[16, Eq. 3.383.10] to write

$\frac{w_{2}(x)}{w_{1}(x)}=\frac{1}{\Gamma(\nu+1)}\int_{-\infty}^{0}\frac{(-t)^{\nu}e^{t}}{x-t}dt.$

For $\nu=-n+v\leq-1$ with $n\in\mathbb{N}$ and $-1<v\leq 0$, the system is
close to a Nikishin system in the sense that the ratio $w_{2}/w_{1}$ will be
equal to $x^{-n}$ times the Stieltjes transform of a gamma density on
$(-\infty,0]$, plus a rational function with a pole of order $n$ at $0$. To
show this, we can apply [16, Eq. 8.356.5] to obtain

$\frac{w_{2}(x)}{w_{1}(x)}=\frac{(-v)_{n}}{x^{n}}e^{x}E_{v+1}(x)+\frac{1}{x^{n}}\sum_{j=0}^{n-1}\frac{x^{j}}{\Gamma(j-v+1)},$

where $(z)_{n}=z(z+1)\dots(z+n-1)$ denotes the Pochhammer symbol, and we can
use the result in the previous case. Note that if $\nu\in\mathbb{Z}_{\leq-1}$
(so $v=0$) the Stieltjes transform disappears and only the rational function
remains, which means that the system $(w_{1},w_{2})$ is not perfect as stated
before.

Another property of the system of weights is that it satisfies the matrix
Pearson equation

$\left[x\begin{pmatrix}w_{1}\\\
w_{2}\end{pmatrix}\right]^{\prime}=\begin{pmatrix}\alpha+1-x&0\\\
-1&\alpha+\nu+1\end{pmatrix}\begin{pmatrix}w_{1}\\\ w_{2}\end{pmatrix}.$

In Section 4, we will extend the results about the type I and type II multiple
orthogonal polynomials to certain mixed type multiple orthogonal polynomials.
First, we will consider mixed type orthogonality with respect to the two pairs
of weights $(e^{-x},E_{\nu+1}(x))$ and $(x^{\alpha},x^{\beta})$ on
$[0,\infty)$. The resulting mixed type multiple orthogonal polynomials are of
interest because they resemble the rational approximants used by Apéry in his
proof of irrationality of $\zeta(3)$. We believe this new family of mixed type
multiple orthogonal polynomials may be relevant in the rational approximation
of Euler’s constant $\gamma$. Afterwards, we will investigate the dual setting
in which we will reverse the order of the weights. In both settings, we will
assume that $\alpha,\beta,\nu,\alpha+\nu,\beta+\nu>-1$ for the same reasons as
before. Additionally, we will exclude $\alpha-\beta\in\mathbb{Z}$ in order for
the pair of weights $(x^{\alpha},x^{\beta})$ to be linearly independent over
the polynomials. It is known that in that case the pair of weights forms an
AT-system on $(0,\infty)$. In fact, under the additional restriction that
$0<\alpha-\beta<1$, it is even a Nikishin system because of the relation (see
[8, Eq. 5.12.3])

$x^{\beta-\alpha}=\frac{1}{\Gamma(1-\alpha+\beta)\Gamma(\alpha-\beta)}\int_{-\infty}^{0}\frac{(-t)^{\beta-\alpha}}{x-t}dt.$

##  2 Type I multiple orthogonal polynomials

A type I multiple orthogonal polynomial for the two weight functions
$(x^{\alpha}e^{-x},x^{\alpha}E_{\nu+1}(x))$ that corresponds to the multi-
index $(n,m)$ is a vector of polynomials $(A_{n,m},B_{n,m})$ with $\deg
A_{n,m}=n-1$ and $\deg B_{n,m}=m-1$ for which the type I function

$F_{n,m}(x)=A_{n,m}(x)e^{-x}+B_{n,m}(x)E_{\nu+1}(x)$ (1)

satisfies the orthogonality conditions

$\int_{0}^{\infty}F_{n,m}(x)x^{\alpha+k}dx=0,\quad k=0,\dots,n+m-2.$

Such polynomials (having maximal degree) exist for every multi-index because
of the perfectness of the system and they are unique up to a non-zero scalar
multiplication. We may normalize the type I function by considering the
normalization condition

$\int_{0}^{\infty}F_{n,m}(x)x^{\alpha+n+m-1}dx=1.$ (2)

The type I function for which this holds will be denoted by $F_{n,m}^{(I)}$.

###  2.1 Convolution and Rodrigues formula

We will be able to identify the type I function, and hence the underlying type
I multiple orthogonal polynomials, for multi-indices with $n+1\geq m$ by
careful analysis of its Mellin transform. The Mellin transform of a locally
Lebesgue integrable function $f:(0,\infty)\to\mathbb{C}$ is defined by the
relation

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

whenever the integral exists.

######  Lemma 2.1.

Suppose that $n+1\geq m$. The Mellin transform of the normalized type I
function $F_{n,m}^{(I)}$ is given by $\mathcal{F}_{n,m}\cdot\hat{F}_{n,m}(s)$,
where

$\hat{F}_{n,m}(s)=\frac{\Gamma(s)}{(s+\nu)_{m}}(\alpha+1-s)_{n+m-1},\quad
s>\max\\{0,-\nu\\},$ (3)

and

$\mathcal{F}_{n,m}=(-1)^{n+m-1}\frac{(\alpha+\nu+1)_{n+2m-1}}{(n+m-1)!(\alpha+1)_{n+m-1}(\alpha+\nu+1)_{n+m-1}\Gamma(\alpha+1)}.$
(4)

###### Proof.

Take a type I function $F_{n,m}$ and expand the polynomials $A_{n,m}(x)$ and
$B_{n,m}(x)$ as

$A_{n,m}(x)=\sum_{k=0}^{n-1}a_{k}x^{k},\quad
B_{n,m}(x)=\sum_{l=0}^{m-1}b_{l}x^{l}.$

The Mellin transform $\hat{F}_{n,m}(s)$ of
$F_{n,m}(x)=A_{n,m}(x)e^{-x}+B_{n,m}(x)E_{\nu+1}(x)$ then becomes

$\hat{F}_{n,m}(s)=\sum_{k=0}^{n-1}a_{k}\int_{0}^{\infty}e^{-x}x^{k+s-1}dx+\sum_{l=0}^{m-1}b_{l}\int_{0}^{\infty}E_{\nu+1}(x)x^{l+s-1}dx,$

which may be written as

$\sum_{k=0}^{n-1}a_{k}\Gamma(k+s)+\sum_{l=0}^{m-1}b_{l}\frac{\Gamma(l+s)}{l+s+\nu}.$

We can take a common factor out of the sums and get

$\Gamma(s)\left(\sum_{k=0}^{n-1}a_{k}(s)_{k}+\sum_{l=0}^{m-1}b_{l}\frac{(s)_{l}}{l+s+\nu}\right)$

Writing the terms of the second sum on a common denominator then gives

$\Gamma(s)\left(p_{n-1}(s)+\frac{q_{2m-2}(s)}{(s+\nu)_{m}}\right)$

in terms of polynomials $p_{n-1}(s)$ and $q_{2m-2}(s)$ with their subscript as
maximal degree. Finally, we take out the denominator $(s+\nu)_{m}$ to obtain

$\hat{F}_{n,m}(s)=\frac{\Gamma(s)}{(s+\nu)_{m}}\left(p_{n-1}(s)(s+\nu)_{m}+q_{2m-2}(s)\right).$

Observe that the degree of the polynomial between the parentheses is less than
or equal to $\max\\{n+m-1,2m-2\\}$. The restriction $m-1\leq n$ now implies
that this maximum is always equal to $n+m-1$. On the other hand, we know that
$\hat{F}_{n,m}(s)$ has zeros at the $n+m-1$ values
$s=\alpha+1,\dots,\alpha+n+m-1$ because of the orthogonality conditions that
the type I function $F_{n,m}$ satisfies. The polynomial between the brackets
is therefore necessarily a non-zero scalar multiple of $(\alpha+1-s)_{n+m-1}$.
Note that $\hat{F}_{n,m}(\alpha+n+m)\neq 0$. The normalization condition (2)
corresponds to taking the scalar
$\mathcal{F}_{n,m}=1/\hat{F}_{n,m}(\alpha+n+m)$, which readily leads to the
stated expression in (4). ∎

We will compare the above to the Mellin transform of a Laguerre polynomial
multiplied by a gamma density, which will be computed below. We will denote by
$L_{n}^{(\alpha)}(x)$ the Laguerre polynomial of degree $n$ with leading
coefficient $(-1)^{n}/n!$ (as in [33, Chapter 5]). It satisfies the
orthogonality conditions

$\int_{0}^{\infty}L_{n}^{(\alpha)}(x)e^{-x}x^{k+\alpha}dx=0,\quad
k=0,\dots,n-1,$

which will also be implied by its Mellin transform.

######  Lemma 2.2.

The Mellin transform of $L_{n}^{(\alpha)}(x)e^{-x}$ is given by

$\int_{0}^{\infty}L_{n}^{(\alpha)}(x)e^{-x}x^{s-1}\,dx=\Gamma(s)\frac{(\alpha+1-s)_{n}}{n!},\quad
s>0.$

###### Proof.

We may use the Rodrigues formula for Laguerre polynomials (see [33, Eq.
(5.1.5)])

$L_{n}^{(\alpha)}(x)=\frac{x^{-\alpha}e^{x}}{n!}\frac{d^{n}}{dx^{n}}\left[x^{n+\alpha}e^{-x}\right]$
(5)

and perform integration by parts $n$ times to get

$\int_{0}^{\infty}L_{n}^{(\alpha)}(x)e^{-x}x^{s-1}dx=\frac{(-1)^{n}}{n!}\int_{0}^{\infty}x^{\alpha+n}e^{-x}\left(\frac{d^{n}}{dx^{n}}x^{s-\alpha-1}\right)dx.$

After working out the derivatives, this becomes

$\int_{0}^{\infty}L_{n}^{(\alpha)}(x)e^{-x}x^{s-1}dx=\frac{(\alpha+1-s)_{n}}{n!}\int_{0}^{\infty}x^{s-1}e^{-x}dx.$

The remaining integral is the gamma function $\Gamma(s)$. ∎

The type I function can now be determined explicitly; it arises as the Mellin
convolution of a Laguerre polynomial multiplied by a gamma density and a beta
density. The Mellin convolution of two appropriate functions $f$ and $g$ is
given by

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

Its key property is that its Mellin transform is multiplicative in terms of
the Mellin transforms of the input functions

$\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).$

######  Theorem 2.3.

Suppose that $n+1\geq m$. The normalized type I function $F_{n,m}^{(I)}$ is
given by $(n+m-1)!\mathcal{F}_{n,m}F_{n,m}/(m-1)!$, where $\mathcal{F}_{n,m}$
is the constant in (4) and

$F_{n,m}(x)=\int_{x}^{\infty}L_{n+m-1}^{(\alpha)}(t)e^{-t}\Bigl{(}1-\frac{x}{t}\Bigr{)}^{m-1}\Bigl{(}\frac{x}{t}\Bigr{)}^{\nu}\,\frac{dt}{t}.$
(6)

###### Proof.

The Mellin transform of the normalized type I function $F_{n,m}^{(I)}$, see
Lemma 2.1, arises as a non-zero scalar multiple of the product of the Mellin
transform of a Laguerre polynomial multiplied with a gamma density

$\int_{0}^{\infty}L_{n+m-1}^{(\alpha)}(x)e^{-x}x^{s-1}dx=\Gamma(s)\frac{(\alpha+1-s)_{n+m-1}}{(n+m-1)!},\quad
s>0,$

see Lemma 2.2, and of a beta density

$\int_{0}^{\infty}(1-x)^{m-1}x^{\nu+s-1}\chi_{[0,1]}(x)dx=B(s+\nu,m)=\frac{(m-1)!}{(s+\nu)_{m}},\quad
s>-\nu.$

The multiplicative property of the Mellin convolution and the uniqueness of
the Mellin transform then imply that $F_{n,m}^{(I)}$ is given by the Mellin
convolution of the functions$(n+m-1)!L_{n+m-1}^{(\alpha)}(x)e^{-x}$ and
$(1-x)^{m-1}x^{\nu}\chi_{[0,1]}(x)/(m-1)!$, which gives the integral in (6). ∎

By working out the convolution formula for the type I function in (6), we may
obtain explicit expressions for the underlying type I polynomials $A_{n,m}$
and $B_{n,m}$. However, since this is rather tedious and not very insightful,
we won’t do it here. Later in Proposition 2.7, we will be able to derive an
expression for $B_{n,m}$ in a more elegant way.

We can use the convolution formula for the type I function to derive a
Rodrigues formula for it. The relevant weight functions are those of the
associated Nikishin system, i.e. $w_{1}$ and $w_{2}/w_{1}$.

######  Theorem 2.4.

Suppose that $n+1\geq m$. The type I function in (6) is given by the Rodrigues
formula

$F_{n,m}(x)=\frac{x^{-\alpha}}{(n+m-1)!(\nu+1)_{m-1}}\frac{d^{n+m-1}}{dx^{n+m-1}}\left[x^{n+m-1}w_{1}(x)\frac{d^{m-1}}{dx^{m-1}}\left[x^{m-1}\frac{w_{2}(x)}{w_{1}(x)}\right]\right].$
(7)

###### Proof.

The result will follow from the Rodrigues formula for the underlying Laguerre
polynomial. We first perform the change of variables $t\mapsto xt$ in (6) to
get

$F_{n,m}(x)=\int_{1}^{\infty}L_{n+m-1}^{(\alpha)}(xt)e^{-xt}\left(1-\frac{1}{t}\right)^{m-1}\frac{dt}{t^{\nu+1}}.$

From the Rodrigues formula $\eqref{Laguerre_Rodr}$ for
$L_{n+m-1}^{(\alpha)}(x)$ and the chain rule, it follows that

$L_{n+m-1}^{(\alpha)}(xt)e^{-xt}=\frac{(xt)^{-\alpha}}{(n+m-1)!}\frac{d^{n+m-1}}{dx^{n+m-1}}\left[(xt)^{n+m-1+\alpha}e^{-xt}\right]t^{-n-m+1},$

or after rearranging some factors depending on $t$,

$L_{n+m-1}^{(\alpha)}(xt)e^{-xt}=\frac{x^{-\alpha}}{(n+m-1)!}\frac{d^{n+m-1}}{dx^{n+m-1}}\left[x^{n+m-1+\alpha}e^{-xt}\right].$

We can use this to take $n+m-1$ derivatives outside of the integral

$F_{n,m}(x)=\frac{x^{-\alpha}}{(n+m-1)!}\frac{d^{n+m-1}}{dx^{n+m-1}}\left[x^{n+m-1+\alpha}\int_{1}^{\infty}e^{-xt}\left(1-\frac{1}{t}\right)^{m-1}\frac{dt}{t^{\nu+1}}\right]$

(the integrand is smooth for $x\in\mathbb{R}$ and $t\geq 1$, so an application
of the Leibniz integral rule is justified). The remaining integral is related
to the confluent hypergeometric function

$U(a,b,x)=\frac{1}{\Gamma(a)}\int_{0}^{\infty}e^{-xt}t^{a-1}(1+t)^{b-a-1}\,dt,$

which is defined for $x,a>0$ and $b\in\mathbb{R}$ (see [8, Eq. 13.4.4]), via
the change of variables $t\mapsto 1+t$

$\int_{1}^{\infty}e^{-xt}\left(1-\frac{1}{t}\right)^{m-1}\frac{dt}{t^{\nu+1}}=(m-1)!e^{-x}U(m,-\nu+1;x).$

The latter can be written in terms of

$U(1,-\nu+1;x)=e^{x}\int_{1}^{\infty}e^{-xt}t^{-\nu-1}\,dt=e^{x}E_{\nu+1}(x)$

by making use of the differentiation formula (see [8, Eq. 13.3.24])

$\frac{d^{n}}{dz^{n}}\left[x^{n+a-1}U(a,b;z)\right]=(a)_{n}(a-b+1)_{n}x^{a-1}U(a+n,b;z).$

Indeed, the special case

$\frac{d^{m-1}}{dx^{m-1}}\left[x^{m-1}U(1,-\nu+1;x)\right]=(m-1)!(\nu+1)_{m-1}U(m,-\nu+1;x)$

leads to the desired result. ∎

###  2.2 Hypergeometric series

The Mellin transform of the type I function can be used to derive a
representation in terms of certain generalized hypergeometric series. The
generalized hypergeometric function ${}_{p}F_{q}$ is a function that depends
on two tuples $(a_{1},\dots,a_{p})\in\mathbb{C}^{p}$ and
$(b_{1},\dots,b_{q})\in(\mathbb{C}\backslash\mathbb{Z}_{\leq 0})^{q}$ and is
defined, for $p\leq q+1$, by the power series

${}_{p}F_{q}\left(\begin{array}[]{c}a_{1},\dots,a_{p}\\\
b_{1},\dots,b_{q}\end{array};x\right)=\sum_{k=0}^{\infty}\frac{(a_{1})_{k}\dots(a_{p})_{k}}{(b_{1})_{k}\dots(b_{q})_{k}}\frac{x^{k}}{k!},$

see, e.g., [8, Chapter 16] for some basic properties.

######  Theorem 2.5.

Suppose that $n+1\geq m$ and $\nu\not\in\mathbb{Z}$. The type I function with
Mellin transform (3) arises as the following combination of ${}_{2}F_{2}$
hypergeometric series

$\displaystyle F_{n,m}(x)=$
$\displaystyle\frac{(\alpha+1)_{n+m-1}}{(\nu)_{m}}{}_{2}F_{2}\left(\begin{array}[]{c}n+m+\alpha,-m-\nu+1\\\
\alpha+1,-\nu+1\end{array};-x\right)$
$\displaystyle+\frac{(\alpha+\nu+1)_{n+m-1}\Gamma(-\nu)}{(m-1)!}{}_{2}F_{2}\left(\begin{array}[]{c}-m+1,n+m+\alpha+\nu\\\
\nu+1,\alpha+\nu+1\end{array};-x\right)x^{\nu}.$

###### Proof.

It will be more convenient to work with the Mellin transform of the complete
type I function $x^{\alpha}F_{n,m}(x)$, which is

$\int_{0}^{\infty}F_{n,m}(x)x^{\alpha+s-1}dx=\frac{\Gamma(s+\alpha)}{(s+\alpha+\nu)_{m}}(1-s)_{n+m-1},\quad
s>-\min\\{\alpha,\alpha+\nu\\}.$

Recall that the gamma function has an analytic continuation to
$\mathbb{C}\backslash\mathbb{Z}_{\leq 0}$ that has simple poles at all
$k\in\mathbb{Z}_{\leq 0}$ with residue $(-1)^{k}/k!$. We can therefore extend
the above function in $s$ analytically to the complex plane except for the
simple poles $-\alpha,-\alpha-1,\dots$ and $-\alpha-\nu,\dots,-\alpha-\nu-m+1$
(here we use that $\nu\not\in\mathbb{Z}$). The described poles all lie in an
interval $(-\infty,c)$ for some $c\in\mathbb{R}$ that satisfies
$-\min\\{\alpha,\alpha+\nu,0\\}<c<1$ (here we use that
$\alpha,\alpha+\nu>-1$). The original function can then be recovered via the
inverse Mellin transform

$x^{\alpha}F_{n,m}(x)=\int_{c-i\infty}^{c+i\infty}\frac{\Gamma(s+\alpha)}{(s+\alpha+\nu)_{m}}(1-s)_{n+m-1}x^{-s}ds$

by making use of the residue theorem. For this we can consider a
counterclockwise oriented rectangle that connects the points $c\pm iN$ and
$(c-M)\pm iN$, add the residues of the poles in this region and let
$N,M\to\infty$. The residue of the integrand at $s=-\alpha-k$ with
$k\in\mathbb{Z}_{\geq 0}$ is

$\frac{(-1)^{k}}{k!}\frac{(k+\alpha+1)_{n+m-1}}{(-k+\nu)_{m}}x^{\alpha+k}=\frac{(-1)^{k}}{k!}\frac{(n+m+\alpha)_{k}(\alpha+1)_{n+m-1}(-m-\nu+1)_{k}}{(\alpha+1)_{k}(-\nu+1)_{k}(\nu)_{m}}x^{\alpha+k}$

and the residue at $s=-\alpha-\nu-k$ with $k\in\\{0,\dots,m-1\\}$ is

$\displaystyle\frac{\Gamma(-k-\nu)}{(-k)_{k}(1)_{m-k-1}}$
$\displaystyle(k+\alpha+\nu+1)_{n+m-1}x^{\alpha+\nu+k}$
$\displaystyle=\frac{(-1)^{k}}{k!}\frac{\Gamma(-\nu)(n+m+\alpha+\nu)_{k}(\alpha+\nu+1)_{n+m-1}(-m+1)_{k}}{(\nu+1)_{k}(\alpha+\nu+1)_{k}(m-1)!}x^{\alpha+\nu+k}.$

These are exactly the terms in the stated expression for $F_{n,m}$. ∎

The fact that part of the above expression is a polynomial times $x^{\nu}$ is
to be expected, because of the following expansions for the weight functions
(see [8, Eq. 8.19.10]) of which $F_{n,m}$ is a polynomial combination:

$e^{-x}={}_{0}F_{0}\left(\begin{array}[]{c}-\\\ -\end{array};-x\right),\quad
E_{\nu+1}(x)=\Gamma(-\nu)x^{\nu}+\frac{1}{\nu}\,{}_{1}F_{1}\left(\begin{array}[]{c}-\nu\\\
-\nu+1\end{array};-x\right).$

It also suggests that this polynomial will be related to $B_{n,m}(x)$. We will
make this rigorous in what follows.

The connection to certain ${}_{2}F_{2}$ hypergeometric series allows us to
obtain many properties of the type I function. For example, we can use it to
determine a differential equation.

######  Proposition 2.6.

Suppose that $n+1\geq m$ and $\nu\not\in\mathbb{Z}$. Every type I function
$y=F_{n,m}(x)$ satisfies the differential equation

$x^{2}y^{\prime\prime\prime}+x(\alpha-\nu+3+x)y^{\prime\prime}+((\alpha+1)(-\nu+1)+x(n+\alpha-\nu+2))y^{\prime}+(n+m+\alpha)(-m-\nu+1)y=0.$

###### Proof.

In general, it is known that a ${}_{p}F_{q}$ hypergeometric series satisfies a
linear differential equation of order $\max\\{p,q+1\\}$, see [8, Eq. 16.8.3]
or [28, §16.8 (ii)]. For the ${}_{2}F_{2}$ hypergeometric series

${}_{2}F_{2}\left(\begin{array}[]{c}a_{1},a_{2}\\\
b_{1},b_{2}\end{array};-x\right),$

this differential equation becomes

$x^{2}y^{\prime\prime\prime}+x(b_{1}+b_{2}+1+x)y^{\prime\prime}+(b_{1}b_{2}+x(a_{1}+a_{2}+1))y^{\prime}+a_{1}a_{2}y=0.$
(8)

There are two other fundamental solutions and they are given by

$x^{1-b_{1}}{}_{2}F_{2}\left(\begin{array}[]{c}a_{1}-b_{1}+1,a_{2}-b_{1}+1\\\
-b_{1}+2,b_{2}-b_{1}+1\end{array};-x\right),\quad
x^{1-b_{2}}{}_{2}F_{2}\left(\begin{array}[]{c}a_{1}-b_{2}+1,a_{2}-b_{2}+1\\\
b_{1}-b_{2}+1,-b_{2}+2\end{array};-x\right),$

see [8, Eq. 16.8.6]. The stated differential equation is the differential
equation that corresponds to the hypergeometric series in the first term of
the formula in Theorem 2.5. In this case, the two other solutions are of the
form

$x^{-\alpha}{}_{2}F_{2}\left(\begin{array}[]{c}n+m,-m-\alpha-\nu+1\\\
-\alpha+1,-\alpha-\nu+1\end{array};-x\right),\quad
x^{\nu}{}_{2}F_{2}\left(\begin{array}[]{c}n+m+\alpha+\nu,-m+1\\\
\alpha+\nu+1,\nu+1\end{array};-x\right).$

The type I function $F_{n,m}(x)$ thus arises as a linear combination of two
solutions of a linear differential equation and is therefore a solution itself
as well. ∎

We can use the differential equation for the type I function to obtain an
explicit expression for the type I polynomial $B_{n,m}(x)$.

######  Proposition 2.7.

Suppose that $n+1\geq m$ and $\nu\not\in\mathbb{Z}$. Then

$B_{n,m}(x)=B_{n,m}(0)\cdot{}_{2}F_{2}\left(\begin{array}[]{c}-m+1,n+m+\alpha+\nu\\\
\nu+1,\alpha+\nu+1\end{array};-x\right).$

###### Proof.

Consider the differential equation in Proposition 2.6, but take the
derivatives of $B_{n,m}(x)E_{\nu+1}(x)$ in
$F_{n,m}(x)=A_{n,m}(x)e^{-x}+B_{n,m}(x)E_{\nu+1}(x)$ via
$(B_{n,m}(x)x^{\nu})\cdot\Gamma_{-\nu}(x)$ and use that
$(\Gamma_{-\nu}(x))^{\prime}=-x^{-\nu-1}e^{-x}$. Since the functions $e^{-x}$
and $\Gamma_{-\nu}(x)$ are linearly independent over the polynomials (compare
for example their Mellin transforms), $B_{n,m}(x)x^{\nu}$ itself must be a
solution to the differential equation as well. Only one of the fundamental
solutions listed in the proof of Proposition 2.6 is of the form $x^{\nu}$
times a polynomial, hence we must have the stated result. ∎

###  2.3 Recurrence relations

It follows from the general theory of multiple orthogonal polynomials that
type I (and type II) multiple orthogonal polynomials satisfy certain nearest
neighbor recurrence relations, see, e.g., [35]. For the normalized type I
functions here, these are of the form

$\displaystyle xF_{n,m}^{(I)}(x)$ $\displaystyle=$ $\displaystyle
F_{n-1,m}^{(I)}(x)+c_{n-1,m}F_{n,m}^{(I)}(x)+a_{n,m}F_{n+1,m}^{(I)}(x)+b_{n,m}F_{n,m+1}^{(I)}(x),$
(9) $\displaystyle xF_{n,m}^{(I)}(x)$ $\displaystyle=$ $\displaystyle
F_{n,m-1}^{(I)}(x)+d_{n,m-1}F_{n,m}^{(I)}(x)+a_{n,m}F_{n+1,m}^{(I)}(x)+b_{n,m}F_{n,m+1}^{(I)}(x).$
(10)

We will determine the coefficients in two ways. Our first proof will be based
on the three-term recurrence relation for Laguerre polynomials, which is (see
[33])

$-xL_{n}^{(\alpha)}(x)=(n+1)L_{n+1}^{(\alpha)}(x)-(2n+\alpha+1)L_{n}^{(\alpha)}(x)+(n+\alpha)L_{n-1}^{(\alpha)}(x).$
(11)

In the second proof, we will only use the Mellin transforms of the type I
functions. It will be more straightforward but less constructive.

The first proof requires the following connection formula. A posteriori, it
can also be obtained by subtracting the two nearest neighbor recurrence
relations (9)–(10).

######  Proposition 2.8.

One has

$F_{n,m-1}^{(I)}(x)-F_{n-1,m}^{(I)}(x)=\frac{(n+m-1)(n+m+\alpha-1)(n+m+\alpha+\nu-1)}{(n+2m+\alpha+\nu-2)(n+2m+\alpha+\nu-1)}F_{n,m}^{(I)}(x).$

###### Proof.

In terms of the type I functions with Mellin transforms (3), we have

$\hat{F}_{n-1,m}(s)=\frac{1}{\alpha+n+m-s-1}\hat{F}_{n,m}(s)$

and

$\hat{F}_{n,m-1}(s)=\frac{s+\nu+m-1}{\alpha+n+m-s-1}\hat{F}_{n,m}(s).$

Hence

$(\alpha+\nu+n+2m-2)\hat{F}_{n-1,m}(s)-\hat{F}_{n,m-1}(s)=\hat{F}_{n,m}(s),$

which for the normalized functions becomes

$(n+2m+\alpha+\nu-2)\frac{\mathcal{F}_{n,m}}{\mathcal{F}_{n-1,m}}\hat{F}_{n-1,m}^{(I)}(x)-\frac{\mathcal{F}_{n,m}}{\mathcal{F}_{n,m-1}}\hat{F}_{n,m-1}^{(I)}(x)=\hat{F}_{n,m}^{(I)}(x).$

If we plug in the values of the ratios

$\displaystyle\frac{\mathcal{F}_{n,m}}{\mathcal{F}_{n-1,m}}$ $\displaystyle=$
$\displaystyle-\frac{(n+2m+\alpha+\nu-1)}{(n+m-1)(n+m+\alpha-1)(n+m+\alpha+\nu-1)},$
$\displaystyle\frac{\mathcal{F}_{n,m}}{\mathcal{F}_{n,m-1}}$ $\displaystyle=$
$\displaystyle-\frac{(n+2m+\alpha+\nu-2)(n+2m+\alpha+\nu-1)}{(n+m-1)(n+m+\alpha-1)(n+m+\alpha+\nu-1)},$

the desired result follows after applying the inverse Mellin transform. ∎

We can now derive the nearest neighbor recurrence relations for the normalized
type I functions.

######  Theorem 2.9.

The coefficients in the nearest neighbor recurrence relations (9)–(10) are
given by

$\displaystyle a_{n,m}$
$\displaystyle=\frac{(n+m)(n+m+\alpha)(n+m+\alpha+\nu)}{n+2m+\alpha+\nu},$
$\displaystyle b_{n,m}$
$\displaystyle=-\frac{(n+m)(n+m+\alpha)(n+m+\alpha+\nu)m(m+\nu)(m+\alpha+\nu)}{(n+2m+\alpha+\nu-1)(n+2m+\alpha+\nu)^{2}(n+2m+\alpha+\nu+1)},$
$\displaystyle c_{n-1,m}$ $\displaystyle=a_{n,m}-a_{n-1,m},$ $\displaystyle
d_{n,m-1}$ $\displaystyle=a_{n,m}-a_{n,m-1}.$

###### Proof.

Consider the type I functions in (6). We can pass the multiplication by $x$
through the convolution as follows

$xF_{n,m}(x)=\int_{x}^{\infty}\left[tL_{n+m-1}^{(\alpha)}(t)\right]e^{-t}\left[\frac{x}{t}\left(1-\frac{x}{t}\right)^{m-1}\right]\left(\frac{x}{t}\right)^{\nu}\frac{dt}{t}.$

If we then combine the three-term recurrence relation (11) for Laguerre
polynomials and the recurrence relation $x(1-x)^{m-1}=(1-x)^{m-1}-(1-x)^{m}$,
we find

$\displaystyle xF_{n,m}(x)=$ $\displaystyle(n+m)[F_{n,m+1}(x)-F_{n+1,m}]$
$\displaystyle-\,(2n+2m+\alpha-1)[F_{n-1,m+1}(x)-F_{n,m}(x)]$
$\displaystyle+\,(n+m+\alpha-1)[F_{n-2,m+1}-F_{n-1,m}(x)],$

For the normalized type I functions, this becomes

$\displaystyle xF_{n,m}^{(I)}(x)=$
$\displaystyle\left[\frac{m\mathcal{F}_{n,m}}{\mathcal{F}_{n,m+1}}F_{n,m+1}^{(I)}(x)-\frac{\mathcal{F}_{n,m}}{\mathcal{F}_{n+1,m}}F_{n+1,m}^{(I)}(x)\right]$
$\displaystyle-\
(2n+2m+\alpha-1)\left[\frac{m\mathcal{F}_{n,m}}{\mathcal{F}_{n-1,m+1}}F_{n-1,m+1}^{(I)}(x)-F_{n,m}^{(I)}(x)\right]$
$\displaystyle+\
(n+m-1)(n+m+\alpha-1)\left[\frac{m\mathcal{F}_{n,m}}{\mathcal{F}_{n-2,m+1}}F_{n-2,m+1}^{(I)}(x)-\frac{\mathcal{F}_{n,m}}{\mathcal{F}_{n-1,m}}F_{n-1,m}^{(I)}(x)\right].$

In order to obtain the nearest neighbor recurrence relations, we need to get
rid of the terms $F_{n-1,m+1}^{(I)}(x)$ and $F_{n-2,m+1}^{(I)}(x)$. We can do
this by making use of the connection formula. One has

$F_{n-2,m+1}^{(I)}(x)=F_{n-1,m}^{(I)}(x)-C_{n-1,m+1}F_{n-1,m+1}^{(I)}(x),$
$F_{n-1,m+1}^{(I)}(x)=F_{n,m}^{(I)}(x)-C_{n,m+1}F_{n,m+1}^{(I)}(x).$

The coefficient of the function $F_{n-1,m}^{(I)}(x)$ then becomes

$(n+m-1)(n+m+\alpha+\nu-1)\left[\frac{m\mathcal{F}_{n,m}}{\mathcal{F}_{n-2,m+1}}-\frac{\mathcal{F}_{n,m}}{\mathcal{F}_{n-1,m}}\right]=1.$

The other coefficients become, using the same notation as in (9),

$\displaystyle a_{n,m}=$
$\displaystyle-\frac{\mathcal{F}_{n,m}}{\mathcal{F}_{n+1,m}},$ $\displaystyle
b_{n,m}=$
$\displaystyle\,\frac{m\mathcal{F}_{n,m}}{\mathcal{F}_{n,m+1}}+(2n+2m+\alpha-1)\frac{m\mathcal{F}_{n,m}}{\mathcal{F}_{n-1,m+1}}C_{n,m+1}$
$\displaystyle+(n+m-1)(n+m+\alpha-1)\frac{m\mathcal{F}_{n,m}}{\mathcal{F}_{n-2,m+1}}C_{n-1,m+1}C_{n,m+1},$
$\displaystyle c_{n-1,m}=$ $\displaystyle-\
(2n+2m+\alpha-1)\left[\frac{m\mathcal{F}_{n,m}}{\mathcal{F}_{n-1,m+1}}-1\right]$
$\displaystyle+(n+m-1)(n+m+\alpha-1)\frac{m\mathcal{F}_{n,m}}{\mathcal{F}_{n-2,m+1}}C_{n-1,m+1}C_{n-1,m+1}.$

After some computations, we obtain the expressions described in the theorem.
For (10), we can use the connection formula to replace $F_{n-1,m}^{(I)}(x)$ in
(9) by $F_{n,m-1}^{(I)}(x)$ and an extra term containing $F_{n,m}^{(I)}(x)$. ∎

Alternative proof of Theorem 2.9. Consider the type I functions with Mellin
transform (3). We may take the Mellin transform of the linear form

$c_{1}F_{n-1,m}(x)+c_{2}F_{n,m}(x)+c_{3}F_{n+1,m}(x)+c_{4}F_{n,m+1}(x)$

and show that for suitable values of the coefficients it will be equal to the
Mellin transform of $xF_{n,m}(x)$. By converting to the normalized type I
functions, we should then end up with the described recurrence coefficients.
After working out the described Mellin transforms, we would like to have

$\displaystyle c_{1}(s+\nu+m)+c_{2}(\alpha-s+n+m-1)(s+\nu+m)$
$\displaystyle+c_{3}(\alpha-s+n+m-1)(\alpha-s+n+m)(s+\nu+m)$
$\displaystyle+c_{4}(\alpha-s+n+m-1)(\alpha-s+n+m)$
$\displaystyle=s(s+\nu)(\alpha-s).$

By comparing the coefficients of $s$ on both sides of the equation, we obtain
four equations for the four unknown $c_{j}$. We can solve this system of
equations via a computer algebra system. We find

$\displaystyle c_{1}$
$\displaystyle=-\frac{(n+m-1)(n+m+\alpha+\nu-1)(n+m+\alpha-1)}{n+2m+\alpha+\nu-1},$
$\displaystyle c_{2}$
$\displaystyle=2n+m+\alpha-1+\frac{(m+n)m(m+\nu)}{n+2m+\alpha+\nu}-\frac{(m+n-1)m(m+\nu)}{n+2m+\alpha+\nu-1},$
$\displaystyle c_{3}$ $\displaystyle=-1,$ $\displaystyle c_{4}$
$\displaystyle=\frac{m(m+\alpha+\nu)(m+\nu)}{(\alpha+\nu+n+2m)(n+2m+\alpha+\nu-1)}.$

The associated recurrence relation for the normalized functions then has
coefficients

$\displaystyle c_{1}^{(I)}$
$\displaystyle=c_{1}\mathcal{F}_{n,m}/\mathcal{F}_{n-1,m}=-c_{1}\frac{(n+2m+\alpha+\nu-1)}{(n+m-1)(n+m+\alpha-1)(n+m+\alpha+\nu-1)}=1,$
$\displaystyle c_{2}^{(I)}$ $\displaystyle=c_{2}=a_{n,m}-a_{n-1,m},$
$\displaystyle c_{3}^{(I)}$
$\displaystyle=c_{3}\mathcal{F}_{n,m}/\mathcal{F}_{n+1,m}=-c_{3}\frac{(n+m)(n+m+\alpha)(n+m+\alpha+\nu)}{n+2m+\alpha}=a_{n,m},$
$\displaystyle c_{4}^{(I)}$
$\displaystyle=c_{4}\mathcal{F}_{n,m}/\mathcal{F}_{n,m+1}=-c_{4}\frac{(n+m)(n+m+\alpha)(n+m+\alpha+\nu)}{(n+2m+\alpha+\nu)(n+2m+\alpha+\nu+1)}=b_{n,m}.$

In order to obtain the formula for $c_{2}$, we used the fact that

$a_{n,m}=(n+m)(n+\alpha)+\frac{(n+m)m(m+\nu)}{n+2m+\alpha+\nu},$

which follows from a partial decomposition of $a_{n,m}$ as a function in
$\alpha$. ∎

##  3 Type II multiple orthogonal polynomials

A type II multiple orthogonal polynomial for the two weight functions
$(x^{\alpha}e^{-x},x^{\alpha}E_{\nu+1}(x))$ that corresponds to the multi-
index $(n,m)$ is a polynomial $P_{n,m}$ with $\deg P_{n,m}=n+m$ that satisfies
the orthogonality relations

$\displaystyle\int_{0}^{\infty}P_{n,m}(x)e^{-x}x^{k+\alpha}dx=0,\quad
k=0,\dots,n-1,$
$\displaystyle\int_{0}^{\infty}P_{n,m}(x)E_{\nu+1}(x)x^{k+\alpha}dx=0,\quad
k=0,\dots,m-1.$

Such polynomials (having maximal degree) exist for every multi-index because
of the perfectness of the system and they are unique up to a non-zero scalar
multiplication. The normalized type II polynomial, which is monic, will be
denoted by $P_{n,m}^{(II)}$. Similarly as before, we will only consider multi-
indices $(n,m)$ with $n+1\geq m$.

###  3.1 Differentiation formula

The type II polynomials will also be connected to the Laguerre polynomials;
now the connection will come in the form of a differential formula. In order
to obtain this formula, we require the following reduction property of the
type II polynomials. Note that we explicitly wrote down the dependence of the
type II polynomial $P_{n,m}=P_{n,m}^{\alpha,\nu}$ on the parameters.

######  Lemma 3.1.

Suppose that $n\geq m+1$. Then

$\frac{d}{dx}\left[x^{\alpha+\nu}P_{n,m}^{\alpha,\nu}(x)\right]\propto
x^{\alpha+\nu-1}P_{n-1,m+1}^{\alpha,\nu-1}(x),$

where the symbol $\propto$ means equality up to a non-zero scalar
multiplication.

###### Proof.

After working out the derivative, it is clear that the left hand side is of
the form $x^{\alpha+\nu-1}$ times a polynomial of degree $n+m$. It then
remains to verify that it satisfies the required orthogonality conditions,
i.e.,

$\int_{0}^{\infty}\frac{d}{dx}\left[x^{\alpha+\nu}P_{n,m}^{\alpha,\nu}(x)\right]e^{-x}x^{k-\nu+1}dx=0,\quad
k=0,\dots,n-2,$
$\int_{0}^{\infty}\frac{d}{dx}\left[x^{\alpha+\nu}P_{n,m}^{\alpha,\nu}(x)\right]E_{\nu}(x)x^{k-\nu+1}dx=0,\quad
k=0,\dots,m.$

We can rewrite both integrals via integration by parts. The first one becomes

$\left[P_{n,m}^{\alpha,\nu}(x)e^{-x}x^{k+\alpha+1}\right]_{0}^{\infty}-\int_{0}^{\infty}x^{\alpha+\nu}P_{n,m}^{\alpha,\nu}(x)\frac{d}{dx}\left[e^{-x}x^{k-\nu+1}\right]dx.$

The restriction $\alpha>-1$ implies that the first part will vanish. After
working out the derivative, we thus get

$-\int_{0}^{\infty}P_{n,m}^{\alpha,\nu}(x)e^{-x}x^{k+\alpha}(-x+(k-\nu+1))dx.$

This vanishes for $k=0,\dots,n-2$ due to the orthogonality of
$P_{n,m}^{\alpha,\nu}$ with respect to the first weight. The other integral
becomes

$\left[P_{n,m}^{\alpha,\nu}(x)E_{\nu}(x)x^{k+\alpha+1}\right]_{0}^{\infty}-\int_{0}^{\infty}x^{\alpha+\nu}P_{n,m}^{\alpha,\nu}(x)\frac{d}{dx}\left[E_{\nu}(x)x^{k-\nu+1}\right]dx,$

where the first part vanishes again because $\alpha>-1$. It is convenient to
use the relation $E_{\nu}(x)=x^{\nu-1}\Gamma_{-\nu+1}(x)$ when computing the
derivative

$\frac{d}{dx}\left[E_{\nu}(x)x^{k-\nu+1}\right]=-x^{k-\nu}e^{-x}+kx^{k-1}\Gamma_{-\nu+1}(x)=-x^{k-\nu}(e^{-x}-kE_{\nu}(x)).$

We can use the recurrence relation $\nu E_{\nu+1}(x)+xE_{\nu}(x)=e^{-x}$, see
[8, Eq. 8.19.12], to write this in terms of $E_{\nu+1}(x)$. The integral then
becomes

$\int_{0}^{\infty}P_{n,m}^{\alpha,\nu}(x)x^{k+\alpha-1}((x-k)e^{-x}+k\nu
E_{\nu+1}(x))dx.$

For $k=0$, this is equal to

$\int_{0}^{\infty}P_{n,m}^{\alpha,\nu}(x)x^{\alpha}e^{-x}dx=0,$

while for $k\geq 1$, we have

$\int_{0}^{\infty}P_{n,m}^{\alpha,\nu}(x)x^{k+\alpha-1}(x-k)e^{-x}dx=0,\qquad\text{
if }k\leq n-1$

and

$\int_{0}^{\infty}P_{n,m}^{\alpha,\nu}(x)x^{k+\alpha-1}E_{\nu+1}(x)dx=0,\qquad\text{
if }k\leq m.$

Since $m\leq n-1$, the desired orthogonality conditions with respect to the
second weight are satisfied as well. ∎

The previous property allows us to express the type II polynomial in terms of
a Laguerre polynomial.

######  Theorem 3.2.

Suppose that $n+1\geq m$. The monic type II multiple orthogonal polynomial
$P_{n,m}^{(II)}$ is given by $\mathcal{P}_{n,m}\cdot P_{n,m}$, where

$P_{n,m}(x)=x^{-\alpha-\nu}\frac{d^{m}}{dx^{m}}\left[x^{m+\alpha+\nu}L_{n+m}^{(\alpha)}(x)\right].$
(12)

and

$\mathcal{P}_{n,m}=(-1)^{n+m}(n+m)!\frac{(\alpha+\nu+1)_{n+m}}{(\alpha+\nu+1)_{n+2m}}.$

###### Proof.

We have already shown that

$x^{\alpha+\nu}P_{n,m}^{\alpha,\nu}(x)\propto\frac{d}{dx}\left[x^{\alpha+\nu+1}P_{n+1,m-1}^{\alpha,\nu+1}(x)\right],$

for $n+1\geq m$, so if we extract $m-1$ additional derivatives, we get

$x^{\alpha}P_{n,m}^{\alpha,\nu}(x)\propto\frac{d^{m}}{dx^{m}}\left[x^{\alpha+m}P_{n+m,0}^{\alpha,\nu+m}(x)\right].$

The remaining polynomial is only orthogonal with respect to the first weight
$x^{\alpha}e^{-x}$, hence we must have $P_{n+m,0}^{\alpha,\nu+m}(x)\propto
L_{n+m}^{(\alpha)}(x)$. The leading coefficient of this Laguerre polynomial is
$(-1)^{n+m}/(n+m)!$, so the leading coefficient of the monic type II
polynomial is $\mathcal{P}_{n,m}(-1)^{n+m}(n+m)!(n+m+\alpha+\nu+1)_{m}$. This
readily leads to the stated expression for the constant $\mathcal{P}_{n,m}$. ∎

Observe that if we would replace the Laguerre polynomial in this formula by
its Rodrigues formula (5), we would obtain a Rodrigues formula for the type II
polynomial as well.

###  3.2 Hypergeometric series

In what follows, we will show that the type II polynomials are hypergeometric
polynomials. This property will be inherited from the underlying Laguerre
polynomial, which is a hypergeometric polynomial as well. More precisely, we
have (see [33, Eq. (5.3.3)])

$L_{n}^{(\alpha)}(x)=\frac{(\alpha+1)_{n}}{n!}\,{}_{1}F_{1}\left(\begin{array}[]{c}-n\\\
\alpha+1\end{array};x\right).$

######  Theorem 3.3.

Suppose that $n+1\geq m$. The type II polynomial in (12) arises as the
generalized hypergeometric series

$P_{n,m}(x)=\frac{(\alpha+1)_{n+m}(\alpha+\nu+1)_{m}}{(n+m)!}\,{}_{2}F_{2}\left(\begin{array}[]{c}-n-m,m+\alpha+\nu+1\\\
\alpha+1,\alpha+\nu+1\end{array};x\right).$ (13)

###### Proof.

Use the previous lemma and the following differentiation formula for the
generalized hypergeometric series

$\frac{d^{n}}{dz^{n}}\left[z^{\gamma}{}_{p}F_{q}\left(\begin{array}[]{c}a_{1},\dots,a_{p}\\\
b_{1},\dots,b_{q}\end{array};z\right)\right]=(\gamma-n+1)_{n}z^{\gamma-n}{}_{p+1}F_{q+1}\left(\begin{array}[]{c}\gamma+1,a_{1},\dots,a_{p}\\\
\gamma+1-n,b_{1},\dots,b_{q}\end{array};z\right),$

see [8, Eq. 16.3.2]. ∎

Note that this hypergeometric series is similar to the hypergeometric series
that is used in the remaining fundamental solution of the differential
equation in (the proof of) Proposition 2.6.

The previous result allows us to determine a differential equation for
$P_{n,m}$.

######  Proposition 3.4.

Every type II polynomial $y=P_{n,m}(x)$ satisfies the third order differential
equation

$x^{2}y^{\prime\prime\prime}+x(2\alpha+\nu+3-x)y^{\prime\prime}+[(\alpha+1)(\alpha+\nu+1)+x(n-\alpha-\nu-2)]y^{\prime}+(n+m)(m+\alpha+\nu+1)y=0.$
(14)

###### Proof.

This is the differential equation associated to the ${}_{2}F_{2}$
hypergeometric series in (13). In general, the hypergeometric series

${}_{2}F_{2}\left(\begin{array}[]{c}a_{1},a_{2}\\\
b_{1},b_{2}\end{array};x\right)$

satisfies the differential equation

$x^{2}y^{\prime\prime\prime}+x(b_{1}+b_{2}+1-x)y^{\prime\prime}+(b_{1}b_{2}-x(a_{1}+a_{2}+1))y^{\prime}-a_{1}a_{2}y=0,$

see [8, Eq. 16.8.3], [28, §16.8 (ii)] or use (8) in a negative argument. ∎

###  3.3 Recurrence relations

The nearest neighbor recurrence relations for the monic type II polynomials
are similar to those for the normalized type I functions (9)–(10); one uses
the same coefficients up to a shift

$\displaystyle xP_{n,m}^{(II)}(x)$ $\displaystyle=$ $\displaystyle
P_{n+1,m}^{(II)}(x)+c_{n,m}P_{n,m}^{(II)}(x)+a_{n,m}P_{n-1,m}^{(II)}(x)+b_{n,m}P_{n,m-1}^{(II)}(x),$
(15) $\displaystyle xP_{n,m}^{(II)}(x)$ $\displaystyle=$ $\displaystyle
P_{n,m+1}^{(II)}(x)+d_{n,m}P_{n,m}^{(II)}(x)+a_{n,m}P_{n-1,m}^{(II)}(x)+b_{n,m}P_{n,m-1}^{(II)}(x).$
(16)

Below we will describe an explicit method to derive these recurrence relations
from the three-term recurrence relation for the underlying Laguerre
polynomials. We again require a connection formula.

######  Lemma 3.5.

One has

$P_{n,m+1}^{(II)}(x)-P_{n+1,m}^{(II)}(x)=\frac{(n+m+\alpha+1)(n+m+\alpha+\nu+1)}{(n+2m+\alpha+\nu+1)(n+2m+\alpha+\nu+2)}P_{n,m}^{(II)}(x).$

###### Proof.

Consider the hypergeometric series representations of the monic type II
polynomials in (13) and compare the coefficients of each monomial $x^{k}$ on
both sides of the equation. If we denote the coefficient of $x^{k}$ in a
polynomial $P(x)$ by $P[k]$, we have

$P_{n,m+1}^{(II)}[k]=\frac{(\alpha+n+m+1)(\alpha+\nu+n+m+1)(\alpha+\nu+m+k+1)}{(\alpha+\nu+n+2m+1)(\alpha+\nu+n+2m+2)(-n-m+k-1)}P_{n,m}^{(II)}[k],$
$P_{n+1,m}^{(II)}[k]=\frac{(\alpha+n+m+1)(\alpha+\nu+n+m+1)}{(\alpha+\nu+n+2m+1)(-n-m+k-1)}P_{n,m}^{(II)}[k].$

Now use

$\frac{\alpha+\nu+m+k+1}{\alpha+\nu+n+2m+2}-1=\frac{-n-m+k-1}{\alpha+\nu+n+2m+2}$

to get the stated result. ∎

Alternative proof of (15)–(16). We may compute one of the derivatives in the
differentiation formula (12) via

$P_{n,m}(x)=x^{-\alpha-\nu}\frac{d^{m}}{dx^{m}}\left[x^{m+\alpha+\nu-1}\cdot
xL_{n+m}^{(\alpha)}(x)\right],$

so that after an application of the Leibniz rule, we get

$P_{n,m}(x)=m!\sum_{k=0}^{m}\binom{m+\alpha+\nu-1}{m-k}\frac{x^{k}}{k!}\frac{d^{k}}{dx^{k}}\left[xL_{n+m}^{(\alpha)}(x)\right].$

The three-term recurrence relation (11) for Laguerre polynomials and the
identity

$\binom{m+\alpha+\nu-1}{m-k}=\binom{m+\alpha+\nu}{m-k}-\binom{m-1+\alpha+\nu}{m-1-k}$

then lead to the recurrence relation

$\displaystyle P_{n,m}(x)=$
$\displaystyle(n+m+1)[mP_{n+2,m-1}(x)-P_{n+1,m}(x)]$
$\displaystyle-(2n+2m+\alpha+1)[mP_{n+1,m-1}(x)-P_{n,m}(x)]$
$\displaystyle+(n+m+\alpha)[mP_{n,m-1}(x)P_{n-1,m}(x)].$

For the monic type II polynomials, this becomes

$\displaystyle P_{n,m}^{(II)}(x)=$
$\displaystyle\left[\frac{m\mathcal{P}_{n,m}}{\mathcal{P}_{n+2,m-1}}P_{n+2,m-1}^{(II)}(x)-\frac{\mathcal{P}_{n,m}}{\mathcal{P}_{n+1,m}}P_{n+1,m}^{(II)}(x)\right]$
$\displaystyle-(2n+2m+\alpha+1)\left[\frac{m\mathcal{P}_{n,m}}{\mathcal{P}_{n+1,m-1}}P_{n+1,m-1}^{(II)}(x)-P_{n,m}^{(II)}(x)\right]$
$\displaystyle+(n+m)(n+m+\alpha)\left[\frac{m\mathcal{P}_{n,m}}{\mathcal{P}_{n,m-1}}P_{n,m-1}^{(II)}(x)-\frac{\mathcal{P}_{n,m}}{\mathcal{P}_{n-1,m}}P_{n-1,m}^{(II)}(x)\right].$

We may use the connection formula (Lemma 3.5) to convert this to the nearest
neighbor recurrence relations. ∎

###  3.4 Asymptotic behavior of the zeros

In what follows, we will investigate the asymptotic behavior of the zeros of
the type II polynomial $P_{\vec{n}}=P_{n_{1},n_{2}}$, with multi-index
$\vec{n}=(n_{1},n_{2})$, as $n=n_{1}+n_{2}\to\infty$. We will now use a multi-
index $(n_{1},n_{2})$ instead of $(n,m)$ to avoid confusion with the degree
$n=n_{1}+n_{2}$ which tends to infinity. The limit will be taken along multi-
indices $\vec{n}$ with $(n_{1}/n,n_{2}/n)\to(q_{1},q_{2})$ as $n\to\infty$. We
necessarily have $0\leq q_{1},q_{2}\leq 1$ and $q_{1}+q_{2}=1$. The
restriction that $n_{1}+1\geq n_{2}$ leads to the condition $q_{1}\geq q_{2}$.

Afterwards, we will compare the results to what is known for the monic type II
multiple Laguerre polynomials of the first kind
$L_{\vec{n}}^{(II\mid\alpha,\beta)}(x)$. These polynomials satisfy the
following orthogonality conditions

$\displaystyle\int_{0}^{\infty}L_{\vec{n}}^{(II\mid\alpha,\beta)}(x)e^{-x}x^{k+\alpha}dx=0,\quad
k=0,\dots,n_{1}-1,$
$\displaystyle\int_{0}^{\infty}L_{\vec{n}}^{(II\mid\alpha,\beta)}(x)e^{-x}x^{k+\beta}dx=0,\quad
k=0,\dots,n_{2}-1.$

####  3.4.1 Smallest zeros

We will focus first on the behavior of the smallest zeros of $P_{\vec{n}}$,
which we may derive from its Mehler-Heine asymptotics. The Mehler-Heine
asymptotics of several families of type II multiple orthogonal polynomials,
including the multiple Laguerre polynomials of the first kind, were determined
in [34]. We will employ the same strategy used there.

######  Proposition 3.6.

For the type II polynomial (13), one has

$\lim_{n\to\infty}\frac{n!}{(\alpha+1)_{n}(\alpha+\nu+1)_{n_{2}}}P_{\vec{n}}(x/n^{2})={}_{0}F_{2}\left(\begin{array}[]{c}-\\\
\alpha+1,\alpha+\nu+1\end{array};-q_{2}x\right).$

As a consequence, if we denote the $k$-th zero of $P_{\vec{n}}(x)$ and of this
hypergeometric series by $p_{k,\vec{n}}$ and $x_{k}$ respectively, we have
$\lim_{n\to\infty}n^{2}p_{k,\vec{n}}=x_{k}$.

###### Proof.

The first result follows immediately from the explicit formula (13) and

$\lim_{n\to\infty}\frac{(-n)_{k}}{n^{k}}=(-1)^{k},\quad\lim_{n\to\infty}\frac{(n_{2}+\alpha+1)_{k}}{n^{k}}=q_{2}^{k}.$

In order to obtain the second one, we can then apply Hurwitz’ theorem. ∎

The associated result for the type II multiple Laguerre polynomials of the
first kind is

$\lim_{n\to\infty}c_{\vec{n}}^{(\alpha,\beta)}L_{\vec{n}}^{(II\mid\alpha,\beta)}(x/n^{2})={}_{0}F_{2}\left(\begin{array}[]{c}-\\\
\alpha+1,\beta+1\end{array};-q_{1}q_{2}x\right),$

see [34, Theorem 3]. Thus, if we denote the $k$-th zero of
$L_{\vec{n}}^{(II\mid\alpha,\beta)}$ by $l_{\vec{n},k}^{(\alpha,\beta)}$ and
apply Hurwitz’ theorem again, we can conclude that

$\lim_{n\to\infty}\frac{p_{\vec{n},k}}{l_{\vec{n},k}^{(\alpha,\alpha+\nu)}}=q_{1}.$

This shows that the smallest zeros of $P_{\vec{n}}$ will be closer to $0$ than
the smallest zeros of $L_{\vec{n}}^{(II\mid\alpha,\alpha+\nu)}$ and that this
effect strengthens for smaller values of $q_{1}$, i.e., if we consider more
orthogonality conditions with respect to the exponential integral weight.

####  3.4.2 Asymptotic zero distribution and largest zero

We will now determine the asymptotic zero distribution of the polynomials
$P_{\vec{n}}$ as $n\to\infty$. In [6], one was able to determine the
asymptotic zero distribution for several families of type II multiple
orthogonal polynomials, including the multiple Laguerre polynomials of the
first kind, for multi-indices on the step-line by investigating an underlying
four-term recurrence relation. Essential to their approach was that the
coefficients in this recurrence relations satisfied a particular asymptotic
property. After using Lemma 3.5 to write the nearest neighbor recurrence
relations (15)–(16) for $P_{\vec{n}}$ as a four-term recurrence relation on
the step-line, it can be verified that the coefficients do not satisfy the
desired property. The asymptotic zero distribution of the polynomials
$P_{\vec{n}}$ therefore doesn’t fit in the class of distributions found in
this article and we will require a more hands on approach to determine it
(which was necessary anyway if we don’t want to restrict to multi-indices on
the step-line).

Since the pair $(w_{1},w_{2})$ forms a Nikishin system, and thus an AT-system,
of weights supported on $[0,\infty)$ and $(-\infty,0]$, $P_{\vec{n}}$ will
have $n$ real simple zeros in $[0,\infty)$ (see, e.g., [27, Chapter 4, §4]).
This interval is unbounded, so we can expect that the (largest) zeros of
$P_{\vec{n}}$ will flow away to $\infty$. In order to proceed, we need to find
an appropriate scaling $\tilde{P}_{\vec{n}}(x)=P_{\vec{n}}(n^{\gamma}x)$ such
that the zeros are in a compact interval $[0,c]$ with $c=c_{q_{1},q_{2}}>0$. A
suitable value of $\gamma$ is indicated by the coefficients of the nearest
neighbor recurrence relations (15)–(16); it is the (smallest) value of
$\gamma$ for which the recurrence relations for the monic scaled polynomial
has bounded coefficients. As $n\to\infty$, we have

$a_{\vec{n}}=\mathcal{O}(n^{2}),\quad b_{\vec{n}}=\mathcal{O}(n^{2}),\quad
c_{\vec{n}}=\mathcal{O}(n),\quad d_{\vec{n}}=\mathcal{O}(n),$

which suggests taking $\gamma=1$.

The right scaling can be derived more rigorously via the following lemma. It
also shows that we can assume that $c\leq 4$. Additionally, the proof gives an
explicit method to prove the fact that $P_{\vec{n}}$ has $n$ simple real zeros
in $[0,\infty)$.

######  Lemma 3.7.

For the largest zero $p_{n,\vec{n}}$ of $P_{\vec{n}}$, we have
$\limsup_{n\to\infty}x_{n,\vec{n}}/n\leq 4$.

###### Proof.

We will show this by comparing the largest zero of $P_{\vec{n}}$ with the
largest zero of the underlying Laguerre polynomial $L_{n}^{(\alpha)}$ in the
differentiation formula (12). The idea will be that the operator

$P(x)\mapsto
x^{-\alpha-\nu}\frac{d^{m}}{dx^{m}}\left[x^{m+\alpha+\nu}P(x)\right]$

will pull the zeros of $L_{n}^{(\alpha)}(x)$ towards $0$. It is known from the
general theory of orthogonal polynomials that $L_{n}^{(\alpha)}(x)$ has $n$
real simple zeros in $[0,\infty)$. If we apply the operator

$P(x)\mapsto x^{-\alpha-\nu}\frac{d}{dx}\left[x^{\alpha+\nu}P(x)\right]$

once to $x^{m}L_{n}^{(\alpha)}(x)$, we get a polynomial with a zero of order
$m-1$ at $0$ and of which the other $n$ zeros are real, simple, in
$[0,\infty)$ and interlace with the zeros of the input polynomial
$L_{n}^{(\alpha)}(x)$ and $0$ (due to Rolle’s theorem). In particular, all of
these zeros are smaller than the largest zero of $L_{n}^{(\alpha)}(x)$.
Applying the operator $m-1$ additional times and repeating this argument, we
recover $P_{\vec{n}}$ together with the conclusion that it has $n$ simple real
zeros in $[0,\infty)$ all smaller than the largest zero of
$L_{n}^{(\alpha)}(x)$. It is known that the largest zero of the Laguerre
polynomial $L_{n}^{(\alpha)}(x)$ behaves as $4n$ as $n\to\infty$ (see, e.g.,
[33, Thm. 6.31.2]). Hence we get the estimate in the lemma. ∎

In what follows, we will use some of the techniques from [23] to derive the
asymptotic zero distribution of $\tilde{P}_{\vec{n}}(x)=P_{\vec{n}}(nx)$.
Denote the zeros of $\tilde{P}_{\vec{n}}(x)$ by
$x_{1,\vec{n}}<\dots<x_{n,\vec{n}}$. Their distribution is given by the
probability measure

$\mu_{\vec{n}}=\frac{1}{n}\sum_{j=1}^{n}\delta_{x_{j,\vec{n}}}$

with Stieltjes transform

$S_{\vec{n}}(z)=\int_{0}^{c}\frac{d\mu_{\vec{n}}(x)}{z-x}=\frac{1}{n}\frac{\tilde{P}^{\prime}_{\vec{n}}(z)}{\tilde{P}_{\vec{n}}(z)}.$

According to Helly’s selection principle, the sequence of probability measures
$(\mu_{\vec{n}})_{n\in\mathbb{N}}$ on $[0,c]$ has a subsequence
$(\mu_{\vec{n}_{k}})_{k\in\mathbb{N}}$ that converges weakly to a probability
measure $\mu=\mu_{q_{1},q_{2}}$ on $[0,c]$. If this limit is independent of
the subsequence, it is called the asymptotic zero distribution of
$\tilde{P}_{\vec{n}}(x)$.

Let $(\mu_{\vec{n}_{k}})_{k\in\mathbb{N}}$ be any subsequence of
$(\mu_{\vec{n}})_{n\in\mathbb{N}}$ that converges weakly to $\mu$. Then it
follows from the Grommer-Hamburger theorem (see, e.g., [14]) that
$(S_{\vec{n}_{k}}(z))_{k\in\mathbb{N}}$ converges uniformly to the Stieltjes
transform

$S(z)=\int_{0}^{c}\frac{d\mu(x)}{z-x}$

on compact subsets of $\mathbb{C}\backslash[0,c]$. Moreover, $zS(z)\to 1$ as
$z\to\infty$.

We can proceed via the (third order) differential equation for
$P_{\vec{n}}(x)$; we will use it to show that $S(z)$ satisfies a third order
algebraic equation. Since this equation will be independent of the selected
subsequence, we will be able to conclude that $S(z)$ is the Stieltjes
transform of the asymptotic zero distribution of $\tilde{P}_{\vec{n}}(x)$.
Moreover, by solving the equation and by looking for the solution(s) with
$zS(z)\to 1$ as $z\to\infty$, we will be able to obtain an explicit expression
for it.

######  Lemma 3.8.

The Stieltjes transform $S(z)$ of $\mu$ satisfies the cubic equation

$z^{2}S^{3}-z^{2}S^{2}+q_{1}zS+q_{2}=0,\quad z\in\mathbb{C}\backslash[0,c].$

###### Proof.

We will start by transforming the differential equation (14) for
$P_{\vec{n}}(z)$ into a differential equation for $\tilde{P}_{\vec{n}}(z)$.
Only the asymptotic structure as $n\to\infty$ will be important. By making use
of the relation $\tilde{P}^{(k)}_{\vec{n}}(z)=n^{k}P^{(k)}_{\vec{n}}(nz)$, we
find that

$z^{2}\tilde{P}^{\prime\prime\prime}_{\vec{n}}(z)/n+z(-z+\mathcal{O}(1))\tilde{P}^{\prime\prime}_{\vec{n}}(z)+[z(n_{1}+\mathcal{O}(1))+\mathcal{O}(1)]\tilde{P}^{\prime}_{\vec{n}}(z)+n(n_{2}+\mathcal{O}(1))P_{\vec{n}}(z)=0,$

as $n\to\infty$. From the identity

$\tilde{P}^{\prime}_{\vec{n}}(z)=n\tilde{P}_{\vec{n}}(z)S_{\vec{n}}(z)$

then follows that

$\tilde{P}^{\prime\prime}_{\vec{n}}(z)=n^{2}\tilde{P}_{\vec{n}}(z)(S_{\vec{n}}(z)^{2}+n^{-1}S_{\vec{n}}^{\prime}(z))$

and similarly that

$\tilde{P}^{\prime\prime\prime}_{\vec{n}}(z)=n^{3}\tilde{P}_{\vec{n}}(z)(S_{\vec{n}}(z)^{3}+3n^{-1}S_{\vec{n}}(z)S_{\vec{n}}^{\prime}(z)+n^{-2}S_{\vec{n}}^{\prime\prime}(z)).$

By passing to the given subsequence $(n_{k})_{k\in\mathbb{N}}$, we may assume
that $S_{\vec{n}}$ converges uniformly to $S$ on compact subsets of
$\mathbb{C}\backslash[0,c]$ as $n\to\infty$. In that case, the same is true
for all the derivatives. We may therefore write

$\tilde{P}^{(k)}_{\vec{n}}(z)=n^{k}\tilde{P}_{\vec{n}}(z)(S_{\vec{n}}^{k}(z)+\mathcal{O}(n^{-1})),\quad
n\to\infty,$

which holds uniformly on compact subsets of $\mathbb{C}\backslash[0,c]$. If we
plug this in the above equation for $\tilde{P}_{\vec{n}}(z)$, divide by
$n^{2}\tilde{P}_{\vec{n}}(z)$ and take the limit $n\to\infty$, we obtain the
desired algebraic equation for $S(z)$. ∎

Observe that in the Laguerre case $(q_{1},q_{2})=(1,0)$, the cubic equation
corresponds to the quadratic equation

$zS^{2}-zS+1=0,$

which is, as expected, the algebraic equation for the Stieltjes transform of
the asymptotic zero distribution of the scaled Laguerre polynomials
$L_{n}^{(\alpha)}(nx)$ (the Marchenko-Pastur distribution).

######  Theorem 3.9.

The density of the asymptotic zero distribution of $P_{\vec{n}}(nx)$ as
$n\to\infty$ is given by

$u(x;q_{1},q_{2})=\frac{1}{2^{4/3}3^{1/2}\pi}\frac{\sqrt[3]{r_{+}(x;q_{1},q_{2})}-\sqrt[3]{r_{-}(x;q_{1},q_{2})}}{\sqrt[3]{x^{2}}},\quad
x\in[0,c],$

where

$r_{\pm}(x;q_{1},q_{2})=2x^{2}-9q_{1}x-27q_{2}\pm
3\sqrt{3}\sqrt{-(1+q_{2})^{2}x^{2}+(4q_{1}^{3}+18q_{1}q_{2})x+27q_{2}^{2}}$

and

$c_{q_{1},q_{2}}=\frac{(2q_{1}^{3}+9q_{1}q_{2})+2(q_{2}^{2}+q_{2}+1)^{3/2}}{(1+q_{2})^{2}}.$

###### Proof.

We can solve the cubic equation (3.8) explicitly by making use of Cardano’s
formula. The three solutions are given by

$S_{k}(z)=\frac{1}{3}+\zeta^{k}T_{+}(z)^{1/3}+\bar{\zeta}^{k}T_{-}(z)^{1/3},\quad
k=0,1,2,$

in terms of the third root of unity $\zeta=e^{2\pi i/3}$ and the functions

$T_{\pm}(z)=R(z)\pm(R(z)^{2}+Q(z)^{3})^{1/2},\quad
R(z)=\frac{-9q_{1}/z-27q_{2}/z^{2}+2}{54},\quad Q(z)=\frac{3q_{1}/z-1}{9}.$

The square and cubic root may be taken along any branch cut; swapping to
another branch of the cubic root would just permute the solutions $S_{k}(z)$,
while changing the branch of the square root would simply interchange the
roles of $T_{\pm}(z)$. We are only interested in the solution(s) with
$zS(z)\to 1$ for $z\to\infty$. By making use of a computer algebra system, we
can find

$zS_{0}(z)\to\infty,\quad zS_{1}(z)\to\frac{q_{1}-q_{2}-1}{2},\quad
zS_{2}(z)\to\frac{q_{1}+q_{2}+1}{2}=1.$

Consequently, $S_{2}(z)$ is the Stieltjes transform of the probability measure
that we are looking for. We can recover the underlying density $u(x)$ by
making use of the Stieltjes-Perron inversion formula (and
$\zeta^{2}=-(1+\sqrt{3}i)/2$)

$u(x)=-\lim_{y\to
0^{+}}\frac{1}{\pi}\text{Im}(S_{2}(x+iy))=\frac{\sqrt{3}}{6\pi}\left(\frac{r_{+}(x)}{2x^{2}}\right)^{1/3}-\frac{\sqrt{3}}{6\pi}\left(\frac{r_{-}(x)}{2x^{2}}\right)^{1/3}.$

Note that the complex roots have become real roots. For the cubic root, this
means that $\sqrt[3]{x}=\text{sgn}(x)\sqrt[3]{|x|}$. It is justified to take
the real square root, because the function in the argument is positive for
$x\in[0,c]$. The constant $c$ corresponds to the largest zero of $u(x)$. The
zeros of $u(x)$ are determined by the equation

$-(1+q_{2})^{2}x^{2}+(4q_{1}^{3}+18q_{1}q_{2})x+27q_{2}^{2}=0,$

which has solutions

$x_{\pm}=\frac{(2q_{1}^{3}+9q_{1}q_{2})\pm
2(q_{2}^{2}+q_{2}+1)^{3/2}}{(1+q_{2})^{2}}.$

We have $x_{-}<x_{+}$ and thus $c=x_{+}$. ∎

In the Laguerre case $(q_{1},q_{2})=(1,0)$, the above expression becomes

$u(x;1,0)=\frac{1}{2\pi}\frac{\sqrt{4-x}}{\sqrt{x}},\quad x\in[0,4],$

which is the expected Marchenko-Pastur distribution. By making use of the
basic identity $a-b=(a^{3}-b^{3})/(a^{2}+ab+b^{2})$, the connection between
the two distributions can be made more clear:

$u(x;q_{1},q_{2})=\frac{3}{2^{1/3}\pi}\frac{\sqrt{-(1+q_{2})^{2}x^{2}+(4q_{1}^{3}+18q_{1}q_{2})x+27q_{2}^{2}}}{\sqrt[3]{x^{2}}(\sqrt[3]{r_{+}(x;q_{1},q_{2})^{2}}+\sqrt[3]{r_{-}(x;q_{1},q_{2})^{2}}+\sqrt[3]{4x}(x-3q_{1}))}.$

In Fig. 1 we plotted the density $u(x;q_{1},q_{2})$ for several values of
$q_{1}\in[\frac{1}{2},1]$.

Figure 1: The density $u(x;q_{1},q_{2})$ for $q_{1}=1$ (solid),
$q_{1}=\frac{3}{4}$ (dash) and $q_{1}=\frac{1}{2}$ (dots)

The figure shows that as $q_{1}$ decreases, or equivalently $q_{2}$ increases,
the zeros of the polynomials $P_{\vec{n}}$ are pulled closer towards $0$.
Typically for Nikishin systems of two weights, the zeros of the type II
polynomials are contained in the support of the first weight and are attracted
to the support of the second weight. By increasing the contribution of the
second weight, which had support $(-\infty,0]$, we thus observe that this
effect becomes stronger.

The behavior of the density $u(x)$ near its endpoints is

$u(x)\sim\frac{\sqrt{3}}{2\pi}q_{2}^{1/3}x^{-2/3},\quad x\to 0,$
$u(x)=\mathcal{O}((c-x)^{1/2}),\quad x\to c.$

Near the right endpoint, we have the usual square root behavior that also
occurs for the asymptotic zero distribution of the Laguerre polynomials and
the multiple Laguerre polynomials of the first kind on the step-line. Near the
left endpoint, we have the same behavior as the latter but with a different
constant.

We end this section by noting that the asymptotic zero distribution $\nu$ of
the scaled type I polynomial $\tilde{B}_{\vec{n}}(x)=B_{\vec{n}}(-nx)$ is the
same as the asymptotic zero distribution $\mu$ of $\tilde{P}_{\vec{n}}(x)$, so
that the asymptotic zero distribution of $B_{\vec{n}}(nx)$ has support
$[-c,0]$ and is the reflection of $\mu$ over $0$. This follows from the fact
that the Stieltjes transform $S(z;\nu)$ satisfies the same algebraic equation
as $S(z;\mu)$ because the coefficients in the hypergeometric differential
equations for $\tilde{B}_{\vec{n}}(x)$ and $\tilde{P}_{\vec{n}}(x)$ are
asymptotically equivalent. It also explains the solution $S_{1}(z)$ in the
proof of Theorem 3.9, which had a finite limit $zS(z)$ as $z\to\infty$ as
well; the change of variables $z\mapsto-z$ corresponds to letting
$q_{1}\mapsto-q_{1}$ so that now $-zS_{1}(-z)\to 1$ as $z\to\infty$.

###  3.5 Connection to random matrix theory

The multiple orthogonal polynomials that were investigated in the previous
sections appear naturally in the study of the squared singular values of the
product of a truncated Haar distributed random unitary matrix with a Ginibre
matrix. A (complex) Ginibre matrix is a random matrix whose entries are
independent standard complex Gaussian random variables, see [15]. There is a
similar connection between the multiple orthogonal polynomials associated with
the $K$-Bessel functions in [37] and the product of two Ginibre matrices as
demonstrated in [22]. In the latter, it is also shown how products of more
than two Ginibre matrices are related to multiple orthogonal polynomials
associated with certain Meijer-$G$ functions. Products of several Ginibre
matrices with a truncated unitary matrix and products that only consist of
truncated unitary matrices were investigated in [21] and [20] respectively. A
potential connection to multiple orthogonal polynomials however was not
mentioned here.

It is known that the joint probability density of the squared singular values
of a size $(n+a)\times n$ upper left corner of a random unitary matrix of size
$(n+a+b)\times(n+a+b)$ (assuming that $b\geq n$) is given by a polynomial
ensemble with probability density

$\propto\Delta_{n}(x)\det\left[x_{j}^{a+k-1}(1-x_{j})^{b-n}\right]_{j,k=1}^{n},\quad\text{
all }x_{j}\in[0,1].$

Here $\Delta_{n}(x)$ denotes the Vandermonde determinant
$\Delta_{n}(x)=\prod_{1\leq i<j\leq n}(x_{j}-x_{i})$. On the other hand, the
joint probability density of the squared singular values of a size
$(n+a)\times n$ Ginibre matrix is determined by a polynomial ensemble with
probability density

$\propto\Delta_{n}(x)\det\left[x_{j}^{a+k-1}e^{-x_{j}}\right]_{j,k=1}^{n},\quad\text{
all }x_{j}>0.$

Observe that in these ensembles a beta and gamma density is used respectively.

Suppose that $n=n_{1}+n_{2}$ and consider the product $TG$ of a size
$(n+a+b)\times(n+a)$ matrix $T$ that arises by truncating a size
$(n+n_{2}+a+b)\times(n+n_{2}+a+b)$ random unitary matrix with a Ginibre matrix
$G$ of size $(n+a)\times n$. It was proven in [20] that the joint probability
density of the squared singular values of $TG$ is then given by a polynomial
ensemble with probability density

$\propto\Delta_{n}(x)\det\left[w_{k-1}(x_{j})\right]_{j,k=1}^{n},\quad\text{
all }x_{j}>0,$ (17)

in terms of weights $w_{k}$ that arise as the Mellin convolution of a beta and
gamma density

$w_{k}(x)=\int_{0}^{1}t^{a+b}(1-t)^{n_{2}-1}(x/t)^{a+k}e^{-x/t}\frac{dt}{t},\quad
x>0.$

These weights are connected to the exponential integral in the following way

$w_{k}(x)=x^{a+k}\sum_{j=0}^{n_{2}-1}(-1)^{j}\binom{n_{2}-1}{k}E_{b+j-k}(x)$

(use the binomial theorem and the change of variables $t\mapsto 1/t$).

Typically one is interested in the correlation kernel

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

of such polynomial ensembles. The kernel may be built out of monic polynomials
$P_{k}$ of degree $k$ and functions $Q_{l}$ in the linear span of
$w_{0},\dots,w_{l}$ 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.$

These are satisfied under the following conditions

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

It is not guaranteed that such a construction is possible, however if it
exists, 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_{0},\dots,w_{n-1}\\}$ respectively. Particular interest goes
to the monic polynomial $P_{n}$ of degree $n$ that extends the above
biorthogonality to $k=n$, because it turns out to be the average
characteristic polynomial $\mathbb{E}[\prod_{j=1}^{n}(x-x_{j})]$ over the
particles in the associated polynomial ensemble.

The result below shows the connection between the functions that satisfy
(18)–(19) and the multiple orthogonal polynomials associated with the
exponential integral that were studied in the previous sections.

######  Proposition 3.10.

The functions $Q_{l}$ that satisfy (19) are related to the type I functions in
(6) in the following way

$Q_{l}(x)\propto x^{a}F_{l-n_{2}+1,n_{2}}^{a,b}(x),\quad 2n_{2}-2\leq l\leq
n-1.$

The relation between the polynomials $P_{k}$ that satisfy (18) and the type II
polynomials in (12) is

$P_{k}(x)\propto P_{k-n_{2},n_{2}}^{a,b}(x),\quad 2n_{2}-1\leq k\leq n.$

###### Proof.

The weight $w_{j}$ has Mellin transform

$\int_{0}^{\infty}w_{j}(x)x^{s-1}dx=(n_{2}-1)!\frac{\Gamma(s+a+j)}{(s+a+b)_{n_{2}}}.$

Using the same argument as in the proof of Lemma 2.1, it can be shown that
$Q_{l}$ has a Mellin transform of the form

$\int_{0}^{\infty}Q_{l}(x)x^{s-1}dx\propto\frac{\Gamma(s+a)}{(s+a+b)_{n_{2}}}(1-s)_{l}.$

The first result then follows by comparing this to the Mellin transform (3) of
the type I function $F_{l-n_{2}+1,n_{2}}^{a,b}$. Similarly as in [21, §4.2],
we may derive an integral representation for $P_{k}$ of the form

$P_{k}(x)\propto\int_{\Sigma_{k}}\frac{1}{(-t)_{k+1}}\frac{(t+a+b+1)_{n_{2}}}{\Gamma(t+a+1)}x^{t}dt,$

where $\Sigma_{k}$ is a closed contour in the domain $\text{Re}(t)>-1$ that
encircles the interval $[0,k]$ once in the positive direction. We may apply
the residue theorem in a similar fashion as in the proof of Proposition 2.5 to
evaluate this contour integral as a ${}_{2}F_{2}$ hypergeometric series. After
comparing with the ${}_{2}F_{2}$ hypergeometric series (13) for
$P_{k-n_{2},n_{2}}^{a,b}(x)$, we obtain the second result. ∎

######  Corollary 3.11.

The average characteristic polynomial $P_{n}$ over the particles in the
polynomial ensemble (17) is given by the monic type II polynomial
$P_{n_{1},n_{2}}^{a,b}$.

##  4 Mixed type multiple orthogonal polynomials

We will investigate two families of mixed type multiple orthogonal polynomials
that involve pairs of weights $(e^{-x},E_{\nu+1}(x))$ and
$(x^{\alpha},x^{\beta})$; the orthogonality conditions will be a mixture of
type I and type II orthogonality conditions. A similar consideration also
appears in Apéry’s construction of rational approximants for $\zeta(3)$, see
Beukers [5]. Mixed type multiple orthogonal polynomials also appear in the
work of Sorokin [31] and in the study of non-intersecting Brownian motions
[7]. We will assume that the parameters satisfy the conditions
$\alpha,\beta,\nu,\alpha+\nu,\beta+\nu>-1$ and
$\alpha-\beta\not\in\mathbb{Z}$.

###  4.1 First family

We will first investigate the mixed type multiple orthogonal polynomials with
respect to the two pairs of weights $(e^{-x},E_{\nu+1}(x))$ and
$(x^{\alpha},x^{\beta})$ on $[0,\infty)$. They will depend on two multi-
indices $\vec{n}=(n_{1},n_{2})$ and $\vec{m}=(m_{1},m_{2})$ that will
determine the degrees of the polynomials and the number of orthogonality
conditions respectively. The multi-indices are connected through the relation
$\left|\vec{n}\right|-1=\left|\vec{m}\right|$. We are looking for polynomials
$A_{\vec{n},\vec{m}}$ and $B_{\vec{n},\vec{m}}$ with $\deg
A_{\vec{n},\vec{m}}\leq n_{1}-1$ and $\deg B_{\vec{n},\vec{m}}\leq n_{2}-1$,
for which the function

$F_{\vec{n},\vec{m}}(x)=A_{\vec{n},\vec{m}}(x)e^{-x}+B_{\vec{n},\vec{m}}(x)E_{\nu+1}(x)$
(20)

satisfies the orthogonality conditions

$\displaystyle\int_{0}^{\infty}F_{\vec{n},\vec{m}}(x)x^{k+\alpha}\,dx$
$\displaystyle=0,\quad k=0,\dots,m_{1}-1,$ (21)
$\displaystyle\int_{0}^{\infty}F_{\vec{n},\vec{m}}(x)x^{k+\beta}\,dx$
$\displaystyle=0,\quad k=0,\dots,m_{2}-1.$ (22)

The polynomials $A_{\vec{n},\vec{m}}$ and $B_{\vec{n},\vec{m}}$ have
$\left|\vec{n}\right|$ coefficients and the system of equations (21)–(22)
corresponds to a homogeneous linear system of $\left|\vec{m}\right|$
equations. Due to the condition $\left|\vec{m}\right|=\left|\vec{n}\right|-1$,
existence of a non-zero solution is guaranteed. If the matrix of the system
has maximal rank, the coefficients of $A_{\vec{n},\vec{m}}$ and
$B_{\vec{n},\vec{m}}$ are uniquely determined up to a scalar multiplication.

The (degenerate) case $m_{2}=0$ corresponds to the setting of the type I
polynomials as before. As a consequence, we will encounter some similar
properties. On the other hand, the (degenerate) case $n_{2}=0$ corresponds to
the setting of the type II multiple Laguerre polynomial of the first kind.
This suggests that, instead of a Laguerre polynomial $L_{n}^{(\alpha)}(x)$,
now a type II multiple Laguerre polynomial of the first kind
$L_{\vec{n}}^{(II\mid\alpha,\beta)}(x)$ may play a role. Explicit formulas for
these polynomials are well-known; a Rodrigues formula can for example be found
in [17, §23.4.1]

$L_{\vec{n}}^{(II\mid\alpha,\beta)}(x)=(-1)^{\left|\vec{n}\right|}e^{x}x^{-\alpha}\frac{d^{n_{1}}}{dx^{n_{1}}}\left[x^{n_{1}+\alpha-\beta}\frac{d^{n_{2}}}{dx^{n_{2}}}\left[x^{n_{2}+\beta}e^{-x}\right]\right].$
(23)

In what follows, we will assume that the multi-index $\vec{n}$ satisfies
$n_{1}+1\geq n_{2}$. This is the natural extension of the restriction we had
before; the case $m_{2}=0$ should correspond to the type I setting.

####  4.1.1 Convolution and Rodrigues formula

Our approach will be similar as in Section 2.1.

######  Lemma 4.1.

There is a mixed type function $F_{\vec{n},\vec{m}}$ with Mellin transform

$\hat{F}_{\vec{n},\vec{m}}(s)=\frac{\Gamma(s)}{(s+\nu)_{n_{2}}}(\alpha+1-s)_{m_{1}}(\beta+1-s)_{m_{2}},\quad
s>\max\\{0,-\nu\\}.$ (24)

###### Proof.

The Mellin transform of a mixed type function may be computed using the same
approach as in the proof of Lemma 2.1. In this way, we can obtain

$\hat{F}_{\vec{n},\vec{m}}(s)=\frac{\Gamma(s)}{(s+\nu)_{n_{2}}}\left(p_{n_{1}-1}(s)(s+\nu)_{n_{2}}+q_{2n_{2}-2}(s)\right)$

in terms of polynomials with their subscript as maximal degree. The degree of
the polynomial between the parentheses is less than or equal to
$\max\\{n_{1}+n_{2}-1,2n_{2}-2\\}$, which is equal to
$n_{1}+n_{2}-1=m_{1}+m_{2}$ under the assumption that $n_{1}\geq n_{2}-1$.
Since the orthogonality conditions imply that $\hat{F}_{\vec{n},\vec{m}}(s)$
has zeros at the $m_{1}+m_{2}$ values $s=\alpha+1,\dots,\alpha+m_{1}$ and
$s=\beta+1,\dots,\beta+m_{2}$, the remaining polynomial must be a non-zero
scalar multiple of $(\alpha+1-s)_{m_{1}}(\beta+1-s)_{m_{2}}$. ∎

We will compare this Mellin transform to the Mellin transform of a type II
Laguerre polynomial of the first kind multiplied by a gamma density, which we
will compute below.

######  Lemma 4.2.

The Mellin transform of $L_{\vec{n}}^{(II\mid\alpha,\beta)}(x)e^{-x}$ is given
by

$\int_{0}^{\infty}L_{\vec{n}}^{(II\mid\alpha,\beta)}(x)e^{-x}x^{s}dx=(-1)^{\left|\vec{n}\right|}\Gamma(s)(\alpha+1-s)_{n_{1}}(\beta+1-s)_{n_{2}},\quad
s>0.$

###### Proof.

We may use the Rodrigues formula (23) and perform integration by parts $n_{1}$
times to get

$\displaystyle\int_{0}^{\infty}L_{\vec{n}}^{(II\mid\alpha,\beta)}(x)e^{-x}x^{s-1}dx$
$\displaystyle=(-1)^{\left|\vec{n}\right|+n_{1}}\int_{0}^{\infty}\left(x^{n_{1}+\alpha-\beta}\frac{d^{n_{2}}}{dx^{n_{2}}}x^{n_{2}+\beta}e^{-x}\right)\left(\frac{d^{n_{1}}}{dx^{n_{1}}}x^{s-\alpha-1}\right)dx$
$\displaystyle=(-1)^{\left|\vec{n}\right|}(\alpha+1-s)_{n_{1}}\int_{0}^{\infty}\left(\frac{d^{n_{2}}}{dx^{n_{2}}}x^{n_{2}+\beta}e^{-x}\right)x^{-\beta+s-1}dx.$

The remaining integral is the Mellin transform of
$n_{2}!L_{n_{2}}^{(\beta)}(x)e^{-x}$ and was computed in Lemma 2.2. ∎

The previous result allows us to write the mixed type function as a Mellin
convolution in which a type II Laguerre polynomial of the first kind appears.

######  Theorem 4.3.

The mixed type function with Mellin transform (24) is given by

$F_{\vec{n},\vec{m}}(x)=\int_{x}^{\infty}L_{\vec{m}}^{(II\mid\alpha,\beta)}(t)e^{-t}\left(1-\frac{x}{t}\right)^{n_{2}-1}\left(\frac{x}{t}\right)^{\nu}\,\frac{dt}{t}$
(25)

multiplied by $(-1)^{\left|\vec{n}\right|}/(n_{2}-1)!$.

###### Proof.

The Mellin transform (24) is $(-1)^{\left|\vec{n}\right|}/(n_{2}-1)!$ times
the product of the Mellin transform of a type II Laguerre polynomial
multiplied with a gamma density

$\int_{0}^{\infty}L_{\vec{m}}^{(II\mid\alpha,\beta)}(x)e^{-x}x^{s-1}\,dx=(-1)^{\left|\vec{n}\right|}\Gamma(s)(\alpha+1-s)_{m_{1}}(\beta+1-s)_{m_{2}},$

(see Lemma 4.2) and of a beta density

$\int_{0}^{\infty}(1-x)^{n_{2}-1}x^{s+\nu-1}\chi_{[0,1]}(x)\,dx=B(s+\nu,n_{2})=\frac{(n_{2}-1)!}{(s+\nu)_{n_{2}}}.$

Hence we can use the multiplicative property of the Mellin convolution and the
uniqueness of the Mellin transform to obtain the above. ∎

The structure of the mixed type function in (25) resembles the function that
is used to generate the rational approximants in Apéry’s proof of the
irrationality of $\zeta(3)$. There one uses a Mellin convolution of a
(multiple) Jacobi-Piñeiro polynomial and a beta density, see [29]. Instead of
a Jacobi-Piñeiro polynomial, we use a multiple Laguerre polynomial multiplied
with a gamma density.

The convolution formula allows us to derive a Rodrigues formula for the mixed
type function.

######  Theorem 4.4.

The mixed type function in (25) is given by the Rodrigues formula

$F_{\vec{n},\vec{m}}(x)=\frac{(-1)^{\left|\vec{n}\right|}}{(\nu+1)_{n_{2}-1}}x^{-\alpha}\frac{d^{m_{1}}}{dx^{m_{1}}}\left[x^{m_{1}+\alpha-\beta}\frac{d^{m_{2}}}{dx^{m_{2}}}\left[x^{m_{2}+\beta}e^{-x}\frac{d^{n_{2}-1}}{dx^{n_{2}-1}}\left[x^{n_{2}-1}e^{x}E_{\nu+1}(x)\right]\right]\right].$

###### Proof.

The proof will be similar to the proof of the Rodrigues formula (7) for the
type I function; now we will use the Rodrigues formula for the underlying
multiple Laguerre polynomial. We again start with the change of variables
$t\mapsto xt$ in order to put the dependence on $x$ in the multiple Laguerre
polynomial

$F_{\vec{n},\vec{m}}(x)=\int_{1}^{\infty}L_{\vec{m}}^{(II\mid\alpha,\beta)}(xt)e^{-xt}\left(1-\frac{1}{t}\right)^{n_{2}-1}\,\frac{dt}{t^{\nu+1}}.$

The Rodrigues formula (5) for $L_{\vec{m}}^{(II\mid\alpha,\beta)}(x)$ and the
chain rule then imply that

$L_{\vec{m}}^{(II\mid\alpha,\beta)}(xt)e^{-xt}=(-1)^{\left|\vec{n}\right|}(xt)^{-\alpha}\frac{d^{m_{1}}}{dx^{m_{1}}}\left[(xt)^{m_{1}+\alpha-\beta}\frac{d^{m_{2}}}{dx^{m_{2}}}\left[(xt)^{m_{2}+\beta}e^{-xt}\right]\cdot
t^{-m_{2}}\right]\cdot t^{-m_{1}},$

or after rearranging some factors depending on $t$,

$L_{\vec{m}}^{(II\mid\alpha,\beta)}(xt)e^{-xt}=(-1)^{\left|\vec{n}\right|}x^{-\alpha}\frac{d^{m_{1}}}{dx^{m_{1}}}\left[x^{m_{1}+\alpha-\beta}\frac{d^{m_{2}}}{dx^{m_{2}}}\left[x^{m_{2}+\beta}e^{-xt}\right]\right].$

We can then take all the derivatives outside of the integral to get

$F_{\vec{n},\vec{m}}(x)=(-1)^{\left|\vec{n}\right|}x^{-\alpha}\frac{d^{m_{1}}}{dx^{m_{1}}}\left[x^{m_{1}+\alpha-\beta}\frac{d^{m_{2}}}{dx^{m_{2}}}\left[x^{m_{2}+\beta}\int_{1}^{\infty}e^{-xt}\left(1-\frac{1}{t}\right)^{n_{2}-1}\,\frac{dt}{t^{\nu+1}}\right]\right].$

The remaining integral can be written in the desired way using the same
argument as in the proof of formula (7). ∎

####  4.1.2 Hypergeometric series

Similarly as with the type I function, we can derive a representation of the
mixed type function in terms of certain generalized hypergeometric series.

######  Proposition 4.5.

Suppose that $\nu\not\in\mathbb{Z}$. The mixed type function with Mellin
transform (24) arises as the following linear combination of ${}_{3}F_{3}$
hypergeometric series

$\displaystyle F_{\vec{n},\vec{m}}(x)=$
$\displaystyle\mathcal{A}_{\vec{n},\vec{m}}\
{}_{3}F_{3}\left(\begin{array}[]{c}m_{1}+\alpha+1,m_{2}+\beta+1,-n_{2}-\nu+1\\\
\alpha+1,\beta+1,-\nu+1\end{array};-x\right)$
$\displaystyle+\mathcal{B}_{\vec{n},\vec{m}}\
{}_{3}F_{3}\left(\begin{array}[]{c}-n_{2}+1,m_{1}+\alpha+\nu+1,m_{2}+\beta+\nu+1\\\
\nu+1,\alpha+\nu+1,\beta+\nu+1\end{array};-x\right)x^{\nu},$

where

$\mathcal{A}_{\vec{n},\vec{m}}=\frac{(\alpha+1)_{m_{1}}(\beta+1)_{m_{2}}}{(\nu)_{n_{2}}},\quad\mathcal{B}_{\vec{n},\vec{m}}=\frac{(\alpha+\nu+1)_{m_{1}}(\beta+\nu+1)_{m_{2}}\Gamma(-\nu)}{(n_{2}-1)!}.$

###### Proof.

We will use the same strategy as in the proof of Theorem 2.5 and take a closer
look at the Mellin transform of $x^{\alpha}F_{\vec{n},\vec{m}}(x)$, which is

$\int_{0}^{\infty}F_{\vec{n},\vec{m}}(x)x^{\alpha+s-1}dx=\frac{\Gamma(s+\alpha)}{(s+\alpha+\nu)_{n_{2}}}(1-s)_{m_{1}}(\beta-\alpha+1-s)_{m_{2}},\quad
s>-\min\\{\alpha,\alpha+\nu\\},$

(starting from $x^{\beta}F_{\vec{n},\vec{m}}(x)$ will lead to the same
result). We can extend this function to a function that is analytic in the
complex plane except for the simple poles $-\alpha,-\alpha-1,\dots$ and
$-\alpha-\nu,\dots,-\alpha-\nu-n_{2}+1$. All these poles lie in an interval
$(-\infty,c)$ for some $c\in\mathbb{R}$ with
$-\min\\{\alpha,\alpha+\nu,0\\}<c<1$. The original function can then be
recovered by evaluating the inverse Mellin transform via the residue theorem

$\displaystyle x^{\alpha}F_{\vec{n},\vec{m}}(x)$
$\displaystyle=\int_{c-i\infty}^{c+\infty}\frac{\Gamma(s+\alpha)}{(s+\alpha+\nu)_{n_{2}}}(1-s)_{m_{1}}(\beta-\alpha+1-s)_{m_{2}}x^{-s}ds$
$\displaystyle=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{k!}\frac{(k+\alpha+1)_{m_{1}}(k+\beta+1)_{m_{2}}}{(-k+\nu)_{n_{2}}}x^{\alpha+k}$
$\displaystyle\quad+\sum_{k=0}^{n_{2}-1}\frac{\Gamma(-k-\nu)}{(-k)_{k}(1)_{n_{2}-k}}(k+\alpha+\nu+1)_{m_{1}}(k+\beta+\nu+1)_{m_{2}}x^{\alpha+\nu+k}.$

After some elementary computations, we obtain the stated result. ∎

The two terms in the above expression are again solutions of the same
hypergeometric differential equation (now of the fourth order). The two other
solutions are

$x^{-\alpha}{}_{3}F_{3}\left(\begin{array}[]{c}-n_{2}-\alpha-\nu+1,m_{1}+1,m_{2}+\beta-\alpha+1\\\
-\alpha-\nu+1,-\alpha+1,\beta-\alpha+1\end{array};-x\right),$
$x^{-\beta}{}_{3}F_{3}\left(\begin{array}[]{c}-n_{2}-\alpha-\nu+1,m_{1}+\alpha-\beta+1,m_{2}+1\\\
-\beta-\nu+1,\alpha-\beta+1,-\beta+1\end{array};-x\right).$

Using the same argument as in Proposition 2.7, we then must have

$B_{\vec{n},\vec{m}}(x)=B_{\vec{n},\vec{m}}(0)\cdot{}_{3}F_{3}\left(\begin{array}[]{c}-n_{2}+1,m_{1}+\alpha+\nu+1,m_{2}+\beta+\nu+1\\\
\nu+1,\alpha+\nu+1,\beta+\nu+1\end{array};-x\right).$

###  4.2 Second family

In this section, we will investigate the mixed type multiple orthogonal
polynomials with respect to the two pairs of weights $(x^{\alpha},x^{\beta})$
and $(e^{-x},E_{\nu+1}(x))$ on $[0,\infty)$. By reversing the order of the
weights, we end up in a setting that is dual to the one before. Now, we are
looking for polynomials $C_{\vec{n},\vec{m}}$ and $D_{\vec{n},\vec{m}}$ with
$\deg C_{\vec{n},\vec{m}}\leq n_{1}-1$ and $\deg D_{\vec{n},\vec{m}}\leq
n_{2}-1$, for which the function

$G_{\vec{n},\vec{m}}(x)=C_{\vec{n},\vec{m}}(x)x^{\alpha}+D_{\vec{n},\vec{m}}(x)x^{\beta}$
(26)

satisfies the orthogonality conditions

$\displaystyle\int_{0}^{\infty}G_{\vec{n},\vec{m}}(x)x^{k}e^{-x}\,dx$
$\displaystyle=0,\qquad k=0,\dots,m_{1}-1,$
$\displaystyle\int_{0}^{\infty}G_{\vec{n},\vec{m}}(x)x^{k}E_{\nu+1}(x)\,dx$
$\displaystyle=0,\qquad k=0,\dots,m_{2}-1.$

The (degenerate) case $n_{2}=0$ corresponds to the setting of the type II
polynomials of before. This will again lead to some similar properties. The
(degenerate) case $m_{2}=0$ corresponds to the setting of the type I multiple
Laguerre polynomials of the first kind. The associated type I functions will
be denoted by $L_{n,m}^{(I\mid\alpha,\beta)}(x)$ and are of the form

$L_{n,m}^{(I\mid\alpha,\beta)}(x)=C_{n,m}(x)x^{\alpha}+D_{n,m}(x)x^{\beta},$
(27)

where $C_{n,m}$ and $D_{n,m}$ are polynomials with $\deg C_{n,m}=n-1$ and
$\deg D_{n,m}=m-1$. They satisfy the orthogonality conditions

$\int_{0}^{\infty}L_{n,m}^{(I,\alpha,\beta)}(x)x^{k}e^{-x}\,dx=0,\quad
k=0,\dots,n+m-2,$

and are normalized such that

$\int_{0}^{\infty}L_{n,m}^{(I,\alpha,\beta)}(x)x^{n+m-1}e^{-x}\,dx=1.$

We don’t believe any explicit formulas for the type I multiple Laguerre
polynomials of the first kind have been written down before, so we will handle
that case first.

In the derivation of the mixed type polynomials, we will assume that the
multi-index $\vec{m}$ satisfies $m_{1}+1\geq m_{2}$. This is the natural
extension of the restriction we had before.

####  4.2.1 Type I multiple Laguerre polynomials of the first kind

######  Proposition 4.6.

The polynomials $C_{n,m}(x)$ and $D_{n,m}(x)$ in (27) are given by

$\displaystyle
C_{n+1,m+1}(x)=\frac{1}{Z_{n,m}}\sum_{k=0}^{n}\binom{n+\alpha-\beta}{n-k}\binom{m-\alpha+\beta}{k}\frac{x^{k}}{\Gamma(k+\alpha+1)},$
(28) $\displaystyle
D_{n+1,m+1}(x)=-\frac{1}{Z_{n,m}}\sum_{k=0}^{m}\binom{m-\alpha+\beta}{m-k}\binom{n+\alpha-\beta}{k}\frac{x^{k}}{\Gamma(k+\beta+1)},$
(29)

where
$Z_{n,m}=(-1)^{n+m}(\alpha-\beta)(\alpha-\beta+1)_{n}(\beta-\alpha+1)_{m}$ is
a normalizing constant. In particular,
$D_{n,m}^{\alpha,\beta}(x)=C_{m,n}^{\beta,\alpha}(x)$.

###### Proof.

We will prove this inductively on the sum $n+m$. For $n+m$ equal to $0$ or
$1$, we can verify the formula by hand or by means of a computer algebra
system. Now, let $(n,m)$ be a multi-index with $n+m\geq 1$ and suppose that
the formulas (28)–(29) hold for all multi-indices that add up to a smaller
amount. We will prove that they are also valid for $(n,m)$ by making use of
the two nearest neighbor recurrence relations

$xL_{n,m}^{(I\mid\alpha,\beta)}(x)=L_{n-1,m}^{(I\mid\alpha,\beta)}(x)+c_{n-1,m}L_{n,m}^{(I\mid\alpha,\beta)}(x)+a_{n,m}L_{n+1,m}^{(I\mid\alpha,\beta)}(x)+b_{n,m}L_{n,m+1}^{(I\mid\alpha,\beta)}(x),$
$xL_{n,m}^{(I\mid\alpha,\beta)}(x)=L_{n,m-1}^{(I\mid\alpha,\beta)}(x)+d_{n,m-1}L_{n,m}^{(I\mid\alpha,\beta)}(x)+a_{n,m}L_{n+1,m}^{(I\mid\alpha,\beta)}(x)+b_{n,m}L_{n,m+1}^{(I\mid\alpha,\beta)}(x).$

If we subtract them, we find another (simpler) recurrence relation

$(d_{n,m-1}-c_{n-1,m})L_{n,m}^{(I\mid\alpha,\beta)}(x)=L_{n-1,m}^{(I\mid\alpha,\beta)}(x)-L_{n,m-1}^{(I\mid\alpha,\beta)}(x).$

We know that the corresponding type II multiple orthogonal polynomials
$L_{n,m}^{(II\mid\alpha,\beta)}$, which are multiple Laguerre polynomials of
the first kind, satisfy nearest neighbor recurrence relations with similar
coefficients

$xL_{n,m}^{(II\mid\alpha,\beta)}(x)=L_{n+1,m}^{(II\mid\alpha,\beta)}(x)+c_{n,m}L_{n,m}^{(II\mid\alpha,\beta)}(x)+a_{n,m}L_{n-1,m}^{(II\mid\alpha,\beta)}(x)+b_{n,m}L_{n,m-1}^{(II\mid\alpha,\beta)}(x),$
$xL_{n,m}^{(II\mid\alpha,\beta)}(x)=L_{n,m-1}^{(II\mid\alpha,\beta)}(x)+d_{n,m}L_{n,m}^{(II\mid\alpha,\beta)}(x)+a_{n,m}L_{n-1,m}^{(II\mid\alpha,\beta)}(x)+b_{n,m}L_{n,m-1}^{(II\mid\alpha,\beta)}(x).$

The coefficients in the nearest neighbor recurrence relations for the multiple
Laguerre polynomials of the first kind were determined in [35]. In particular,

$c_{n,m}=2n+m+\alpha+1,\quad d_{n,m}=n+2m+\beta+1.$

Let us denote the coefficient of the monomial $x^{k}$ in a polynomial $P(x)$
by $P[k]$. If we compare the coefficient of $x^{k}$ in the $C$-polynomial on
each side of the equation

$(-n+m-\alpha+\beta)L_{n,m}^{(I\mid\alpha,\beta)}(x)=L_{n-1,m}^{(I\mid\alpha,\beta)}(x)-L_{n,m-1}^{(I\mid\alpha,\beta)}(x),$

we get

$(-n+m-\alpha+\beta)C_{n,m}[k]=C_{n-1,m}[k]-C_{n,m-1}[k].$

We may compute this difference via the induction hypothesis; it is equal to

$\displaystyle\frac{1}{Z_{n-1,m}}\binom{n-1+\alpha-\beta}{n-1-k}\binom{m-\alpha+\beta}{k}-\frac{1}{Z_{n,m-1}}\binom{n+\alpha-\beta}{n-k}\binom{m-1-\alpha+\beta}{k}$
$\displaystyle=\frac{1}{Z_{n,m}}\binom{n+\alpha-\beta}{n-k}\binom{m-\alpha+\beta}{k}\left(\frac{Z_{n,m}}{Z_{n-1,m}}\frac{n-k}{n+\alpha-\beta}-\frac{Z_{n,m}}{Z_{n,m-1}}\frac{m-k-\alpha+\beta}{m-\alpha+\beta}\right)$
$\displaystyle=\frac{1}{Z_{n,m}}\binom{n+\alpha-\beta}{n-k}\binom{m-\alpha+\beta}{k}(-n+m-\alpha+\beta).$

This is precisely what we wanted to show for $C_{n,m}[k]$. We have therefore
proven (28). The formula for $D_{n,m}(x)$ in (29) then follows from the
symmetry of the problem. ∎

Some elementary computations show that the type I Laguerre polynomials of the
first kind are the following hypergeometric polynomials

$C_{n,m}(x)=\frac{(-1)^{n+m+1}}{\Gamma(\alpha+1)(\beta-\alpha)_{m}(n-1)!}{}_{2}F_{2}\left(\begin{array}[]{c}-n+1,-m+\alpha-\beta+1\\\
\alpha+1,\alpha-\beta+1\end{array};x\right),$
$D_{n,m}(x)=\frac{(-1)^{n+m+1}}{\Gamma(\beta+1)(\alpha-\beta)_{n}(m-1)!}{}_{2}F_{2}\left(\begin{array}[]{c}-m+1,-n-\alpha+\beta+1\\\
\beta+1,-\alpha+\beta+1\end{array};x\right).$

####  4.2.2 Differentiation formula

We will be able to determine an explicit expression for the mixed type
polynomials by making use of the lemma below, which is the analogue of the
reduction property in Lemma 3.1. In this lemma, we explicitly wrote down the
dependence of $G_{\vec{n},\vec{m}}=G_{\vec{n},\vec{m}}^{\alpha,\beta,\nu}$ on
the underlying parameters.

######  Lemma 4.7.

One has,

$\frac{d}{dx}\left[x^{\nu}G_{\vec{n},\vec{m}}^{\alpha,\beta,\nu}(x)\right]\propto
x^{\nu-1}G_{\vec{n},(m_{1}-1,m_{2}+1)}^{\alpha,\beta,\nu-1}(x).$

###### Proof.

We will first verify that the left hand side has the right structure; it is
given by

$((\alpha+\nu)C_{\vec{n},\vec{m}}(x)+xC_{\vec{n},\vec{m}}^{\prime}(x))x^{\alpha+\nu-1}+((\beta+\nu)D_{\vec{n},\vec{m}}(x)+xD_{\vec{n},\vec{m}}^{\prime}(x))x^{\beta+\nu-1},$

where the polynomials between the brackets are still of degree $n_{1}$ and
$n_{2}$ respectively. It then remains to show that the required orthogonality
conditions are satisfied, i.e.,

$\int_{0}^{\infty}\frac{d}{dx}\left[x^{\nu}G_{\vec{n},\vec{m}}^{\alpha,\beta,\nu}(x)\right]e^{-x}x^{k-\nu+1}\,dx=0,\quad
k=0,\dots,m_{1}-2,$
$\int_{0}^{\infty}\frac{d}{dx}\left[x^{\nu}G_{\vec{n},\vec{m}}^{\alpha,\beta,\nu}(x)\right]E_{\nu}(x)x^{k-\nu+1}\,dx=0,\quad
k=0,\dots,m_{2}.$

We can work out both integrals via integration by parts. The first one becomes

$-\int_{0}^{\infty}G_{\vec{n},\vec{m}}^{\alpha,\beta,\nu}(x)e^{-x}x^{k}(-x+(k-\nu+1))\,dx,$

and indeed vanishes for $k=0,\dots,m_{1}-2$. Using the same argument as in the
proof of Lemma 3.1, the second integral becomes

$\int_{0}^{\infty}G_{\vec{n},\vec{m}}^{\alpha,\beta,\nu}(x)x^{k-1}((x-k)e^{-x}+k\nu
E_{\nu+1}(x))\,dx.$

For $k=0$, this is

$\int_{0}^{\infty}G_{\vec{n},\vec{m}}(x)e^{-x}\,dx=0,$

while for $k>0$, the integral vanishes if $k\leq m_{1}-1$ and $k\leq m_{2}$.
We thus also obtain the desired orthogonality conditions with respect to the
second weight, since $m_{1}\geq m_{2}+1$ by assumption. ∎

We can now express the mixed type function in terms of the type I function of
the type I multiple Laguerre polynomials of the first kind.

######  Theorem 4.8.

There is a mixed type function $G_{\vec{n},\vec{m}}$ given by

$G_{\vec{n},\vec{m}}(x)=x^{-\nu}\frac{d^{m_{2}}}{dx^{m_{2}}}\left[x^{m_{2}+\nu}L_{\vec{n}}^{(I\mid\alpha,\beta)}(x)\right].$
(30)

###### Proof.

It follows from the previous lemma that

$G_{\vec{n},\vec{m}}^{\alpha,\beta,\nu}(x)\propto
x^{-\nu}\frac{d}{dx}\left[x^{\nu+1}G_{\vec{n},(m_{1}+1,m_{2}-1)}^{\alpha,\beta,\nu+1}(x)\right],$

for $m_{1}+1\geq m_{2}$. After extracting $m_{2}-1$ additional derivatives, we
thus get

$G_{\vec{n},\vec{m}}^{\alpha,\beta,\nu}(x)\propto
x^{-\nu}\frac{d^{m_{2}}}{dx^{m_{2}}}\left[x^{m_{2}+\nu}G_{\vec{n},(m_{1}+m_{2},0)}^{\alpha,\beta,\nu+m_{2}}(x)\right].$

The multi-index $(m_{1}+m_{2},0)$ corresponds to the type I setting with
$m_{1}+m_{2}=n_{1}+n_{2}-1$ orthogonality conditions and therefore
$G_{\vec{n},(m_{1}+m_{2},0)}^{\alpha,\beta,\nu+m_{2}}(x)\propto
L_{\vec{n}}^{(I\mid\alpha,\beta)}(x)$. ∎

####  4.2.3 Hypergeometric series

The mixed type polynomials can be represented as ${}_{3}F_{3}$ hypergeometric
polynomials in the following way.

######  Proposition 4.9.

The mixed type polynomials $C_{\vec{n},\vec{m}}(x)$ and
$D_{\vec{n},\vec{m}}(x)$ in the mixed type function (30) arise as the
hypergeometric series

$C_{\vec{n},\vec{m}}(x)=\mathcal{C}_{\vec{n},\vec{m}}\,{}_{3}F_{3}\left(\begin{array}[]{c}-n_{1}+1,-n_{2}+\alpha-\beta+1,m_{2}+\alpha+\nu+1\\\
\alpha+1,\alpha-\beta+1,\alpha+\nu+1\end{array};x\right),$
$D_{\vec{n},\vec{m}}(x)=\mathcal{D}_{\vec{n},\vec{m}}\,{}_{3}F_{3}\left(\begin{array}[]{c}-n_{2}+1,-n_{1}-\alpha+\beta+1,m_{2}+\beta+\nu+1\\\
\beta+1,-\alpha+\beta+1,\beta+\nu+1\end{array};x\right),$

where

$\mathcal{C}_{\vec{n},\vec{m}}=\frac{(-1)^{n_{1}+n_{2}+1}(\alpha+\nu+1)_{m_{2}}}{\Gamma(\alpha+1)(\beta-\alpha)_{n_{2}}(n_{1}-1)!},\quad\mathcal{D}_{\vec{n},\vec{m}}=\frac{(-1)^{n_{1}+n_{2}+1}(\beta+\nu+1)_{m_{2}}}{\Gamma(\beta+1)(\alpha-\beta)_{n_{1}}(n_{2}-1)!}.$

###### Proof.

This follows from the hypergeometric series representations for the type I
multiple Laguerre polynomials of the first kind and the differentiation
properties of the generalized hypergeometric series. ∎

##  5 Related multiple orthogonal polynomials

The type I and II multiple orthogonal polynomials that were investigated here
are related to the two-parameter subfamily of the Jacobi-Piñeiro polynomials
in, e.g., [29], and the multiple orthogonal polynomials associated with
$K$-Bessel functions in [37], with confluent hypergeometric functions in [24]
and with Gauss’ hypergeometric functions in [25]. In the latter, one considers
multiple orthogonal polynomials on the step-line associated with two weights
$\mathcal{W}(x;a,b,c,d)$ and $\mathcal{W}(x;a,b+1,c+1,d)$. The weights are
defined in terms of Gauss’ hypergeometric function and have moments of the
form

$\int_{0}^{\infty}\mathcal{W}(x;a,b,c,d)x^{k}dx=\frac{(a+1)_{k}(b+1)_{k}}{(c+1)_{k}(d+1)_{k}}.$

One assumes that $a,b,c,d\in\mathbb{R}_{>-1}$, $c+1,d>a$ and $c,d>b$ for
regularity reasons.

The connection to the multiple orthogonal polynomials in this article is given
in Figure 2 below. It is based on the formulas for the Mellin transforms of
the type I functions (where we might have had to rebalance the underlying
weights) as given in Figure 3.

Figure 2: Connection to other weights.

The vertical order is according to the type of hypergeometric function that is
used in the expression for the type I functions (and the type II polynomials)
which can be obtained by applying the inverse Mellin transform. Moreover, the
asymptotic zero distribution of the type II multiple orthogonal polynomials on
the step-line is the same on each row. Moving horizontally in the figure, one
obtains formulas that are valid for a larger amount of multi-indices.
Classical orthogonal polynomials appear first as degenerate cases of the
multiple orthogonal polynomials in the second column.

Figure 3 suggests that in a setting with a weight that has moments of the form

$\frac{(a_{1}+1)_{k}\dots(a_{r}+1)_{k}}{(b_{1}+1)_{k}\dots(b_{r}+1)_{k}},$

the Mellin transforms of the type I functions on the step-line will be given
by

$\frac{\Gamma(s+a_{1})\dots\Gamma(s+a_{r})}{\Gamma(s+b_{1}+n_{1})\dots\Gamma(s+b_{r}+n_{r})}(1-s)_{n_{1}+\dots
n_{r}-1}$

and that we have to consider $r$ weights in total to attain this. For an
appropriate choice of weights, it can be verified that this is indeed the
case. In that case, the type I function arises as an $r$-fold Mellin
convolution in which Jacobi polynomials and beta densities appear. An
application of the inverse Mellin transform then allows us to obtain a formula
for the type I function in terms of $r$ hypergeometric functions of type
${}_{r+1}F_{r}$. From this, we may derive expressions for the type II
polynomials as well.

Figure 3: Underlying Mellin transforms.

By taking suitable limits, it is also possible to consider weights with fewer
Pochhammer symbols in the denominator. The associated multiple orthogonal
polynomials will appear in the study of the singular values of (mixed)
products of Ginibre matrices and truncated unitary matrices. The multiple
orthogonal polynomials associated with weights with no Pochhammer in the
denominator were already investigated in [22], where they appeared in the
study of products of Ginibre matrices. Weight functions for which the moments
are ratios of Pochhammer symbols also appear naturally in the work of Sokal
[30], which makes a connection between multiple orthogonal polynomials,
production matrices and branched continued fractions.

We may extend the multiple orthogonal polynomials in this article to multiple
orthogonal polynomials with respect to more weights by considering additional
exponential integral weights $(w_{1},\dots,w_{r})$, with
$w_{j}(x)=x^{\alpha}E_{\nu_{j}+1}(x)$, that have parameters $\nu_{1}=-1$ (so
that $w_{1}(x)=x^{\alpha}e^{-x}$) and $\nu_{i}-\nu_{j}\notin\mathbb{Z}$
whenever $i\neq j$. This system of weights will again be a Nikishin system
under some additional restrictions on the parameters.

## References

  * [1] C. Álvarez-Fernández, U. Fidalgo Prieto, M. Mañas, Multiple orthogonal polynomials of mixed type: Gauss-Borel factorization and the multi-component 2D Toda hierarchy, Adv. Math. 227 (2011), no. 4, 1451–1525.
  * [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] A. I. Aptekarev, A. Branquinho, W. Van Assche, Multiple orthogonal polynomials for classical weights, Trans. Amer. Math. Soc. 355 (2003), no. 10, 3887–3914.
  * [5] 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.
  * [6] E. Coussement, J. Coussement, W. Van Assche, Asymptotic zero distribution for a class of multiple orthogonal polynomials, Trans. Amer. Math. Soc. 360 (2008), no. 10, 5571–5588.
  * [7] E. Daems, A.B.J. Kuijlaars, Multiple orthogonal polynomials of mixed type and non-intersecting Brownian motions, J. Approx. Theory 146 (2007), no. 1, 91–114.
  * [8] NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/, Release 1.1.3 of 2021-09-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] U. Fidalgo Prieto, A. López García, G. López Lagomasino, V.N. Sorokin, Mixed type multiple orthogonal polynomials for two Nikishin systems, Constr. Approx. 32 (2010), no. 2, 255–306.
  * [12] U. Fidalgo, S. Medina Peralta, J. Mínguez Ceniceros, Mixed type multiple orthogonal polynomials: perfectness and interlacing properties of zeros, Linear Algebra Appl. 438 (2013), no. 3, 1229–1239.
  * [13] W. Gautschi, Polynomials orthogonal with respect to exponential integrals, Numer. Algor. 70 (2015), 215-–226.
  * [14] J.S. Geronimo, T.P. Hill, Necessary and sufficient condition that the limit of Stieltjes transforms is a Stieltjes transform, J. Approx. Theory 121 (2003), 54–60.
  * [15] J. Ginibre, Statistical ensembles of complex, quaternion, and real matrices, J. Math. Phys. 6 (1965), 440–449.
  * [16] I.S. Gradshteyn, I.M. Ryzhik, Table of Integrals, Series and Products, Academic Press, New York, 1980 (corrected and enlarged edition).
  * [17] M.E.H. Ismail, Classical and Quantum Orthogonal Polynomials in One Variable, Encyclopedia of Mathematics and its Applications 98, Cambridge University Press, 2008.
  * [18] A.B.J. Kuijlaars, Multiple orthogonal polynomials in random matrix theory, Proceedings of the International Congress of Mathematicians. Volume III, Hindustan Book Agency, New Delhi, 2010, pp. 1417–1432.
  * [19] 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, 2010, pp. 155–176
  * [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, A. Loureiro, Multiple orthogonal polynomials associated with confluent hypergeometric functions, J. Approx. Theory 260 (2020), 105484.
  * [25] H. Lima, A. Loureiro, Multiple orthogonal polynomials with respect to Gauss’ hypergeometric function, Stud. in Appl. Math. 148 (2021), no. 1., 154–185.
  * [26] K. Mahler, Perfect systems, Compos. Math. 19 (1968), 95–166.
  * [27] E.M. Nikishin, V.N. Sorokin, Rational Approximations and Orthogonality, Translations of Mathematical Monographs 92, Amer. Math. Soc., Providence RI, 1991.
  * [28] F.W.J. Olver, D.W. Lozier, R.F. Boisvert, C.W. Clark, NIST Handbook of Mathematical Functions, Cambridge University Press, 2010.
  * [29] 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.
  * [30] A.D. Sokal, Multiple orthogonal polynomials, $d$-orthogonal polynomials, production matrices, and branched continued fractions, arXiv:2204.11528 [math.CA], to appear in Trans. Amer. Math. Soc.
  * [31] V.N. Sorokin, Simultaneous approximation of several linear forms. (Russian. English summary), Vestnik Moskov. Univ. Ser. I Mat. Mekh. 101(1983), no. 1, 44–47;
  * [32] V.N. Sorokin, Simultaneous Padé approximation for functions of Stieltjes type, Sib. Mat. Zh. 31, no. 5 (1990), 128-–137; English translation in Sib. Math. J. 31 (1990), no. 5, 809–817.
  * [33] G. Szegő, Orthogonal Polynomials, Amer. Math. Soc. Colloq. Publ. 23, Amer. Math. Soc., Providience RI, 1939 (fourth edition 1975).
  * [34] W. Van Assche, Mehler-Heine asymptotics for multiple orthogonal polynomials, Proc. Amer. Math. Soc. 145 (2017), 303–314.
  * [35] W. Van Assche, Nearest neighbor recurrence relations for multiple orthogonal polynomials, J. Approx. Theory 163 (2011), 1427–1448.
  * [36] 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.
  * [37] W. Van Assche, S. Yakubovich Multiple orthogonal polynomials associated with Macdonald functions, Integral Transf. Spec. Funct. 9 (2000), no. 3, 229–244.