A Riemann–Hilbert Approach to the Heun Equation

A Riemann–Hilbert Approach to the Heun Equation††This paper is a contribution
to the Special Issue on Painlevé Equations and Applications in Memory of
Andrei Kapaev. The full collection is available at
https://www.emis.de/journals/SIGMA/Kapaev.html

Boris DUBROVIN † and Andrei KAPAEV ‡

B. Dubrovin and A. Kapaev

† SISSA, Via Bonomea 265, 34136, Trieste, Italy<EMAIL_ADDRESS>
‡ Deceased

Received February 07, 2018, in final form August 15, 2018; Published online
September 07, 2018

We describe the close connection between the linear system for the sixth
Painlevé equation and the general Heun equation, formulate the Riemann–Hilbert
problem for the Heun functions and show how, in the case of reducible
monodromy, the Riemann–Hilbert formalism can be used to construct explicit
polynomial solutions of the Heun equation.

Heun polynomials; Riemann–Hilbert problem; Painlevé equations

34M03; 34M05; 34M35; 34M55; 57M50

##  1 Introduction

General Heun equation (GHE) [14] is the 2nd order linear ODE with four
distinct Fuchsian singularities depending on $6$ arbitrary complex parameters.
Without loss of generality, three of the singular points can be placed at $0$,
$1$ and $\infty$ while the position $a$ of the fourth singularity remains not
fixed. The canonical form of the general Heun equation (GHE) reads

$\displaystyle\frac{{\rm d}^{2}y}{{\rm
d}z^{2}}+\left(\frac{\gamma}{z}+\frac{\kappa}{z-1}+\frac{\epsilon}{z-a}\right)\frac{{\rm
d}y}{{\rm d}z}+\frac{\alpha\beta z-q}{z(z-1)(z-a)}y=0,$ (1.1)
$\displaystyle\gamma+\kappa+\epsilon=\alpha+\beta+1,\qquad a\neq 0,1,$

where the parameters $\alpha$, $\beta$, $\gamma$, $\kappa$, $\epsilon$
determine the characteristic exponents at the singular points,

$\displaystyle z=0\colon\quad$ $\displaystyle\\{0,1-\gamma\\},$ $\displaystyle
z=1\colon\quad$ $\displaystyle\\{0,1-\kappa\\},$ $\displaystyle
z=a\colon\quad$ $\displaystyle\\{0,1-\epsilon\\},$ $\displaystyle
z=\infty\colon\quad$ $\displaystyle\\{\alpha,\beta\\},$

while the remaining accessory parameter $q$ depends on global monodromy
properties of solutions to (1.1).

The general Heun equation together with its confluent and transformed
(trigonometric and elliptic) counterparts like Mathieu, spheroidal wave,
Leitner–Meixner, Lamé and Coulomb spheroidal equations finds numerous
applications in quantum and high energy physics, general relativity,
astrophysics, molecular physics, crystalline materials, 3d wave in atmosphere,
Bethe ansatz systems etc., see [1, 15, 25, 27] for a comprehensive but not
exhaustive list of references.

Importance of the Heun equation (1.1) is due to the fact that any Fuchsian
second order linear ODE with $4$ singular points can be reduced to (1.1) by
elementary transformations, while a trigonometric or elliptic change of the
independent variable yield linear ODEs with periodic and double periodic
coefficients, see [4, 23, 29].

The GHE is the classical example of the Fuchsian ODE that does not admit any
continuous isomonodromy deformation. To overcome the difficulty in the
construction of the isomonodromy problem, R. Fuchs [7, 8] added to four
conventional Fuchsian singularities one apparent singularity that is presented
in the equation but is absent in the solution. Position $y$ of this apparent
singularity is changed together with the position of the Fuchsian singularity
$x$ according to the second order nonlinear ODE known now as the sixth
Painlevé equation P${}_{\text{\rm\scriptsize VI}}$. In [5], it was observed
that under certain assumptions the apparent singularity disappears at the
critical values and movable poles of $y$, and the linear Fuchsian ODE turns
into the general Heun equation.

In the present paper, instead of the scalar second order differential equation
with apparent fifth singular point we shall use its first order $2\times 2$
matrix version with four Fuchsian singular points. We will give a detailed
proof that the general Heun equation appears at the poles of
P${}_{\text{\rm\scriptsize VI}}$, formulate the Riemann–Hilbert (RH) problem
for the general Heun equation and explore some implications of the RH problem
scheme to the Heun transcendents. In fact, we show how one can obtain the Heun
polynomials (i.e., polynomial solutions of (1.1)) within the suggested
Riemann–Hilbert formalism.

##  2 Reduction of the linear differential system for
P$\boldsymbol{{}_{\text{\rm\bf\scriptsize VI}}}$
to the general Heun equation (GHE)

###  2.1 Isomonodromy deformations of a Fuchsian linear ODE
with four singularities

The modern theory of the isomonodromy deformation was developed in the
pioneering work of M. Jimbo, T. Miwa, and K. Ueno [18, 19], although its
origin goes back to the classical papers of R. Fuchs [7], R. Garnier [9, 10],
and L. Schlesinger [24]. We shall briefly outline the theory in the case of
the $2\times 2$ matrix Fuchsian ODE with four singular points. The reader can
find more details in the works [18, 19] (see also Part 1 of the monograph
[5]).

The generic first order $2\times 2$ matrix Fuchsian ODE with four singular
points can be written in the form

$\displaystyle\frac{{\rm d}\Psi}{{\rm d}\lambda}=A(\lambda)\Psi,$ (2.1)

with the coefficient matrix

$\displaystyle
A(\lambda)=\frac{A_{1}}{\lambda}+\frac{A_{2}}{\lambda-x}+\frac{A_{3}}{\lambda-1},\qquad\sum_{j=1}^{3}A_{j}=-\delta\sigma_{3},\qquad\delta\neq
0,$ (2.2)

where $\sigma_{3}=\left(\begin{smallmatrix}1&0\\\
0&-1\end{smallmatrix}\right)$. Below we assume that

$\displaystyle\operatorname{Tr}A(\lambda)\equiv 0,$

and denote by

$\displaystyle\pm\alpha_{j},\qquad j=1,2,3,\qquad
2\alpha_{j}\notin{\mathbb{Z}}$

the eigenvalues of the matrix residues $A_{j}$.

The deformation,

$\displaystyle A(\lambda)\equiv A(\lambda,x),$

with respect to the position of the singularity $\lambda=x$ is isomonodromic
iff $\Psi(\lambda)\equiv\Psi(\lambda,x)$ satisfies an auxiliary linear ODE
with respect to this variable

$\displaystyle\frac{\partial\Psi}{\partial x}=B(\lambda)\Psi,\qquad
B(\lambda)=-\frac{A_{2}}{\lambda-x}.$

In [19], it is shown that the unique zero $y\equiv y(x)$ of the entry
$A_{12}(\lambda)$( which is a rational function whose numerator is linear in
$\lambda$) satisfies the classical sixth Painlevé equation
P${}_{\text{\rm\scriptsize VI}}$,

$\displaystyle
y_{xx}=\frac{1}{2}\left(\frac{1}{y}+\frac{1}{y-1}+\frac{1}{y-x}\right)y_{x}^{2}-\left(\frac{1}{x}+\frac{1}{x-1}+\frac{1}{y-x}\right)y_{x}$
$\displaystyle\hphantom{y_{xx}=}{}+\frac{y(y-1)(y-x)}{x^{2}(x-1)^{2}}\left\\{\alpha_{0}+\beta_{0}\frac{x}{y^{2}}+\gamma_{0}\frac{x-1}{(y-1)^{2}}+\delta_{0}\frac{x(x-1)}{(y-x)^{2}}\right\\},$

where

$\displaystyle\alpha_{0}=2\left(\delta-\frac{1}{2}\right)^{2},\qquad\beta_{0}=-2\alpha_{1}^{2},\qquad\gamma_{0}=2\alpha_{3}^{2},\qquad\delta_{0}=-2\left(\alpha_{2}^{2}-\frac{1}{4}\right).$

Suitable parameterization of the coefficient matrix $A(\lambda)$ of the
Fuchsian equation (2.2) and appearance of the sixth Painlevé equation
P${}_{\text{\rm\scriptsize VI}}$ satisfied by this zero $y(x)$ are explained
in more detail in Appendix A.1.

In [5, p. 86], it was observed that, at the critical values $y=0,1,x$ and
movable poles $y=\infty$, the linear matrix equation (2.2) with the
parametrization (A.1) becomes equivalent to the GHE. In the following
subsections, we describe the way in which the Heun equation emerges at the
movable poles of the Painlevé functions with more details than in [5].

###  2.2 Movable poles of P$\boldsymbol{{}_{\text{\rm\bf\scriptsize VI}}}$

If $\delta\neq\frac{1}{2}$, then equation P${}_{\text{\rm\scriptsize VI}}$
admits a 2-parameter family of solutions with the following leading terms [11]
of the Laurent expansion

$\displaystyle
y(x)=c_{-1}(x-a)^{-1}+c_{0}+c_{1}(x-a)+c_{2}(x-a)^{2}+{\mathcal{O}}\big{(}(x-a)^{3}\big{)},$
(2.3)

where $a\in{\mathbb{C}}\backslash\\{0,1\\}$, $c_{0}\in{\mathbb{C}}$ are
arbitrary, while all other coefficients are determined recursively by $a$,
$c_{0}$, $\sigma=\pm 1$ and the local monodromies $\alpha_{j}$, $j=1,2,3$,

$\displaystyle
c_{-1}=\sigma\frac{a(a-1)}{2\big{(}\delta-\frac{1}{2}\big{)}},\qquad\sigma\in\\{1,-1\\},$
$\displaystyle
c_{1}=1-\frac{c_{0}-\frac{1}{3}}{a}-\frac{c_{0}-\frac{2}{3}}{a-1}+\frac{\sigma}{3}\left(\delta-\frac{1}{2}\right)\left\\{1-\frac{1}{a}\big{(}6c_{0}^{2}-4c_{0}+1\big{)}+\frac{1}{a-1}\big{(}6c_{0}^{2}-8c_{0}+3\big{)}\right\\}$
$\displaystyle\hphantom{c_{1}=}{}+\frac{\sigma}{2\big{(}\delta-\frac{1}{2}\big{)}}\left\\{1-\frac{2}{3}\left(\alpha_{2}^{2}-\frac{1}{4}\right)+\frac{2}{3a}\left(\alpha_{3}^{2}-\frac{1}{4}\right)-\frac{2}{3(a-1)}\left(\alpha_{1}^{2}-\frac{1}{4}\right)\right\\},$

etc.

If $\delta=\frac{1}{2}$, then the movable poles of solutions to
P${}_{\text{\rm\scriptsize VI}}$ are double [11],

$\displaystyle
y(x)=c_{-2}(x-a)^{-2}+c_{-1}(x-a)^{-1}+c_{0}+{\mathcal{O}}(x-a),$ (2.4)

where $a\in{\mathbb{C}}\backslash\\{0,1\\}$,
$c_{-2}\in{\mathbb{C}}\backslash\\{0\\}$ are arbitrary, and all other
coefficients are determined recursively by $a$ and $c_{-2}$ and the local
monodromies $\alpha_{j}$, $j=1,2,3$,

$\displaystyle c_{-1}=\frac{2a-1}{a(a-1)}c_{-2},$ $\displaystyle
c_{0}=\frac{1}{3}(a+1)$
$\displaystyle\hphantom{c_{0}=}{}+\frac{c_{-2}}{12a^{2}(a-1)^{2}}\left(12a(a-1)+1-4a\alpha_{1}^{2}+4(a-1)\alpha_{3}^{2}-4a(a-1)\left(\alpha_{2}^{2}-\frac{1}{4}\right)\right),$

etc.

###  2.3 The coefficient matrix $\boldsymbol{A(\lambda)}$ at the movable poles
of P$\boldsymbol{{}_{\text{\rm\bf\scriptsize VI}}}$

In this subsection, our concern is the behavior of $A(\lambda)$ at the poles
of $y(x)$. We shall use the parameterization of $A(\lambda)$ given in (A.1).
As we will see, the coefficient matrix is continuous at the simple poles with
positive $\sigma$ (that is, $\delta\neq\frac{1}{2}$, $\sigma=+1$) and is
singular at the poles of any other kind. In the latter case, the linear ODE
can be regularized by a suitable Schlesinger transformation [19], and all
three resulting regular linear ODEs are equivalent to the GHE.

######  Theorem 2.1.

At any pole of a solution to P${}_{\text{\rm\scriptsize VI}}$, the associated
linear ODE (2.1) is equivalent $($in some cases after a suitable
regularization$)$ to GHE (1.1). Moreover, the pole position becomes the
position of the fourth singularity in GHE while the free parameter of the
Laurent expansion of the Painlevé function determines the accessory parameter
in GHE.

######  Remark 2.2.

The part of this statement concerning the relation of the free parameter of
the Laurent expansion of the Painlevé function and the accessory parameter in
GHE, in the case of $\delta\neq\frac{1}{2}$, has been already established in
[20]. In [20] the authors are using a very different approach based on the
discovered in [20] remarkable connection of the classical conformal blocks and
the sixth Painlevé equation. The authors of [20] also make use of the striking
fact (first observed in [26]) that the Heun equation can be thought of as the
quantization of the classical Hamiltonian of P${}_{\text{\rm\scriptsize VI}}$.

In Sections 2.3.1–2.3.5, we give the detailed proof of Theorem 2.1 considering
each case individually.

####  2.3.1 $\boldsymbol{\delta\neq 0,\frac{1}{2}}$, and
$\boldsymbol{\sigma=+1}$ (regular case)

In this case, the coefficient matrix remains continuous. In more details,
using (2.3) along with (A.2), one finds

$\displaystyle\kappa=\kappa_{0}(x-a)+{\mathcal{O}}\big{(}(x-a)^{2}\big{)},\qquad\kappa_{0}={\rm
const},\qquad x\to a$

(see equation (A.1) in Appendix A for the definition of the functions
$\kappa=\kappa(x)$ and $\tilde{\kappa}=\tilde{\kappa}(x)$). This zero cancels
the simple pole of $y(x)$, so the $(1,2)$-entry of the matrix $A(\lambda)$
remains bounded. Furthermore, the direct computer-aided computation yields the
continuity of all other entries as well,

$\displaystyle
A(\lambda)=\frac{a_{3}\lambda^{2}+b_{3}\lambda+c_{3}}{\lambda(\lambda-1)(\lambda-a)}\sigma_{3}+\frac{c_{+}}{\lambda(\lambda-1)(\lambda-a)}\sigma_{+}$
$\displaystyle\hphantom{A(\lambda)=}{}+\frac{b_{-}\lambda+c_{-}}{\lambda(\lambda-1)(\lambda-a)}\sigma_{-}+{\mathcal{O}}(x-a),\qquad
x\to a.$ (2.5)

Here and below $\sigma_{+}=\left(\begin{smallmatrix}0&1\\\
0&0\end{smallmatrix}\right)$, $\sigma_{-}=\left(\begin{smallmatrix}0&0\\\
1&0\end{smallmatrix}\right)$. The coefficients $a_{3}$, $b_{3}$, $c_{3}$,
$c_{+}$, $b_{-}$, $c_{-}$ are expressed in terms of the local monodromies
$\delta$, $\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$, the position of the pole
$a$ and the coefficient $c_{0}$ of the Laurent series (2.3). Complete details
can be found in Appendix A.2.

####  2.3.2 $\boldsymbol{\delta\neq 0,\frac{1}{2},1}$ and
$\boldsymbol{\sigma=-1}$ (generic singular case)

Now, the function $\kappa$ develops a pole as $x\to a$,

$\displaystyle\kappa=\kappa_{0}(x-a)^{-1}+{\mathcal{O}}(1),\qquad\kappa_{0}={\rm
const},\qquad x\to a,$

and the matrix $A(\lambda)$ becomes singular. To overcome this difficulty, one
can apply a suitable Schlesinger transformation [19]

$\displaystyle A(\lambda)\mapsto
R(\lambda)A(\lambda)R(\lambda)^{-1}+R_{\lambda}(\lambda)R(\lambda)^{-1},$

that regularizes the equation (2.2) as $x\to a$.

Assume first that $\delta\neq 1$. The needed Schlesinger transformation shifts
the formal monodromy at infinity $\delta$ by $-1$, i.e.,
$\delta\mapsto\delta-1$; it is given explicitly by

$\displaystyle
R_{0}(\lambda)=\begin{pmatrix}\lambda+\dfrac{1+x-2p-y(2\delta+1)}{2(\delta-1)}&-\dfrac{\kappa}{2\delta-1}\vspace{1mm}\\\
\dfrac{2\delta-1}{\kappa}&0\end{pmatrix},\qquad\delta\neq\frac{1}{2},1.$

It is straightforward to check that the transformed matrix $\hat{A}$ remains
regular at the pole $x=a$, and the limiting coefficient matrix is as follows

$\displaystyle\hat{A}=\big{(}R_{0}AR_{0}^{-1}+(R_{0})_{\lambda}R_{0}^{-1}\big{)}\bigr{|}_{x=a}$
$\displaystyle\hphantom{\hat{A}}{}=\frac{\hat{a}_{3}\lambda^{2}+\hat{b}_{3}\lambda+\hat{c}_{3}}{\lambda(\lambda-1)(\lambda-a)}\sigma_{3}+\frac{\hat{b}_{+}\lambda+\hat{c}_{+}}{\lambda(\lambda-1)(\lambda-a)}\sigma_{+}+\frac{\hat{c}_{-}}{\lambda(\lambda-1)(\lambda-a)}\sigma_{-},$
(2.6)

where expressions for the constant parameters $\hat{a}_{3}$, $\hat{b}_{3}$,
$\hat{c}_{3}$, $\hat{b}_{+}$, $\hat{c}_{+}$, $\hat{c}_{-}$ in terms of the
parameters $a$, $c_{0}$ and the local monodromies can be found in Appendix
A.3.

####  2.3.3 $\boldsymbol{\delta=1}$, $\boldsymbol{\sigma=-1}$ (the first
special singular case)

Let us proceed to the case $\delta=1$ and $\sigma=-1$. We choose the
Schlesinger transformation that changes the formal monodromy at infinity and
at the origin by one half

$\displaystyle\delta=1\mapsto\delta-\frac{1}{2}=\frac{1}{2},\qquad\alpha_{1}\mapsto\alpha_{1}-\frac{1}{2},$

and is given by the gauge matrix $R_{1}(\lambda)$,

$\displaystyle
R_{1}(\lambda)=\frac{1}{\sqrt{\lambda}}\begin{pmatrix}\lambda+g&-\kappa\\\
-\dfrac{g}{\kappa}&1\end{pmatrix},\qquad g=-p-y+\frac{z+\alpha_{1}x}{y}.$

The transformed coefficient matrix $\check{A}$ is regular at the pole $x=a$,
and its limiting value is as follows

$\displaystyle\check{A}=\big{(}R_{1}AR_{1}^{-1}+{R_{1}}_{\lambda}R_{1}^{-1}\big{)}\bigr{|}_{x=a}$
$\displaystyle\hphantom{\check{A}}{}=\frac{\check{a}_{3}\lambda^{2}+\check{b}_{3}\lambda+\check{c}_{3}}{\lambda(\lambda-1)(\lambda-a)}\sigma_{3}+\frac{\check{b}_{+}\lambda+\check{c}_{+}}{\lambda(\lambda-1)(\lambda-a)}\sigma_{+}+\frac{\check{b}_{-}}{(\lambda-1)(\lambda-a)}\sigma_{-},$
(2.7)

where the explicit expressions for the constant coefficients are given in
Appendix A.4.

####  2.3.4 $\boldsymbol{\delta=\frac{1}{2}}$

If $\delta=\frac{1}{2}$ then the Laurent expansion with the double pole (2.4)
implies that the coefficient matrix $A(\lambda)$ is singular at $x=a$. The
chosen regularizing Schlesinger transformation shifts $\delta\mapsto\delta+1$,

$\displaystyle R_{2}(\lambda)=\begin{pmatrix}0&-\dfrac{2}{\tilde{\kappa}}\\\
\dfrac{\tilde{\kappa}}{2}&\lambda+g_{2}\end{pmatrix},\qquad
g_{2}=-\frac{1}{3}(2p+2y-2\tilde{y}+x+1).$

The transformed matrix,
$\tilde{A}=R_{2}AR_{2}^{-1}+{R_{2}}_{\lambda}R_{2}^{-1}$, is regular at the
pole $x=a$ and, at this point, takes the value

$\displaystyle\tilde{A}(\lambda)=\frac{\tilde{a}_{3}\lambda^{2}+\tilde{b}_{3}\lambda+\tilde{c}_{3}}{\lambda(\lambda-1)(\lambda-a)}\sigma_{3}+\frac{\tilde{c}_{+}}{\lambda(\lambda-1)(\lambda-a)}\sigma_{+}+\frac{\tilde{b}_{-}\lambda+\tilde{c}_{-}}{\lambda(\lambda-1)(\lambda-a)}\sigma_{-},$
(2.8)

where the coefficients $\tilde{a}_{3}$, $\tilde{b}_{3}$, $\tilde{c}_{3}$,
$\tilde{c}_{+}$, $\tilde{b}_{-}$, $\tilde{c}_{-}$ are presented in Appendix
A.5.

####  2.3.5 GHE from the linear ODEs at the poles of
P$\boldsymbol{{}_{\text{\rm\bf\scriptsize VI}}}$

In this subsection, we show that all the linear matrix ODEs corresponding to
the poles of the sixth Painlevé function are equivalent to the GHE.

Observe that the coefficient matrices of the form (2.5) corresponding to
$\delta\neq 0,\frac{1}{2},1$ and $\sigma=+1$ as well as the regularized
coefficient matrix (2.8) corresponding to $\delta=\frac{1}{2}$ coincide with
each other modulo notations. Similarly, mutatis mutandis, the matrix (2.6) for
$\delta\neq 0,\frac{1}{2},1$, $\sigma=-1$ is the $\sigma_{1}$-conjugate of the
previous coefficient matrices. The coefficient matrix (2.7) for $\delta=1$,
$\sigma=-1$, is an inessential modification of (2.6). It is enough then to
consider the cases (2.5) and (2.7).

Consider first the coefficient matrix (2.5).

The first order matrix equation for the function $\Psi(\lambda)$ is always
equivalent to the second order Fuchsian ODE for the entry $\Psi_{1*}(\lambda)$
of the first row of the matrix function $\Psi(\lambda)$. However, extra
(apparent) singularities might appear in the process of excluding the entry
$\Psi_{2*}(\lambda)$. In the case of (2.5), however, the rational function
representing the $12$ – entry of matrix (2.5) does not have $\lambda$ in its
numerator. Hence, when the entry $\Psi_{2*}(\lambda)$ is excluded from the
system, no apparent singularities appear. Therefore, in the case (2.5), the
entry $\Psi_{1*}(\lambda)$ of the first row of the matrix function
$\Psi(\lambda)$ satisfies a linear 2nd order Fuchsian ODE with $4$ singular
points without any apparent singularity and therefore is equivalent to GHE. It
is, in fact, straightforward to check that the function

$\displaystyle
u(\lambda)=\lambda^{\alpha_{1}}(\lambda-a)^{\alpha_{2}}(\lambda-1)^{\alpha_{3}}\Psi_{1*}(\lambda)$
(2.9)

satisfies the general Heun equation in its canonical form (1.1)

$\displaystyle
u^{\prime\prime}+\left(\frac{1-2\alpha_{1}}{\lambda}+\frac{1-2\alpha_{2}}{\lambda-a}+\frac{1-2\alpha_{3}}{\lambda-1}\right)u^{\prime}+\frac{\mu\lambda+\nu}{\lambda(\lambda-a)(\lambda-1)}u=0,$
$\displaystyle\mu=(\alpha_{1}+\alpha_{2}+\alpha_{3}-\delta)(\alpha_{1}+\alpha_{2}+\alpha_{3}+\delta-2),$
(2.10)
$\displaystyle\nu=\alpha_{1}+\alpha_{2}-(\alpha_{1}+\alpha_{2})^{2}+\alpha_{3}^{2}-\delta^{2}+a(\alpha_{1}+\alpha_{3}-(\alpha_{1}+\alpha_{3})^{2}+\alpha_{2}^{2}-\delta^{2})+b_{3}(2\delta-1).$

Observe that the expression for the accessory parameter $\nu$ in (2.10)
besides the pole position and the local monodromies, involves the coefficient
$b_{3}$ in the parameterization of the entry $A_{11}$. Thus, taking into
account formula for $b_{3}$ in (A.3), we see that the accessory parameter
$\nu$ is determined by the free coefficient $c_{0}$ (or $c_{-2}$ in the case
(2.8 – see (A.4))) in the Laurent expansion of the sixth Painlevé
transcendent.

For the coefficient matrix (2.6), corresponding to $\delta\neq 1$, $\sigma=-1$
a similar statement is valid for the entries of the second row of
$\Psi(\lambda)$,

$\displaystyle
v(\lambda)=\lambda^{\alpha_{1}}(\lambda-a)^{\alpha_{2}}(\lambda-1)^{\alpha_{3}}\Psi_{2*}(\lambda),$
$\displaystyle
v^{\prime\prime}+\left(\frac{1-2\alpha_{1}}{\lambda}+\frac{1-2\alpha_{2}}{\lambda-a}+\frac{1-2\alpha_{3}}{\lambda-1}\right)v^{\prime}+\frac{\hat{\mu}\lambda+\hat{\nu}}{\lambda(\lambda-a)(\lambda-1)}v=0,$
$\displaystyle\hat{\mu}=(\alpha_{1}+\alpha_{2}+\alpha_{3}-\delta-1)(\alpha_{1}+\alpha_{2}+\alpha_{3}+\delta-1),$
$\displaystyle\hat{\nu}=\alpha_{1}+\alpha_{2}-(\alpha_{1}+\alpha_{2})^{2}+\alpha_{3}^{2}-(\delta-1)^{2}$
$\displaystyle\hphantom{\hat{\nu}=}{}+a\big{(}\alpha_{1}+\alpha_{3}-(\alpha_{1}+\alpha_{3})^{2}+\alpha_{2}^{2}-(\delta-1)^{2}\big{)}+\hat{b}_{3}(2\delta-1).$

In the case (2.7) corresponding to $\delta=1$, $\sigma=-1$, the function

$\displaystyle\check{v}(\lambda)=\lambda^{\alpha_{1}-\frac{1}{2}}(\lambda-a)^{\alpha_{2}}(\lambda-1)^{\alpha_{3}}\Psi_{2*}(\lambda)$

satisfies the following Heun equation

$\displaystyle\check{v}^{\prime\prime}+\left(\frac{1-2\alpha_{1}}{\lambda}+\frac{1-2\alpha_{2}}{\lambda-a}+\frac{1-2\alpha_{3}}{\lambda-1}\right)\check{v}^{\prime}+\frac{\check{\mu}\lambda+\check{\nu}}{\lambda(\lambda-1)(\lambda-a)}\check{v}=0,$
$\displaystyle\check{\mu}=(\alpha_{1}+\alpha_{2}+\alpha_{3}-2)(\alpha_{1}+\alpha_{2}+\alpha_{3}),$
$\displaystyle\check{\nu}=b_{3}-\frac{1}{2}(a+1)+\alpha_{1}+\alpha_{2}-(\alpha_{1}+\alpha_{2})^{2}+\alpha_{3}^{2}+a\big{(}\alpha_{1}+\alpha_{3}-(\alpha_{1}+\alpha_{3})^{2}+\alpha_{2}^{2}\big{)}.$

This completes the proof of Theorem 2.1.

##  3 Riemann–Hilbert problem approach to the Heun equation

Main result of this section is the formulation of the RH problem for the
general Heun functions in the generic case

$\displaystyle 2\alpha_{1},2\alpha_{2},2\alpha_{3},2\delta\notin{\mathbb{Z}}.$

We shall start, following closely references [5, 17], with the standard
definition of the monodromy data for Fuchsian system (2.1), (2.2) and with the
related Riemann–Hilbert problem for the sixth Painlevé equation.

###  3.1 Monodromy data

$\lambda_{1}$$E_{1}$$M_{1}$$\infty$$\lambda_{2}$$E_{2}$$M_{2}$$\lambda_{3}$$E_{3}$$M_{3}$$M_{\infty}$$\lambda_{0}$
Figure 1: The jump contour $\gamma$ for the RH Problem 3.1.

Let $\lambda_{1}=0$ and $\lambda_{3}=1$. Then fix a point
$\lambda_{2}=x\in{\mathbb{C}}\backslash\\{0,1,\infty\\}$, choose a base point
$\lambda_{0}\in{\mathbb{C}}\backslash\\{0,1,x,\infty\\}$ and cut the complex
plane along the segments

$\displaystyle[\lambda_{0},\lambda_{1}]\cup[\lambda_{0},\lambda_{2}]\cup[\lambda_{0},\lambda_{3}]\cup[\lambda_{0},\infty].$

Encircle the points $\lambda_{1}=0$, $\lambda_{2}=x$ and $\lambda_{3}=1$ using
non-intersecting circles ${\mathcal{C}}_{j}$, $j=1,2,3$. Denote $\gamma$ the
graph

$\displaystyle\gamma=[\lambda_{0},\lambda_{1}]\cup[\lambda_{0},\lambda_{2}]\cup[\lambda_{0},\lambda_{3}]\cup[\lambda_{0},\infty]\cup\cup_{j=1}^{3}{\mathcal{C}}_{j}$

and orient it as in Fig. 1. Denote also $D_{j}$, $j=1,2,3$ the interiors of
the circles ${\mathcal{C}}_{j}$ and $D_{\infty}$ the domain

$\displaystyle
D_{\infty}=\bigl{(}{\mathbb{C}}\backslash\gamma\bigr{)}\backslash\bigl{(}\cup_{j=1}^{3}D_{j}\bigr{)}$

The domains $D_{j}$ are assigned to the principal branches of the Frobenius
(canonical) solutions to (2.1) at the Fuchsian singular points defined by the
conditions

$\displaystyle\Psi_{j}(\lambda)=T_{j}(I+{\mathcal{O}}(\lambda-\lambda_{j}))(\lambda-\lambda_{j})^{\alpha_{j}\sigma_{3}},$
$\displaystyle\lambda\to\lambda_{j},\qquad\lambda\in D_{j},\qquad
j=1,2,3,\qquad\det T_{j}=1,$
$\displaystyle\Psi_{\infty}=\big{(}I+{\mathcal{O}}\big{(}\lambda^{-1}\big{)}\big{)}\lambda^{-\delta\sigma_{3}},\qquad\lambda\to\infty,\qquad\lambda\in
D_{\infty}.$ (3.1)

Here, the branches of $(\lambda-\lambda_{j})^{\alpha_{j}}$, $j=1,2,3$, and
$\lambda^{-\delta}$ are fixed by the condition

$\displaystyle\arg(\lambda-\lambda_{j})\to\pi,\qquad\arg\lambda\to\pi,\qquad\mbox{as}\quad\lambda\to-\infty.$

Given a pair of the characteristic exponents $(\alpha_{j},-\alpha_{j})$, the
matrix of eigenvectors $T_{j}\in{\rm SL}(2,{\mathbb{C}})$ is determined up to
a right diagonal factor. In contrast, the Frobenius solution
$\Psi_{\infty}(\lambda)$ is normalized and therefore, as soon as the pair
$(-\delta,\delta)$ is fixed, it is determined uniquely.

The matrices of the local monodromy are defined as the branch matrices of the
Frobenius solutions

$\displaystyle\Psi_{j}\big{(}\lambda_{j}+(\lambda-\lambda_{j}){\rm
e}^{2\pi{\rm i}}\big{)}=\Psi_{j}(\lambda){\rm e}^{2\pi{\rm
i}\alpha_{j}\sigma_{3}},\qquad\Psi_{\infty}\big{(}{\rm e}^{2\pi{\rm
i}}\lambda\big{)}=\Psi_{\infty}(\lambda){\rm e}^{-2\pi{\rm
i}\delta\sigma_{3}}.$

Introduce also the connection matrices between Frobenius solutions at infinity
and at the finite singularities

$\displaystyle\Psi_{j}(\lambda)=\Psi_{\infty}(\lambda)E_{j}.$

Similar to $\Psi_{j}(\lambda)$, the connection matrices $E_{j}$ are determined
modulo arbitrary right diagonal factors. In contrast, the monodromy matrices
$M_{j}$,

$\displaystyle M_{j}=E_{j}{\rm e}^{2\pi{\rm
i}\alpha_{j}\sigma_{3}}E_{j}^{-1},$ (3.2)

are determined uniquely. The monodromy matrices are the branching matrices of
the solution $\Psi_{\infty}(\lambda)$ at the singular points $\lambda_{j}$,
$j=1,2,3$; namely, one has that

$\displaystyle\Psi_{\infty}\big{(}\lambda_{j}+(\lambda-\lambda_{j}){\rm
e}^{2\pi{\rm i}}\big{)}=\Psi_{\infty}(\lambda)M_{j},\qquad j=1,2,3.$ (3.3)

Together with the matrix

$\displaystyle M_{\infty}:={\rm e}^{-2\pi{\rm i}\delta\sigma_{3}}$ (3.4)

they generate the monodromy group of equation (2.1)

$\displaystyle{\mathbb{M}}=\langle M_{1},M_{2},M_{3},M_{\infty}\rangle,\qquad
M_{j}\in{\rm SL}(2,{\mathbb{C}}),$

and are subject of one (cyclic) constraint

$\displaystyle M_{1}M_{2}M_{3}=M_{\infty}.$ (3.5)

Given the local monodromies, each of the monodromy matrices $M_{j}$ depends on
2 parameters. The total set of the 6 parameters determining the monodromy
matrices $M_{j}$, $j=1,2,3,\infty$, is subject to a system of 3 scalar
constraints. Thus the parameter set of the monodromy data involves generically
$3$ parameters. One of these parameters corresponds to the constant factor
$\kappa_{0}$ determining the auxiliary function $\kappa$ – see (A.1), (A.2).
This is a reflection of the possible conjugation of $A(\lambda)$ by a constant
diagonal matrix – the action which does not affect the zero of
$A_{12}(\lambda)$, i.e., the P${}_{\text{\rm\scriptsize VI}}$ function $y(x)$.
Neglecting this auxiliary parameter, the space of essential monodromy data,
${\mathcal{M}}$, is invariant with respect to an overall conjugation by a
diagonal matrix and can be identified with an algebraic variety – the
monodromy surface, of dimension 2 (see below equation (3.6)). At the same
time, the full space of monodromy data, can be represented as

$\displaystyle{\mathfrak{M}}={\mathcal{M}}\times{\mathbb{C}},\qquad\dim{\mathcal{M}}=2.$

In [17], M. Jimbo has proposed a parameterization of the 2-dimensional
monodromy surface ${\mathcal{M}}$ by the trace coordinates invariant with
respect to the overall diagonal conjugation. Namely, letting

$\displaystyle a_{j}=\operatorname{Tr}M_{j}=2\cos(2\pi\alpha_{j}),\qquad
j=1,2,3,\infty,\qquad\alpha_{\infty}=\delta,$ $\displaystyle
t_{ij}=\operatorname{Tr}(M_{i}M_{j})=2\cos(2\pi\sigma_{ij}),\qquad i,j=1,2,3,$

one finds the relation between all these parameters for a 2-dimensional
surface called the Fricke cubic

$\displaystyle
t_{12}t_{23}t_{31}+t_{12}^{2}+t_{23}^{2}+t_{31}^{2}-(a_{1}a_{2}+a_{3}a_{\infty})t_{12}-(a_{2}a_{3}+a_{1}a_{\infty})t_{23}-(a_{3}a_{1}+a_{2}a_{\infty})t_{31}$
$\displaystyle\qquad{}+a_{1}^{2}+a_{2}^{2}+a_{3}^{2}+a_{\infty}^{2}+a_{1}a_{2}a_{3}a_{\infty}-4=0.$
(3.6)

According to [16], apart from the singular points of the surface (3.6), the
monodromy matrices can be written explicitly in terms of the variables
$t_{ij}$. Exact formulae can be found in [17] and [16].

Each point of the surface $\mathcal{M}$ represents an isomonodromic family of
equations (2.1) which in turns generates a solution $y(x)$ of the Painlevé VI
equation. Hence, the P${}_{\text{\rm\scriptsize VI}}$ transcendents can be
parameterized by the points of $\mathcal{M}$. In fact, at generic points of
the Fricke cubic one can use any pair of the parameters $t_{ij}$ or
$\sigma_{ij}$ to parameterize the set of the corresponding Painlevé functions.
For instance, one can choose,

$\displaystyle t:=t_{12}=\operatorname{Tr}(M_{1}M_{2}),\qquad
s:=t_{1,3}=\operatorname{Tr}(M_{1}M_{3}),$ (3.7)

so that we have the parameterization of the P${}_{\text{\rm\scriptsize VI}}$
functions by the pair $(t,s)$,

$\displaystyle y\equiv y(x;t,s).$ (3.8)

It also should be mentioned that some of the physically important solutions,
e.g., the so-called classical solutions to P${}_{\text{\rm\scriptsize VI}}$,
correspond to non-generic points of the monodromy data set and for their
parameterization one can use the full monodromy space ${\mathfrak{M}}$. We
refer to [13] for more detail on this issue.

###  3.2 Riemann–Hilbert problem for P$\boldsymbol{{}_{\text{\rm\bf\scriptsize
VI}}}$

The inverse monodromy problem, i.e., the problem of reconstruction of the
function $\Psi$, and hence of the corresponding Painlevé function $y(x)$, from
their monodromy data is formulated as a Riemann–Hilbert (RH) problem. The
direct and inverse monodromy problems associated with the equation
P${}_{\text{\rm\scriptsize VI}}$ were studied by several authors. We mention
here the pioneering paper [17], and subsequent papers [2, 3, 13].

We shall now formulate precisely the Riemann–Hilbert problem corresponding to
the inverse monodromy problem for Fuchsian $2\times 2$ system (2.1).

######  Riemann–Hilbert Problem 3.1.

Given $x$, $\delta$, $\alpha_{j}$, $2\delta,2\alpha_{j}\notin{\mathbb{Z}}$,
$j=1,2,3$, the oriented graph $\gamma$ shown in Fig. 1 $($with
$\lambda_{1}=0$, $\lambda_{2}=x$, and $\lambda_{3}=1)$ and the jump matrices
$M_{j}$, $E_{j}$, $j=1,2,3$, all assigned to the branches of $\gamma$ and
satisfying the conditions (3.2), (3.4), (3.5), find a piecewise holomorphic
$2\times 2$ matrix function $\Psi(\lambda)$ with the following properties:

  1. $1)$

$\|\Psi(\lambda)\lambda^{\delta\sigma_{3}}-I\|\leq C|\lambda|^{-1}$ where $C$
is a constant and $\lambda\to\infty$,

  2. $2)$

$\|\Psi(\lambda)(\lambda-\lambda_{j})^{-\alpha_{j}\sigma_{3}}\|\leq C_{j}$ as
$\lambda\to\lambda_{j}$, where $C_{j}$ are some constants, $j=1,2,3$,

  3. $3)$

$\|\Psi(\lambda)\|\leq C_{0}$, where $C_{0}$ is a constant and $\lambda$
approaches any nodal point of the graph $\gamma$,

  4. $4)$

across the piecewise oriented contour $\gamma$, the discontinuity condition
holds,

$\displaystyle\Psi_{+}(\lambda)=\Psi_{-}(\lambda)G(\lambda),\qquad\lambda\in\gamma,$

where $\Psi_{+}(\lambda)$ and $\Psi_{-}(\lambda)$ are the left and right
limits of $\Psi(\lambda)$ as $\lambda$ transversally approaches the contour
$\gamma$, and $G(\lambda)$ is the piecewise constant matrix defined on
$\gamma$, see Fig. 1.

######  Proposition 3.2.

If a solution to the RH Problem 3.1 exists it is unique.

######  Proposition 3.3.

Having the canonical solutions $\Psi_{j}(\lambda)$ of (2.1), the equations

$\displaystyle\Psi(\lambda)=\Psi_{j}(\lambda),\qquad\lambda\in D_{j},\qquad
j=1,2,3,\infty,$ (3.9)

define the function which solves the RH Problem 3.1 whose data are determined
by the corresponding monodromy data.

######  Proposition 3.4.

Conversely, if for given $x$, $\delta$, $\alpha$, and matrices $M_{j}$,
$E_{j}$ the RH Problem 3.1 is solvable then the function $\Psi(\lambda)$
satisfies the Fuchsian system (2.1) whose Frobenius solutions are determined
by the solution $\Psi(\lambda)$ of the RH problem according to the equations
(3.9) $($read backwards$)$ and whose monodromy data coincide with the given RH
data.

The proofs of these propositions are standard, see, e.g., [5].

Assuming that the RH Problem 3.1 is solvable, the Painlevé function can be
extracted from the asymptotics of its solution at infinity. Indeed,
introducing the matrices
${{\mathcal{E}}_{k}}=\\{\delta_{ik}\delta_{jk}\\}_{i,j=1,2}$, $k=1,2$, by
straightforward computations we find the asymptotics of
$\Psi_{\infty}(\lambda)$,

$\displaystyle\Psi_{\infty}(\lambda)=\left(I+\frac{1}{\lambda}\psi_{1}+\frac{1}{\lambda^{2}}\psi_{2}+{\mathcal{O}}\left(\frac{1}{\lambda^{3}}\right)\right){\rm
e}^{\frac{1}{\lambda}d_{1}\sigma_{3}+\frac{1}{\lambda^{2}}(d_{21}{\mathcal{E}}_{1}+d_{22}{\mathcal{E}}_{2})}\lambda^{-\delta\sigma_{3}},\qquad\lambda\to\infty,$
(3.10)

where the coefficient matrices $\psi_{1}$ and $\psi_{2}$ are off-diagonal. The
expressions of $\psi_{k}$, $d_{kl}$ in terms of the coefficients of
$A(\lambda)$ (see (A.1)) can be found in Appendix A.6.

Using (A.5), namely, the $\sigma_{+}$-components $(\psi_{1})_{+}$ and
$(\psi_{2})_{+}$ of $\psi_{1}$ and $\psi_{2}$ respectively, and the scalars
$d_{1}$, $d_{21}$, $d_{22}$, we find

$\displaystyle\kappa=(2\delta-1)(\psi_{1})_{+},\qquad
p=-\delta(x+1)+2\delta\frac{(\delta-1)(\psi_{2})_{+}}{\big{(}\delta-\frac{1}{2}\big{)}(\psi_{1})_{+}}+d_{1}\frac{\delta+\frac{1}{2}}{\delta-\frac{1}{2}},$
$\displaystyle
y=x+1-\frac{(\delta-1)(\psi_{2})_{+}}{\big{(}\delta-\frac{1}{2}\big{)}(\psi_{1})_{+}}-\frac{d_{1}}{\delta-\frac{1}{2}},$
$\displaystyle z=-d_{21}+d_{22}-2\delta(d_{21}+d_{22})-\delta x$
$\displaystyle\hphantom{z=}{}-\left(\frac{\delta(\delta-1)(\psi_{2})_{+}}{\big{(}\delta-\frac{1}{2}\big{)}(\psi_{1})_{+}}+\frac{d_{1}}{2\big{(}\delta-\frac{1}{2}\big{)}}-\delta(x+1)\right)\left(\frac{(\delta-1)(\psi_{2})_{+}}{\big{(}\delta-\frac{1}{2}\big{)}(\psi_{1})_{+}}+\frac{d_{1}}{\delta-\frac{1}{2}}\right).$
(3.11)

###  3.3 Riemann–Hilbert problem for the Heun function

Main result of this section states that the RH problem for the Heun function
coincides with that for P${}_{\text{\rm\scriptsize VI}}$ supplemented by the
additional condition of triangularity of the sub-leading term of the
asymptotic expansion of $\Psi(\lambda)$ as $\lambda\to\infty$.

Again, we consider the non-resonant case $2\alpha_{1}$, $2\alpha_{2}$,
$2\alpha_{3}$, $2\delta\notin{\mathbb{Z}}$.

####  3.3.1 Limiting equation (2.2) and the $\boldsymbol{\Psi}$-function
at the pole $\boldsymbol{x=a}$ of $\boldsymbol{y(x)}$ as
$\boldsymbol{\delta\neq\frac{1}{2}}$ and $\boldsymbol{\sigma=+1}$

Consider the Laurent expansion (2.3) with $\sigma=+1$ and the corresponding
coefficient matrix (2.5). For $x$ in a punctured neighborhood of the pole
$x=a$, there exists an isomonodromy family of $\Psi$-functions. Furthermore,
the continuity of the coefficient matrix with respect to $x$ implies the
continuity of the Frobenius solutions (3.1) at $x=a$ as well.

On the other hand, the form of the coefficient matrix (2.5) with
$a_{3}=-\delta$ implies the following asymptotics of the solution to the
linear ODE $\Psi_{\lambda}=A\Psi$,

$\displaystyle\Psi(\lambda)=\left(I+\frac{1}{\lambda}\psi_{1}+\frac{1}{\lambda^{2}}\psi_{2}+{\mathcal{O}}\left(\frac{1}{\lambda^{3}}\right)\right){\rm
e}^{\frac{1}{\lambda}d_{1}\sigma_{3}+\frac{1}{\lambda^{2}}d_{2}\sigma_{3}}\lambda^{-\delta\sigma_{3}},\qquad\lambda\to\infty,$
(3.12)

where the main difference from the asymptotic parameters in (3.10), (A.5) is
the lower-triangular structure of the coefficient $\psi_{1}$,

$\displaystyle\psi_{1}=-\frac{b_{-}}{2\delta+1}\sigma_{-},\qquad
d_{1}=-b_{3}+\delta(a+1),$
$\displaystyle\psi_{2}=\frac{c_{+}}{2(\delta-1)}\sigma_{+}-\frac{b_{-}(a+1+2b_{3})+c_{-}(2\delta+1)}{4(\delta+1)\big{(}\delta+\frac{1}{2}\big{)}}\sigma_{-},$
$\displaystyle
d_{2}=\frac{1}{2}\bigl{(}-b_{3}(a+1)-c_{3}+\delta\big{(}1+a+a^{2}\big{)}\bigr{)}.$

We point out that the lower-triangular structure of $\psi_{1}$ implies the
lower-triangular structure of the
${\mathcal{O}}\big{(}\lambda^{-1}\big{)}$-term in the expansion (3.12).

All other principal analytic properties of the limiting function
$\Psi(\lambda)$ including the leading order asymptotics at the singular points
and the monodromy properties coincide with those of the function
$\Psi(\lambda)$ at the regular points of the Painlevé transcendent located in
a sufficiently small neighborhood of the pole $x=a$.

Thus we have shown that the RH problem for the function $\Psi(\lambda)$ at the
pole of the sixth Painlevé transcendent $y(x)$ as $\delta\neq\frac{1}{2}$ and
$\sigma=+1$, and therefore the RH problem for a solution of the general Heun
equation (2.10), coincides with the RH Problem 3.1 supplemented by the
condition of the lower triangularity of the coefficient $\psi_{1}$ at
infinity:

######  Riemann–Hilbert Problem 3.5.

Given $x$, $\delta$, $\alpha_{j}$, $2\delta,2\alpha_{j}\notin{\mathbb{Z}}$,
$j=1,2,3$, the oriented graph $\gamma$ shown in Fig. 1 $($with
$\lambda_{1}=0$, $\lambda_{2}=a$, and $\lambda_{3}=1)$ and the jump matrices
$M_{j}$, $E_{j}$, $j=1,2,3$, all assigned to the branches of $\gamma$ and
satisfying the conditions (3.2), (3.4), (3.5), find a piecewise holomorphic
$2\times 2$ matrix function $\Psi(\lambda)$ with the following properties:

  1. $1)$

$\|\Psi(\lambda)\lambda^{\delta\sigma_{3}}-I\|\leq C|\lambda|^{-1}$, $C={\rm
const}$, as $\lambda\to\infty$,

  2. $2)$

$\lim\limits_{\lambda\to\infty}\lambda(\Psi(\lambda)\lambda^{\delta\sigma_{3}}-I)_{12}=0$,

  3. $3)$

$\|\Psi(\lambda)(\lambda-\lambda_{j})^{-\alpha_{j}\sigma_{3}}\|\leq C_{j}$,
$C_{j}={\rm const}$, as $\lambda\to\lambda_{j}$, where $\lambda_{1}=0$,
$\lambda_{2}=a$, $\lambda_{3}=1$,

  4. $4)$

$\|\Psi(\lambda)\|\leq C_{0}$, $C_{0}={\rm const}$, as $\lambda$ approaches
any nodal point of the graph $\gamma$,

  5. $5)$

across the oriented contour $\gamma$, the discontinuity condition holds,

$\displaystyle\Psi_{+}(\lambda)=\Psi_{-}(\lambda)G(\lambda),\qquad\lambda\in\gamma,$

where $\Psi_{+}(\lambda)$ and $\Psi_{-}(\lambda)$ are the left and right
continuous limits of $\Psi(\lambda)$ as $\lambda$ approaches the contour
$\gamma$, and $G(\lambda)$ is the piecewise constant matrix defined on
$\gamma$, see Fig. 1.

Let us show that, conversely, this RH problem leads to the structure (2.5) of
the coefficient matrix $A(\lambda)$. Let $\Psi(\lambda)$ be a unique solution
of the RH Problem 3.5. First, $\det\Psi(\lambda)\equiv 1$ since this
determinant is piecewise holomorphic, continuous across the graph $\gamma$,
bounded at $\lambda_{j}$, $j=1,2,3$, and at the nodes of the graph $\gamma$
and approaches the unit as $\lambda\to\infty$. Consider now the function
$A(\lambda)=\Psi_{\lambda}\Psi^{-1}$. It is piecewise holomorphic, continuous
across $\gamma$, has simple poles at $\lambda=\lambda_{j}$, $j=1,2,3$, and
$\lambda=\infty$, is bounded at the nodes of $\gamma$ and therefore it is
rational.

The conditions (1) and (2) of the RH Problem 3.5 imply the asymptotics

$\displaystyle\Psi(\lambda)=\bigl{(}I+\psi_{1}\lambda^{-1}+{\mathcal{O}}\big{(}\lambda^{-2}\big{)}\bigr{)}\lambda^{-\delta\sigma_{3}},\qquad\lambda\to\infty,$

where $(\psi_{1})_{12}=0$. Therefore

$\displaystyle
A(\lambda)=-\frac{\delta}{\lambda}\sigma_{3}-\frac{1}{\lambda^{2}}\bigl{(}\psi_{1}+[\psi_{1},\sigma_{3}]\bigr{)}+{\mathcal{O}}\big{(}\lambda^{-3}\big{)},\qquad\lambda\to\infty,$

i.e.,

$\displaystyle(A(\lambda))_{12}={\mathcal{O}}\big{(}\lambda^{-3}\big{)}.$

All these properties of $A(\lambda)$ imply its structure given in (2.5). Thus,
using (2.9), the first row of $\Psi(\lambda)$ determines a fundamental system
of solutions to GHE (2.10).

We have proved the following

######  Proposition 3.6.

Solution of the RH Problem 3.5, if it exists, determines a fundamental system
of solutions to GHE (2.10) with the prescribed monodromy properties.

The accessory parameter $\nu$ in (2.10) is also determined via the RH Problem
3.5. Indeed, $\nu$ can be extracted from the asymptotics of $\Psi(\lambda)$ at
infinity. Namely, the parameter $d_{1}$, i.e., the diagonal part of the term
${\mathcal{O}}\big{(}\lambda^{-1}\big{)}$ of the asymptotic expansion of
$\Psi\lambda^{\delta\sigma_{3}}$, determines the coefficient $b_{3}$ in the
coefficient matrix (2.5) and hence the free coefficient $c_{0}$ in the Laurent
expansion (2.3) and the accessory parameter $\nu$,

$\displaystyle b_{3}=-d_{1}+\delta(a+1),\qquad
c_{0}=\frac{b_{3}+a-1}{2\delta-1}=\frac{-d_{1}+a(\delta+1)+\delta-1}{2\delta-1},$
$\displaystyle\nu=-d_{1}(2\delta-1)+\alpha_{1}+\alpha_{2}-(\alpha_{1}+\alpha_{2})^{2}+\alpha_{3}^{2}-\delta^{2}$
$\displaystyle\hphantom{\nu=}{}+a\big{(}\alpha_{1}+\alpha_{3}-(\alpha_{1}+\alpha_{3})^{2}+\alpha_{2}^{2}-\delta^{2}\big{)}+\delta(2\delta-1)(a+1).$
(3.13)

######  Remark 3.7.

The condition (2) of the RH Problem 3.5 can be in fact thought of as an
addition to cyclic relation (3.5) restriction on the monodromy matrices
$\\{M_{j}\\}$ which would also involve the point $a$. This restriction can be
formulated in terms of the sixth Painlevé transcendent in two different but
equivalent ways. Firstly, introducing on the monodromy surface $\mathcal{M}$
coordinates $t$ and $s$ (see (3.7)), let $y(x)\equiv y(x;t,s)$ be the
corresponding sixth Painlevé function (cf. (3.8)). Then the point $a$ must be
one of the $(\sigma=+1)$ poles of $y(x;t,s)$. Alternatively, assuming the
position $a$ of the pole of $y(x)$ to be a free parameter, one can
parameterize $y(x)$ by the pair111This parameterization of $y(x)$ is more
subtle than by the monodromy pair $(t,s)$. Indeed, P${}_{\text{\rm\scriptsize
VI}}$ function might have infinitely many poles and therefore one might have
an infinite discrete set of the pairs $(t,a_{n})$ corresponding to the same
P${}_{\text{\rm\scriptsize VI}}$ function $y(x)\equiv y(x;t,s)$. Hence, one
has to be careful describing the global properties of the parameterization,
$y(x)\equiv y(x;t,a)$. Also, the function $s(t,a)$ implicitly defined by
equation (3.14) below might have infinitely many branches. This means that for
each pair $(t,a)$ there might be an infinite discrete set of the values of the
second monodromy coordinate, $s$, and hence an infinite discrete set of the
P${}_{\text{\rm\scriptsize VI}}$ functions $y(x)$ with the same values of $t$
and $a$. The discussion of these important issues is, however, beyond the
scope of this work. $(t,a)$, $y(x)\equiv y(x;t,a)$. Then, the second monodromy
data $s$ becomes the function of $(t,a)$ which can be described implicitly as
follows. Note that together with $s$, the coefficient $c_{0}$ in the Laurent
expansion (2.3) and the accessory parameter $\nu$ also become the functions of
$t$ and $a$,

$\displaystyle s\equiv s(t,a),\qquad c_{0}\equiv
c_{0}(t,a),\qquad\nu\equiv\nu(t,a).$

At the same time, the RH Problem 3.1 determines the solution $\Psi(\lambda)$
and hence all the objects related to it, specifically the coefficients of its
expansion (3.10) at $\lambda=\infty$, as the functions of the point on the
monodromy surface, that is as the functions of $t$ and $s$. In particular, we
have that $d_{1}=d_{1}(x;t,s)$. Using now the second equation in (3.13), the
function $s(t,a)$ can be defined implicitly via the equation

$\displaystyle
c_{0}(t,a)=\frac{-d_{1}(a;t,s)+a(\delta+1)+\delta-1}{2\delta-1}.$ (3.14)

######  Remark 3.8.

Excluding the parameter $d_{1}$ from the second and third equations in (3.13),
we obtain the formula relating the accessory parameter $\nu(t,a)$ and the free
coefficient $c_{0}(t,a)$ in the Laurent expansion (2.3),

$\displaystyle\nu(t,a)=(1-2\delta)^{2}c_{0}(t,a)+3\delta^{2}(1+a)$
$\displaystyle\hphantom{\nu(t,a)=}{}+a\big{(}\alpha_{1}+\alpha_{3}-(\alpha_{1}+\alpha_{3})^{2}+\alpha_{2}^{2}\big{)}+\alpha_{1}+\alpha_{2}-(\alpha_{1}+\alpha_{2})^{2}+\alpha_{3}^{2}.$
(3.15)

As it has already been mentioned before, this relation has been already found
in [20] using a heuristic technique based on the remarkable connection (also
discovered in [20]) between the classical conformal blocks and the sixth
Painlevé equation.

######  Remark 3.9.

In this section we considered the regular case, i.e., $\delta\neq
0,\frac{1}{2}$ and $\sigma=+1$ only. Our arguments, however, can be easily
extended to the generic singular case, i.e., $\delta\neq 0,\frac{1}{2},1$ and
$\sigma=-1$, and to the special singular cases, i.e., $\delta=1$ and
$\sigma=-1$ and $\delta=\frac{1}{2}$. One only needs, before producing the
relevant analogs of the RH problem 3.5, to make the preliminary Schlesinger
transformations of the RH Problem 3.1 with the gauge matrices $R(\lambda)$
discussed in details in Sections 2.3.2, 2.3.3 and 2.3.4.

######  Remark 3.10.

In this paper we are dealing with the poles of $y(x)$. Similar results
concerning the reduction of the RH Problem 3.1 to a Riemann–Hilbert problem
for Heun equation can be obtained for two other critical values of $y(x)$,
i.e., when $a$ is either a zero of $y(x)$ or $y(a)=1$.

We believe that the Riemann–Hilbert technique we are developing here can be
used to study effectively the Heun functions. In particular, the relation
(3.15) allows one to get nontrivial information about the accessory parameter
$\nu(t,a)$. In particular, using the known connection formulae for the pole
distributions of the P${}_{\text{\rm\scriptsize VI}}$ equation [12] (obtained
with the help of the isomonodromy RH Problem 3.1; see also [13] for complete
list of the asymptotic connection formulae and the history of the question),
one can obtain the asymptotic expansion of the Laurent coefficient
$c_{0}(t,a)$ for the either large values of $a$ or for small values of $a$ or
for the values of $a$ close to 1. This in turn would yield the explicit
formulae for the asymptotic behavior of the accessory parameter $\nu(t,a)$ as
$a\sim\infty$, $a\sim 0$ and $a\sim 1$. This question should be addressed in
details in the future work on this subject. In the rest of this paper, we will
demonstrate the usefulness of the RH Problem 3.5 in the study of another issue
related to the Heun equations which is the construction of its explicit
solutions.

###  3.4 Example: reducible monodromy, generalized Jacobi
and Heun polynomials

Consider the case of reducible monodromy when all the monodromy matrices are
upper triangular. Reducible monodromy for P${}_{\text{\rm\scriptsize VI}}$
system was studied, e.g., in [22] where it was shown that this class of
Painlevé functions contains classical and all rational solutions to
P${}_{\text{\rm\scriptsize VI}}$.

In [5], the relevant function $\Psi(\lambda)$ was constructed explicitly for
$0\leq\operatorname{Re}\alpha_{j}<\frac{1}{2}$, $j=1,2,3$, and
$\alpha_{1}+\alpha_{2}+\alpha_{3}+\delta=0$. Below, we use the approach of [5]
assuming that

$\displaystyle\operatorname{Re}\alpha_{j}\in\big{[}0,\tfrac{1}{2}\big{)},\qquad\alpha_{1}+\alpha_{2}+\alpha_{3}+\delta=-n,\qquad
n\in{\mathbb{N}}.$ (3.16)

Let all the monodromy matrices $M_{j}$, $j=1,2,3$, be upper triangular and
thus each of them depends on one free parameter $s_{j}$,

$\displaystyle M_{j}=\begin{pmatrix}{\rm e}^{-2\pi{\rm i}\alpha_{j}}&s_{j}\\\
0&{\rm e}^{2\pi{\rm i}\alpha_{j}}\end{pmatrix},\qquad s_{j}\neq 0,\qquad
j=1,2,3.$ (3.17)

The cyclic relation $M_{1}M_{2}M_{3}=M_{\infty}$ implies that

$\displaystyle\alpha_{1}+\alpha_{2}+\alpha_{3}+\delta=-n\in{\mathbb{Z}},\qquad
s_{1}{\rm e}^{2\pi{\rm i}\alpha_{2}}+s_{2}{\rm e}^{-2\pi{\rm
i}\alpha_{1}}+s_{3}{\rm e}^{2\pi{\rm i}\delta}=0.$

Thus the set of reducible monodromy data form a 2-dimensional linear space.

Following [5], we first simplify the RH jump graph replacing $\gamma$ by the
broken line
$[\lambda_{1},\lambda_{2}]\cup[\lambda_{2},\lambda_{3}]\cup[\lambda_{3},\infty)$,
see Fig. 2.

$\lambda_{1}$$M_{1}$$\lambda_{2}$$M_{1}M_{2}$$\lambda_{3}$$M_{\infty}$$\infty$
Figure 2: Simplified jump contour for RH Problem 3.1.

On the plane, make a cut along the broken line
$[\lambda_{1},\lambda_{2}]\cup[\lambda_{2},\lambda_{3}]\cup[\lambda_{3},\infty)$,
define the function

$\displaystyle
f(\lambda)=(\lambda-\lambda_{1})^{\alpha_{1}}(\lambda-\lambda_{2})^{\alpha_{2}}(\lambda-\lambda_{3})^{\alpha_{3}}.$

Although until we reach Proposition 3.14 it is not really important, we remind
that

$\displaystyle\lambda_{1}=0,\quad\lambda_{2}=x,\qquad\mbox{and}\qquad\lambda_{3}=1.$

Observe the following properties of $f(\lambda)$:

$\displaystyle
f(\lambda)=\lambda^{-\delta-n}\big{(}1+{\mathcal{O}}\big{(}\lambda^{-1}\big{)}\big{)},\qquad\lambda\to\infty,$
$\displaystyle f_{+}(\lambda)=f_{-}(\lambda){\rm e}^{-2\pi{\rm
i}\alpha_{1}},\qquad\lambda\in(\lambda_{1},\lambda_{2}),$ $\displaystyle
f_{+}(\lambda)=f_{-}(\lambda){\rm e}^{-2\pi{\rm
i}(\alpha_{1}+\alpha_{2})},\qquad\lambda\in(\lambda_{2},\lambda_{3}),$
$\displaystyle f_{+}(\lambda)=f_{-}(\lambda){\rm e}^{2\pi{\rm
i}\delta},\qquad\lambda\in(\lambda_{3},\infty).$

Let us represent the solution $\Psi(\lambda)$ as the product

$\displaystyle\Psi(\lambda)=\Phi(\lambda)f^{\sigma_{3}}(\lambda).$

The function $\Phi(\lambda)$ thus has the following properties:

######  Riemann–Hilbert Problem 3.11.

  1. $1)$

$\Phi(\lambda)=\big{(}I+{\mathcal{O}}\big{(}\lambda^{-1}\big{)}\big{)}\lambda^{n\sigma_{3}}$
as $\lambda\to\infty$,

  2. $2)$

$\|\Phi(\lambda)f^{\sigma_{3}}(\lambda)E_{j}f^{-\sigma_{3}}(\lambda)\|\leq C$
as $\lambda\to\lambda_{j}$, $j=1,2,3$,

  3. $3)$

$\Phi_{+}(\lambda)=\Phi_{-}(\lambda)G_{f}(\lambda)$,
$\lambda\in(\lambda_{1},\lambda_{2})\cup(\lambda_{2},\lambda_{3})$, moreover

$\displaystyle G_{f}(\lambda)=I+g(\lambda)\sigma_{+},\qquad
g(\lambda)=\begin{cases}s_{1}f_{+}(\lambda)f_{-}(\lambda),&\lambda\in(\lambda_{1},\lambda_{2}),\\\
-s_{3}{\rm e}^{2\pi{\rm
i}\delta}f_{+}(\lambda)f_{-}(\lambda),&\lambda\in(\lambda_{2},\lambda_{3}).\end{cases}$
(3.18)

######  Proposition 3.12.

Solution to the RH Problem 3.11 is given by polynomials orthogonal with
respect to the weight function $g(\lambda)$ on
$(\lambda_{1},\lambda_{2})\cup(\lambda_{2},\lambda_{3})$ and by their Cauchy
integrals.

###### Proof.

First of all, by the conventional arguments, if a solution to this problem
exists, it is unique.

Next, upper triangularity of all jump matrices, see condition (3), means that
the first column of $\Phi(\lambda)$ is single-valued and continuous across the
broken line $(\lambda_{1},\lambda_{2})\cup(\lambda_{2},\lambda_{3})$.
Condition (2) means that the first column is bounded at $\lambda=\lambda_{j}$,
$j=1,2,3$. Then condition (1) yields that the entries of the first column are
some polynomials of degree $n$ and $n-1$, respectively.

Denote the relevant monic polynomials of degree $n$ and $n-1$ by
$\pi_{n}(\lambda)$ and $\pi_{n-1}(\lambda)$, respectively. Consider an
auxiliary matrix function $Y(\lambda)$ (cf. [6])

$\displaystyle
Y(\lambda)=\begin{pmatrix}\pi_{n}(\lambda)&\displaystyle\frac{1}{2\pi{\rm
i}}\int_{\ell}\pi_{n}(\zeta)g(\zeta)\frac{{\rm
d}\zeta}{\zeta-\lambda}\vspace{1mm}\\\
c_{n-1}\pi_{n-1}(\lambda)&\displaystyle\frac{c_{n-1}}{2\pi{\rm
i}}\int_{\ell}\pi_{n-1}(\zeta)g(\zeta)\frac{{\rm
d}\zeta}{\zeta-\lambda}\end{pmatrix},$

where $\ell=(\lambda_{1},\lambda_{2})\cup(\lambda_{2},\lambda_{3})$ and the
constant $c_{n-1}$ is defined by

$\displaystyle c_{n-1}=-\frac{2\pi{\rm
i}}{\int_{\ell}\pi_{n-1}^{2}(\zeta)g(\zeta){\rm d}\zeta}.$

As it is easy to see, $Y(\lambda)$ satisfies the jump condition (3).

At the points $\lambda=\lambda_{j}$, $j=1,2,3$, the function $Y(\lambda)$ has
the algebraic branching, $(\lambda-\lambda_{j})^{2\alpha_{j}}$, and therefore,
taking into account assumption (3.16), condition (2) is satisfied.

Finally, consider the series expansion of the entries of the second column of
$Y(\lambda)$ as $|\lambda|$ is large enough,

$\displaystyle\frac{1}{2\pi{\rm i}}\int_{\ell}\pi_{m}(\zeta)g(\zeta)\frac{{\rm
d}\zeta}{\zeta-\lambda}=-\frac{1}{2\pi{\rm
i}}\sum_{k=0}^{\infty}\lambda^{-k-1}\int_{\ell}\pi_{m}(\zeta)\zeta^{k}g(\zeta){\rm
d}\zeta,$

$m=n$ or $m=n-1$. Then condition (1) means that the polynomials
$p_{n}(\lambda)$ and $p_{n-1}(\lambda)$ are both orthogonal to all lower
degree monomials $\lambda^{k}$ for $k=0,1,\dots,n-1$ and $k=0,1,\dots,n-2$,
respectively,

$\displaystyle\int_{\ell}\pi_{n}(\zeta)\zeta^{k}g(\zeta){\rm d}\zeta=0,\qquad
k=0,1,\dots,n-1,$
$\displaystyle\int_{\ell}\pi_{n-1}(\zeta)\zeta^{k}g(\zeta){\rm
d}\zeta=0,\qquad k=0,1,\dots,n-2.$

The choice of $c_{n-1}$ made above implies the normalization of $Y(\lambda)$
at infinity

$\displaystyle
Y(\lambda)=\big{(}I+{\mathcal{O}}\big{(}\lambda^{-1}\big{)}\big{)}\lambda^{n\sigma_{3}}$

that complies with the condition (1). Thus the explicitly constructed function
$Y(\lambda)$ solves the RH Problem 3.11. Since its solution is unique, proof
of proposition is completed. ∎

The weight function $g(\lambda)$ in (3.18) can be understood as a
generalization of the hypergeometric weight, thus we call the polynomials
$p_{n}(\lambda)$ and $p_{n-1}(\lambda)$ determined by the RH Problem 3.11 the
generalized Jacobi polynomials.

Now, we are going to construct the generalized Jacobi polynomials explicitly
and relate them to the polynomial solutions of the Heun equation, i.e., to the
Heun polynomials.

To this end, we look for solution to the RH Problem 3.11 in the form

$\displaystyle\Phi(\lambda)=R(\lambda)\begin{pmatrix}1&\phi(\lambda)\\\
0&1\end{pmatrix},$ (3.19)

where $R(\lambda)$ is a matrix-valued polynomial. Using the jump properties of
$\Phi(\lambda)$, we find the jump properties of the scalar function
$\phi(\lambda)$,

$\displaystyle\phi_{+}(\lambda)-\phi_{-}(\lambda)=g(\lambda),\qquad\lambda\in(\lambda_{1},\lambda_{2})\cup(\lambda_{2},\lambda_{3}).$

One of the solutions to this scalar jump problem is given explicitly

$\displaystyle\phi(\lambda)=\frac{1}{2\pi{\rm
i}}\int_{\lambda_{1}}^{\lambda_{3}}\frac{g(\zeta)}{\zeta-\lambda}{\rm
d}\zeta.$ (3.20)

Observe the behavior of $\phi(\lambda)$ (3.20) at the singularities

$\displaystyle\phi(\lambda)={\mathcal{O}}\big{(}\lambda^{-1}\big{)},\qquad\lambda\to\infty,$
$\displaystyle\phi(\lambda)={\mathcal{O}}\big{(}(\lambda-\lambda_{j})^{2\alpha_{j}}\big{)}+{\mathcal{O}}(1),\qquad\lambda\to\lambda_{j},\qquad
j=1,2,3.$

Below, we use the coefficients $\phi_{k}$ of the expansion of $\phi(\lambda)$
(3.20) near infinity

$\displaystyle\phi(\lambda)=\frac{1}{2\pi{\rm
i}}\int_{\lambda_{1}}^{\lambda_{3}}\frac{g(\zeta)}{\zeta-\lambda}{\rm
d}\zeta=\sum_{k=1}^{\infty}\phi_{k}\lambda^{-k},$

closely related to the moments of the weight function $g(\lambda)$,

$\displaystyle\phi_{k}=-\frac{1}{2\pi{\rm
i}}\int_{\lambda_{1}}^{\lambda_{3}}g(\zeta)\zeta^{k-1}{\rm d}\zeta,\qquad
k\in{\mathbb{N}}.$ (3.21)

The left factor $R(\lambda)$ is a polynomial matrix of the Schlesinger
transformation at infinity [19] and for $n\geq 0$ it can be found explicitly.
For instance,

$\displaystyle\text{if}\ n=0\colon\quad
R(\lambda)=R_{0}(\lambda)=I,\qquad\pi_{0}(\lambda)\equiv
1,\qquad\pi_{-1}(\lambda)\equiv 0,$ $\displaystyle\text{if}\quad
n=1\colon\quad
R(\lambda)=R_{1}(\lambda)=\begin{pmatrix}\lambda-\dfrac{\phi_{2}}{\phi_{1}}&-\phi_{1}\vspace{1mm}\\\
\dfrac{1}{\phi_{1}}&0\end{pmatrix},\qquad\pi_{1}(\lambda)=\lambda-\frac{\phi_{2}}{\phi_{1}}.$
(3.22)

In particular, relations (3.22) imply that the RH Problem 3.11 for $n=0$ is
always solvable, cf. [5], while for $n=1$, this problem is solvable if
$\phi_{1}\neq 0$.

To formulate the result for any fixed $n\geq 2$, introduce the following
explicit form of the polynomial matrix $R(\lambda)=R_{n}(\lambda)$,

$\displaystyle
R_{n}(\lambda)=\begin{pmatrix}\displaystyle\sum_{k=0}^{n}p_{k}^{(n)}\lambda^{k}&\displaystyle\sum_{k=0}^{n-1}q_{k}^{(n)}\lambda^{k}\vspace{1mm}\\\
\displaystyle\sum_{k=0}^{n-1}r_{k}^{(n)}\lambda^{k}&\displaystyle\sum_{k=0}^{n-2}s_{k}^{(n)}\lambda^{k}\end{pmatrix},\qquad
p_{n}^{(n)}=1,$ (3.23)

and the column vectors of the coefficients of the polynomial entries of
$R_{n}(\lambda)$,

$\displaystyle{\bf p}_{n}=\begin{pmatrix}p_{0}^{(n)}\\\ \vdots\\\
p_{n-1}^{(n)}\end{pmatrix},\qquad{\bf q}_{n}=\begin{pmatrix}q_{0}^{(n)}\\\
\vdots\\\ q_{n-1}^{(n)}\end{pmatrix},\qquad{\bf
r}_{n}=\begin{pmatrix}r_{0}^{(n)}\\\ \vdots\\\
r_{n-1}^{(n)}\end{pmatrix},\qquad{\bf s}_{n}=\begin{pmatrix}s_{0}^{(n)}\\\
\vdots\\\ s_{n-2}^{(n)}\end{pmatrix}.$

We also need the Hankel matrix
${\mathcal{H}}_{n}=\\{\phi_{i+j-1}\\}_{i,j=1,n}$ of the moments (3.21) along
with its determinant $\triangle_{n}$,

$\displaystyle{\mathcal{H}}_{n}=\begin{pmatrix}\phi_{1}&\phi_{2}&\cdots&\phi_{n}\\\
\phi_{2}&\phi_{3}&\cdots&\phi_{n+1}\\\ \vdots&\vdots&\ddots&\vdots\\\
\phi_{n}&\phi_{n+1}&\cdots&\phi_{2n-1}\\\
\end{pmatrix},\qquad\triangle_{n}=\det{\mathcal{H}}_{n}.$

Finally define the sequence of coefficients $f_{k}$ for asymptotic expansion
of $f(\lambda)=\prod\limits_{j=1,2,3}(\lambda-\lambda_{j})^{\alpha_{j}}$ at
infinity

$\displaystyle f(\lambda)=\lambda^{-\delta-n}{\rm
e}^{\sum\limits_{k=1}^{\infty}f_{k}\lambda^{-k}},\qquad
f_{k}=-\frac{1}{k}\sum_{j=1}^{3}\alpha_{j}\lambda_{j}^{k},\qquad\lambda\to\infty.$

Then we have the following

######  Proposition 3.13.

RH Problem 3.11 is solvable if and only if $\triangle_{n}\neq 0$. Its solution
has the form (3.19) where the coefficients of the polynomial entries of the
matrix $R(\lambda)$ (3.23) are given by

$\displaystyle{\bf p}_{n}=-{\mathcal{H}}_{n}^{-1}\begin{pmatrix}\phi_{n+1}\\\
\vdots\\\ \phi_{2n}\end{pmatrix},\qquad p_{n}^{(n)}=1,$ $\displaystyle
q_{k}^{(n)}=-\sum_{m=k+1}^{n}\phi_{m-k}p_{m}^{(n)},\qquad k=0,\dots,n-1,$
(3.24) $\displaystyle{\bf r}_{n}={\mathcal{H}}_{n}^{-1}\begin{pmatrix}0\\\
\vdots\\\ 0\\\ 1\end{pmatrix},\qquad r_{n-1}^{(n)}=c_{n-1},\qquad
s_{k}^{(n)}=-\sum_{m=k+1}^{n-1}\phi_{m-k}r_{m}^{(n)},\qquad k=0,\dots,n-2.$

###### Proof.

Requiring that the product (3.19) has the canonical asymptotics

$\displaystyle R_{n}(\lambda)\begin{pmatrix}1&\phi(\lambda)\\\
0&1\end{pmatrix}=\begin{pmatrix}\lambda^{n}+{\mathcal{O}}\big{(}\lambda^{n-1}\big{)}&{\mathcal{O}}\big{(}\lambda^{-n-1}\big{)}\\\
{\mathcal{O}}\big{(}\lambda^{n-1}\big{)}&\lambda^{-n}+{\mathcal{O}}\big{(}\lambda^{-n-1}\big{)}\end{pmatrix},$

we find the following expansions in the second column of the product

$\displaystyle(12)\colon\quad\sum_{k=0}^{n-1}q_{k}^{(n)}\lambda^{k}+\sum_{m=0}^{n}\sum_{l=1}^{\infty}\phi_{l}p_{m}^{(n)}\lambda^{m-l}$
$\displaystyle\qquad\quad\
{}=\sum_{k=0}^{n-1}\lambda^{k}\left(q_{k}^{(n)}+\sum_{m=k+1}^{n}\phi_{m-k}p_{m}^{(n)}\right)+\sum_{k=1}^{\infty}\lambda^{-k}\sum_{m=0}^{n}\phi_{k+m}p_{m}^{(n)},$
$\displaystyle(22)\colon\quad\sum_{k=0}^{n-2}s_{k}^{(n)}\lambda^{k}+\sum_{m=0}^{n-1}\sum_{l=1}^{\infty}\phi_{l}r_{m}^{(n)}\lambda^{m-l}$
$\displaystyle\qquad\quad\
{}=\sum_{k=0}^{n-2}\lambda^{k}\left(s_{k}^{(n)}+\sum_{m=k+1}^{n-1}\phi_{m-k}r_{m}^{(n)}\right)+\sum_{k=1}^{\infty}\lambda^{-k}\sum_{m=0}^{n-1}\phi_{k+m}r_{m}^{(n)}.$
(3.25)

Thus expressions for $q_{k}^{(n)}$ and $s_{k}^{(n)}$ in (3.26) in terms of
$p_{m}^{(n)}$ and $r_{m}^{(n)}$ come from (3.25) at the non-negative degrees
of $\lambda$.

Evaluating the prescribed terms at all orders from $\lambda^{-1}$ to
$\lambda^{-n}$, we find equations for the vector coefficients ${\bf r}_{n}$
and ${\bf p}_{n}$. Combining them into the matrix form, it follows

$\displaystyle{\mathcal{H}}_{n}{\bf
p}_{n}+(\phi_{n+1},\phi_{n+2},\dots,\phi_{2n})^{\rm
T}=0,\qquad{\mathcal{H}}_{n}{\bf r}_{n}=(0,\dots,0,1)^{\rm T}.$ (3.26)

If $\triangle_{n}\neq 0$, the moment matrix ${\mathcal{H}}_{n}$ is invertible,
and the coefficient vectors ${\bf p}_{n}$, ${\bf r}_{n}$ are computed by
(3.24).

If $\triangle_{n}=0$ then for the equations (3.26) on the vectors ${\bf
p}_{n}$ and ${\bf r}_{n}$, there are two alternatives. The first one implies
that $R(\lambda)$ does not exist. The second alternative implies an infinite
number of $R(\lambda)$ and therefore contradicts the uniqueness of solution to
the RH Problem 3.11.

This completes the proof. ∎

In the next proposition, we find the classical sixth Painlevé functions
determined by the explicitly solvable RH Problem 3.11. For instance, in the
simplest cases,

$\displaystyle\text{if}\ n=0\colon\quad
y=x+1-\frac{(\delta-1)\phi_{2}}{\big{(}\delta-\frac{1}{2}\big{)}\phi_{1}}-\frac{f_{1}}{\delta-\frac{1}{2}},$
$\displaystyle\text{if}\ n=1\colon\quad
y=x+1-\frac{\delta-1}{\delta-\frac{1}{2}}\frac{\left|\begin{matrix}\phi_{1}&\phi_{2}\\\
\phi_{3}&\phi_{4}\end{matrix}\right|}{\left|\begin{matrix}\phi_{1}&\phi_{2}\\\
\phi_{2}&\phi_{3}\end{matrix}\right|}+\frac{\delta}{\delta-\frac{1}{2}}\frac{\phi_{2}}{\phi_{1}}-\frac{f_{1}}{\delta-\frac{1}{2}}.$

######  Proposition 3.14.

If $\triangle_{n}\triangle_{n+1}\neq 0$ then the RH Problem 3.11 determines
the classical solution to P${}_{\text{\rm\scriptsize VI}}$ corresponding to
the $\Psi$ function with the prescribed reducible monodromy data in terms of
the moment functions $\phi_{k}$, $k=1,\dots,\phi_{2n+2}$ as follows

$\displaystyle
y=x+1-\frac{f_{1}}{\delta-\frac{1}{2}}-\frac{\delta-1}{\delta-\frac{1}{2}}\frac{\triangle_{n}}{\triangle_{n+1}}\phi_{2n+2}+\frac{\delta-1}{\delta-\frac{1}{2}}\frac{\triangle_{n}}{\triangle_{n+1}}\begin{pmatrix}\phi_{n+2}&\dots&\phi_{2n+1}\end{pmatrix}{\mathcal{H}}_{n}^{-1}\begin{pmatrix}\phi_{n+1}\\\
\vdots\\\ \phi_{2n}\end{pmatrix}$
$\displaystyle\hphantom{y=}{}+\frac{\delta}{\delta-\frac{1}{2}}\begin{pmatrix}\phi_{n+1}&\dots&\phi_{2n}\end{pmatrix}{\mathcal{H}}_{n}^{-1}\begin{pmatrix}0\\\
\vdots\\\ 0\\\ 1\end{pmatrix}.$ (3.27)

###### Proof.

According to (3.11), in order to find the Painlevé function $y$ we have to
compute two first coefficients $(\psi_{1})_{+}$ and $(\psi_{2})_{+}$ along
with the coefficient $d_{1}$ of the asymptotic expansion of $\Psi(\lambda)$ at
$\lambda\to\infty$, see (3.10). Using (3.25) and (3.26) again, we find

$\displaystyle\Psi(\lambda)=R_{n}(\lambda)\begin{pmatrix}1&\phi(\lambda)\\\
0&1\end{pmatrix}f^{\sigma_{3}}(\lambda)$
$\displaystyle\hphantom{\Psi(\lambda)}{}=\begin{pmatrix}\displaystyle\sum_{k=0}^{n}\lambda^{-k}p_{n-k}^{(n)}&\displaystyle\sum_{k=1}^{\infty}\lambda^{-k}\sum_{m=0}^{n}\phi_{k+n+m}p_{m}^{(n)}\vspace{1mm}\\\
\displaystyle\sum_{k=1}^{n}\lambda^{-k}r_{n-k}^{(n)}&\displaystyle
1+\sum_{k=1}^{\infty}\lambda^{-k}\sum_{m=0}^{n-1}\phi_{k+n+m}r_{m}^{(n)}\end{pmatrix}{\rm
e}^{\sum\limits_{k=1}^{\infty}\lambda^{-k}f_{k}\sigma_{3}}\lambda^{-\delta\sigma_{3}}$
$\displaystyle\hphantom{\Psi(\lambda)}{}=\begin{pmatrix}1&\dfrac{1}{\lambda}(\psi_{1})_{+}+\dfrac{1}{\lambda^{2}}(\psi_{2})_{+}+{\mathcal{O}}\left(\dfrac{1}{\lambda^{3}}\right)\vspace{1mm}\\\
\dfrac{1}{\lambda}r_{n-1}^{(n)}+{\mathcal{O}}\left(\dfrac{1}{\lambda^{2}}\right)&1\end{pmatrix}$
$\displaystyle\hphantom{\Psi(\lambda)=}{}\times\left(I+{\mathcal{O}}\left(\frac{1}{\lambda^{2}}\right)\right){\rm
e}^{\frac{1}{\lambda}(f_{1}+p_{n-1}^{(n)})\sigma_{3}}\lambda^{-\delta\sigma_{3}},$

where

$\displaystyle(\psi_{1})_{+}=\sum_{m=0}^{n}\phi_{n+m+1}p_{m}^{(n)}=-\begin{pmatrix}\phi_{n+1}&\phi_{n+2}&\cdots&\phi_{2n}\end{pmatrix}{\mathcal{H}}_{n}^{-1}\begin{pmatrix}\phi_{n+1}\\\
\phi_{n+2}\\\ \vdots\\\
\phi_{2n}\end{pmatrix}+\phi_{2n+1}=\frac{\triangle_{n+1}}{\triangle_{n}},$
$\displaystyle(\psi_{2})_{+}=\sum_{m=0}^{n}\phi_{n+m+2}p_{m}^{(n)}-\sum_{m=0}^{n-1}\phi_{n+m+1}r_{m}^{(n)}\sum_{k=0}^{n}\phi_{n+k+1}p_{k}^{(n)}$
$\displaystyle\hphantom{(\psi_{2})_{+}}{}=\phi_{2n+2}+\begin{pmatrix}\phi_{n+2}&\phi_{n+3}&\dots&\phi_{2n+1}\end{pmatrix}{\bf
p}_{n}-\phi_{2n+1}\begin{pmatrix}\phi_{n+1}&\phi_{n+2}&\dots&\phi_{2n}\end{pmatrix}{\bf
r}_{n}$
$\displaystyle\hphantom{(\psi_{2})_{+}=}{}-\begin{pmatrix}\phi_{n+1}&\phi_{n+2}&\dots&\phi_{2n}\end{pmatrix}{\bf
r}_{n}\cdot\begin{pmatrix}\phi_{n+1}&\phi_{n+2}&\dots&\phi_{2n}\end{pmatrix}{\bf
p}_{n}$
$\displaystyle\hphantom{(\psi_{2})_{+}}{}=\phi_{2n+2}-\begin{pmatrix}\phi_{n+2}&\phi_{n+3}&\dots&\phi_{2n+1}\end{pmatrix}{\mathcal{H}}_{n}^{-1}\begin{pmatrix}\phi_{n+1}\\\
\vdots\\\ \phi_{2n}\end{pmatrix}$
$\displaystyle\hphantom{(\psi_{2})_{+}=}{}-\frac{\triangle_{n+1}}{\triangle_{n}}\begin{pmatrix}\phi_{n+1}&\phi_{n+2}&\dots&\phi_{2n}\end{pmatrix}{\mathcal{H}}_{n}^{-1}\begin{pmatrix}0\\\
\vdots\\\ 0\\\ 1\end{pmatrix},$ $\displaystyle
p_{n-1}^{(n)}=-\begin{pmatrix}\phi_{n+1}&\dots&\phi_{2n}\end{pmatrix}{\mathcal{H}}_{n}^{-1}\begin{pmatrix}0\\\
\vdots\\\ 0\\\ 1\end{pmatrix}.$ (3.28)

Finally, elementary manipulations yield (3.27). ∎

######  Proposition 3.15.

If $\triangle_{n}\neq 0$ and $\triangle_{n+1}=0$, then the RH Problem 3.11
determines a polynomial solution of the Heun equation $($the Heun
polynomial$)$ of degree $n$.

###### Proof.

The coefficient $(\psi_{1})_{+}$ is computed in (3.28),

$\displaystyle(\psi_{1})_{+}=\frac{\triangle_{n+1}}{\triangle_{n}}.$ (3.29)

The RH Problem 3.11 transforms to the Heun RH Problem 3.5 if
$(\psi_{1})_{+}=0$ – condition (2) of the RH Problem 3.5. In virtue of (3.29)
this happens only if $\triangle_{n+1}=0$. ∎

The equation

$\displaystyle\triangle_{n+1}=0$ (3.30)

is an equation on the variable $\lambda_{2}=x$. This equation determines those
$a=x$ for which the Heun equation with the given $\alpha_{j}$, $\delta$
satisfying (3.16) admits polynomial solutions given by $\Phi_{11}(\lambda)$.
Alternatively, one can keep $x=a$ free, then (3.30) would be equation on the
parameters $s_{1}$ and $s_{3}$ of the monodromy matrices (3.17). This in turn
would yield a restriction on the accessory parameter of the Heun equation
which would guarantee the existence of its polynomial solutions for given $a$.

######  Remark 3.16.

Determinantal conditions on the accessory parameter yielding existence of the
polynomial solution of the Heun equation has been known. It would be of course
interesting to compare those with the condition which comes from (3.30).

######  Remark 3.17.

There is a vast literature devoted to the exact solutions of the Heun equation
– see work [21] and references therein. We also want to mention the works [28]
and [30]. In [28], the indicated at the beginning of this paper relation of
the Heun equation and the Lame equation has been used to obtain the solvable
Heun equations and their polynomial solutions from Lame polynomials
corresponding to the finite-gap Lame potentials. In [30], the hypergeometric
type integral representations for the solutions of the Heun equation are found
in the case when one of the singularities $\\{0,1,t,\infty\\}$ is apparent. It
would be important to understand these explicit solutions within the
Riemann–Hilbert formalism which we are developing in this article.

##  Appendix A Parameterizations of the coefficient matrix for
P$\boldsymbol{{}_{\text{\rm\bf\scriptsize VI}}}$ and GHE

###  A.1 Parametrization of $\boldsymbol{A(\lambda)}$ for
P$\boldsymbol{{}_{\text{\rm\bf\scriptsize VI}}}$

Introducing ${\mathfrak{s}}{\mathfrak{l}}_{2}({\mathbb{C}})$ generators
$\sigma_{3}=\left(\begin{smallmatrix}1&0\\\ 0&-1\end{smallmatrix}\right)$,
$\sigma_{+}=\left(\begin{smallmatrix}0&1\\\ 0&0\end{smallmatrix}\right)$,
$\sigma_{-}=\left(\begin{smallmatrix}0&0\\\ 1&0\end{smallmatrix}\right)$, and
taking into account the simplifying assumption that the rational coefficient
matrix $A(\lambda)$ is diagonal at infinity, consider the parameterization
similar to the one used in [19]

$\displaystyle
A(\lambda)=\frac{-\delta(\lambda-y)^{2}+p(\lambda-y)+z}{\lambda(\lambda-1)(\lambda-x)}\sigma_{3}+\frac{\kappa(\lambda-y)}{\lambda(\lambda-1)(\lambda-x)}\sigma_{+}+\frac{\tilde{\kappa}(\lambda-\tilde{y})}{\lambda(\lambda-1)(\lambda-x)}\sigma_{-}.$
(A.1)

The residue matrices $A_{j}$ in (2.2) have the eigenvalues $\pm\alpha_{j}$,
$j=1,2,3$ respectively.

As soon as $\delta\neq 0$, the parameters $p$, $\tilde{y}$ and
$\tilde{\kappa}$ can be expressed in terms of the parameters $x$, $y$, $z$,
$\kappa$, $\delta$ and the eigenvalues $\alpha_{j}$, $j=1,2,3$,

$\displaystyle p=\frac{\alpha_{1}^{2}x}{2\delta
y}-\frac{\alpha_{3}^{2}(x-1)}{2\delta(y-1)}+\frac{\alpha_{2}^{2}x(x-1)}{2\delta(y-x)}-\frac{1}{2}\delta(3y-x-1)-\frac{z^{2}}{2\delta
y(y-1)(y-x)},$
$\displaystyle\tilde{y}=(y-1)\big{(}\alpha_{1}^{2}x^{2}-\big{(}\delta
y^{2}+py-z\big{)}^{2}\big{)}\bigl{\\{}\alpha_{1}^{2}x^{2}(y-1)-\alpha_{3}^{2}(x-1)^{2}y+(\delta-p)^{2}y-4\delta^{2}y^{2}$
$\displaystyle\hphantom{\tilde{y}=}{}+6\delta
py^{2}-p^{2}y^{2}+6\delta^{2}y^{3}-4\delta py^{3}-3\delta^{2}y^{4}-2\delta
yz+2\delta y^{2}z+z^{2}\bigr{\\}}^{-1},$
$\displaystyle\tilde{\kappa}=\frac{\alpha_{1}^{2}x^{2}-\big{(}\delta
y^{2}+py-z\big{)}^{2}}{\kappa y\tilde{y}}.$

Considering $y$, $z$ and $\kappa$ as the differentiable functions of $x$, they
satisfy the system

$\displaystyle\frac{{\rm d}(\ln\kappa)}{{\rm
d}x}=(2\delta-1)\frac{y-x}{x(x-1)},$ $\displaystyle\frac{{\rm d}y}{{\rm
d}x}=\frac{y^{2}-y+2z}{x(x-1)},$ $\displaystyle\frac{{\rm d}z}{{\rm
d}x}=\frac{1}{x(x-1)y(y-1)(y-x)}\bigl{\\{}z^{2}\big{(}x-2y-2xy+3y^{2}\big{)}+zy(y-1)(y-x)(y+x-1)$
$\displaystyle\hphantom{\frac{{\rm d}z}{{\rm
d}x}=}{}-\alpha_{1}^{2}x(y-1)^{2}(y-x)^{2}+\alpha_{3}^{2}(x-1)y^{2}(y-x)^{2}-\alpha_{2}^{2}x(x-1)y^{2}(y-1)^{2}$
$\displaystyle\hphantom{\frac{{\rm d}z}{{\rm
d}x}=}{}+\delta(\delta-1)y^{2}(y-1)^{2}(y-x)^{2}\bigr{\\}}.$ (A.2)

Finally, eliminating $z$, one arrives at the second order ODE for $y$

$\displaystyle
y_{xx}=\frac{1}{2}\left(\frac{1}{y}+\frac{1}{y-1}+\frac{1}{y-x}\right)y_{x}^{2}-\left(\frac{1}{x}+\frac{1}{x-1}+\frac{1}{y-x}\right)y_{x}$
$\displaystyle\hphantom{y_{xx}=}{}+\frac{y(y-1)(y-x)}{x^{2}(x-1)^{2}}\left\\{2\left(\delta-\frac{1}{2}\right)^{2}-2\alpha_{1}^{2}\frac{x}{y^{2}}+2\alpha_{3}^{2}\frac{x-1}{(y-1)^{2}}-2\left(\alpha_{2}^{2}-\frac{1}{4}\right)\frac{x(x-1)}{(y-x)^{2}}\right\\},$

which is the classical equation P${}_{\text{\rm\scriptsize VI}}$.

###  A.2 Parameterization of $\boldsymbol{A(\lambda)}$ at the pole of
$\boldsymbol{y(x)}$ as $\boldsymbol{\delta\neq 0,\frac{1}{2}}$, and
$\boldsymbol{\sigma=+1}$

The coefficients of the matrix $A(\lambda)$ (2.5) are expressed explicitly in
terms of the pole position $a$, free parameter $c_{0}$ of the Laurent
expansion at the pole and the local monodromies $\alpha_{j}$, $j=1,2,3$, as
follows

$\displaystyle a_{3}=-\delta,\qquad b_{3}=c_{0}(2\delta-1)+1-a,$
$\displaystyle
c_{3}=\frac{1}{2\delta}\bigl{[}(\delta-1-c_{0}(2\delta-1))^{2}-a\alpha_{1}^{2}+(a-1)\alpha_{3}^{2}-a(a-1)\alpha_{2}^{2}$
$\displaystyle\hphantom{c_{3}=}{}+a\big{(}\delta^{2}-2-2c_{0}\big{(}2\delta^{2}+\delta-1\big{)}\big{)}+a^{2}(\delta+1)^{2}\bigr{]},$
$\displaystyle c_{+}=-\kappa_{0}\frac{a(a-1)}{2\delta-1},$ $\displaystyle
b_{-}=\frac{2\delta-1}{\kappa_{0}\delta
a(a-1)}\bigl{\\{}a^{3}\big{(}\alpha_{2}^{2}-(1+\delta)^{2}\big{)}+a^{2}\bigl{[}3-2\alpha_{2}^{2}-3c_{0}+\alpha_{2}^{2}c_{0}+2\delta+\alpha_{2}^{2}\delta$
$\displaystyle\hphantom{b_{-}=}{}+2c_{0}\delta-2\alpha_{2}^{2}c_{0}\delta-2\delta^{2}+7c_{0}\delta^{2}-\delta^{3}+2c_{0}\delta^{3}+\alpha_{1}^{2}(1+\delta)-\alpha_{3}^{2}(1+\delta)\bigr{]}$
$\displaystyle\hphantom{b_{-}=}{}-a\bigl{[}\alpha_{1}^{2}(1-\delta+c_{0}(2\delta-1)-\alpha_{3}^{2}(2+\delta+c_{0}(2\delta-1))-\alpha_{2}^{2}(1-\delta+c_{0}(2\delta-1))$
$\displaystyle\hphantom{b_{-}=}{}+(1-\delta+c_{0}(2\delta-1))\big{[}3+\delta-\delta^{2}+c_{0}\big{(}{-}3+4\delta+4\delta^{2}\big{)}\big{]}\bigr{]}$
$\displaystyle\hphantom{b_{-}=}{}+(1+c_{0}(2\delta-1))\big{[}{-}\alpha_{3}^{2}+(-1+c_{0}+\delta-2c_{0}\delta)^{2}\big{]}\bigr{\\}},$
$\displaystyle
c_{-}=\frac{2\delta-1}{4\kappa_{0}\delta^{2}a(a-1)}\bigl{\\{}a^{4}\bigl{(}\alpha_{2}^{2}-(\delta+1)^{2}\bigr{)}^{2}-2a^{3}\bigl{(}\alpha_{2}^{2}-(\delta+1)^{2}\bigr{)}\bigl{(}{-}2-\alpha_{1}^{2}+\alpha_{3}^{2}+\alpha_{2}^{2}+\delta^{2}$
$\displaystyle\hphantom{c_{-}=}{}+2c_{0}\big{(}1-\delta-2\delta^{2}\big{)}\bigr{)}+2a\bigl{(}\alpha_{3}^{2}-(\delta-1+c_{0}(1-2\delta))^{2}\bigr{)}\bigl{(}2+\alpha_{1}^{2}-\alpha_{3}^{2}-\alpha_{2}^{2}-\delta^{2}$
$\displaystyle\hphantom{c_{-}=}{}+2c_{0}\big{(}{-}1+\delta+2\delta^{2}\big{)}\bigr{)}+\bigl{(}\alpha_{3}^{2}-(\delta-1+c_{0}(1-2\delta))^{2}\bigr{)}^{2}\bigr{\\}}.$
(A.3)

###  A.3 Parameterization of $\boldsymbol{A(\lambda)}$ at the pole of
$\boldsymbol{y(x)}$ as $\boldsymbol{\delta\neq 0,\frac{1}{2},1}$, and
$\boldsymbol{\sigma=-1}$

The constant parameters in (2.6) are determined by the parameters of the
Laurent expansion as follows

$\displaystyle\hat{a}_{3}=-\delta+1,\qquad\hat{b}_{3}=a-1+c_{0}(2\delta-1),$
$\displaystyle\hat{c}_{3}=\frac{1}{2(\delta-1)}\bigl{\\{}a^{2}\bigl{[}(\delta-2)^{2}-\alpha_{2}^{2}\bigr{]}+a\bigl{[}{-}\alpha_{1}^{2}+\alpha_{3}^{2}+\alpha_{2}^{2}+(\delta-1)^{2}-2-c_{0}\big{(}4(\delta-1)^{2}-2\delta\big{)}\bigr{]}$
$\displaystyle\hphantom{\hat{c}_{3}=}{}-\alpha_{3}^{2}+(\delta-
c_{0}(2\delta-1))^{2}\bigr{\\}},$
$\displaystyle\hat{b}_{+}=\frac{\kappa_{0}}{a(a-1)(\delta-1)(2\delta-1)}\bigl{\\{}a^{3}\big{(}{-}\alpha_{2}^{2}+(\delta-2)^{2}\big{)}+a^{2}\bigl{[}\alpha_{2}^{2}-2+c_{0}\big{(}\alpha_{2}^{2}-8\big{)}$
$\displaystyle\hphantom{\hat{b}_{+}=}{}+(\delta-2)\big{(}\alpha_{1}^{2}-\alpha_{3}^{2}\big{)}+\delta\big{(}\alpha_{2}^{2}-5\big{)}+2c_{0}\delta\big{(}11-\alpha_{2}^{2}\big{)}+c_{0}\delta^{2}(2\delta-13)-\delta^{2}(\delta-5)\bigr{]}$
$\displaystyle\hphantom{\hat{b}_{+}=}{}-(c_{0}(2\delta-1)-1)\big{(}\alpha_{3}^{2}-(c_{0}+\delta-2c_{0}\delta)^{2}\big{)}-a\bigl{[}c_{0}^{2}(1-2\delta)^{2}(2\delta-5)$
$\displaystyle\hphantom{\hat{b}_{+}=}{}-\delta\big{(}3+\alpha_{1}^{2}-\alpha_{2}^{2}+\delta-\delta^{2}\big{)}+c_{0}(2\delta-1)\big{(}3+\alpha_{1}^{2}-\alpha_{2}^{2}+6\delta-3\delta^{2}\big{)}$
$\displaystyle\hphantom{\hat{b}_{+}=}{}+\alpha_{3}^{2}(3-\delta-
c_{0}(2\delta-1))\bigr{]}\bigr{\\}},$
$\displaystyle\hat{c}_{+}=\frac{\kappa_{0}}{4a(a-1)(\delta-1)^{2}(2\delta-1)}\bigl{\\{}a^{4}\big{(}\alpha_{2}^{2}-(\delta-2)^{2}\big{)}^{2}$
$\displaystyle\hphantom{\hat{c}_{+}=}{}-2a^{3}\big{(}\alpha_{2}^{2}-(\delta-2)^{2}\big{)}\bigl{[}{-}1-2\delta+\delta^{2}-\alpha_{1}^{2}+\alpha_{3}^{2}+\alpha_{2}^{2}+c_{0}\big{(}{-}4+10\delta-4\delta^{2}\big{)}\bigr{]}$
$\displaystyle\hphantom{\hat{c}_{+}=}{}+a^{2}\bigl{[}\alpha_{1}^{4}+\alpha_{3}^{4}+\big{(}\alpha_{2}^{2}-1\big{)}^{2}+8c_{0}\big{(}1-\alpha_{2}^{2}\big{)}+c_{0}^{2}\big{(}24-2\alpha_{2}^{2}\big{)}+4\delta\big{(}1-\alpha_{2}^{2}\big{)}$
$\displaystyle\hphantom{\hat{c}_{+}=}{}+c_{0}\delta\big{(}12+16\alpha_{2}^{2}\big{)}+c_{0}^{2}\delta\big{(}{-}120+8\alpha_{2}^{2}\big{)}+\delta^{2}\bigl{(}10-88c_{0}+c_{0}^{2}\big{(}198-8\alpha_{2}^{2}\big{)}\bigr{)}$
$\displaystyle\hphantom{\hat{c}_{+}=}{}+\delta^{3}\bigl{(}{-}12+72c_{0}-120c_{0}^{2}\bigr{)}+\delta^{4}\big{(}3-16c_{0}+24c_{0}^{2}\big{)}-2\alpha_{1}^{2}\big{(}1+\alpha_{3}^{2}+\alpha_{2}^{2}-4c_{0}$
$\displaystyle\hphantom{\hat{c}_{+}=}{}-6\delta+10c_{0}\delta+3\delta^{2}-4c_{0}\delta^{2}\big{)}+2\alpha_{3}^{2}\big{(}{-}5+2\alpha_{2}^{2}+2\delta-2c_{0}\big{(}2-5\delta+2\delta^{2}\big{)}\big{)}\bigr{]}$
$\displaystyle\hphantom{\hat{c}_{+}=}{}+2a\big{(}\alpha_{3}^{2}-(c_{0}+\delta-2c_{0}\delta)^{2}\big{)}\bigl{[}1+2\delta-\delta^{2}+\alpha_{1}^{2}-\alpha_{3}^{2}-\alpha_{2}^{2}+c_{0}\big{(}4-10\delta+4\delta^{2}\big{)}\bigr{]}$
$\displaystyle\hphantom{\hat{c}_{+}=}{}+\big{(}\alpha_{3}^{2}-(c_{0}+\delta-2c_{0}\delta)^{2}\big{)}^{2}\bigr{\\}},$
$\displaystyle\hat{c}_{-}=-\frac{a(a-1)(2\delta-1)}{\kappa_{0}}.$

###  A.4 Parameterization of $\boldsymbol{A(\lambda)}$ at the pole of
$\boldsymbol{y(x)}$ as $\boldsymbol{\delta=1}$ and $\boldsymbol{\sigma=-1}$

The parameters in the coefficient matrix (2.7) are given by

$\displaystyle\check{a}_{3}=-\frac{1}{2}=-\delta+\frac{1}{2},\qquad\check{b}_{3}=c_{0}+\frac{3}{2}a-\frac{1}{2},\qquad\check{c}_{3}=a\left(\alpha_{1}-\frac{1}{2}\right),$
$\displaystyle\check{b}_{+}=\frac{\kappa_{0}}{a(a-1)}\bigl{(}a^{2}\big{(}1-\alpha_{2}^{2}\big{)}-a\big{(}2+\alpha_{1}^{2}-\alpha_{3}^{2}-\alpha_{2}^{2}-2c_{0}\big{)}-\alpha_{3}^{2}+(c_{0}-1)^{2}\bigr{)},$
$\displaystyle\check{c}_{+}=\kappa_{0}\left({-}\alpha_{3}^{2}+\alpha_{2}^{2}+\frac{a+1}{a-1}\alpha_{1}^{2}+2\alpha_{1}\left(\frac{c_{0}}{a-1}+1\right)\right),\qquad\check{b}_{-}=-\frac{a(a-1)}{\kappa_{0}}.$

###  A.5 Parameterization of $\boldsymbol{A(\lambda)}$ at the pole of
$\boldsymbol{y(x)}$ as $\boldsymbol{\delta=\frac{1}{2}}$

The coefficients of $A(\lambda)$ (2.8) in terms of the monodromy exponents,
the pole position $a\neq 0,1$ and the arbitrary coefficient $c_{-2}\neq 0$ of
the Laurent expansion (2.4) are as follows

$\displaystyle\tilde{a}_{3}=-\frac{3}{2}=-\delta-1,\qquad\tilde{b}_{3}=a+1-\frac{a^{2}(a-1)^{2}}{2c_{-2}},$
$\displaystyle\tilde{c}_{3}=\frac{a^{4}(a-1)^{4}}{12c_{-2}^{2}}+\frac{a^{2}(a-1)^{2}(a+1)}{6c_{-2}}$
$\displaystyle\hphantom{\tilde{c}_{3}=}{}+\frac{1}{12}\bigr{(}a^{2}\big{(}1-4\alpha_{2}^{2}\big{)}+a\big{(}{-}7-4\alpha_{1}^{2}+4\alpha_{3}^{2}+4\alpha_{2}^{2}\big{)}+1-4\alpha_{3}^{2}\bigr{)},$
$\displaystyle\tilde{c}_{+}=-\frac{\kappa_{0}c_{-2}}{a^{2}(a-1)^{2}},$
$\displaystyle\tilde{b}_{-}=\frac{a^{8}(a-1)^{8}}{12\kappa_{0}c_{-2}^{4}}-\frac{(a-1)^{4}a^{4}}{12\kappa_{0}c_{-2}^{2}}\bigl{(}3\big{(}a^{2}-a+1\big{)}+4a\alpha_{1}^{2}-4(a-1)\alpha_{3}^{2}+4a(a-1)\alpha_{2}^{2}\bigr{)}$
$\displaystyle\hphantom{\tilde{b}_{-}=}{}+\frac{(a-1)^{2}a^{2}}{12\kappa_{0}c_{-2}}\bigl{(}-2a^{3}+3a^{2}+3a-2-4a(a+1)\alpha_{1}^{2}+4\big{(}a^{2}-3a+2\big{)}\alpha_{3}^{2}$
$\displaystyle\hphantom{\tilde{b}_{-}=}{}+4a\big{(}2a^{2}-3a+1\big{)}\alpha_{2}^{2}\bigr{)},$
$\displaystyle\tilde{c}_{-}=\frac{a^{10}(a-1)^{10}}{144\kappa_{0}c_{-2}^{5}}+\frac{a^{8}(a-1)^{8}(a+1)}{36\kappa_{0}c_{-2}^{4}}$
$\displaystyle\hphantom{\tilde{c}_{-}=}{}-\frac{a^{6}(a-1)^{6}}{72\kappa_{0}c_{-2}^{3}}\bigl{(}{-}3\big{(}a^{2}-a+1\big{)}+4a\alpha_{1}^{2}-4(a-1)\alpha_{3}^{2}+4a(a-1)\alpha_{2}^{2}\bigr{)}$
$\displaystyle\hphantom{\tilde{c}_{-}=}{}-\frac{a^{4}(a-1)^{4}(a+1)}{36\kappa_{0}c_{-2}^{2}}\bigl{(}{-}a^{2}+7a-1+4a\alpha_{1}^{2}-4(a-1)\alpha_{3}^{2}+4a(a-1)\alpha_{2}^{2}\bigr{)}$
$\displaystyle\hphantom{\tilde{c}_{-}=}{}+\frac{a^{2}(a-1)^{2}}{144\kappa_{0}c_{-2}}\bigl{(}\big{(}1-4\alpha_{3}^{2}\big{)}^{2}+a^{4}\big{(}1-4\alpha_{2}^{2}\big{)}^{2}+2a\big{(}{-}1+4\alpha_{3}^{2}\big{)}\big{(}7+4\alpha_{1}^{2}-4\alpha_{3}^{2}-4\alpha_{2}^{2}\big{)}$
$\displaystyle\hphantom{\tilde{c}_{-}=}{}+2a^{3}\big{(}7+4\alpha_{1}^{2}-4\alpha_{3}^{2}-4\alpha_{2}^{2}\big{)}\big{(}{-}1+4\alpha_{2}^{2}\big{)}+a^{2}\big{(}51+16\alpha_{1}^{4}+16\alpha_{3}^{4}-64\alpha_{2}^{2}+16\alpha_{2}^{4}$
$\displaystyle\hphantom{\tilde{c}_{-}=}{}+64\alpha_{3}^{2}\big{(}{-}1+\alpha_{2}^{2}\big{)}-8\alpha_{1}^{2}\big{(}11+4\alpha_{3}^{2}+4\alpha_{2}^{2}\big{)}\big{)}\bigr{)}.$
(A.4)

###  A.6 Asymptotic coefficients of $\boldsymbol{\Psi_{\infty}(\lambda)}$ as
$\boldsymbol{\lambda\to\infty}$

The asymptotic parameters in (3.10) are expressed as follows

$\displaystyle\delta\neq\pm\frac{1}{2},\pm
1\colon\quad\psi_{1}=\frac{\kappa}{2\delta-1}\sigma_{+}-\frac{\tilde{\kappa}}{2\delta+1}\sigma_{-},\qquad
d_{1}=-p-\delta(2y-x-1),$
$\displaystyle\psi_{2}=\kappa\frac{2p+(2\delta+1)y-x-1}{4(\delta-1)\big{(}\delta-\frac{1}{2}\big{)}}\sigma_{+}-\tilde{\kappa}\frac{2p-(2\delta+1)\tilde{y}+x+1+4\delta
y}{4(\delta+1)\big{(}\delta+\frac{1}{2}\big{)}}\sigma_{-},$ $\displaystyle
d_{21}=\frac{\kappa\tilde{\kappa}}{2(1+2\delta)}+\frac{1}{2}\bigl{(}p(y-x-1)+\delta\big{(}(y-x-1)^{2}-x\big{)}-z\bigr{)},$
$\displaystyle
d_{22}=\frac{\kappa\tilde{\kappa}}{2(1-2\delta)}-\frac{1}{2}\bigl{(}p(y-x-1)+\delta\big{(}(y-x-1)^{2}-x\big{)}-z\bigr{)}.$
(A.5)

### Acknowledgments

A.K. was supported by the project SPbGU 11.38.215.2014. He also thanks the of
SISSA where this project was originated. Many thanks to the anonymous referees
for their suggestions towards improving the manuscript.

## References

  * [1] Arscott F.M., Darai A., Curvilinear coordinate systems in which the Helmholtz equation separates, IMA J. Appl. Math. 27 (1981), 33–70.
  * [2] Boalch P., From Klein to Painlevé via Fourier, Laplace and Jimbo, Proc. London Math. Soc. 90 (2005), 167–208, math.AG/0308221.
  * [3] Dubrovin B., Mazzocco M., Monodromy of certain Painlevé-VI transcendents and reflection groups, Invent. Math. 141 (2000), 55–147, math.AG/9806056.
  * [4] Erdélyi A., Magnus W., Oberhettinger F., Tricomi F.G., Higher transcendental functions, McGraw-Hill, New York, 1953.
  * [5] Fokas A.S., Its A.R., Kapaev A.A., Novokshenov V.Yu., Painlevé transcendents. The Riemann–Hilbert approach, Mathematical Surveys and Monographs, Vol. 128, Amer. Math. Soc., Providence, RI, 2006.
  * [6] Fokas A.S., Its A.R., Kitaev A.V., The isomonodromy approach to matrix models in $2$D quantum gravity, Comm. Math. Phys. 147 (1992), 395–430.
  * [7] Fuchs R., Sur quelques équations différentielles linéaires du second ordre, C. R. Acad. Sci. Paris 141 (1906), 555–558.
  * [8] Fuchs R., Über lineare homogene Differentialgleichungen zweiter Ordnung mit drei im Endlichen gelegenen wesentlich singulären Stellen, Math. Ann. 63 (1907), 301–321.
  * [9] Garnier R., Sur des équations différentielles du troisième ordre dont l’intégrale générale est uniforme et sur une classe d’équations nouvelles d’ordre supérieur dont l’intégrale générale a ses points critiques fixes, Ann. Sci. École Norm. Sup. 29 (1912), 1–126.
  * [10] Garnier R., Étude de l’intégrale générale de l’équation VI de M. Painlevé dans le voisinage de ses singularités transcendantes, Ann. Sci. École Norm. Sup. 34 (1917), 239–353.
  * [11] Gromak V.I., Laine I., Shimomura S., Painlevé differential equations in the complex plane, De Gruyter Studies in Mathematics, Vol. 28, Walter de Gruyter & Co., Berlin, 2002.
  * [12] Guzzetti D., Pole distribution of PVI transcendents close to a critical point, Phys. D 241 (2012), 2188–2203, arXiv:1104.5066.
  * [13] Guzzetti D., Tabulation of Painlevé 6 transcendents, Nonlinearity 25 (2012), 3235–3276, arXiv:1108.3401.
  * [14] Heun K., Zur Theorie der Riemann’schen Functionen zweiter Ordnung mit vier Verzweigungspunkten, Math. Ann. 33 (1888), 161–179.
  * [15] Hortaçsu M., Heun functions and their uses in physics, in Proceedings of the 13th Regional Conference on Mathematical Physics (October 27–31, 2010, Antalya, Turkey), Editors U. Camci, I. Semi, World Sci. Publ., Hackensack, NJ, 2013, 23–39, arXiv:1101.0471.
  * [16] Iwasaki K., An area-preserving action of the modular group on cubic surfaces and the Painlevé VI equation, Comm. Math. Phys. 242 (2003), 185–219.
  * [17] Jimbo M., Monodromy problem and the boundary condition for some Painlevé equations, Publ. Res. Inst. Math. Sci. 18 (1982), 1137–1161.
  * [18] Jimbo M., Miwa T., Ueno K., Monodromy preserving deformation of linear ordinary differential equations with rational coefficients. I. General theory and $\tau$-function, Phys. D 2 (1981), 306–352.
  * [19] Jimbo M., Miwa T., Monodromy preserving deformation of linear ordinary differential equations with rational coefficients. II, Phys. D 2 (1981), 407–448.
  * [20] Litvinov A., Lukyanov S., Nekrasov N., Zamolodchikov A., Classical conformal blocks and Painlevé VI, J. High Energy Phys. 2014 (2014), no. 7, 144, 20 pages, 1309.4700.
  * [21] Maier R.S., The 192 solutions of the Heun equation, Math. Comp. 76 (2007), 811–843, math.CA/0408317.
  * [22] Mazzocco M., Rational solutions of the Painlevé VI equation, J. Phys. A: Math. Gen. 34 (2001), 2281–2294, nlin.SI/0007036.
  * [23] Ronveaux A. (Editor), Heun’s differential equations, Oxford University Press, Oxford – New York – Tokyo, 1995.
  * [24] Schlesinger L., Über eine Klasse von Differentialsystemen beliebiger Ordnung mit festen kritischen Punkten, J. Reine Angew. Math. 141 (1912), 96–145.
  * [25] Schmidt D., Wolf G., A method of generating integral relations by the simultaneous separability of generalized Schrödinger equations, SIAM J. Math. Anal. 10 (1979), 823–838.
  * [26] Slavyanov S.Yu., Painlevé equations as classical analogues of Heun equations, J. Phys. A: Math. Gen. 29 (1996), 7329–7335.
  * [27] Slavyanov S.Yu., Lay W., Special functions: a unified theory based on singularities, Oxford Mathematical Monographs, Oxford University Press, Oxford, 2000.
  * [28] Smirnov A.O., Finite-gap solutions of the Fuchsian equations, Lett. Math. Phys. 76 (2006), 297–316, math.CA/0310465.
  * [29] Stieltjes T.J., Sur certains polynômes: qui vérifient une équation différentielle linéaire du second ordre et sur la theorie des fonctions de Lamé, Acta Math. 6 (1885), 321–326.
  * [30] Takemura K., Integral transformation of Heun’s equation and some applications, J. Math. Soc. Japan 69 (2017), 849–891, arXiv:1008.4007.