Differential equations for the recurrence coefficients limits for multiple orthogonal polynomials from a Nevai class

TL;DR

This paper investigates the limiting behavior of recurrence coefficients for multiple orthogonal polynomials in a Nevai class, describing their limits via differential equations and illustrating results numerically for Angelesco systems.

Contribution

It introduces a novel approach to characterize the limits of recurrence coefficients using differential equations for multiple orthogonal polynomials in a Nevai class.

Findings

Limits of recurrence coefficients are described by PDEs or ODEs.

Numerical illustrations are provided for Angelesco systems.

The approach advances understanding of asymptotic properties of multiple orthogonal polynomials.

Abstract

A limiting property of the nearest-neighbor recurrence coefficients for multiple orthogonal polynomials from a Nevai class is investigated. Namely, assuming that the nearest-neighbor coefficients have a limit along rays of the lattice, we describe it in terms of the solution of a system of partial differential equations. In the case of two orthogonality measures the differential equation becomes ordinary. For Angelesco systems, the result is illustrated numerically.

Equations

$$x{P}_n(x)=P_{n+1}(x)+b_n{P}_n(x)+a_{n-1}{P}_{n-1}(x)$$

$${\mathcal{J}}=\left(\begin{array}[]{cccc}b_{0}&\sqrt{a_{0}}&{0}&\\ \sqrt{a_{0}}&b_{1}&\sqrt{a_{1}}&\ddots\\ {0}&\sqrt{a_{1}}&b_{2}&\ddots\\ &\ddots&\ddots&\ddots\end{array}\right).$$

$$\int P_{\vec{n}}(x)x^{k}d{\mu }_i=0,k\in \{0,\dots ,n_i-1\},i\in \{1,\dots ,d\}.$$

$$z{P}_{\vec{n}}(z)=P_{\vec{n}+{\vec{e}}_j}(z)+b_{\vec{n},j}{P}_{\vec{n}}(z)+\sum _{i=1}^d{a}_{\vec{n},i}{P}_{\vec{n}-{\vec{e}}_i}(z),$$

$$b_{\vec{n}+{\vec{e}}_j,i}-b_{\vec{n},i}=b_{\vec{n}+{\vec{e}}_i,j}-b_{\vec{n},j},i<j,$$

$$\det \left(\begin{array}{cc}{b}_{\vec{n}+{\vec{e}}_j,i}& b_{\vec{n},i}\\ b_{\vec{n}+{\vec{e}}_i,j}& b_{\vec{n},j}\end{array}\right)=\sum _{k=1}^d{a}_{\vec{n}+{\vec{e}}_j,k}-\sum _{k=1}^d{a}_{\vec{n}+{\vec{e}}_i,k},i<j,$$

$$\frac{a_{\vec{n},i}}{a_{\vec{n}+{\vec{e}}_j,i}}=\frac{b_{\vec{n}-{\vec{e}}_i,j}-b_{\vec{n}-{\vec{e}}_i,i}}{b_{\vec{n},j}-b_{\vec{n},i}},i\ne j.$$

$${\nabla }_j{b}_{\vec{n},i}={\nabla }_i{b}_{\vec{n},j},i<j,$$

$$b_{\vec{n},j}{\nabla }_i{b}_{\vec{n},j}-b_{\vec{n},i}{\nabla }_j{b}_{\vec{n},i}=⟨({\overrightarrow{\nabla }}_j-{\overrightarrow{\nabla }}_i),{\vec{a}}_{\vec{n}}⟩,i<j,$$

$$({\nabla }_i\ln )a_{\vec{n},j}=({\nabla }_j\ln )(b_{\vec{n}-{\vec{e}}_j,i}-b_{\vec{n}-{\vec{e}}_j,j}),i\ne j,$$

$${\nabla }_j{b}_{\vec{n},i}:=b_{\vec{n}+{\vec{e}}_j,i}-b_{\vec{n},i},{\overrightarrow{\nabla }}_i:=({\nabla }_i,\cdots ,{\nabla }_i),({\nabla }_i\ln )a_{\vec{n},j}:=\left(\frac{a_{\vec{n}+{\vec{e}}_i,j}}{a_{\vec{n},j}}-1\right).$$

$$a_{\vec{n},j}=0,\text{mbox}whenever{n}_j=0,$$

$$\text{supp}({\mu }_i)=[{\alpha }_i,{\beta }_i],\text{with }{\alpha }_i<{\beta }_i<{\alpha }_{i+1}\text{for all }i.$$

$$P_{\vec{n}}(z)=:\prod _{i=1}^d\prod _{l=1}^{n_i}(z-x_{\vec{n},i,l}),x_{\vec{n},i,l}\in [{\alpha }_i,{\beta }_i].$$

$$n_i=t_i|\vec{n}|+o\left(|\vec{n}|\right),|\vec{t}|=1$$

$$\omega (x):=\underset{\mathcal{N}}{\lim }\frac{1}{|\vec{n}|}\sum _{i=1}^d\sum _{l=1}^{n_i}\delta (x-x_{\vec{n},i,l}),$$

$$M_{\vec{t}}\big{(}\{\alpha_{i},\beta_{i}\}_{1}^{d}\big{)}:=\big{\{}\vec{\nu}=(\nu_{1},\ldots,\nu_{d}):~{}\nu_{i}\in M_{t_{i}}(\alpha_{i},\beta_{i}),~{}i\in\{1,\ldots,d\}\big{\}},$$

$$I[\>\vec{\omega}\>]=\min_{\nu\in M_{\vec{t}}(\{\alpha_{i},\beta_{i}\}_{1}^{d})}I[\>\vec{\nu}\>],\qquad I[\>\vec{\nu}\>]:=\sum_{i=1}^{d}\bigg{(}2I[\nu_{i}]+\sum_{k\neq i}I[\nu_{i},\nu_{k}]\bigg{)},$$

$$\text{supp}({\omega }_i)=[{\alpha }_{\vec{t},i},{\beta }_{\vec{t},i}]\subseteq [{\alpha }_i,{\beta }_i],i\in \{1,\dots ,d\}.$$

$$d{\mu }_i(x)={\rho }_i(x)dx,$$

$$\underset{{\mathcal{N}}_{\vec{t}}}{\lim }{a}_{\vec{n},i}=A_{\vec{t},i}\text{and}\underset{{\mathcal{N}}_{\vec{t}}}{\lim }{b}_{\vec{n},i}=B_{\vec{t},i},i\in \{1,\dots ,d\},$$

$$\underset{{\mathcal{N}}_{\vec{t}}}{\lim }{a}_{\vec{n},i}\text{and}\underset{{\mathcal{N}}_{\vec{t}}}{\lim }{b}_{\vec{n},i}$$

$$S_{d-1}:=\{\vec{s}\in [0,1{]}^{d-1}:|\vec{s}|\le 1\}.$$

$$A_j(\vec{s})=\underset{\mathcal{N}}{\lim }{a}_{\vec{n},j},$$

$$B_j(\vec{s})=\underset{\mathcal{N}}{\lim }{b}_{\vec{n},j},$$

$$A_j^{(m)}(\vec{s})=a_{\vec{n},j}.$$

$$A_j^{(m)}(\vec{u}):=\sum _{k=1}^{2^{d-2}}\left[A_j^{(m)}(P^{(k)})\prod _{l=1}^{d-1}\frac{q_l^{(k)}-u_l}{q_l^{(k)}-p_l^{(k)}}+A_j^{(m)}(Q^{(k)})\prod _{l=1}^{d-1}\frac{p_l^{(k)}-u_l}{p_l^{(k)}-q_l^{(k)}}\right].$$

$$|A_j^{(m)}(\vec{s})-A_j(\vec{s})|$$

$$|B_j^{(m)}(\vec{s})-B_j(\vec{s})|$$

$$\nabla B_i(\vec{s})\cdot \left({\vec{\delta }}_j-\vec{s}\right)=\nabla B_j(\vec{s})\cdot \left({\vec{\delta }}_i-\vec{s}\right),i<j,$$

$$B_j(\vec{s})\nabla B_j(\vec{s})\cdot \left({\vec{\delta }}_i-\vec{s}\right)-B_i(\vec{s})\nabla B_i(\vec{s})\cdot \left({\vec{\delta }}_j-\vec{s}\right)=\left(\sum _{l=1}^d\nabla A_l\right)\cdot \left({\vec{\delta }}_j-{\vec{\delta }}_i\right),i<j,$$

$$A_j(\vec{s})\nabla \left(B_i(\vec{s})-B_j(\vec{s})\right)\cdot \left(\vec{s}-{\vec{\delta }}_j\right)+(B_i(\vec{s})-B_j(\vec{s}))\nabla A_j(\vec{s})\cdot \left({\vec{\delta }}_i-\vec{s}\right)=0,i\ne j.$$

$$\left|\frac{\partial }{\partial s_k}{A}_j^{(m)}(\vec{s})-\frac{\partial }{\partial s_k}{A}_j(\vec{s})\right|$$

$$\left|\frac{\partial }{\partial s_k}{B}_j^{(m)}(\vec{s})-\frac{\partial }{\partial s_k}{B}_j(\vec{s})\right|$$

$$\left|\frac{\partial }{\partial s_1}{A}_j(\vec{s})-\frac{\partial }{\partial s_1}{A}_j(\vec{u})\right|\le \frac{\epsilon }{4}$$

Full text

In the case of two orthogonality measures the differential equation becomes ordinary. For Angelesco systems, the result is illustrated numerically.

1.  Introduction

**1.1.  **Orthogonal polynomials on the real line and the Jacobi matrices

Given a probability measure $\mu$ on ${\mathbb{R}}$ with infinite support, the sequence of its monic orthogonal polynomials $\{P_{k}\}_{k=0}^{\infty}$ satisfies the well-known three-term recurrence relation

[TABLE]

with $P_{-1}=0$, $P_{0}=1$, where the recurrence coefficients $\{a_{k},b_{k}\}_{k=0}^{\infty}$ satisfy $a_{k}>0$, $b_{k}\in{\mathbb{R}}$.

The corresponding Jacobi matrix is defined to be

[TABLE]

Assuming $\{a_{k}\}_{k=0}^{\infty}$ and $\{b_{k}\}_{k=0}^{\infty}$ are bounded, the spectral measure of ${\mathcal{J}}$ with respect to $(1,0,0,\ldots)^{T}$ coincides with the orthogonality measure $\mu$. Favard’s theorem establishes a one-to-one correspondence between all $\mu$ with compact infinite support and all such bounded self-adjoint Jacobi matrices ${\mathcal{J}}$.

We say that a probability measure on ${\mathbb{R}}$ belongs to the Nevai class $N(a,b)$ if its Jacobi coefficients (in (1.1)) satisfy $b_{n}\to b$ and $a_{n}\to a$ as $n\to\infty$.

Weyl’s theorem on compact perturbations implies that any measure in $N(a,b)$ has $\sigma_{ess}(\mu)=[b-2\sqrt{a},b+2\sqrt{a}]$. For the converse direction, we have the Denisov–Rakhmanov theorem stating that if $\sigma_{ess}(\mu)=[\alpha,\beta]$ and $\frac{d\mu}{dx}>0$ a.e. on $[\alpha,\beta]$ then \mu\in N(\big{(}\tfrac{\beta-\alpha}{4}\big{)}^{2},\tfrac{\alpha+\beta}{2}).

See, e.g., [14] for more details from the theory of orthogonal polynomials.

**1.2.  **Multiple orthogonal polynomials and the nearest neighbor recurrence relations

Let us now describe multiple orthogonality situation with respect to the vector-measure $\vec{\mu}:=\{\mu_{i}\}_{i=1}^{d}$ on ${\mathbb{R}}$. For the rest of the paper we will use the notation $|\vec{v}|:=v_{1}+...+v_{d}$ for any vector-valued object $\vec{v}=(v_{j})_{j=1}^{d}$.

For any $\vec{n}=(n_{1},\ldots,n_{d})\in\mathbb{Z}^{d}_{+}$, let $P_{\vec{n}}$ be the monic polynomial of smallest degree which satisfies

[TABLE]

The polynomial $P_{\vec{n}}(x)$ is called the type II multiple orthogonal polynomial (MOP). Obviously, $P_{\vec{n}}$ is uniquely determined and $\deg P_{\vec{n}}\leq|\vec{n}|$. When $\deg P_{\vec{n}}=|\vec{n}|$ the multi-index $\vec{n}$ is said to be normal. If all multi-indices of the lattice $\mathbb{Z}^{d}_{+}$ are normal then the system of measures $\{\mu_{i}\}_{i=1}^{d}$ is called perfect. It is known [15, 16], that (similarly to the case with one measure) MOPs for the perfect systems satisfy the following nearest neighbor recurrence relations (NNRR)

[TABLE]

where $\vec{e}_{j}$ is the $j$-th standard basis vector of ${\mathbb{R}}^{d}$. Here we have $d$ recurrence relations for $j=1,\ldots,d$. Thus for each $\vec{n}\in\mathbb{Z}^{d}_{+}$ we have two sets of the coefficients for NNRR, namely $\{b_{\vec{n},j}\}_{j=1}^{d}$ and $\{a_{\vec{n},i}\}_{i=1}^{d}$. Note that for each fixed $j$, $\{a_{k\vec{e}_{j},j}\}_{k=1}^{\infty}$ and $\{b_{k\vec{e}_{j},j}\}_{k=1}^{\infty}$ are the $\{a_{k}\}_{k=0}^{\infty}$ and $\{b_{k}\}_{k=0}^{\infty}$ from the usual three-term recurrence (1.1) for the measure $\mu_{j}$.

In order to define by means of (1.4) the polynomials $\{P_{\vec{n}}(z)\}$ in unique way the NNRR coefficients cannot be taken arbitrary. As was shown in [16],the recurrence coefficients must satisfy the compatibility conditions (CC):

[TABLE]

It is not hard to see that these $2d(d-1)$ equalities can be rewritten as

[TABLE]

where we denote

[TABLE]

The system of difference equations (1.8)–(1.10) together with the marginal conditions

[TABLE]

is also called Discrete Integrable System (DIS) for details see [3]. The boundary problem for DIS (1.8)–(1.10) in $\mathbb{Z}^{d}_{+}$ means the following. Given the boundary data: coefficients of the $d$-collections of the three-terms recurrence relations, corresponding to usual orthogonal polynomials with respect to each $\{\mu_{i}\}_{i=1}^{d}$ measure. Then solving equations (1.8)–(1.10) we have to find all NNRR coefficients $\{b_{\vec{n},j}\}_{j=1}^{d}$ and $\{a_{\vec{n},i}\}_{i=1}^{d}$.

**1.3.  **Zero asymptotics and limits of the recurrence coefficients

Our goal is to investigate the asymptotic behavior of the recurrence coefficients \big{\{}a_{\vec{n},i},b_{\vec{n},i}\big{\}} as $|\>\vec{n}\>|$ grows. This behavior is intimately connected to the asymptotic zero distribution of multiple orthogonal polynomials $P_{\vec{n}}$. To state the problem, we need to place some restrictions on the way $|\>\vec{n}\>|$ approaches infinity as well as the measures $\mu_{i}$. At the same time we have to be in the class of the perfect systems to keep NNRR.

The important example of a perfect system of measures $\{\mu_{i}\}$ is the so-called Angelesco system defined by 111If supports of measures are intervals with nonintersecting interiors then system $\{\mu_{i}\}$ is perfect as well.

[TABLE]

Multiple orthogonal polynomial with respect to Angelesco system has the form:

[TABLE]

Moreover, we restrict our attention to sequences of multi-indices such that

[TABLE]

for some $\vec{t}\in(0,1)^{d}$. We denote $\lim_{\mathcal{N}}$ to be the limit as $|\vec{n}|\to\infty$ along the sequence of multi-indices satisfying (1.13). Asymptotic zero distribution for $P_{\vec{n}}(z)$ (or limiting zero counting measure):

[TABLE]

for Angelesco systems (1.12) with $\mu^{\prime}_{i}>0$ a.e. on $[\alpha_{i},\beta_{i}]$ in the regime (1.13) was obtained by Gonchar and Rakhmanov [10]. To state their result we fix $\vec{t}$ as in (1.13), and denote

[TABLE]

where $M_{t}(\alpha,\beta)$ is the set of positive Borel measures of mass $t$ supported on $[\alpha,\beta]$.

Theorem 1** ([10]).**

1)There exists the unique vector of measures \vec{\omega}\in M_{\vec{t}}\big{(}\{\alpha_{i},\beta_{i}\}_{1}^{d}\big{)}\,:

[TABLE]

where $I[\nu_{i}]:=I[\nu_{i},\nu_{i}]$ and $I[\nu_{i},\nu_{k}]:=-\int\int\log|z-x|{d}\nu_{i}(x){d}\nu_{k}(z)$.

2)* Moreover, for the limiting counting measure (1.14) it holds: $\omega=|\vec{\omega}|$.*

An important feature of the case $d>1$ (in comparison with the classic $d=1$) is the fact that measures $\omega_{i}$ might no longer be supported on the whole intervals $[\alpha_{i},\beta_{i}]$ (the so-called pushing effect), but in general it holds that

[TABLE]

Namely the supports of the extremal measures (not the supports of the multiple orthogonality measures 222For $d=1$ both of these notions coincide. ) define the recurrence coefficients limits.

To describe the asymptotics of the recurrence coefficients, we shall need a $(d+1)$-sheeted compact Riemann surface, say $\mathfrak{R}_{\vec{t}}$, that we realize in the following way. Take $d+1$ copies of $\overline{\mathbb{C}}$. Cut one of them along the union \bigcup_{i=1}^{d}\big{[}\alpha_{\vec{t},i},\beta_{\vec{t},i}\big{]}, which henceforth is denoted by $\mathfrak{R}_{\vec{t}}^{(0)}$. Each of the remaining copies are cut along only one interval \big{[}\alpha_{\vec{t},i},\beta_{\vec{t},i}\big{]} so that no two copies have the same cut and we denote them by $\mathfrak{R}_{\vec{t}}^{(i)}$. To form $\mathfrak{R}_{\vec{t}}$, take $\mathfrak{R}_{\vec{t}}^{(i)}$ and glue the banks of the cut \big{[}\alpha_{\vec{t},i},\beta_{\vec{t},i}\big{]} crosswise to the banks of the corresponding cut on $\mathfrak{R}_{\vec{t}}^{(0)}$. It can be easily verified that thus constructed Riemann surface has genus 0. Denote by $\pi$ the natural projection from $\mathfrak{R}_{\vec{t}}$ to $\overline{\mathbb{C}}$. We also shall employ the notations z for a point on $\mathfrak{R}_{\vec{t}}$ and $z^{(i)}$ for a point on $\mathfrak{R}_{\vec{t}}^{(i)}$ with $\pi(\textbf{z})=\pi(z^{(i)})=z$.

Since $\mathfrak{R}_{\vec{t}}$ has genus zero, one can arbitrarily prescribe zero/pole multisets of rational functions on $\mathfrak{R}_{\vec{t}}$ as long as the multisets have the same cardinality. Hence, we define $\Upsilon_{i}$, $i\in\{1,\ldots,d\}$, to be the rational function on $\mathfrak{R}_{\vec{t}}$ with a simple zero at $\infty^{(0)}$, a simple pole at $\infty^{(i)}$, and otherwise non-vanishing and finite. We normalize it so that $\Upsilon_{i}(z^{(i)})/z\to 1$ as $z\to\infty$. Then the following theorem holds.

Theorem 2** ([2]).**

Let $\{\mu_{i}\}_{i=1}^{d}$ be a system of measures satisfying (1.12) and such that

[TABLE]

where $\rho_{i}$ is holomorphic and non-vanishing in some neighborhood of $[\alpha_{i},\beta_{i}]$. Further, let $\mathcal{N}_{\vec{t}}=\{\>\vec{n}\>\}$ be a sequence of multi-indices as in (1.13) for some $\vec{t}\in(0,1)^{d}$. Then the recurrence coefficients \big{\{}a_{\vec{n},j},b_{\vec{n},j}\big{\}} given by (1.4) and (1.3) satisfy

[TABLE]

where $A_{\vec{t},i}$ and $B_{\vec{t},i}$ are constants: \textbf{z}^{2}\Upsilon_{i}(z^{(0)})=A_{\vec{t},i}(z+B_{\vec{t},i})+\mathcal{O}\big{(}z^{-1}\big{)} as $z\to\infty$.

Remarks. 1) We note that Theorem 2 is valid for $d=1$ as well.

2) It is not too difficult to extend the proof (from [10]) of Theorem 1 to include the case of touching intervals.

3) We also can affirm (at least for $d=2$) that Theorem 2 remains valid for the case of touching intervals (technicalities can be taken from [7]) and for weight functions (1.17) with singularities of the types: Jacobi and Fisher-Hartwig weights [18]. $\blacksquare$

Let us make the following definition by analogy with the scalar case (see Section 1.1).

Definition**.**

We say that a perfect system of measures $\{\mu_{i}\}_{i=1}^{d}$ belongs to the multiple Nevai class if for each $i\in\{1,\ldots,d\}$ the limits

[TABLE]

*exist along each sequence (1.13) for any $\vec{t}\in[0,1]^{d}$, $|\vec{t}|=1$. *

Perfect systems from multiple Nevai class appear naturally in various contexts [1, 4, 6, 11, 17], e.g., in random matrix theory [8]. Note that if a system of measures belongs to a multiple Nevai class, then the recurrence along the step-line has asymptotically periodic recurrence coefficients.

Notice that Theorem 2 can be viewed as a partial analogue of the Denisov–Rakhmanov theorem, and Angelesco systems from Theorem 2 belong to the multiple Nevai class. It is an interesting open problem to generalize this analogue of Denisov–Rakhmanov result to more general measures (i.e. to Angelesco systems with $\mu_{j}>0$ a.e. on $\sigma_{ess}(\mu_{j})$).

The organization of the paper is as follows. In Section 2 we state and prove our main result: a conditional theorem on partial differential equations for the limiting value (in the regime (1.13)) of the NNRR coefficients. In Section 3 we discuss the special case of two $d=2$ orthogonality measures when our partial differential equations become ordinary differential equations. In Section 4, using a parametrization of $\mathfrak{R}_{\vec{t}}$ from [13], we give a constructive procedure for determination of limits in (1.18). Finally, in Section 5 we present numeric illustrations.

2.  Differential equations for the limits of NNRR coefficients

**2.1.  **Construction of the approximating functions

For the rest of the paper, let us denote

[TABLE]

Assume that $\{\mu_{j}\}_{j=1}^{d}$ form a perfect system from the multiple Nevai class.

This means that there exist $S_{d-1}\to{\mathbb{R}}$ functions $A_{j}(\vec{s}),B_{j}(\vec{s})$ ($1\leq j\leq d$) defined via

[TABLE]

where $\lim_{\mathcal{N}}$ notation is defined in Section 1.3 with $\vec{t}=\{\vec{s},1-|\vec{s}|\}$ (that is, $\vec{s}$ consists of the first $d-1$ coordinates of $\vec{t}$ which defines the direction of the approach to infinity).

In this paper we investigate the possibility of describing functions $\{A_{j},B_{j}\}_{j=1}^{d}$ through differential equations. This is done in Theorem 3 below.

Before stating the main result, let us introduce the families of approximations $A_{j}^{(m)}$ and $B_{j}^{(m)}$ of the limiting functions $A_{j}$ and $B_{j}$.

Fix $m\in{\mathbb{Z}}_{+}$ and $1\leq j\leq d$. We take all the coefficients $\{a_{\vec{n},j}\}$ with $|\vec{n}|=m$ and form an approximating function $A^{(m)}(\vec{s})$ as follows. First, for any $\vec{n}$ with $|\vec{n}|=m$, define $\vec{s}\in S_{d-1}$ via $s_{j}:=\tfrac{n_{j}}{m}$ ($1\leq j\leq d-1$) and let

[TABLE]

For points in $\tfrac{1}{m}{\mathbb{Z}}_{+}^{d-1}$ that are not in $S_{d-1}$ we can choose $A^{(m)}$ to be zero. Then we can extend $A^{(m)}_{j}$ to the rest of the simplex $S_{d-1}$ via the multilinear interpolation which can be written as follows. Choose a cube $K$ of side length $\tfrac{1}{m}$ with vertices in $\tfrac{1}{m}{\mathbb{Z}}_{+}^{d-1}$; let us denote them $\{P^{(k)},Q^{(k)}\}_{k=1}^{2^{d-2}}$, where for each $k$, vertices $P^{(k)}$ and $Q^{(k)}$ are opposite of each other. If $P^{(k)}=(p^{(k)}_{1},\ldots,p^{(k)}_{d-1})$ and $Q^{(k)}=(q^{(k)}_{1},\ldots,q^{(k)}_{d-1})$ then we let

[TABLE]

for $\vec{u}\in K$.

The main features of this multilinear interpolation function (2.4) that are important to us are:

1. The right-hand side of (2.4) agrees with the left-hand side of (2.4) when $\vec{u}\in\{P^{(k)},Q^{(k)}\}_{k=1}^{2^{d-2}}$, so that the function is well defined at the vertices of our cubes;

2. For $\vec{u}$ belonging to any face of a cube $K$, the expression (2.4) reduces to the multilinear interpolation of one dimension lower over the vertices of that face. As a result,  (2.4) on a face of a cube $K$ will agree with (2.4) defined through another cube sharing the same face. So the function $A^{(m)}$ is well-defined on $S_{d-1}$. Moreover, it is continuous on $S_{d-1}$ and is differentiable on the interiors of each of the cubes $K$;

3. In each of the $d-1$ variables $u_{l}$, the function $A_{j}^{(m)}$ is linear within each of the cubes $K$. This will be used in the proof of Theorem 4 below;

4. Partial derivatives of the right-hand side of (2.4) are linear functions along each path parallel to the coordinate axes. In particular, it implies that the maxima and minima over $K$ of partial derivatives of $A^{(m)}_{j}$ are attained at $\{P^{(k)},Q^{(k)}\}_{k=1}^{2^{d-2}}$. This will be used in the proof of Lemma 1 below.

We can do the same construction with coefficients $b_{\vec{n},j}$ to form the multilinear approximations $B^{(m)}_{j}:S_{d-1}\to{\mathbb{R}}$ for functions $B_{j}$.

Notice that (2.2)–(2.3) implies pointwise convergence $A^{(m)}_{j}$ and $B^{(m)}_{j}$ on $S_{d-1}$ to $A_{j}$ and $B_{j}$, respectively, as $m\to\infty$.

**2.2.  **The main theorem

For the rest of the paper we assume that the functions $A_{j}$ and $B_{j}$ ($1\leq j\leq d$) are piecewise continuously differentiable on $S_{d-1}$ in the following sense. We suppose that $S_{d-1}$ can be decomposed into a finite union of closed sets $\{D_{i}\}$ such that:

(i) $A_{j}$ and $B_{j}$ are differentiable on the interior $\operatorname{Int}(D_{i})$;

(ii) Each of the partial derivatives of $A_{j}$ and $B_{j}$ are continuous $\operatorname{Int}(D_{i})$ and can be continuously extended to $D_{i}$.

Note that the latter condition means that each of the partial derivatives of $A_{j}$ and $B_{j}$ is uniformly continuous on $\operatorname{Int}(D_{i})$, a fact that we use in the proof of Lemma 1.

We also assume that sets $D_{i}$ are not pathological, in particular, the closure of $\operatorname{Int}(D_{i})$ is assumed to be $D_{i}$.

Recall that $\{\vec{e}_{j}\}_{j=1}^{d}$ is the standard basis of ${\mathbb{R}}^{d}$. For the notational convenience, let us denote $\vec{\delta}_{j}$ ($1\leq j\leq d-1$) to be the $j$-th standard basis vector in ${\mathbb{R}}^{d-1}$, while $\vec{\delta}_{d}$ to be the zero vector in ${\mathbb{R}}^{d-1}$.

Theorem 3**.**

Assume that we have a perfect system $\{\mu_{j}\}_{j=1}^{d}$ from the multiple Nevai class satisfying the conditions

• (i)

$A_{j}$* and $B_{j}$ are piecewise continuously differentiable on $S_{d-1}$ for each $1\leq j\leq d$;*

• (ii)

For each $1\leq j\leq d$, we have uniform convergence:

[TABLE]

as $m\to\infty$, where sequences $o(\tfrac{1}{m})$ are independent of $\vec{s}\in S_{d-1}$.

Then the limiting functions $A_{j}$ and $B_{j}$, $1\leq j\leq d$, satisfy the following system of $2d(d-1)$ differential equations:

[TABLE]

In the system (2.7)–(2.9), $\vec{u}\cdot\vec{v}$ stands for the standard inner product in ${\mathbb{R}}^{d-1}$, and $\nabla$ for the gradient operator for a function of $d-1$ variables.

Remarks.

1) Condition (i) is fulfilled for Angelesco systems from Theorem 2. This follows from smoothness of the dependence of the residues of $\Upsilon$ on $\vec{t}$. We show it explicitly for $d=2$ in the last section. As for (ii), (2.5)–(2.6) holds uniformly on compacts of $\operatorname{Int}(S_{d-1})$ (this follows from the proof of Theorem 2). Whether this can be extended to the whole $S_{d-1}$ is still unknown.

2) Since the system $\{\mu_{j}\}_{j=1}^{d}$ is in the multiple Nevai class determined by the functions $\{A_{j},B_{j}\}_{j=1}^{d}$, each of the measures $\mu_{j}$ is in the Nevai class, in particular its essential support is an interval. These intervals (together with (1.11)) allow one to establish boundary conditions for the functions $\{A_{j},B_{j}\}_{j=1}^{d}$. We do this explicitly for $d=2$ in the next section.

**2.3.  **Convergence of the derivatives

In order to prove Theorem 3, we will need to control the derivatives of our approximation functions. This is the purpose of the following lemma.

Lemma 1**.**

Suppose (i)–(ii) of Theorem 3 hold. Then for $1\leq k\leq d-1$ and any point $\vec{s}_{0}$ in $\operatorname{Int}(D_{i})$, there exists a neighbourhood $U(\vec{s}_{0})\subset\operatorname{Int}(D_{i})$ such that

[TABLE]

for all $\vec{s}\in U(\vec{s}_{0})$ as $m\to\infty$, where $o(1)$ is independent of $\vec{s}\in U(\vec{s}_{0})$.

Remark.

Partial derivatives of $A^{(m)}_{j}$ and $B^{(m)}_{j}$ have jump discontinuities along each side of the $\tfrac{1}{m}{\mathbb{Z}}_{+}^{d-1}$ cubes (see Section 2.1). At a point of discontinuity, we interpret $\frac{\partial}{\partial s_{k}}A^{(m)}_{j}(\vec{s})$ and $\frac{\partial}{\partial s_{k}}B^{(m)}_{j}(\vec{s})$ in (2.10) and (2.11) as one of the limiting values of these functions from the inside of one of the cubes.

Proof.

Fix $j$. Let us prove (2.10) for $k=1$.

Choose $M_{1}\in{\mathbb{N}}$ large enough so that a cube with side length $\tfrac{2}{M_{1}}$ centered at $\vec{s}_{0}$ belongs to $D_{i}$. Let $U(\vec{s}_{0})$ be the cube centered at $\vec{s}_{0}$ of side length $\tfrac{1}{M_{1}}$.

Let $\varepsilon>0$ be arbitrary. By the discussion in the beginning of the section, $\frac{\partial}{\partial s_{1}}A_{j}$ is uniformly continuous on $D_{i}$. We can therefore find $M_{2}\in{\mathbb{N}}$ so that

[TABLE]

for all $\vec{s}$ and $\vec{u}$ in $D_{i}$ satisfying $||\vec{s}-\vec{u}||\leq\tfrac{1}{M_{2}}$. By  (2.5) we can find $M_{3}\in{\mathbb{N}}$ so that

[TABLE]

for all $\vec{s}\in S_{d-1}$ and $m\geq M_{3}$. Now let $M=\max\{M_{1},M_{2},M_{3}\}$.

For any $\vec{s}$ in $U(\vec{s}_{0})$ and any $m\geq M$, choose a cube $K(m)$ of side length $\tfrac{1}{m}$ containing $\vec{s}$ whose vertices are at $\tfrac{1}{m}{\mathbb{Z}}_{+}^{d-1}$ (as in Section 2.1). By the construction, $K$ belongs to $D_{i}$, and  (2.12) and (2.13) hold for our $m$.

Let us first show that the inequality (2.10) holds for the case when $\vec{s}$ is a vertex of $K$. If $\vec{s}+\tfrac{1}{m}\vec{e}_{1}$ is also a vertex of $K$ (arguments for $\vec{s}-\tfrac{1}{m}\vec{e}_{1}$ are identical), then by the discussion after (2.4),

[TABLE]

for some $\vec{\theta}\in(\vec{s},\vec{s}+\tfrac{1}{m}\vec{e}_{1})$. Here we used (2.13) twice and the Mean Value Theorem. The last expression is $\leq 3\varepsilon/4$ by (2.12).

Now if $\vec{s}$ is not a vertex of $K$, then by the discussion after (2.4), there are vertices $\vec{z}_{1}$ and $\vec{z}_{2}$ of $K$ such that $\frac{\partial}{\partial s_{1}}A^{(m)}_{j}(\vec{z}_{1})\leq\frac{\partial}{\partial s_{1}}A^{(m)}_{j}(\vec{s})\leq\frac{\partial}{\partial s_{1}}A^{(m)}_{j}(\vec{z}_{2})$. By (2.12), $\frac{\partial}{\partial s_{1}}A_{j}(\vec{z}_{2})-\tfrac{\varepsilon}{4}\leq\frac{\partial}{\partial s_{1}}A_{j}(\vec{s})\leq\frac{\partial}{\partial s_{1}}A_{j}(\vec{z}_{1})+\tfrac{\varepsilon}{4}$. Combining these two inequalities together with the estimate at the vertices, we get $\left|\frac{\partial}{\partial s_{1}}A^{(m)}_{j}(\vec{s})-\frac{\partial}{\partial s_{1}}A_{j}(\vec{s})\right|\leq\varepsilon$. ∎

**2.4.  **Proof of Theorem 3

Let $\vec{s}\in S_{d-1}$ belongs to the interior of some $D_{i}$. Choose a neighbourhood $U\subset D_{i}$ of $\vec{s}$ as in Lemma 1. We can assume $\overline{U}\subset D_{i}$ (just shrink $U$ if needed). Let a sequence of multi-indices $\vec{n}$ be given satisfying (1.13) with $\vec{t}=\{\vec{s},1-|\vec{s}|\}$, and as a result (2.2), (2.3) also. For each such $\vec{n}$, let $m=|\vec{n}|$ and define $\vec{s}\,^{(m)}\in S_{d-1}$ with $s^{(m)}_{j}=\frac{n_{j}}{m}$. Then $\vec{s}\,^{(m)}\to\vec{s}$. For each $m$ let $K_{m}$ be a cube of side length $\tfrac{1}{m+1}$ containing $\vec{s}\,^{(m)}$ whose vertices are at $\tfrac{1}{m+1}{\mathbb{Z}}_{+}^{d-1}$ (as in Section 2.1). We consider $m$ large enough so that each $K_{m}$ belongs to $U$.

Let $1\leq i\leq d-1$. Notice that by Taylor’s theorem

[TABLE]

where on the last step we used (ii) of Theorem 4 and Lemma 1. However the $o(\tfrac{1}{m})$ error term in (2.15) is dependent on $\vec{s}\,^{(m)}$ and can in principle be non-uniform. To justify uniformity in (2.17) we proceed as follows. We start with (2.14), and note that $\tfrac{m}{m+1}\vec{s}\,^{(m)}+\tfrac{1}{m+1}\vec{\delta}_{i}=\vec{s}\,^{(m)}+\sum_{l=1}^{d-1}\vec{\Delta}^{(m,i)}_{l}$ where $\vec{\Delta}^{(m,i)}_{i}=\frac{1-s_{i}^{(m)}}{m+1}\vec{\delta}_{i}$ and $\vec{\Delta}^{(m,i)}_{l}=-\frac{s_{l}^{(m)}}{m+1}\vec{\delta}_{l}$ for $l\neq i$. These $\vec{\Delta}^{(m,i)}_{l}$ are just the increment $\tfrac{1}{m+1}\left(\vec{\delta}_{i}-\vec{s}\,^{(m)}\right)$ from $\vec{s}\,^{(m)}$ to $\tfrac{m}{m+1}\vec{s}\,^{(m)}+\tfrac{1}{m+1}\vec{\delta}_{i}$ separated in coordinates, and $\sum_{l=1}^{d-1}\vec{\Delta}^{(m,i)}_{l}=\tfrac{1}{m+1}\left(\vec{\delta}_{i}-\vec{s}\,^{(m)}\right)$. Now recall that the multilinear approximation function $A_{j}^{(m)}$ (2.4) is linear along coordinate axes, so applying this for each of the $d-1$ increment we get:

[TABLE]

where on the last step we used (ii) of Theorem 4 and Lemma 1 (notice that now $o(\tfrac{1}{m})$ is uniform!). Now for any $p$, $\nabla A_{j}\left(\vec{s}\,^{(m)}+\sum_{l=1}^{p-1}\vec{\Delta}^{(m,i)}_{l}\right)=\nabla A_{j}\left(\vec{s}\,^{(m)}\right)+o(1)$ (with uniform $o(1)$), since $\vec{\Delta}^{(m,i)}_{l}=o(1)$ for each $l$ and $\nabla A_{j}$ is continuous and therefore uniformly continuous on $\overline{U}$. Plugging this into the last equation and using $\vec{\Delta}^{(m,i)}_{p}=o(\tfrac{1}{m})$ implies (2.17) with uniform $o(\tfrac{1}{m})$.

Similar arguments give us for $1\leq i\leq d-1$,

[TABLE]

with uniform $o(\tfrac{1}{m})$. For $i=d$, we get the following expressions instead:

[TABLE]

with uniform $o(\tfrac{1}{m})$. Notice that these expressions for $a_{\vec{n}\pm\vec{e}_{d},j}$ agree with the expressions for $a_{\vec{n}\pm\vec{e}_{i},j}$ (with $i\leq d-1$) if we adopt our notation $\vec{\delta}_{d}=\vec{0}\in{\mathbb{R}}^{d-1}$.

Analogous equalities hold for the $b$-coefficients and the corresponding $B_{j}$ functions.

Let us now plug these equalities into (1.8). For any $1\leq i<j\leq d$, we get:

[TABLE]

Now multiply by $m+1$, and take limit as $m\to\infty$. Using continuity of $\nabla B_{j}$’s, we obtain  (2.7).

Similar computations show that  (1.9) leads to  (2.8) (for any $1\leq i<j\leq d$) and (1.10) produces  (2.9) (for any $i\neq j$, $1\leq i,j\leq d$).

3.  $d=2$ case: system of ordinary differential equations

**3.1.  **The main theorem: $d=2$

In the case $d=2$, we have four functions $A_{1},A_{2},B_{1},B_{2}$ of one variable $s\in[0,1]$, and the corresponding differential system takes the form stated below.

Theorem 4**.**

1)* Assume that we have a perfect system $\mu_{1},\mu_{2}$ from the multiple Nevai class satisfying the conditions*

• (i)

$A_{j}$* and $B_{j}$ are piecewise continuously differentiable on $[0,1]$ for each $1\leq j\leq 2$;*

• (ii)

For each $1\leq j\leq 2$, we have uniform convergence:

[TABLE]

as $m\to\infty$, where sequences $o(\tfrac{1}{m})$ are independent of $s\in[0,1]$.

Then the limiting functions $A_{j}$ and $B_{j}$, $1\leq j\leq 2$, satisfy the following system of ordinary differential equations:

[TABLE]

where

[TABLE]

2)* Suppose an Angelesco system satisfies conditions of Theorem 2 and (ii). Then there exist $c_{1},c_{2}\in(0,1)$ such that the functions $A_{1},A_{2},B_{1},B_{2}:[0,1]\to{\mathbb{R}}$ are smooth on $[0,c_{1})$ and $(c_{2},1]$, and satisfy the system of differential equations*

[TABLE]

with initial/boundary conditions

[TABLE]

where

[TABLE]

and (3.4). Moreover, $A_{1},A_{2},B_{1},B_{2}$ are constant on the interval $[c_{1},c_{2}]$.

Remarks. 1) We note that general (and conditional) part 1) of Theorem 4 admits presence inside $[0,1]$ of a subdomain, where $A_{1},A_{2},B_{1},B_{2}$ are constant. For Angelesco systems it is a generic situation which happens when “pushing” is not active, see [10].

2) We note that known information about support of zero counting measure of MOP for Angelesco system (see [10]) allows us to identify the subdomain where $A_{1},A_{2},B_{1},B_{2}$ are constant, i.e. interval $[c_{1},c_{2}]$. Then it is possible, using BC (3.6) and (3.7) to solve the ODE system (3.5) on $[0,c_{1}]$ and $[c_{2},1]$.

**3.2.  **Proof of Theorem 4

Taking $d=2$ in  (2.7) (with $i=1,j=2$),  (2.8) (with $i=1,j=2$), and  (2.9) (with $i=2,j=1$; then $i=1,j=2$) gives us four ODE’s:

[TABLE]

Let us simplify this system. First of all, let

[TABLE]

Using (3.8) and (3.9), we get $B_{1}^{\prime}=\frac{A_{1}^{\prime}+A_{2}^{\prime}}{sB}$, $B_{2}^{\prime}=-\frac{A_{1}^{\prime}+A_{2}^{\prime}}{(1-s)B}$, so $B^{\prime}=B_{2}^{\prime}-B_{1}^{\prime}=-\frac{A_{1}^{\prime}+A_{2}^{\prime}}{s(1-s)B}$. This equation together (3.10) and (3.11) established (3.3). Part 1) of the Theorem 4 is proved.

Let us divide interval $[0,1]$ into two disjoint sets:

[TABLE]

From [10] we know that: $I_{1}$ consists of one point if $\Delta_{1}$ and $\Delta_{2}$ are touching, and otherwise $I_{1}$ is an interval $[c_{1},c_{2}]$ inside $(0,1)$.

For $s\in I_{2}$, the determinant of the matrix in (3.3) must be zero, i.e.,

[TABLE]

which implies

[TABLE]

on the set where $B(s)\neq 0$. This means that

[TABLE]

Plugging this into the third equation of  (3.3), we get

[TABLE]

which simplifies to

[TABLE]

The first two equations in (3.3) can be solved for $\tfrac{B^{\prime}(s)}{B(s)}$ giving us

[TABLE]

So our new system of two ODE’s is

[TABLE]

for $s\in I_{2}$.

It is not hard to notice from (3.13) that $A_{1}$ and $A_{2}$ have double zeros at [math] and $1$, respectively. So let

[TABLE]

Then our system (3.3) becomes:

[TABLE]

Using $B(s)^{2}=C_{1}(s)+C_{2}(s)$, we can eliminate $B$:

[TABLE]

Finally, let us deal with the boundary conditions for all of our functions.

Since our system $\{\mu_{j}\}_{j=1}^{2}$ is from a multiple Nevai class, we also have that $\mu_{1}$ and $\mu_{2}$ are in the (scalar) Nevai class $N({A_{1}(1)},B_{1}(1))$ and $N({A_{2}(0)},B_{2}(0))$, respectively. Since $\text{\rm{supp}}(\mu_{j})=[\alpha_{j},\beta_{j}]$, Weyl’s theorem (see Section 1.1) gives us:

[TABLE]

The marginal conditions (1.11) give us

[TABLE]

We also need the other two boundary conditions

[TABLE]

which can be obtained from Section 4 below.

This means that $B$ has boundary values

[TABLE]

Finally, to get the boundary values for $C_{1}(t)$ and $C_{2}(t)$, we recall that on the region $I_{2}$ where $B\neq 0$, we have

[TABLE]

Taking $s\to 0$, we therefore get $A_{1}(0)=A_{1}^{\prime}(0)=0$ and

[TABLE]

which implies $A_{1}^{\prime\prime}(0)=2B(0)^{2}-2A_{2}(0)$. Similarly, $s\to 1$ gives us $A_{2}(1)=A_{2}^{\prime}(1)=0$ and

[TABLE]

which implies $A_{2}^{\prime\prime}(1)=2B(1)^{2}-2A_{1}(1)$. Then $C_{1}(0)=\tfrac{1}{2}A_{1}^{\prime\prime}(0)$, $C_{2}(0)=A_{2}(0)$, $C_{1}(1)=A_{1}(1)$, $C_{2}(1)=\tfrac{1}{2}A_{2}^{\prime\prime}(1)$ which result in our boundary conditions (3.6)–(3.7). Part 2) of the theorem is now proved.

4.  Determination of the limits by means of parametrization of $\mathfrak{R}_{\vec{t}}$

In this section we employ an algebraically-geometric approach in order to determine the limits of the NNRR’s coefficients. We restrict the consideration to the case of Angelesco system with two orthogonality measures (we allow the supports to have a common endpoint). Thus in this setting we set

[TABLE]

Our input is the supports measures of orthogonality (1.12)

[TABLE]

Note that using the linear map $y(x)=(x-\beta_{1})/(\beta_{2}-\beta_{1})$, these segments can be transformed to

[TABLE]

where $y(\alpha_{1})=-\alpha$ and $y(\alpha_{2})=\beta$. Thus, without loss of generality, we can use (4.2) as the input.

Our goal is to construct the following procedure: based on Theorem 2, find the limits (2.2), (2.3) via computing the residues of $\Upsilon_{i},\,\,i=1,2$.

In order to reach this goal we have to solve two problems:

Problem 1. For each $s\in(0,1)$, find the segments $[\alpha_{s,i},\beta_{s,i}],\,\,\,i=1,2$, of the support of the extremal vector-measure $\vec{\omega}=(\omega_{1},\omega_{2})$, minimizing the energy functional (1.15).

Problem 2. Using the endpoints $\{\alpha_{s,i},\beta_{s,i}\}_{i=1}^{2}$ as the branch points of the Riemann surface $\mathfrak{R}_{s}$ (which is defined in subsection 1.3), find the limits $\{A_{s,i},B_{s,i}\}_{i=1}^{2}$ by computing the residues of the meromorphic on $\mathfrak{R}_{s}$ functions $\Upsilon_{i},\,\,i=1,2$.

**4.1.  **Parametrization

of $\mathfrak{R}_{s}$ and solution to Problem 2

To solve both problems we use (introduced in [5] and developed in [12], [13]) parametrization of the three-sheeted Riemann surfaces with four branch points.

We fix $s\in(0,1)$ and start with parametrization of $\mathfrak{R}(\alpha,\beta):=\mathfrak{R}_{s}$, where we take (4.2) for the intervals $[\alpha_{s,i},\beta_{s,i}],i=1,2$. We define

[TABLE]

It is not difficult to check that for $\alpha,\beta$ given in (4.2) there exists a unique solution of the equation

[TABLE]

We have the following

Theorem 5** ([5, 12, 13]).**

Riemann surface $\mathfrak{R}(\alpha,\beta)$ can be defined by means of the conformal map of the Riemann sphere $\overline{\mathbb{C}}\ni w\mapsto\textbf{z}(w)\in\mathfrak{R}(\alpha,\beta)$ given by

[TABLE]

where $\pi:\mathfrak{R}(\alpha,\beta)\to\overline{\mathbb{C}}$ is the natural projection.

Let $\tau_{0},\tau_{1},\tau_{2}$ be $\textbf{z}^{-1}(\infty^{(0)}),\textbf{z}^{-1}(\infty^{(1)}),\textbf{z}^{-1}(\infty^{(2)})$, respectively. Substituting $R_{u}$ from (4.3) into (4.5), we obtain

[TABLE]

where $\tau_{0}\equiv\tau_{\alpha\beta}\;$ satisfies

[TABLE]

and $\tau_{1},\tau_{2}$ are roots of the quadratic equation

[TABLE]

Solution of Problem 2 is given by the following corollary of Theorem  5.

Corollary**.**

Let (4.2) be supports (1.16) of extremal measures (1.15) for some fixed $s\in(0,1)$ of Angelesco system (4.1), and let $(u_{\alpha\beta},\tau_{\alpha\beta})$ be the images of transformations (4.4), (4.7). Then for limits (1.18) of the corresponding NNRR coefficients we have

[TABLE]

where parameters $\tau_{0}:=\tau_{\alpha\beta},\tau_{1},\tau_{2}$ are defined in (4.7), (4.8), and

[TABLE]

Formulas for $A_{2},B_{2}$ can be obtained by the swap of indices $1\to 2,2\to 1$.

Proof of this corollary is presented below in subsection 4.4.

**4.2.  **Parametrization of supports and ray directions and solution to Problem 1

Before we start dealing with Problem 1, let us come back to the parametrization (4.3) and consider $(u,\tau)$ on the half-strip $\bigsqcup:=(1,2)\times(1,\infty)$. If we invert map (4.4), (4.7), then we get a smooth diffeomorphism $(\mathcal{A},\mathcal{B}):\;\bigsqcup\ni(u,\tau)\mapsto(\alpha,\beta)\in(0,+\infty)\times(0,1)$:

[TABLE]

which by means of coordinates $(u,\tau)$ parametrize the branch points $\{-\alpha,\beta\}$ of the Riemann surface $\mathfrak{R}(\alpha,\beta)$, i.e., the left endpoints of segments (4.2) of supports of the extremal vector-measure $\overrightarrow{\omega}$.

In [13] there was introduced a parametrization of the direction $(t_{1},t_{2})$, see  (1.13), that corresponds to the masses for the extremal measures $(\omega_{1},\omega_{2})$ which have supports $[-\alpha,0],[\beta,1]$ when the vector equilibrium problem is formulated on $[-\alpha,0],[0,1]$. It is given by the function

[TABLE]

[TABLE]

Now we can deal with Problem 1. Without loss of generality (we make it clear below in subsection 4.3), it is enough to consider the Angelesco system on touching intervals ($\beta=0$):

[TABLE]

Problem 1 can be decomposed into two parts:

Problem 1.1. Given $\alpha$, find $s_{\alpha}\in(0,1)$ such that segments (4.13) are supports of the extremal measure of problem (1.15).

Problem 1.2. For fixed $s\in(s_{\alpha},1)$ find the value of $\beta_{s}$ so that:

[TABLE]

Solution of these problems is given in the following theorem.

Theorem 6** (for proof see [13]).**

Given $\alpha$ in (4.13):

1) Excluding variable $\tau$ from the system of equations

[TABLE]

we get the value of $\theta(\alpha)=:\theta_{\alpha}$. Then the answer to Problem 1.1 is $s_{\alpha}=\displaystyle\frac{1+\theta_{\alpha}}{2}$.

2) For each $s\in(s_{\alpha},1)$, let $\theta=2s-1\in(\theta_{\alpha},1)$. Then the system

[TABLE]

has a unique solution $(\tilde{u},\tilde{\tau})$, and $\beta_{s}:=\mathcal{B}(\tilde{u},\tilde{\tau})$ is the answer to Problem 1.2.

Summarizing, we have for $d=2$ the following Procedure for finding limits (2.2), (2.3) of NNRR coefficients

[TABLE]

for the Angelesco systems of MOPs (1.3), (1.12) defined on intervals (4.13).

1. Solve Problem 1.1: find $\theta_{\alpha},s_{\alpha}=(1+\theta_{\alpha})/2$. To do this, evaluate functions $\mathcal{A}(2,\tau),\Theta(2,\tau)$ by (4.11), (4.12), which determines the value and $\theta_{\alpha}$ according to Theorem 6 1).

2. For each $s\in(s_{\alpha},1)$ solve Problem 1.2: find $\beta_{s}$ from (4.14). To do this, solve the system from Theorem 6 2) for $\theta:=2s-1$ and substitute its solution $(\tilde{u},\tilde{\tau})$ into the function $\mathcal{B}$ to find $\beta_{s}$.

3. For each $s\in(s_{\alpha},1)$ find $A_{j}(s),B_{j}(s),\;j=1,2$. To do this, apply Corollary of Theorem 5 with the supports of the extremal measure being $[-\alpha,0]$ and $[\beta_{s},1]$, i.e., solve equations (4.4), (4.7), (4.8) and substitute the resulting $u_{\alpha\beta},\tau_{\alpha\beta},\tau_{1},\tau_{2}$ into the formulas  (4.9) for $A_{j}(s),B_{j}(s),\;j=1,2$.

To find limits (2.2), (2.3) for $s\in(0,s_{\alpha})$, we do the following:

4. We make reflection with respect to 0 and scaling (by $k=1/\alpha$) to get the system of intervals to the form (4.13). As a result, the new intervals are $[-\widehat{\alpha},0]$, $[0,1]$ with $\widehat{\alpha}=1/\alpha$.

5. We apply the above steps 1, 2, 3 of the Procedure to this new system of intervals to get the limits $\widehat{A}_{j}(s)$, $\widehat{B}_{j}(s)$, $j=1,2$ for $s\in(\widehat{s}_{\widehat{\alpha}},1)$ (note that $\widehat{s}_{\widehat{\alpha}}=1-s_{\alpha}$).

6. Then $A_{j}(s)=\widehat{A}(1-s)/k^{2}$ and $B_{j}(s)=-\widehat{B}(1-s)/k$ for $j=1,2$ and $s\in(0,s_{\alpha})$. Indeed, scaling by $k$ stretches all the $b_{\vec{n},j}$-coefficients by $k$ and all the $a_{\vec{n},j}$-coefficients by $k^{2}$. Reflection multiplies the $b_{\vec{n},j}$-coefficients by $-1$, keeps $a_{\vec{n},j}$’s intact and flips $s$ to $1-s$.

**4.3.  **Remark on Problem 1 for the

measures with non-touching supports

At first we provide an equivalent characterization of the extremal vector-measure $\vec{\omega}=(\omega_{1},\omega_{2})$ of the functional (1.15). We have (see [10]):

[TABLE]

where $V^{\nu}(z)=-\int\log|z-x|\,d\nu(x)$ is log-potential of measure $\nu$.

If we consider the vector potential

[TABLE]

where $\mathbb{A}$ is called the Angelesco matrix of interaction, then from (4.15) we can see that components of $\overrightarrow{W}$ possess the equilibrium property

[TABLE]

Thus the extremal measure $\vec{\omega}$ is also called the equilibrium measure.

Many properties of the equilibrium measure follow from equilibrium relations (4.15) and from the fact that log-potential is a convex function outside of the measure support. For example, for $d=2$ the strict inclusion $\text{\rm{supp}}\,\omega_{i}\subsetneq[\alpha_{i},\beta_{i}]$ may happen only for one component $i=1$ or $i=2$. As another example, if we have for fixed $s$ in (4.15)

[TABLE]

then $\text{\rm{supp}}\,\omega_{2}=[\alpha_{2}^{*},\beta_{2}]$ with $\alpha_{2}<\alpha_{2}^{*}$, and for this $s$ the extremal measure $\vec{\omega}$ is the same as for all Angelesco systems with supports

[TABLE]

Using this property we can reduce the solution of Problem 1 for the Angelesco systems with non-touching supports (4.2) to the case (4.13) considered above. Indeed, for the non-touching case we start with case (4.13) anyway, i.e., with intervals $[-\alpha,0]$ and $[0,1]$ and perform step 1 of the above Procedure: find $s_{\alpha}$. Then we perform a new step:

1.5. Find $s^{\beta}\in(s_{\alpha},1)$ such that for the Angelesco system supported by (4.13) we have333This can be done by executing step 3 of the Procedure for $s>s_{\alpha}$ until (4.16) happens.

[TABLE]

We note, that the obtained $s^{\beta}$ is equal $c_{2}$ from point 2) of Theorem 4:

[TABLE]

Then, performing steps 2 and 3 for $s\in(c_{2},1)$ we obtain $A_{j}(s),B_{j}(s),j=1,2$.

In an analogous way we obtain value of $s=c_{1}<c_{2}$ and $A_{j}(s),B_{j}(s),j=1,2$ for $s\in(0,c_{1})$. At the end we recall that for $s\in(c_{1},c_{2})$ limits $A_{j}(s),B_{j}(s),j=1,2$ remain to be the constants.

**4.4.  **Proof of Corollary of Theorem 5

From Theorem 5 we know that the function $\textbf{z}:\overline{\mathbb{C}}\to\mathfrak{R}_{s}:=\mathfrak{R}(\alpha,\beta)$ is a conformal map, where

[TABLE]

see (4.6). Meromorphic on $\mathfrak{R}_{s}$ function $\Upsilon_{1}$ is defined by its divisor and normalization:

[TABLE]

Our goal is to obtain two terms of the power series expansion of $\Upsilon_{1}(\textbf{z})$ at the point $\infty^{(0)}$, namely to find the coefficients $A_{1},B_{1}$ in

[TABLE]

In coordinates $w\in\overline{\mathbb{C}}$ we have

[TABLE]

where $C_{1}$ is determined from the normalization at the point $\infty^{(1)}$, see (4.18):

[TABLE]

For the coefficients $A_{1},B_{1}$ we have from (4.19)

[TABLE]

Thus substituting (4.20), (4.17) in (4.21) we obtain $C_{1}$ in (4.10):

[TABLE]

Analogously, plugging (4.20), (4.17) into (4.22) for $A_{1}$, we obtain (4.9):

[TABLE]

and plugging (4.20), (4.17), (4.23) into (4.22) for $B_{1}$, we get:

[TABLE]

Using the notation $P(w):=\dfrac{w^{2}(w-\gamma)}{(w-\tau_{1})^{2}(w-\tau_{2})}$, we continue:

[TABLE]

To compute $P^{\prime}(w)$ we use

[TABLE]

This allows us to arrive to (4.9):

[TABLE]

[TABLE]

Corollary of Theorem 5 is proved.

5.  Comparing numerics: Angelesco system $d=2$

**5.1.  **Numerics: two touching intervals

For the Angelesco systems with two intervals we now have three methods of numerically estimating the limits $A_{1}(s),A_{2}(s),B_{1}(s),B_{2}(s)$ ($0\leq s\leq 1$) of the NNRR’s coefficients:

• (i)

by computing $a_{\vec{n},j}$ and $b_{\vec{n},j}$ recursively (through (1.5)–(1.7), see [9]) for large enough $|\vec{n}|$;

• (ii)

through the system of ODE’s in Section 3 (namely, (3.5));

• (iii)

through the algebraically-geometric approach of Section 4.

On Fig. 1 we present the numerics in Wolfram Mathematica for the case $[\alpha_{1},\beta_{1}]=[-2,0]$, $[\alpha_{2},\beta_{2}]=[0,1]$. In (i) $|\vec{n}|$ was taken $1500$ (blue plot); in (ii) the in-built NDSolve Mathematica function was used (orange plot); notice that the ODE for $B_{1}$ in (3.4) has a singular behavior at $s=0$ and the same is true for $B_{2}$ at $s=1$, so one should use

[TABLE]

instead (these follow from  (3.4) and (3.15)); in (iii) the interval $s\in[0,1]$ was divided into $3000$ subintervals (red plot). The three plots are effectively indistinguishable.

**5.2.  **Numerics: two non-touching intervals

On Figure 2 we present the limits $A_{1}(s),A_{2}(s),B_{1}(s),B_{2}(s)$ for an Angelesco system with $[\alpha_{1},\beta_{1}]=[-2,0]$, $[\alpha_{2},\beta_{2}]=[0.25,1]$. The blue plot corresponds to the computation of $a_{\vec{n},j}$ and $b_{\vec{n},j}$ recursively (via  (1.5)–(1.7)) with $|\vec{n}|=1500$; the purple plot corresponds to the numerical approximation of the solution to the system of ODE’s (via (3.5)) with the boundary conditions at $s=0$; the green plot corresponds to the numerical approximation of the solution to the system of ODE’s (via (3.5)) with the boundary conditions at $s=1$. Equivalently, the purple plot corresponds to the coefficients’ limits for the Angelesco system with supports of $\mu_{1}$ and $\mu_{2}$ being $[-2,0.25]$ and $[0.25,1]$, while the green plot corresponds to the supports $[-2,0]$ and $[0,1]$. See Subsection 4.3 for the explanation of this phenomenon. This can also be seen from the fact that (3.6) is independent of $\beta_{1}$ and that  (3.7) is independent of $\alpha_{2}$. Again, the plots effectively overlap (away from the plateau regions).

Acknowledgments

The authors are grateful to the anonymous referees for their corrections and careful proof-reading of the paper. The second author thanks W. Van Assche for the excellent mini-course on multiple orthogonal polynomials at the Summer School on OPSF at the University of Kent, which lead to the idea of the current paper. He also thanks the organizers of the Summer School and W. Van Assche and A. Martínez-Finkelshtein for useful discussions.