Figure 1. Zeros of the multiple Hermite polynomial $\mathfrak{H}_{\bm{n}}$
(indicated by empty circles, all on the real axis) and of $S_{\bm{n},1}$
(filled circles) for $\bm{n}=(35,35)$, $c_{1}=-c_{2}$, for different values of
$c$: $c=1$ (top left), $c=4$ (top right), $c=8$ (bottom left) and $c=16$
(bottom right). The zeros of $S_{\bm{n},1}$ that are on the left semi-axis
apparently interlace with the zeros of $\mathfrak{H}_{\bm{n}}$.

A generalization of multiple Hermite polynomials to the case of polynomials
$P_{\bm{n}}$, $\bm{n}=(n,n)$, satisfying the varying orthogonality conditions

$\int_{-\infty}^{\infty}x^{k}P_{\bm{n}}(x)e^{-n(V(x)\pm cx)}\,dx=0,\quad
k=0,1,\dots,n-1,$

where $V$ is a polynomial of even degree and positive leading coefficient, has
been carried out in [19], again associated to random matrix models with
external source. The limit zero distribution of $P_{\bm{n}}$’s was described
there in terms of a constrained 2-component vector equilibrium problem, with
one of the components on the imaginary axis. The two-cut situation in the
asymptotics of multiple Hermite polynomials, and thus the reduction to the
Angelesco equilibrium, is in this case equivalent to the constraint on the
imaginary axis to be not achieved (not “saturated”). The study in [19] has
been extended in [75] (see also [11] for the quartic case) to a non-symmetric
situation, showing that it can be alternatively characterized in terms of a
3-component critical vector measure with the interaction matrix from Corollary
6.3. The curves outlined by the zeros of $P_{\bm{n}}$ and $S_{\bm{n},1}$ in
Figure 1 are consistent with the support of the three components described in
[75, Theorem C].

###  8.3. Multiple Laguerre polynomials of the first kind

These polynomials are defined by the orthogonality conditions

$\begin{aligned}
&\int_{0}^{\infty}x^{k}\mathfrak{L}_{\bm{n}}(x)x^{\alpha_{1}}e^{-x}\,dx=0,\quad
k=0,1,\ldots,n_{1}-1\\\
&\int_{0}^{\infty}x^{k}\mathfrak{L}_{\bm{n}}(x)x^{\alpha_{2}}e^{-x}\mathrm{~{}d}x=0,\quad
k=0,1,\ldots,n_{2}-1,\end{aligned}\quad\deg\mathfrak{L}_{\bm{n}}\leq
N=n_{1}+n_{2},$

where $\alpha_{1},\alpha_{2}>0$ and $\alpha_{1}-\alpha_{2}\notin\mathbb{Z}$,
under which condition the weights form an AT-system; see [7], [52, section
23.4.1], and [100, Section 3.2]. Not only that, since $x^{\beta}$, for
$\beta<0$ and $x>0$, can be written as the Cauchy integral of a positive
weight on $(-\infty,0)$, coinciding up to a multiplicative constant with
$|x|^{\beta}$, we conclude that this pair of weight forms a Nikishin system
with $\Delta_{1}=\Delta_{2}=[0,+\infty)$ and $[c,d]=[-\infty,0]$, see Section
5.2.2.

Polynomials $\mathfrak{L}_{\bm{n}}$ can be obtained by using the Rodrigues
formula,

$(-1)^{N}e^{-x}\mathfrak{L}_{\bm{n}}(x)=\prod_{j=1}^{2}\left(x^{-\alpha_{j}}\frac{d^{n_{j}}}{dx^{n_{j}}}x^{n_{j}+\alpha_{j}}\right)e^{-x},\quad
N=n_{1}+n_{2},$

from which one can find the explicit expression

$\displaystyle(-1)^{N}e^{-x}\mathfrak{L}_{\bm{n}}(x)$
$\displaystyle=\left(\alpha_{1}+1\right)_{n_{1}}\left(\alpha_{2}+1\right)_{n_{2}}{}_{2}F_{2}\left(\begin{array}[]{c}\alpha_{1}+n_{1}+1,\alpha_{2}+n_{2}+1\\\
\alpha_{1}+1,\alpha_{2}+1\end{array}\bigg{|}-x\right).$

In our notation,

$A_{1}(x)=A_{2}(x)=x,\quad B_{i}(x)=\alpha_{i}-x,\quad\sigma_{i}=0,\quad
i=1,2,$

with $\Delta_{1}=\Delta_{2}=[0,+\infty)$. These polynomials satisfy the
differential equation (7.2), which takes the form

$x^{2}R_{\bm{n}}(x)y^{\prime\prime\prime}(x)+\left[(B_{1}(x)+B_{2}(x)+4)xR_{\bm{n}}(x)-x^{2}R_{\bm{n}}^{\prime}(x)\right]y^{\prime\prime}(x)+E_{\bm{n}}(x)y^{\prime}(x)+F_{\bm{n}}(x)y(x)=0.$
(8.5)

It follows from formula (5.26) that $R_{\bm{n}}(x)/x$ is a constant. Thus, the
equation reduces to

$x^{3}y^{\prime\prime\prime}+x^{2}\left[B_{1}(x)+B_{2}(x)+3\right]y^{\prime\prime}+E_{\bm{n}}y^{\prime}+F_{\bm{n}}y=0,$

where $E_{\bm{n}}$ and $F_{\bm{n}}$ are of degrees at most $3$ and $2$,
respectively. As in Section 8.2, the asymptotics of the corresponding
functions of second kind $q_{\bm{n},1}$ and $q_{\bm{n},2}$ at $\infty$ yields
some constraints on the coefficients of $E_{\bm{n}}$ and $F_{\bm{n}}$, which
unfortunately are not sufficient to determine the polynomials in this case. We
need to make use also of the predicted behavior of the solutions at the
origin.

Notice that $0$ is a regular singular point (a Fuchsian singularity) of (8.5).
The fact that the weights constitute an AT-system on $[0,+\infty)$ implies
also that $\mathfrak{L}_{n_{1},n_{2}}(0)\neq 0$. In consequence,
$E_{\bm{n}}(0)=0$ and $F_{\bm{n}}(0)=0$. Expanding the solutions at the origin
we conclude that the indicial polynomial must vanish at $0$, $-\alpha_{1}$ and
$-\alpha_{2}$. All this additional information allows us to determine
$E_{\bm{n}}$ and $F_{\bm{n}}$:

$\displaystyle E_{\bm{n}}$
$\displaystyle=x(x^{2}+(n_{1}+n_{2}-\alpha_{1}-\alpha_{2}-3)x+(1+\alpha_{1})(1+\alpha_{2})\,,$
$\displaystyle F_{\bm{n}}$
$\displaystyle=x(-(n_{1}+n_{2})x+n_{1}+n_{2}+n_{1}n_{2}+\alpha_{1}n_{2}+\alpha_{2}n_{1})\,.$

Canceling the common factor $x$ in the four coefficients of (8.5) yields the
equation that appeared already in [7, Section 4.3],

$\displaystyle x^{2}y^{\prime\prime\prime}(x)$
$\displaystyle+\left(-2x^{2}+\left(\alpha_{1}+\alpha_{2}-3\right)x\right)y^{\prime\prime}(x)+\left(x^{2}-x\left(\alpha_{1}+\alpha_{2}-n_{1}-n_{2}+3\right)\right.$
$\displaystyle\left.+\left(\alpha_{1}+1\right)\left(\alpha_{2}+1\right)\right)y^{\prime}(x)-\left(x(n_{1}+n_{2})-\left(n_{1}+n_{2}+n_{1}n_{2}+\alpha_{1}n_{2}+\alpha_{2}n_{1}\right)\right)y(x)=0.$

From the conclusions of Section 5.2.2 it follows that the zeros of
$\mathfrak{L}_{\bm{n}}$ are all positive, and those of $S_{\bm{n},1}$ are
negative. By Theorem 5.10, the discrete vector measure
$\vec{\nu}_{1}:=(\nu(\mathfrak{L}_{\bm{n}}),\nu(S_{\bm{n},1}))$ is a critical
vector measure for the energy functional $\mathcal{E}_{\vec{\varphi},a}$, with
$a=-1/2$ and

$\vec{\varphi}(x)=\frac{1}{2}\left(x-(\alpha_{1}+1)\log(x),(\alpha_{1}-\alpha_{2}-1)\log|x|\right).$

Direct computations show that

$\int_{0}^{\infty}x^{k+\alpha_{j}}e^{-x}dx=\begin{cases}\sqrt{\pi}\,2^{-k-1}(2k+1)!!,&j=1,\\\
(k+1)!,&j=2,\end{cases}\quad k=0,1,\dots,$

which allows to find the moments of $w_{i}$, the asymptotic expansion of
$\mathfrak{C}_{w_{i}}[\mathfrak{L}_{\bm{n}}]$ at infinity, and in consequence,
$S_{\bm{n},i}$ (using formula (5.12)).

Let us consider the particular case of $n_{1}=n_{2}=n$, with $\alpha_{1}=1/2$,
$\alpha_{2}=1$. Then, for $\bm{n}=(5,5)$,

$\displaystyle\mathfrak{L}_{(5,5)}(x)=$ $\displaystyle
x^{10}-\frac{165x^{9}}{2}+\frac{5445x^{8}}{2}-\frac{186615x^{7}}{4}+\frac{7224525x^{6}}{16}-\frac{80613225x^{5}}{32}+\frac{127182825x^{4}}{16}$
$\displaystyle-\frac{107120475x^{3}}{8}+10758825x^{2}-\frac{13253625x}{4}+\frac{467775}{2},$

and up to normalization,

$\displaystyle S_{\bm{n},1}(x)$
$\displaystyle=8x^{5}+2720x^{4}+107500x^{3}+1945020x^{2}+46682295x+1425581520,$
$\displaystyle S_{\bm{n},2}(x)$
$\displaystyle=32x^{5}+2960x^{4}+67424x^{3}+1313480x^{2}+37066290x+1173966885.$

All zeros of $\mathfrak{L}_{(5,5)}$ are positive and simple, while each
$S_{\bm{n},j}$ has one negative and two pairs of complex conjugate simple
zeros, all of them simple.

Asymptotics of sequences of (rescaled) multiple Laguerre polynomials of the
first kind

$p_{(n,n)}(t)=\mathfrak{L}_{(n,n)}\left(2nt\right)$

was obtained in [31] and [77]. It was shown that the support of the weak-*
limit of the zero-counting measures $\nu(p_{(n,n)})$ is the interval
$\left[0,27/8\right]$, with the density presenting the usual square root
vanishing at the rightmost endpoint of the support. The expression of the
density was derived from the recurrence relation satisfied by polynomials
$\mathfrak{L}_{(n_{1},n_{2})}$ and no equilibrium problem associated to that
distribution was given.

\begin{overpic}[scale={0.7}]{Fig2} \end{overpic} Figure 2. Zeros of the
multiple Laguerre polynomial of the first kind $\mathfrak{L}_{\bm{n}}$
(indicated by empty circles, all on the positive semiaxis) and of
$S_{\bm{n},1}$ (filled circles, all on the negative semiaxis) for
$\bm{n}=(35,35)$ and $\alpha_{1}=1/2$, $\alpha_{2}=1$. Nine real zeros of
$S_{\bm{n},1}$ (ranging from $-74000$ to $-201.53$, are not represented.

Again, it is interesting to compare these conclusions with results of the
numerical experiments presented in Figure 2, where we take $\bm{n}=(35,35)$,
with $\alpha_{1}=1/2$, and $\alpha_{2}=1$. The largest zero of
$\mathfrak{L}_{(35,35)}$ is $217.597$, which is consistent with the expected
value of

$\frac{27}{8}\times 70=236.25.$

###  8.4. Multiple Laguerre polynomials of the second kind

These polynomials are defined by the orthogonality conditions

$\displaystyle\int_{0}^{\infty}x^{k}\mathcal{L}_{\bm{n}}(x)x^{\alpha}e^{-c_{1}x}\,dx=0,\quad
k=0,1,\ldots,n_{1}-1,$
$\displaystyle\int_{0}^{\infty}x^{k}\mathcal{L}_{\bm{n}}(x)x^{\alpha}e^{-c_{2}x}\,dx=0,\quad
k=0,1,\ldots,n_{2}-1,$

where we assume that $\alpha>0$ and $c_{1},c_{2}>0$ with $c_{1}\neq c_{2}$,
under which condition the weights form an AT-system; see e.g. [21], [63], or
[100, section 3.3].

An explicit expression can be found in [7, Section 3] or [52, Section 23.4]:

$\mathcal{L}_{\bm{n}}(x)=\sum_{k_{1}=0}^{n_{1}}\sum_{k_{2}=0}^{n_{2}}(-1)^{k_{1}+k_{2}}\frac{(k_{1}+k_{2})!}{c_{1}^{k_{1}}c_{2}^{k_{2}}}\binom{n_{1}}{k_{1}}\binom{n_{2}}{k_{2}}\binom{N+\alpha}{k_{1}+k_{2}}x^{N-k_{1}-k_{2}}.$

In our notation,

$A_{1}(x)=A_{2}(x)=x,\quad B_{i}(x)=\alpha-c_{i}x,\quad\sigma_{i}=0,\quad
i=1,2,$

with $\Delta_{1}=\Delta_{2}=[0,+\infty)$. These polynomials also satisfy the
differential equation (8.5), that is,

$x^{2}R_{\bm{n}}(x)y^{\prime\prime\prime}(x)+\left[(-(c_{1}+c_{2})x^{2}+2(\alpha+2)x)R_{\bm{n}}(x)-x^{2}R_{\bm{n}}^{\prime}(x)\right]y^{\prime\prime}+E_{\bm{n}}(x)y^{\prime}(x)+F_{\bm{n}}(x)y(x)=0.$

Formula (5.26) shows that now $R_{\bm{n}}(x)/x$ is a polynomial of degree $1$.
Furthermore, by (5.21) and since $B_{1}-B_{2}$ vanishes at $0$,
$R_{\bm{n}}\mathfrak{L}_{\bm{n}}$ has a double root at $0$. Again, the fact
that the weights constitute an AT-system on $[0,+\infty)$ implies also that
$\mathfrak{L}_{\bm{n}}(0)\neq 0$. Hence, $R_{\bm{n}}(x)/x^{2}$ is a constant,
and the third order differential equation is

$x^{4}y^{\prime\prime\prime}-x^{3}((c_{1}+c_{2})x-\alpha-2)y^{\prime\prime}+E_{\bm{n}}y^{\prime}+F_{\bm{n}}y=0.$

Arguments as described in Sections 8.2 and 8.3, using the asymptotics of the
functions of second kind both at $0$ and at $\infty$ and the fact that the
indicial polynomial vanishes at $0$ and $-\alpha$, yield the expressions for
polynomials $E_{\bm{n}}$ and $F_{\bm{n}}$:

$\displaystyle E_{\bm{n}}$
$\displaystyle=c_{1}c_{2}x^{4}-[(c_{1}+c_{2})(\alpha+1)-n_{1}c_{1}-n_{2}c_{2}]x^{3}+\alpha(\alpha+1)x^{2},$
$\displaystyle F_{\bm{n}}$
$\displaystyle=-c_{1}c_{2}(n_{1}+n_{2})x^{3}+\alpha(c_{1}n_{1}+c_{2}n_{2})x^{2}.$

Canceling the common factor $x^{2}$ in the differential equation yields

$\begin{gathered}x^{2}y^{\prime\prime\prime}(x)-\left(x^{2}\left(c_{1}+c_{2}\right)-2x(\alpha+1)\right)y^{\prime\prime}(x)+\left(x^{2}c_{1}c_{2}-x\left[\left(c_{1}+c_{2}\right)(\alpha+1)-n_{1}c_{1}-n_{2}c_{2}\right]\right.\\\
+\alpha(\alpha+1))y^{\prime}(x)-\left(xc_{1}c_{2}(n_{1}+n_{2})-\alpha\left(n_{1}c_{1}+n_{2}c_{2}\right)\right)y(x)=0,\end{gathered}$

which matches the equation found in [7, section 4.3] and [43, section 5.2].

By Theorem 5.10, the discrete vector measure
$\vec{\nu}_{1}:=\left(\nu(\mathcal{L}_{\bm{n}}),\nu(S_{\bm{n},1})\right)$ is a
critical vector measure for the energy functional
$\mathcal{E}_{\vec{\varphi},a}$, with $a=-1/2$ and

$\vec{\varphi}(x)=\frac{1}{2}\left(c_{1}x-(\alpha+1)\log
x,(c_{2}-c_{1})x\right),\quad z=x+iy.$

Since for $c>0$ and $\alpha>-1$,

$\int_{0}^{\infty}x^{\alpha}e^{-cx}dx=c^{-\alpha-1}\Gamma(\alpha+1),$

we can easily calculate the moments of $w_{i}$, the asymptotic expansion of
$\mathfrak{C}_{w_{i}}[\mathcal{L}_{\bm{n}}]$ at infinity, and in consequence,
$S_{\bm{n},i}$ (using formula (5.12)).

Let us consider the particular case of $n_{1}=n_{2}=n$, with $\alpha=1$,
$c_{1}=1$, and $c_{2}=2$. Then, for $\bm{n}=(5,5)$,

$\displaystyle\mathcal{L}_{(5,5)}(x)=$ $\displaystyle
x^{10}-\frac{165x^{9}}{2}+2750x^{8}-\frac{96525x^{7}}{2}+487575x^{6}-\frac{5831595x^{5}}{2}+10239075x^{4}$
$\displaystyle-20270250x^{3}+20790000x^{2}-9355500x+1247400,$

and up to normalization,

$\displaystyle S_{\bm{n},1}(x)$
$\displaystyle=2x^{5}+25x^{4}+605x^{3}+16580x^{2}+506065x+16197810,$
$\displaystyle S_{\bm{n},2}(x)$
$\displaystyle=4x^{5}-320x^{4}+9975x^{3}-151645x^{2}+1115560x-2967600.$

All zeros of $\mathcal{L}_{(5,5)}$ are positive and simple, $S_{\bm{n},1}$ has
one negative and two pairs of complex conjugate simple zeros, while
$S_{\bm{n},2}$, has one positive and two pairs of complex conjugate roots, all
of them simple.

Asymptotics of sequences of (rescaled) multiple Laguerre polynomials of the
second kind,

$p_{(n,n)}(t)=\mathcal{L}_{(n,n)}\left(nt\right),$

with varying $0<c_{1}<c_{2}$ proportional to $n$, has been studied by Lysov
and Wielonsky in [63] using the Riemann-Hilbert technique and the analysis of
the Riemann surface derived from the differential equation. In particular,
they found that there is a critical value $\kappa\approx 12.11\dots$ such that
for $0<c_{2}/c_{1}<\kappa$, the support of the limit of the zero-counting
measures $\nu(p_{(n,n)})$ is a single interval of the form $[0,d]$, $d>0$,
while for $c_{2}/c_{1}>\kappa$, it is comprised of two real intervals
$[0,a]\cup[b,d]$, with $0<a<b<d$. The expression of the density was also
derived, but no equilibrium problem associated to that distribution was given.
It is interesting to compare these conclusions with results of the numerical
experiments presented in Figure 3, where we take $\bm{n}=(35,35)$. We observe
that for small values of $c_{2}/c_{1}$ the zeros of $S_{\bm{n},1}$ sit on a
curve on the complex plane. However, for large ratios $c_{2}/c_{1}$, zeros of
$\mathcal{L}_{(n,n)}$ split into two groups, and the zeros of $S_{\bm{n},1}$,
all real, approximately interlace with the zeros of $\mathcal{L}_{(n,n)}$ on
the leftmost subinterval. This allows us to conjecture that in this case, the
asymptotic zero distribution can be described again in terms of the Angelesco-
type vector equilibrium, see Section 6.3.

\begin{overpic}[scale={0.55}]{Fig3a} \end{overpic} | \begin{overpic}[scale={0.55}]{Fig3b} \end{overpic}
---|---
\begin{overpic}[scale={0.55}]{Fig3c} \end{overpic} | \begin{overpic}[scale={0.55}]{Fig3d} \end{overpic}
Figure 3. Zeros of the multiple Laguerre polynomial of the second kind
$\mathcal{L}_{\bm{n}}$ (indicated by empty circles, all on the positive
semiaxis) and of $S_{\bm{n},1}$ (filled circles) for $\bm{n}=(35,35)$, and
$(c_{1},c_{2})=(35,70)$ (top left), $(c_{1},c_{2})=(35,140)$ (top right), and
$(c_{1},c_{2})=(35,525)$ (bottom left). Bottom right: zoom of the interval
$(0,0.25)$ for $(c_{1},c_{2})=(35,525)$.

###  8.5. Jacobi-Piñeiro polynomials

The Jacobi-Piñeiro polynomials are multiple orthogonal polynomials associated
with an AT system consisting of Jacobi weights on $[0,1]$ with different
powers at 0 and the same behavior at $1$. They are defined by the
orthogonality conditions

$\displaystyle\int_{0}^{1}x^{k}P_{\bm{n}}(x)x^{\beta_{1}}(1-x)^{\alpha}\,dx=0,\quad
k=0,1,\ldots,n_{1}-1,$
$\displaystyle\int_{0}^{1}x^{k}P_{\bm{n}}(x)x^{\beta_{2}}(1-x)^{\alpha}\,dx=0,\quad
k=0,1,\ldots,n_{2}-1.$

In our notation, $\Delta_{1}=\Delta_{2}=[0,1]$ and
$w_{j}(x)=x^{\beta_{j}}(1-x)^{\alpha}$, $j=1,2$, with
$\alpha,\beta_{1},\beta_{2}>-1$ and $\beta_{1}-\beta_{2}\notin\mathbb{Z}$. We
have

$A_{1}(x)=A_{2}(x)=x(x-1),\quad
B_{i}(x)=(\beta_{i}+\alpha)x-\beta_{i},\quad\sigma_{i}=0,\quad i=1,2.$

For the same reason mentioned at the beginning of Section 8.3, this pair of
weight forms a Nikishin system with $\Delta_{1}=\Delta_{2}=[0,1]$ and
$[c,d]=[-\infty,0]$, see Section 5.2.2.

These polynomials were first studied by Piñeiro [83] when $\alpha=0$. The
general case appears in [81, p. 162]. There is a Rodrigues formula for Jacobi-
Piñeiro polynomials $P_{\bm{n}}$, $\bm{n}=(n_{1},n_{2})$, see [52, Section
23.3.2]: with $N=n_{1}+n_{2}$, and up to a constant factor,

$\displaystyle P_{\bm{n}}(x)$
$\displaystyle=(1-x)^{-\alpha}\prod_{j=1}^{2}\left(x^{-\beta_{j}}\frac{d^{n_{j}}}{dx^{n_{j}}}x^{n_{j}+\beta_{j}}\right)(1-x)^{\alpha+N}.$

There is even an explicit expression [100, Section 3.1]: again, up to a
multiplicative constant,

$P_{\bm{n}}(x)=\sum_{k=0}^{n_{1}}\binom{\beta_{1}+n_{1}}{k}\binom{\alpha+N}{n_{1}-k}\sum_{j=0}^{n_{2}}\binom{\beta_{2}+N-k}{j}\binom{\alpha+k+n_{2}}{n_{2}-j}x^{N-k-j}(x-1)^{k+j}.$

Proposition 5.5 assures that polynomial $R_{\bm{n}}$ has degree at most $3$.
By (5.21), $R_{\bm{n}}P_{\bm{n}}$ has a double root at $1$ (observe that
$B_{1}-B_{2}$ also vanish at $1$) and a simple one at $0$. Since the weights
form an AT-system, $P_{\bm{n}}(0)\neq 0$ and $P_{\bm{n}}(1)\neq 0$, so that,
up to a multiplicative constant, $R_{\bm{n}}(z)=z(z-1)^{2}$ .

By Theorem 5.10, and taking into account the expression for $R_{\bm{n}}$, we
conclude that the discrete vector measure
$\vec{\nu}_{1}:=(\nu(\mathfrak{L}_{\bm{n}}),\nu(S_{\bm{n},1}))$ is a critical
vector measure for the energy functional $\mathcal{E}_{\vec{\varphi},a}$, with
$a=-1/2$ and

$\vec{\varphi}=\left(\frac{\beta_{1}+1}{2}\log\frac{1}{\left|x\right|}+\frac{\alpha+1}{2}\log\frac{1}{\left|x-1\right|},\left(\beta_{1}-\beta_{2}+\frac{1}{2}\right)\log\frac{1}{\left|z\right|}\right),\quad
z=x+iy.$

The third order differential equation can be also obtained following the
arguments used in Sections 8.2–8.4, and making use of the asymptotics at
$\infty$, $0$, and $1$, and of the known roots of the indicial polynomials at
both finite points. For instance, in the case $\alpha=0$ (studied by Piñeiro
in [83]), polynomials $E_{\bm{n}}$ and $F_{\bm{n}}$ (coefficients of
$y^{\prime}$ and $y$, respectively) are

$\displaystyle E_{\bm{n}}$
$\displaystyle=-x(x-1)^{3}((1+\beta_{1})(1+\beta_{2})+x(n_{1}^{2}+n_{2}^{2}+n_{1}n_{2}+n_{1}(1+\beta_{1})+n_{2}(1+\beta_{2})-(2+\beta_{1})(2+\beta_{2}))),$
$\displaystyle F_{\bm{n}}$
$\displaystyle=-x(x-1)^{3}(n_{1}+n_{2})(1+n_{1}+\beta_{1})(1+n_{2}+\beta_{2}).$

Canceling the common factor $x(x-1)^{3}$ we obtain the differential equation

$\displaystyle
x^{2}(x-1)y^{\prime\prime\prime}+x(x(5+\beta_{1}+\beta_{2})-3-\beta_{1}-\beta_{2})y^{\prime\prime}$
$\displaystyle-(x(n_{1}^{2}+n_{2}^{2}+n_{1}n_{2}+n_{1}(1+\beta_{1})+n_{2}(1+\beta_{2})-(2+\beta_{1})(2+\beta_{2}))-(1+\beta_{1})(1+\beta_{2}))y^{\prime}$
$\displaystyle-(n_{1}+n_{2})(1+n_{1}+\beta_{1})(1+n_{2}+\beta_{2})y=0.$

This is the particular case (after canceling the common factor $x-1$) of the
equation derived in [7, Section 4.3].

As for the electrostatic partners, let us consider the example when
$\alpha=\beta_{1}=0$ and $\beta_{2}=-1/2$. Direct computation shows that the
moments of $w_{i}$ are

$\int_{0}^{1}x^{k}w_{j}(x)dx=\begin{cases}(k+1)^{-1},&j=1,\\\
2/(2j+1),&j=2,\end{cases}\quad k=0,1,\dots,$

which allows to find the asymptotic expansion of
$\mathfrak{C}_{w_{i}}[P_{\bm{n}}]$ at infinity, and in consequence,
$S_{\bm{n},i}$ (using formula (5.12)) by means of symbolic computation. For
instance, in the case $\bm{n}=(5,5)$ we obtain that

$\displaystyle\mathcal{L}_{(5,5)}(x)=$ $\displaystyle
x^{10}-\frac{380x^{9}}{87}+\frac{1615x^{8}}{203}-\frac{20672x^{7}}{2639}+\frac{9044x^{6}}{2001}-\frac{5168x^{5}}{3335}+\frac{204x^{4}}{667}$
$\displaystyle-\frac{64x^{3}}{2001}+\frac{x^{2}}{667}-\frac{4x}{182091}+\frac{1}{30045015},$

and up to normalization,

$\displaystyle S_{\bm{n},1}(x)=$ $\displaystyle\
882230895x^{5}+4709406975x^{4}+5720142090x^{3}+8795888965x^{2}$
$\displaystyle+11696347475x+11645469674,$ $\displaystyle S_{\bm{n},2}(x)=$
$\displaystyle\
5192762585x^{5}+313459871725x^{4}+662076961780x^{3}+782465377400x^{2}$
$\displaystyle+1267133219685x+1386883054197.$

All zeros of $P_{(5,5)}$ are positive and simple, both $S_{\bm{n},j}$ have one
negative and two pairs of complex conjugate simple zeros, all of them simple.

Zero asymptotics for sequences of Jacobi-Piñeiro polynomials $P_{(n,n)}$ as
$n\to\infty$ and $\alpha$, $\beta_{j}$’s fixed, was obtained in [31] and [77].
Again, the expression of the density, this time on $[0,1]$, was derived from
the recurrence relation satisfied by polynomials
$\mathfrak{L}_{(n_{1},n_{2})}$ and no equilibrium problem associated to that
distribution was given. For comparison, results of the numerical experiments
are presented in Figure 4, where we take $\bm{n}=(75,75)$, with $\beta_{1}=0$,
$\beta_{2}=-1/2$, and $\alpha=0$. According to our discussion in Section
5.2.2, the zeros of $S_{\bm{n},1}$ are real and negative, and the asymptotic
zero distribution can be described in terms of the Nikishin-type vector
equilibrium, see Section 6.4.

\begin{overpic}[scale={0.6}]{Fig4a} \end{overpic} | \begin{overpic}[scale={0.6}]{Fig4b} \end{overpic}
---|---
Figure 4. Left: zeros of the Jacobi-Piñeiro polynomial $P_{\bm{n}}$ (all on
$[0,1]$) and of $S_{\bm{n},1}$ (filled circles, all negative) for
$\bm{n}=(75,75)$, with $\beta_{1}=0$, $\beta_{2}=-1/2$, and $\alpha=0$;
approximately 19 real zeros of $S_{\bm{n},1}$ (ranging from $-771$ to
$-1.74$), are not represented. Right: the histogram of the zeros of
$P_{\bm{n}}$ and the plot of the asymptotic density, predicted in [31].

###  8.6. Angelesco-Jacobi polynomials

These polynomials, known also as Jacobi-Jacobi polynomials (see [12]) are
Hermite–Padé polynomials $P_{\bm{n}}$, $\bm{n}=(n_{1},n_{2})$, satisfying
orthogonality relations

$\begin{split}&\int_{a}^{0}x^{k}P_{\bm{n}}(x)(x-a)^{\alpha}|x|^{\beta}(1-x)^{\gamma}\,dx=0,\quad
k=0,1,2,\dots,n_{1}-1,\\\
&\int_{0}^{1}x^{k}P_{\bm{n}}(x)(x-a)^{\alpha}x^{\beta}(1-x)^{\gamma}x^{k}\,dx=0,\quad
k=0,1,2,\dots,n_{2}-1,\end{split}$

with $a<0$. In our notation,

$A(x)=x(x-a)(x-1),\quad B(x)=\alpha x(x-1)+\beta(x-a)(x-1)+\gamma x(x-a),$

and

$\displaystyle w_{1}(x)$
$\displaystyle=w_{2}(x)=w(x)=(x-a)^{\alpha}|x|^{\beta}(1-x)^{\gamma},$
$\displaystyle v_{1}(x)$
$\displaystyle=v_{2}(x)=v(x)=(x-a)^{\alpha+1}|x|^{\beta+1}(1-x)^{\gamma+1},$

with $\alpha,\beta,\gamma>-1$ and

$\Delta_{1}=[a,0],\quad\Delta_{2}=[0,1]$

(cf. Example 3.2).

This is an Angelesco system of semiclassical weights of class $\sigma=1$, see
(3.4). As it follows from [12, Theorem 2.1], for $n_{1}=n_{2}=n$,
$\bm{n}=(n,n)$, $n\in\mathbb{N}$, $P_{\bm{n}}$ can be expressed using a
Rodrigues formula,

$P_{\bm{n}}(x)=\frac{1}{w(x)}\left(\frac{d}{dx}\right)^{n}\left[A^{n}(x)w(x)\right],$

or explicitly, see [100, Section 3.5]: up to a constant factor,

$P_{\bm{n}}(x)=\sum_{k=0}^{n}\sum_{j=0}^{n-k}\frac{(-n)_{k+j}(-\alpha-n)_{j}(-\gamma-n)_{k}}{(\beta+1)_{k+j}k!j!}(x-a)^{n-k}(x-1)^{k+j}x^{n-j}.$

Kaliaguin [54] studied the case of $a=-1$, with the particular sub-case of
$B\equiv 0$ or $w(x)\equiv 1$ going back to the work of Appell [3]. The case
of $-1<a<0$ was addressed in [56], but see [12] for further historical
details.

If $B\equiv 0$ (that is, $\alpha=\beta=\gamma=0$), definition (3.13) reduces
to

$S_{\bm{n},1}=A\times\mathfrak{Wrons}\left[P_{\bm{n}},\mathfrak{C}_{w_{1}}[P_{\bm{n}}]\right],\quad\mathfrak{C}_{w_{1}}[P_{\bm{n}}](x)=\int_{a}^{0}\frac{P_{\bm{n}}(t)}{t-x}dt.$

Notice that in particular, $x=1$ is one of the $n+1$ zeros of $S_{\bm{n},1}$,
$n-1$ of which interlace with the zeros of $P_{\bm{n}}$ on $[0,1)$.

\begin{overpic}[scale={0.55}]{Fig5a} \end{overpic} | \begin{overpic}[scale={0.55}]{Fig5b} \end{overpic}
---|---
Figure 5. Appell’s polynomials ($\alpha=\beta=\gamma=0$). Left: graph of
$P_{\bm{n}}$ (dashed line) and $S_{\bm{n},1}$ (thick line) on $[0,1]$ for
$\bm{n}=(6,6)$ in the case. Right: zeros of $P_{\bm{n}}$ (empty circles, all
on $[-1,1]$) and of $S_{\bm{n},1}$ (filled circles, all on $[0,1]$) for
$\bm{n}=(15,15)$.

In Appell’s case, $w\equiv 1$ and $a=-1$,

$P_{\bm{n}}(x)=\left(\frac{d}{dx}\right)^{n}\left[x(x^{2}-1)\right]^{n},$

and

$S_{\bm{n},1}(x)=(-1)^{n+1}S_{\bm{n},2}(-x),\quad\deg S_{\bm{n},1}=n+1.$

These explicit formulas allow to use symbolic computation (e.g. Mathematica)
to find explicit expressions for small values of $n$. For instance, for
$\bm{n}=(6,6)$, and up to normalization,

$\displaystyle P_{\bm{n}}(x)$
$\displaystyle=x^{12}-\frac{44x^{10}}{17}+\frac{165x^{8}}{68}-\frac{220x^{6}}{221}+\frac{75x^{4}}{442}-\frac{2x^{2}}{221}+\frac{1}{18564},$
$\displaystyle S_{\bm{n},1}(x)$
$\displaystyle=(x-1)\left(12288x^{6}-38763x^{5}+47253x^{4}-27822x^{3}+8018x^{2}-991x+33\right),$

and their graphs are plotted in Figure 5. We can clearly observe the
interlacing predicted by Proposition 5.11, which in the limit $n\to\infty$
gives the description in therms of the Angelesco equilibrium problem, as
described in Section 6.3.

However, according to Remark 5.3, the critical configuration for the zeros of
$P_{\bm{n}}$ in this case is not unique: we can use for the second component
the zeros of any linear combination of $S_{\bm{n},1}$ and $S_{\bm{n},2}$. In
Figure 6 we illustrate the behavior of zeros of $S_{\bm{n},1}+tS_{\bm{n},2}$,
for different values of $0<t\leq 1$ (notice that the representation of the
zeros of $S_{\bm{n},1}$ in Figure 5, right, corresponds to $t=0$).

\begin{overpic}[scale={0.55}]{Fig6a} \end{overpic} | \begin{overpic}[scale={0.55}]{Fig6b} \end{overpic}
---|---
\begin{overpic}[scale={0.55}]{Fig6c} \end{overpic} | \begin{overpic}[scale={0.55}]{Fig6d} \end{overpic}
Figure 6. Appell’s polynomials ($\alpha=\beta=\gamma=0$) and $\bm{n}=(35,35)$:
zeros of $P_{\bm{n}}$ (indicated by empty circles, all on $(-1,1)$) and of
$S_{\bm{n},1}+tS_{\bm{n},2}$ (filled circles) for $t=10^{-10}$ (top left),
$t=10^{-5}$ (top right), $t=1$ (bottom left), and $t=10^{5}$ (bottom right).

###  8.7. Multiple orthogonal polynomials for the cubic weight

Although we have not discussed the purely complex weights, we finish our
presentation with the illustrative example of polynomials $P_{\bm{n}}$,
$\bm{n}=(n_{1},n_{2})$, satisfying orthogonality relations

$\begin{split}&\int_{\Delta_{1}}z^{k}P_{\bm{n}}(z)e^{-z^{3}}dz=0,\quad
k=0,1,2,\dots,n_{1}-1,\\\
&\int_{\Delta_{2}}z^{k}P_{\bm{n}}(z)e^{-z^{3}}dz=0,\quad
k=0,1,2,\dots,n_{2}-1,\end{split}$

where $\Delta_{1}$ and $\Delta_{2}$ are contours on the complex plane,
extending to $\infty$ on their two ends along the directions determined by the
angles $-2\pi/3$ and $0$, and $-2\pi/3$ and $2\pi/3$, respectively. They were
introduced in [101] and studied in full generality in [74].

In our notation,

$A_{1}(x)=A_{2}(x)=A(x)=1,\quad B_{1}(x)=B_{2}(x)=B(x)=-3x^{2}.$

Detailed explanation of the algorithm of computation of the zeros of
$P_{\bm{n}}$ was given in [74, Section 9]. Since the moments of the weight
were given explicitly, we use the expansion of $\mathfrak{D}_{w}[P_{\bm{n}}]$
at infinity to calculate the expressions for $S_{\bm{n},1}$ and
$S_{\bm{n},2}$, see Figure 7.

It is interesting to compare the location of the zeros of $S_{\bm{n},1}$ with
those of type I multiple orthogonal polynomials $\mathfrak{A}_{\bm{n}}$ and
$\mathfrak{B}_{\bm{n}}$, defined by the following conditions:

$\deg\mathfrak{A}_{\bm{n}}\leq n-1,\quad\deg\mathfrak{B}_{\bm{n}}\leq m-1,$

and

$\displaystyle\int_{\Delta_{1}}z^{k}\mathfrak{A}_{\bm{n}}(z)e^{-z^{3}}dz+\int_{\Delta_{2}}z^{k}\mathfrak{B}_{\bm{n}}(z)e^{-z^{3}}dz=0,$
$\displaystyle\quad k=0,\dots,N-2,$ (8.6)
$\displaystyle\int_{\Delta_{1}}z^{k}\mathfrak{A}_{\bm{n}}(z)e^{-z^{3}}dz+\int_{\Delta_{2}}z^{k}\mathfrak{B}_{\bm{n}}(z)e^{-z^{3}}dz=1,$
$\displaystyle\quad k=N-1,$

where $N=n+m$. According to Figure 7, the zeros of $\mathfrak{B}_{\bm{n}}$
“interlace” with the zeros of $S_{\bm{n},1}$, which brings up a natural
question of a possible connection of these two polynomials.

\begin{overpic}[scale={0.7}]{Fig7} \end{overpic} Figure 7. Zeros of the
multiple orthogonal polynomial $P_{\bm{n}}$ with respect to the cubic weight
(indicated by empty circles, part of them on the positive semiaxis, forming a
symmetric star) and of $S_{\bm{n},1}$ (filled circles) for $\bm{n}=(25,25)$.
For comparison, zeros of type I MOP $\mathfrak{B}_{\bm{n}}$ (filled squares)
are also represented.

##  Appendix A Properties of the electrostatic partner

######  Proposition A.1.

If $P$ is a polynomial of degree $N\in\mathbb{N}$, $q$ the function of the
second kind (3.9), and $S$ is the electrostatic partner of $P$ defined in
(3.21) then:

  1. a)

if $z_{0}\in\mathbb{C}$ is a zero of $P$ of multiplicity $k\geq 1$ then $S$
also has a root at $z_{0}$ with multiplicity at least $k-1$. If in addition
$A(z_{0})=0$ then the multiplicity of $z_{0}$ in $S$ is at least $k$.

  2. b)

if $z_{0}\in\mathbb{C}\setminus\Delta$ is a zero of $\mathfrak{C}_{w}[P]$ of
multiplicity $k\geq 1$ then $S$ also has a root at $z_{0}$ with multiplicity
at least $k-1$. If in addition $A(z_{0})=0$ then the multiplicity of $z_{0}$
in $S$ is at least $k$.

  3. c)

Let $\Omega$ be a simply-connected domain such that $A(z)\neq 0$ for
$z\in\Omega$, $\mathfrak{C}_{w}[P]$ holomorphic in $\Omega$ and let $w$ be a
holomorphic branch of this function in $\Omega$. If $z_{0}$ is a boundary
point of $\Omega$, then

$\lim_{z\to z_{0}}(z-z_{0})A^{k}(z)w(z)\frac{d^{k}}{dz^{k}}q(z)=0,\qquad
k=0,1,\dots,$ (A.1)

where we take non-tangential limit with $z\in\Omega$.

###### Proof.

Assume $z_{0}\notin\Delta$. If $z_{0}\in\mathbb{C}$ is a zero of $p$ of
multiplicity $k\geq 1$ then the assertion in a) follows directly from the
definition (3.13). Same argument works to prove a) if we assume that
$z_{0}\in\mathbb{C}$ is a zero of $\mathfrak{C}_{w}[P]$ of multiplicity $k\geq
1$.

On the other hand, from the assumption (3.2) it follows that for
$z_{0}\in\Delta$,

$\lim_{z\to z_{0}}(z-z_{0})\,\mathfrak{C}_{w}[P](z)=0,$ (A.2)

which proves b) also in this case.

We turn now to c). If $A(z_{0})\neq 0$ then both $w^{-1}$ and $q$ are bounded
at $z_{0}$ and the assertion is obvious. Hence, let $A(z_{0})=0$.

Denote

$h_{k}(z):=(z-z_{0})A^{k}(z)w(z)\frac{d^{k}}{dz^{k}}q(z),\quad k=0,1,2,\dots$

Obviously, $h_{0}(z)=(z-z_{0})\,\widehat{p}(z)$, and (A.1) for $k=0$ is
consequence of (A.2).

Using the induction in $k$, assume that (A.1) is establised for a certain
$k\geq 0$, i.e.

$\lim_{z\to z_{0}}h_{k}(z)=0,$

where we always take non-tangential limit from $\Omega$. By Lemma A.2 below,

$\lim_{z\to z_{0}}(z-z_{0})h^{\prime}_{k}(z)=0.$ (A.3)

Thus, using (3.3),

$\displaystyle 0$ $\displaystyle=\lim_{z\to
z_{0}}(z-z_{0})h^{\prime}_{k}(z)=\lim_{z\to z_{0}}A(z)h^{\prime}_{k}(z)$
$\displaystyle=\lim_{z\to
z_{0}}\left[h_{k}(z)\left(\frac{A(z)}{z-z_{0}}+kA^{\prime}(z)+A(z)\frac{w^{\prime}(z)}{w(z)}\right)+h_{k+1}(z)\right]$
$\displaystyle=\lim_{z\to
z_{0}}\left[h_{k}(z)\left(\frac{A(z)}{z-z_{0}}+kA^{\prime}(z)+B(z)\right)+h_{k+1}(z)\right]=\lim_{z\to
z_{0}}h_{k+1}(z),$

which proves (A.1). ∎

######  Lemma A.2.

Let $\Omega$ be a domain, $z_{0}\in\partial\Omega$ a boundary point satisfying
the following property: there exists a sector
$\Gamma:=\\{z\in\mathbb{C}:\,|\arg(z-z_{0})-\theta_{0}|<\delta\\}$, with
certain $\theta_{0}\in[0,2\pi)$, $\delta>0$, such that for a sufficiently
small $r>0$, $\Gamma\cap\\{z\in\mathbb{C}:\,|z-z_{0}|<r\\}\subset\Omega$. If
function $f$ is holomorphic in $\Omega$ and

$\lim_{z\to z_{0},\,z\in\Gamma}\frac{f(z)}{z-z_{0}}=0,$

then

$\lim_{z\to z_{0},\,z\in\Gamma}f^{\prime}(z)=0.$

###### Proof.

Without loss of generality we can assume that $z_{0}=0$ and $\theta_{0}=0$.

Each $0<q<1$ defines a sub-sector $\Gamma_{q}$ of $\Gamma$ given by

$\Gamma_{q}:=\\{z\in\Gamma:\,\left|\tan\left(\delta-\arg(z)\right)\right|\geq
q\\}.$

Notice that as $q\to 0$, $\Gamma_{q}$ exhausts $\Gamma$.

Fix a $0<q<1$ and an arbitrarily small $\varepsilon>0$. By assumptions, there
exists $0<r^{\prime}=r^{\prime}(\varepsilon)<r$ such that

$t\in\Gamma,\,|t|<r^{\prime}\quad\Rightarrow\quad|f(t)|\leq\varepsilon|t|.$

Let $z\in\Gamma_{q}$; by construction, the circle
$C_{z}:=\\{t\in\mathbb{C}:\,|t-z|=q|z|\\}\subset\Gamma$. We assume $|z|$ small
enough so that $|t|<r^{\prime}$ for all $t\in C_{z}$. Then

$f^{\prime}(z)=\frac{1}{2\pi i}\oint_{C_{z}}\frac{f(t)}{(t-z)^{2}}dt,$

and

$|f^{\prime}(z)|\leq\frac{1}{q|z|}\max_{t\in
C_{z}}|f(t)|\leq\frac{\varepsilon}{q|z|}\max_{t\in
C_{z}}|t|=\frac{\varepsilon}{q}\frac{|z|+q|z|}{|z|}=\frac{\varepsilon(1+q)}{q},$

which proves the assertion. ∎

##  Appendix B Electrostatic partner in the real case

The construction in Section 3 is carried out in a very general setup. In this
Appendix, we prove a technical result valid in the real case, when $A$ and $B$
have real coefficients,

$\Delta\subset\mathbb{R}\quad\text{ and }\quad w(x)\geq 0\text{ on }\Delta.$
(B.1)

We assume that $P\not\equiv 0$ is a polynomial with real coefficients, and as
before, denote $S=\mathfrak{D}_{w}[P]$.

######  Lemma B.1.

Let $\zeta_{1}<\zeta_{2}$ be two consecutive real zeros of $P$, such that
$(\zeta_{1},\zeta_{2})\cap\Delta=\emptyset$, $S(\zeta_{j})\neq 0$ for $j=1,2$,
and $A$ preserves sign on $(\zeta_{1},\zeta_{2})$. Then
$S\,\mathfrak{C}_{w}[P]$ does change sign in $(\zeta_{1},\zeta_{2})$.

Analogously, if $y_{1}<y_{2}$ are two consecutive real zeros of
$\mathfrak{C}_{w}[P]$, such that $(y_{1},y_{2})\cap\Delta=\emptyset$,
$(AS)(y_{j})\neq 0$ for $j=1,2$, and $A$ preserves sign on $(y_{1},y_{2})$,
then $SP$ changes sign in $(y_{1},y_{2})$.

###### Proof.

Since $S(\zeta_{j})\neq 0$, $j=1,2$, it follows from b) in Proposition A.1
that these are simple zeros of $P$, $A(\zeta_{j})\neq 0$, $j=1,2$, and thus,
$P^{\prime}(\zeta_{1})P^{\prime}(\zeta_{2})<0$. Evaluating in the definition
(3.13), we get

$S(\zeta_{j})=\mathfrak{D}_{w}[P](\zeta_{j})=-A(\zeta_{j})\mathfrak{C}_{w}[P](\zeta_{j})P^{\prime}(\zeta_{j}),\quad
j=1,2,$

so that

$\left(S\,\mathfrak{C}_{w}[P]\right)(\zeta_{j})=-\left(A(\mathfrak{C}_{w}[P])^{2}P^{\prime}\right)(\zeta_{j}),\quad
j=1,2.$

Since $A$ preserves sign on $[\zeta_{1},\zeta_{2}]$, we get that

$\left(S\mathfrak{C}_{w}[P]\right)(\zeta_{1})\left(S\mathfrak{C}_{w}[P]\right)(\zeta_{2})<0.$

This proves the first assertion.

Similarly from (3.13),

$S(y_{j})=\left(A\left(\mathfrak{C}_{w}[P]\right)^{\prime}P\right)(y_{j}),\quad
j=1,2,$

so that

$(SP)(y_{j})=\left(A\left(\mathfrak{C}_{w}[P]\right)^{\prime}P^{2}\right)(y_{j}),\quad
j=1,2.$

If $y_{1}<y_{2}$ do not coincide with the zeros of $AS$, then as (3.22) shows,
these are simple zeros of $\widehat{p}$, so that

$\left(\mathfrak{C}_{w}[P]\right)^{\prime}(y_{1})\left(\mathfrak{C}_{w}[P]\right)^{\prime}(y_{2})<0,$

and we conclude that

$\left(SP\right)(y_{1})\left(SP\right)(y_{2})<0.$

∎

## Acknowledgments

The first author was partially supported by Simons Foundation Collaboration
Grants for Mathematicians (grant 710499). He also acknowledges the support of
the Spanish Government and the European Regional Development Fund (ERDF)
through grant PID2021-124472NB-I00, Junta de Andalucía (research group FQM-229
and Instituto Interuniversitario Carlos I de Física Teórica y Computacional),
and by the University of Almería (Campus de Excelencia Internacional del Mar
CEIMAR) in the early stages of this project.

The second and third authors were partially supported by Spanish Ministerio de
Ciencia, Innovación y Universidades, under grant MTM2015-71352-P.

The third author was additionally supported by Junta de Andalucía (research
group FQM-384), the University of Granada (Research Project ERDF-UGR
A-FQM-246-UGR20), and by the IMAG–Maria de Maeztu grant
CEX2020-001105-M/AEI/10.13039/501100011033.

The authors also grateful to Alexandre Eremenko, who suggested the idea of the
proof of Lemma A.2, and to the anonymous referees, whose outstanding work
helped to considerably improve this manuscript.

## References

  * [1] M. Adler and J. Moser.  On a class of polynomials connected with the Korteweg-de Vries equation.  Comm. Math. Phys., 61(1):1–30, 1978.
  * [2] M. A. Angelesco.  Sur deux extensions des fractions continues algébriques.  C.R. Acad. Sci. Paris, (18):262–263, 1919.
  * [3] P. Appell.  Sur une suite des polynômes ayant toutes leurs racines réelles.  Arch. Math. Phys., 1(3):69–71, 1901.
  * [4] A. I. Aptekarev.  Strong asymptotics of polynomials of simultaneous orthogonality for Nikishin systems.  Mat. Sb., 190(5):3–44, 1999.
  * [5] A. I. Aptekarev.  Asymptotics of Hermite-Padé approximants for a pair of functions with branch points.  Dokl. Akad. Nauk, 422(4):443–445, 2008.  English transl. in Dokl. Math. 78:2 (2008), 717-719.
  * [6] A. I. Aptekarev, P. M. Bleher, and A. B. J. Kuijlaars.  Large $n$ limit of Gaussian random matrices with external source. II.  Comm. Math. Phys., 259(2):367–389, 2005.
  * [7] A. I. Aptekarev, A. Branquinho, and W. Van Assche.  Multiple orthogonal polynomials for classical weights.  Transactions AMS, 355:3887–3914, 2003.
  * [8] A. I. Aptekarev and V. A. Kalyagin.  Analytic properties of functions representable by $P^{(m)}$-fractions with periodic coefficients.  Akad. Nauk SSSR Inst. Prikl. Mat. Preprint, (57):12, 1986.
  * [9] A. I. Aptekarev and A. B. J. Kuijlaars.  Hermite-Padé approximations and multiple orthogonal polynomial ensembles.  Russ. Math. Surv. 66, No. 6, 1133–1199 (2011); translation from Usp. Mat. Nauk 66, No. 6, 123–190 (2011).
  * [10] A. I. Aptekarev and V. G. Lysov.  Systems of Markov functions generated by graphs and the asymptotics of their Hermite-Padé approximants.  Mat. Sb., 201(2):29–78, 2010.  English trans. in Sb. Math. 201:2 (2010), 183–234.
  * [11] A. I. Aptekarev, V. G. Lysov, and D. N. Tulyakov.  Random matrices with an external source and the asymptotics of multiple orthogonal polynomials.  Mat. Sb., 202(2):3–56, 2011.
  * [12] A. I. Aptekarev, F. Marcellán, and I. A. Rocha.  Semiclassical multiple orthogonal polynomials and the properties of Jacobi-Bessel polynomials.  J. Approx. Theory, 90:117–146, 1997.
  * [13] A. I. Aptekarev and H. Stahl.  Asymptotics of Hermite-Padé polynomials.  In Progress in approximation theory (Tampa, FL, 1990), volume 19 of Springer Ser. Comput. Math., pages 127–167. Springer, New York, 1992.
  * [14] A. I. Aptekarev, W. Van Assche, and M. L. Yattselev.  Hermite-Padé approximants for a pair of Cauchy transforms with overlapping symmetric supports.  Comm. Pure Appl. Math., 70(3):444–510, 2017.
  * [15] H. Aref.  Vortices and polynomials.  Fluid Dynam. Res., 39(1-3):5–23, 2007.
  * [16] A. M. Barry, F. Hajir, and P. G. Kevrekidis.  Generating functions, polynomials and vortices with alternating signs in Bose-Einstein condensates.  J. Phys. A, 48(15):155205, 18, 2015.
  * [17] A. B. Bartman.  A new interpretation of the Adler-Moser KdV polynomials: interaction of vortices.  In Nonlinear and turbulent processes in physics, Vol. 3 (Kiev, 1983), pages 1175–1181. Harwood Academic Publ., Chur, 1984.
  * [18] M. Bertola, E. Chavez-Heredia, and T. Grava.  The Stieltjes–Fekete problem and degenerate orthogonal polynomials.  Preprint arXiv:2206.06861, 2022.
  * [19] P. Bleher, S. Delvaux, and A. B. J. Kuijlaars.  Random matrix model with external source and a constrained vector equilibrium problem.  Comm. Pure Appl. Math., 64(1):116–160, 2011.
  * [20] P. Bleher and A. B. J. Kuijlaars.  Large $n$ limit of Gaussian random matrices with external source. I.  Comm. Math. Phys., 252(1-3):43–76, 2004.
  * [21] P. M. Bleher and A. B. J. Kuijlaars.  Integral representations for multiple Hermite and multiple Laguerre polynomials.  Ann. Inst. Fourier (Grenoble), 55(6):2001–2014, 2005.
  * [22] P. M. Bleher and A. B. J. Kuijlaars.  Large $n$ limit of Gaussian random matrices with external source. III. Double scaling limit.  Comm. Math. Phys., 270(2):481–517, 2007.
  * [23] M. Bôcher.  The roots of polynomials that satisfy certain differential equations of the second order.  Bull. Amer. Math. Soc., 4:256–258, 1897.
  * [24] J. W. Brown.  On Angelesco-type polynomials.  Ricerca (Napoli) (3), 24(maggio-agosto):3–7, 1973.
  * [25] J. L. Burchnall and T. W. Chaundy.  A Set of Differential Equations which can be Solved by Polynomials.  Proc. London Math. Soc. (2), 30(6):401–414, 1930.
  * [26] J. Bustamante and G. López Lagomasino.  Hermite–Padé approximation for Nikishin systems of analytic functions.  Russian Acad. Sci. Sb. Math., 77:367–384, 1994.
  * [27] K. Castillo, M. N. de Jesus, and J. Petronilho.  An electrostatic interpretation of the zeros of sieved ultraspherical polynomials.  J. Math. Phys., 61(5):053501, 19, 2020.
  * [28] T. S. Chihara.  On quasi-orthogonal polynomials.  Proc. Amer. Math. Soc., 8:765–767, 1957.
  * [29] T. S. Chihara.  An Introduction to Orthogonal Polynomials.  Gordon and Breach, New York, 1978.
  * [30] P. A. Clarkson.  Vortices and polynomials.  Stud. Appl. Math., 123(1):37–62, 2009.
  * [31] E. Coussement, J. Coussement, and W. Van Assche.  Asymptotic zero distribution for a class of multiple orthogonal polynomials.  Trans. Am. Math. Soc., 360(10):5571–5588, 2008.
  * [32] J. Coussement and W. Van Assche.  Differential equations for multiple orthogonal polynomials with respect to classical weights: raising and lowering operators.  J. Phys. A, 39(13):3311–3318, 2006.
  * [33] P. A. Deift.  Orthogonal polynomials and random matrices: a Riemann-Hilbert approach.  New York University Courant Institute of Mathematical Sciences, New York, 1999.
  * [34] M. V. Demina and N. A. Kudryashov.  Special polynomials and rational solutions of the hierarchy of the second Painlevé equation.  Teoret. Mat. Fiz., 153(1):58–67, 2007.
  * [35] M. V. Demina and N. A. Kudryashov.  Point vortices and classical orthogonal polynomials.  Regul. Chaotic Dyn., 17(5):371–384, 2012.
  * [36] M. V. Demina and N. A. Kudryashov.  Vortices and polynomials: non-uniqueness of the Adler-Moser polynomials for the Tkachenko equation.  J. Phys. A, 45(19):195205, 12, 2012.
  * [37] M. V. Demina and N. A. Kudryashov.  Multi-particle dynamical systems and polynomials.  Regul. Chaotic Dyn., 21(3):351–366, 2016.
  * [38] M. V. Dëmina, N. A. Kudryashov, and D. I. Sinel′shchikov.  The polygonal method for constructing exact solutions of some nonlinear differential equations to describe waves on water.  Zh. Vychisl. Mat. Mat. Fiz., 48(12):2151–2162, 2008.
  * [39] D. K. Dimitrov and B. Shapiro.  Electrostatic problems with a rational constraint and degenerate Lamé equations.  Potential Anal., 52(4):645–659, 2020.
  * [40] D. K. Dimitrov and W. Van Assche.  Lamé differential equations and electrostatics.  Proc. Amer. Math. Soc., 128(12):3621–3628, 2000.  Erratum in Proc. Amer. Math. Soc. 131 (2003), no. 7, 2303.
  * [41] K. Driver and H. Stahl.  Normality in Nikishin systems.  Indag. Math., 5:161ñ–187, 1994.
  * [42] U. Fidalgo Prieto and G. López Lagomasino.  Nikishin systems are perfect.  Constr. Approx., 34:297–356, 2011.
  * [43] G. Filipuk, W. Van Assche, and L. Zhang.  Ladder operators and differential equations for multiple orthogonal polynomials.  J. Phys. A, 46(20):205204, 24, 2013.
  * [44] P. J. Forrester and J. B. Rogers.  Electrostatics and the zeros of the classical polynomials.  SIAM J. Math. Anal., 17(2):461–468, 1986.
  * [45] F. D. Gakhov.  Boundary value problems.  Dover Publications Inc., New York, 1990.  Translated from the Russian, Reprint of the 1966 translation.
  * [46] A. A. Gonchar, E. A. Rakhmanov, and V. N. Sorokin.  Hermite-Padé approximants for systems of Markov-type functions.  Sb. Math., 188:671–696, 1997.
  * [47] F. A. Grünbaum.  Variations on a theme Heine and Stieltjes: An electrostatic interpretation of the zeros of certain polynomials.  99, pages 189–194, 1998.
  * [48] E. Heine.  Handbuch der Kugelfunctionen, volume II.  G. Reimer, Berlin, $2$nd. edition, 1878.
  * [49] M. E. H. Ismail.  An electrostatic model for zeros of general orthogonal polynomials.  Pacific J. Math., 193:355–369, 2000.
  * [50] M. E. H. Ismail.  More on electrostatic models for zeros of orthogonal polynomials.  In Proceedings of the International Conference on Fourier Analysis and Applications (Kuwait, 1998), volume 21, pages 191–204, 2000.
  * [51] M. E. H. Ismail.  Functional equations and electrostatic models for orthogonal polynomials.  In Random matrix models and their applications, volume 40 of Math. Sci. Res. Inst. Publ., pages 225–244. Cambridge Univ. Press, Cambridge, 2001.
  * [52] M. E. H. Ismail.  Classical and quantum orthogonal polynomials in one variable, volume 98 of Encyclopedia of Mathematics and its Applications.  Cambridge University Press, Cambridge, 2005.  With two chapters by Walter Van Assche, With a foreword by Richard A. Askey.
  * [53] M. E. H. Ismail and X.-S. Wang.  On quasi-orthogonal polynomials: their differential equations, discriminants and electrostatics.  J. Math. Anal. Appl., 474(2):1178–1197, 2019.
  * [54] V. A. Kaliaguine.  On a class of polynomials defined by two orthogonality relations.  Math. USSR Sbornik, 38(4):563–580, 1981.
  * [55] V. A. Kaliaguine.  On operators associated with Angelesco systems.  East J. Approx., 1(2):157–170, 1995.
  * [56] V. A. Kaliaguine and A. Ronveaux.  On a system of “classical” polynomials of simultaneous orthogonality.  J. Comput. Appl. Math., 67:207–217, 1996.
  * [57] N. A. Kudryashov and M. V. Demina.  Relations between zeros of special polynomials associated with the Painlevé equations.  Phys. Lett. A, 368(3-4):227–234, 2007.
  * [58] A. B. J. Kuijlaars, A. Martínez-Finkelshtein, and R. Orive.  Orthogonality of Jacobi polynomials with general parameters.  Electron. Trans. Numer. Anal., 19:1–17 (electronic), 2005.
  * [59] A. B. J. Kuijlaars, A. Martínez-Finkelshtein, and F. Wielonsky.  Non-intersecting squared Bessel paths and multiple orthogonal polynomials for modified Bessel weights.  Comm. Math. Phys., 286(1):217–275, 2009.
  * [60] G. López-Lagomasino.  An introduction to multiple orthogonal polynomials and Hermite-Padé approximation.  In Orthogonal polynomials: current trends and applications, volume 22 of SEMA SIMAI Springer Ser., pages 237–271. Springer, Cham, 2021\.
  * [61] G. López Lagomasino and W. Van Assche.  The Riemann–Hilbert analysis for a Nikishin system.  Sbornik: Mathematics, 209(7):1019–1050, 2018.
  * [62] I. Loutsenko.  Equilibrium of charges and differential equations solved by polynomials.  J. Phys. A, 37(4):1309–1321, 2004.
  * [63] V. Lysov and F. Wielonsky.  Strong asymptotics for multiple Laguerre polynomials.  Constr. Approx., 28(1):61–111, 2008.
  * [64] V. G. Lysov.  Mixed type Hermite-Padé approximants for a Nikishin system.  Proc. Steklov Inst. Math., 311:199–213, 2020.
  * [65] A. P. Magnus.  Painlevé–type differential equations for the recurrence coefficients of semi–classical orthogonal polynomials.  J. Comp. Appl. Math., 57:215–237, 1995.
  * [66] K. Mahler.  Perfect systems.  Compositio Math., 19:95–166, 1968.
  * [67] F. Marcellán, A. Martínez-Finkelshtein, and P. Martínez-González.  Electrostatic models for zeros of polynomials: old, new, and some open problems.  J. Comput. Appl. Math., 207(2):258–272, 2007.
  * [68] F. Marcellán and I. A. Rocha.  On semiclassical linear functionals: integral representations.  In Proceedings of the Fourth International Symposium on Orthogonal Polynomials and their Applications (Evian-Les-Bains, 1992), volume 57, pages 239–249, 1995.
  * [69] F. Marcellán and I. A. Rocha.  Complex path integral representation for semiclassical linear functionals.  J. Approx. Theory, 94(1):107–127, 1998.
  * [70] A. Martínez-Finkelshtein and E. A. Rakhmanov.  On asymptotic behavior of Heine-Stieltjes and Van Vleck polynomials.  In Recent trends in orthogonal polynomials and approximation theory, volume 507 of Contemp. Math., pages 209–232. Amer. Math. Soc., Providence, RI, 2010.
  * [71] A. Martínez-Finkelshtein and E. A. Rakhmanov.  Critical measures, quadratic differentials, and weak limits of zeros of Stieltjes polynomials.  Comm. Math. Phys., 302(1):53–111, 2011.
  * [72] A. Martínez-Finkelshtein and E. B. Saff.  Asymptotic properties of Heine-Stieltjes and Van Vleck polynomials.  J. Approx. Theory, 118(1):131–151, 2002.
  * [73] A. Martínez-Finkelshtein and G. L. F. Silva.  Critical measures for vector energy: global structure of trajectories of quadratic differentials.  Adv. Math., 302:1137–1232, 2016.
  * [74] A. Martínez-Finkelshtein and G. L. F. Silva.  Critical measures for vector energy: asymptotics of non-diagonal multiple orthogonal polynomials for a cubic weight.  Adv. Math., 349:246–315, 2019.
  * [75] A. Martínez-Finkelshtein and G. L. F. Silva.  Spectral curves, variational problems and the hermitian matrix model with external source.  Comm. Math. Phys., 383:2163–2242, 2021.
  * [76] E. Mukhin and A. Varchenko.  Multiple orthogonal polynomials and a counterexample to the Gaudin Bethe ansatz conjecture.  Trans. Amer. Math. Soc., 359(11):5383–5418, 2007.
  * [77] T. Neuschel and W. Van Assche.  Asymptotic zero distribution of Jacobi-Piñeiro and multiple Laguerre polynomials.  J. Approx. Theory, 205:114–132, 2016.
  * [78] P. K. Newton.  The $N$-vortex problem, volume 145 of Applied Mathematical Sciences.  Springer-Verlag, New York, 2001.  Analytical techniques.
  * [79] E. M. Nikishin.  On simultaneous Padé approximants.  Math. USSR. Sb., 41:409–425, 1982.
  * [80] E. M. Nikishin.  Asymptotic behavior of linear forms for simultaneous Padé approximants.  Izv. Vyssh. Uchebn. Zaved. Mat., pages 33–41, 84, 1986.
  * [81] E. M. Nikishin and V. N. Sorokin.  Rational approximations and orthogonality.  American Mathematical Society, Providence, RI, 1991.
  * [82] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors.  NIST handbook of mathematical functions.  U.S. Department of Commerce, National Institute of Standards and Technology, Washington, DC; Cambridge University Press, Cambridge, 2010.  With 1 CD-ROM (Windows, Macintosh and UNIX).
  * [83] L. R. Pineiro Diaz.  On simultaneous approximations for some collection of Markov functions.  Mosc. Univ. Math. Bull. 42, No. 2, 52–55 (1987); translation from Vestn. Mosk. Univ., Ser. I 1987, No. 2, 67–70 (1987).
  * [84] E. A. Rakhmanov.  On the asymptotics of Hermite-Padé polynomials for two Markov functions.  Mat. Sb., 202(1):133–140, 2011.
  * [85] E. A. Rakhmanov and S. P. Suetin.  Asymptotic behavior of Hermite-Padé polynomials of the first kind for a pair of functions forming a Nikishin system.  Uspekhi Mat. Nauk, 67(5(407)):177–178, 2012.
  * [86] E. A. Rakhmanov and S. P. Suetin.  Distribution of zeros of Hermite-Padé polynomials for a pair of functions forming a Nikishin system.  Mat. Sb., 204(9):115–160, 2013.
  * [87] A. Ronveaux, editor.  Heun’s differential equations.  The Clarendon Press Oxford University Press, New York, 1995.  With contributions by F. M. Arscott, S. Yu. Slavyanov, D. Schmidt, G. Wolf, P. Maroni and A. Duval.
  * [88] P. Rutka and R. Smarzewski.  Complete solution of the electrostatic equilibrium problem for classical weights.  Appl. Math. Comput., 218(10):6027–6037, 2012.
  * [89] B. Shapiro.  Algebro-geometric aspects of Heine-Stieltjes theory.  J. Lond. Math. Soc. (2), 83(1):36–56, 2011.
  * [90] J. A. Shohat.  A differential equation for orthogonal polynomials.  Duke Math. J., 5:401–417, 1939.
  * [91] B. Simanek.  An electrostatic interpretation of the zeros of paraorthogonal polynomials on the unit circle.  SIAM J. Math. Anal., 48(3):2250–2268, 2016.
  * [92] H. Stahl and V. Totik.  General orthogonal polynomials, volume 43 of Encyclopedia of Mathematics and its Applications.  Cambridge University Press, Cambridge, 1992.
  * [93] S. Steinerberger.  Electrostatic interpretation of zeros of orthogonal polynomials.  Proc. Am. Math. Soc., 146(12):5323–5331, 2018.
  * [94] T. J. Stieltjes.  Sur certains polynômes que vérifient une équation différentielle linéaire du second ordre et sur la teorie des fonctions de Lamé.  Acta Math., 6:321–326, 1885.
  * [95] T. J. Stieltjes.  Sur les polynômes de Jacobi.  C. R. Acad. Sci., Paris, 100:439–440, 1885.
  * [96] T. J. Stieltjes.  Un théorème d’algèbre.  Acta Math., 6(1):319–320, 1885.  Extrait d’une lettre adressée à M. Hermite.
  * [97] G. Szegő.  Orthogonal polynomials.  Amer. Math. Soc., Colloq. Publ., 1975.
  * [98] G. Valent and W. Van Assche.  The impact of Stieltjes’ work on continued fractions and orthogonal polynomials: additional material.  In Proceedings of the International Conference on Orthogonality, Moment Problems and Continued Fractions (Delft, 1994), volume 65, pages 419–447, 1995.
  * [99] W. Van Assche.  Padé and Hermite-Padé approximation and orthogonality.  Surv. Approx. Theory, 2:61–91, 2006.
  * [100] W. Van Assche and E. Coussement.  Some classical multiple orthogonal polynomials.  volume 127, pages 317–347. 2001.  Numerical analysis 2000, Vol. V, Quadrature and orthogonal polynomials.
  * [101] W. Van Assche, G. Filipuk, and L. Zhang.  Multiple orthogonal polynomials associated with an exponential cubic weight.  J. Approx. Theory, 190:1–25, 2015.
  * [102] W. Van Assche and A. Vuerinckx.  Multiple Hermite polynomials and simultaneous Gaussian quadrature.  Electron. Trans. Numer. Anal., 50:182–198, 2018.
  * [103] J. F. van Diejen.  Gradient system for the roots of the Askey-Wilson polynomial.  Proc. Am. Math. Soc., 147(12):5239–5249, 2019.
  * [104] J. F. van Diejen.  Stable equilibria for the roots of the symmetric continuous Hahn and Wilson polynomials.  In Orthogonal polynomials: current trends and applications. Proceedings of the 7th EIBPOA conference, Universidad Carlos III de Madrid, Leganés, Spain, July 3–6, 2018, pages 171–192. Cham: Springer, 2021.
  * [105] E. B. Van Vleck.  On the polynomials of Stieltjes.  Bull. Amer. Math. Soc., 4:426–438, 1898.
  * [106] H. von Helmholtz.  Über Integrale der hydrodynamischen Gleichungen, welche den Wirbelbewegungen entsprechen.  J. Reine Angew. Math., 55:25–55, 1858.  English translation by Tait, P.G., 1867. On integrals of the hydrodynamical equations, which express vortex-motion. Philos. Mag. 33(4), 485–512.