PG course on \SPECIAL FUNCTIONS AND THEIR SYMMETRIES

PG course on
SPECIAL FUNCTIONS AND THEIR SYMMETRIES

Vadim KUZNETSOV

Course Outline

1  Gamma and Beta functions
    1.1  Introduction
    1.2  Gamma function
    1.3  Beta function
    1.4  Other beta integrals
        1.4.1  Second beta integral
        1.4.2  Third (Cauchy's) beta integral
        1.4.3  A complex contour for the beta integral
        1.4.4  The Euler reflection formula
        1.4.5  Double-contour integral
2  Hypergeometric functions
    2.1  Introduction
    2.2  Definition
    2.3  Euler's integral representation
    2.4  Two functional relations
    2.5  Contour integral representations
    2.6  The hypergeometric differential equation
    2.7  The Riemann-Papperitz equation
    2.8  Barnes' contour integral for F(a,b;c;x)
3  Orthogonal polynomials
    3.1  Introduction
    3.2  General orthogonal polynomials
    3.3  Zeros of orthogonal polynomials
    3.4  Gauss quadrature
    3.5  Classical orthogonal polynomials
    3.6  Hermite polynomials
4  Separation of variables and special functions
    4.1  Introduction
    4.2  SoV for the heat equation
    4.3  SoV for a quantum problem
    4.4  SoV and integrability
    4.5  Another SoV for the quantum problem
5  Integrable systems and special functions
    5.1  Introduction
    5.2  Calogero-Sutherland system
    5.3  Integral transform
    5.4  Separated equation
    5.5  Integral representation for Jack polynomials
Index

1  Gamma and Beta functions

1.1  Introduction

This course is about special functions and their properties. Many known functions could be called special. They certainly include elementary functions like exponential and more generally, trigonometric, hyperbolic functions and their inverses, logarithmic functions and poly-logarithms, but the class also expands into transcendental functions like Lamé and Mathieu functions. Usually one deals first with a special function of one variable before going into a study of its multi-variable generalisation, which is not unique and which opens up a link with the theory of integrable systems.

We will restrict ourselves to hypergeometric functions which are usually defined by series representations.

Definition 1 A hypergeometric series is a series ∑n=0 an with an+1/an a rational function of n.

Bessel, Legendre, Jacobi functions, parabolic cylinder functions, 3j- and 6j-symbols arising in quantum mechanics and many more classical special functions are partial cases of the hypergeometric functions. However, the trancendental functions ``lying in the land beyond Bessel'' are out of the hypergeormetric class and, thereby, will not be considered in this course.

Euler, Pfaff and Gauss first introduced and studied hypergeometric series, paying special attention to the cases when a series can be summed into an elementary function. This gives one of the motivations for studying hypergeometric series, i.e. the fact that the elementary functions and several other important functions in mathematics can be expressed in terms of hypergeometric functions.

Hypergeometric functions can also be described as solitions of special differential equations, the hypergeometric differential equations. Riemann was first who exploited this idea and introduced a special symbol to classify hypergeometic functions by singularities and exponents of differential equations they satisfy. In this way we come up with an alternative definition of a hypergeometric function.

Definition 2 A hypergeometric function is a solution of a Fuchsian differential equation which has at most three regular singularities.

Notice that transcendental special functions of the Heun class, so-called Heun functions which are ``beyond Bessel'', are defined as special solutions of a generic linear second order Fuchsian differential equation with four regular singularities.

Of course, when talking about 3 or 4 regular singularities it means that the number of singularities can be less than that, either by trivilising the equation or by merging any two regular singularities in order to obtain an irregular singular point thus leading to corresponding confluent cases of hypergeomeric or Heun functions. So, Bessel and parabolic cylinder functions are special cases of the confluent and double-confluent hypergeometric function, while Lamé and Mathieu functions are special cases of Heun and confluent Heun functions, respectively. When a maximal allowed number of singularities of a differential equation grows it results in a more trancendental special function associated with it. A short introduction to the theory of Fuchsian equations with n regular singularities will be given later on in the course.

In the first decade of the XXth century E.W. Barnes introduced yet another approach to hypergeometric functions based on contour integral representations. Such representations are important because they can be used for derivation of many relations between hypergeometric functions and also for studying their asymptotics.

The whole class of hypergeometric functions is very distinguished comparing to other special functions, because only for this class one can have explicit series and integral representations, contiguous and connection relations, summation and transformation formulas, and many other beautiful equations relating one hypergeometric function with another. This is a class of functions for which one can probably say that any meaningful formula can be written explicitly; it does not say though that it is always easy to find the one. Also for that reason this is the class of functions to start from and to put as a basis of an introductory course in special functions.

The main reason of many applications of hypergeometric functions and special functions in general is their usefulness. Summation formulas find their way in combinatorics; classical orthogonal polynomials give explicit bases in several important Hilbert spaces and lead to constructive Hamonic Analysis with applications in quantum physics and chemistry; q-hypergeometric series are related to elliptic and theta-functions and therefore find their application in integration of systems of non-linear differential equations and in some areas of numeric analysis and discrete mathematics.

For this part of the course the main reference is the recent book by G.E. Andrews, R. Askey and R. Roy ``Special Functions'', Encyclopedia of Mathematics and its Applications 71, Cambridge University Press, 1999. The book by N.M. Temme ``Special functions: an introduction to the classical functions of mathematical physics'', John Wiley & Sons, Inc., 1996, is also recommended as well as the classical reference: E.T. Whittaker and G.N. Watson ``A course of modern analysis'', Cambridge University Press, 1927.

1.2  Gamma function

The Gamma function Γ(x) was discovered by Euler in the late 1720s in an attempt to find an analytical continuation of the factorial function. This function is a cornerstone of the theory of special functions.

Thus Γ(x) is a meromorphic function equal to (x−1)! when x is a positive integer. Euler found its representation as an infinite integral and as a limit of a finite product. Let us derive the latter representation following Euler's generalization of the factorial.

ga01.png
Figure 1: The graph of absolute value of Gamma function of a complex variable z=x+iy from MuPAD computer algebra system. There are seen poles at z=0,−1,−2,−3,−4,….

Let x and n be nonnegative integers. For any a ∈ C define the shifted factorial (a)n by
(a)n=a(a+1)…(a+n−1)        for    n > 0,       (a)0=1.
(1)
Then, obviously,
x!=  (x+n)!

(x+1)n
=  n!(n+1)x

(x+1)n
=  n!nx

(x+1)n
 ·   (n+1)x

nx
.
(2)
Since
limn→∞  (n+1)x

nx
=1,
(3)
we conclude that
x!=limn→∞  n!nx

(x+1)n
 .
(4)
The limit exists for ∀x ∈ C such that x ≠ −1,−2,−3,…. for
 n!nx

(x+1)n
=
 n

n+1

x

 
n

j=1 

1+  x

j

−1

 

1+  1

j

x

 
(5)
and

1+  x

j

−1

 

1+  1

j

x

 
=1+  x(x−1)

2j2
+O
 1

j3

.
(6)

Definition 3 For ∀x ∈ C, x ≠ 0,−1,−2,…, the gamma function Γ(x) is defined by
Γ(x)=limk→∞  k!kx−1

(x)k
 .
(7)

Three immediate consequences are
Γ(1)=1,       Γ(x+1)=xΓ(x)       and    Γ(n+1)=n!.
(8)

>From the definition it follows that the gamma function has poles at zero and the negative integers, but 1/Γ(x) is an entire function with zeros at these points. Every entire function has a product representation.

Theorem 4
 1

Γ(x)
=xeγx

n=1 



1+  x

n

e−x/n

,
(9)
where γ is Euler's constant given by
γ = limn→∞
n

k=1 
 1

k
−log n
.
(10)

PROOF.
 1

Γ(x)
=
limn→∞  x(x+1)…(x+n−1)

n!nx−1
=
limn→∞  x
1+  x

1


1+  x

2


1+  x

n

e−xlog  n
=
limn→∞  xex(1+[ 1/2]+…+[ 1/n]−log n) n

k=1 



1+  x

k

e−x/k

=
xeγx

n=1 



1+  x

n

e−x/n

.
The infinite product in (9) exists because

1+  x

n

e−x/n=
1+  x

n


1−  x

n
+  x2

2n2
 …
= 1−  x2

2n2
+O
 1

n3

.
(11)
[¯]

1.3  Beta function

Definition 5 The beta integral is defined for ℜ x > 0, ℜ y > 0 by
B(x,y)=
1

0 
tx−1(1−t)y−1dt.
(12)
One may also speak of the beta function B(x,y), which is obtained from the integral by analytic continuation.

The beta function can be expessed in terms of gamma functions.

Theorem 6
B(x,y)=  Γ(x)Γ(y)

Γ(x+y)
 .
(13)

PROOF. From the definition of beta integral we have the following contiguous relation between three functions
B(x,y+1)=B(x,y)−B(x+1,y),       ℜ x > 0,   ℜ y > 0.
(14)
However, integration by parts of the integral in the left hand side gives
B(x,y+1)=  y

x
 B(x+1,y).
(15)
Combining the last two we get the functional equation of the form
B(x,y)=  x+y

y
 B(x,y+1).
(16)
Iterating this equation we obtain
B(x,y)=  (x+y)n

(y)n
 B(x,y+n).
(17)
Rewrite this relation as
B(x,y)=  (x+y)n

n!nx+y−1
   n!ny−1

(y)n

n

0 
tx−1
1−  t

n

n+y−1

 
dt.
(18)
As n→∞, we have
B(x,y)=  Γ(y)

Γ(x+y)



0 
tx−1e−tdt.
(19)
Set y=1 to arrive at
 1

x
=
1

0 
tx−1dt=B(x,1) =  Γ(1)

Γ(x+1)



0 
tx−1e−tdt.
(20)
Hence
Γ(x)=


0 
tx−1e−tdt,        ℜ x > 0.
(21)
This is the integral representation for the gamma function, which appears here as a byproduct. Now use it to prove the theorem for ℜ x > 0 and ℜ y > 0 and then use the standard argument of analytic continuation to finish the proof. [¯]

An important corollary is an integral representation for the gamma function which may be taken as its definition for ℜ x > 0.

Corollary 7 For ℜ x > 0
Γ(x)=


0 
tx−1e−tdt.
(22)

Use it to explicitly represent the poles and the analytic continuation of Γ(x):
Γ(x)=
1

0 
tx−1e−tdt+


1 
tx−1e−tdt

=

n=0 
   (−1)n

(n+x)n!
+


1 
tx−1e−tdt.
The integral is an intire function and the sum gives the poles at x=−n, n=0,1,2, … with the residues equal to (−1)n/n!.

Several other useful forms of the beta integral can be derived by a change of variables. For example, take t=sin2θ in (12) to get

π/2

0 
sin2x−1θ cos2y−1θ dθ =  Γ(x)Γ(y)

2Γ(x+y)
 .
Put x=y=1/2. The result is
Γ(1/2)=

 

π
 
.

The substitution t=(u−a)/(b−a) gives

b

a 
(b−u)x−1(u−a)y−1du=(b−a)x+y−1   Γ(x)Γ(y)

Γ(x+y)
 ,
which can be rewritten in the alternative form:

b

a 
   (b−u)x−1

Γ(x)
   (u−a)y−1

Γ(y)
du=  (b−a)x+y−1

Γ(x+y)
 .

1.4  Other beta integrals

There are several kinds of integral representations for the beta function. All of them can be brought to the following form



C 
[l1(t)]p[l2(t)]qdt,
where l1(t) and l2(t) are linear functions of t, and C is an appropriate curve. The representation (12) is called Euler's first beta integral. Fot it, the curve consisits of a line segment connecting the two zeros of l-functions. We introduce now four more beta integrals. For the second beta integral, the curve is a half line joining one zero with infinity while the other zero is not on this line. For the third (Cauchy's) beta integral, it is a line with zeros on opposite sides. For the last two beta integrals, the curve is a complex contour. In the first case, it starts and ends at one zero and encircles the other zero in positive direction. In the second case, the curve is a double loop winding around two zeros, once in a positive direction and the second time in the negative direction.

1.4.1  Second beta integral

Set t=s/(s+1) in (12) to obtain the second beta integral with integration over half line,



0 
   sx−1

(1+s)x+y
 ds=  Γ(x)Γ(y)

Γ(x+y)
 .
(23)

1.4.2  Third (Cauchy's) beta integral

The beta integral due to Cauchy is defined by
C(x,y)=


−∞ 
   dt

(1+it)x(1−it)y
=  π22−x−yΓ(x+y−1)

Γ(x)Γ(y)
 ,       ℜ(x+y) > 1.

PROOF. To prove this, first show that integration by parts gives
C(x,y+1)=−  x

y
 C(x+1,y).
Also,
C(x,y)=


−∞ 
   (−1−it)+2

(1+it)x(1−it)y+1
=2C(x,y)−C(x−1,y+1).
Last two combine to give the functional equation
C(x,y)=  2y

x+y−1
 C(x,y+1) .
Iteration gives
C(x,y)=  22n(x)n(y)n

(x+y−1)2n
 C(x+n,y+n) .
Now,
C(x+n,y+n)=


−∞ 
 dt

(1+t2)n(1+it)x(1−it)y
 .
Set t→ t/√n and let n→ ∞. [¯]

The substitution t=tanθ leads to the integral:

π/2

0 
cosx+y−2θ cos(x−y)θ dθ =  π21−x−yΓ(x+y−1)

Γ(x)Γ(y)
 ,    ℜ(x+y) > 1.

1.4.3  A complex contour for the beta integral

Consider the integral
Ix,y=  1

2πi
 
(1+)

0 
wx−1(w−1)y−1dw,
with ℜ x > 0 and y ∈ C. The contour starts and ends at the origin, and encircles the point 1 in positive direction. The phase of w−1 is zero at positive points larger than 1. When ℜ y > 0 we can deform the contour along (0,1). Then we obtain Ix,y=B(x,y)sin(πy )/π. It follows that
B(x,y)=  1

2isinπy

(1+)

0 
wx−1(w−1)y−1dw.
The integral is defined for any complex value of y. For y=1,2,…, the integral vanishes; this is canceled by the infinite values of the term in front of the integral.

There is a similar contour integral representing the gamma function. Let us first prove Hankel's contour integral for the reciprocal gamma function, which is one of the most beautiful and useful representations for this function. It has the following form:
 1

Γ(z)
=  1

2πi



L 
s−zesds,      z ∈ C.
(24)
The contour of integration L is the Hankel contour that runs from −∞, arg s=−π, encircles the origin in positive direction and terminates at −∞, now with arg s=+π. For this we also use the notation ∫−∞(0+). The multi-valued function s−z is assumed to be real for real values of z and s, s > 0.

A proof of (24) follows immediately from the theory of Laplace transforms: from the well-known integral
 Γ(z)

sz
=


0 
tz−1 e−st dt
(24) follows as a special case of the inversion formula. A direct proof follows from a special choice of the contour L: the negative real axis. When ℜ z < 1 we can pull the contour onto the negative axis, where we have
 1

2πi


0

 
(se−iπ)−ze−sds−


0 
(se+iπ)−ze−sds
=  1

π
 sinπz  Γ(1−z).
Using the reflection formula (cf. next subsection for a proof),
Γ(x)Γ(1−x)=  π

sinπx
 ,
(25)
we see that this is indeed the left hand side of (24). In a final step the principle of analytic continuation is used to show that (24) holds for all finite complex values of z. Namely, both the left- and the right-hand side of (24) are entire functions of z.

Another form of (24) is
Γ(z)=  1

2isinπz
 


L 
sz−1 es ds.

1.4.4  The Euler reflection formula

The Euler reflection formula (25) connects the gamma function with the sine function. In a sense, it shows that 1/Γ(x) is `half of the sine function'. To prove the formula (25), set y=1−x, 0 < x < 1 in (23) to obtain
Γ(x)Γ(1−x)=


0 
 tx−1

1+t
 dt.
To compute the integral, consider the following contour integral



C 
 zx−1

1−z
 dz,
where C consists of two circles about the origin of radii R and ε respectively, which are joined along the negative real axis from −R to −ε. Move along the outer circle in the counterclockwise direction, and along the inner circle in the clockwise direction. By the residue theorem



C 
 zx−1

1−z
 dz=−2πi,
when zx−1 has its principal value. Thus
−2πi=
π

−π 
 iRxeixθ

1−Re
 dθ+
ε

R 
 tx−1eixπ

1+t
 dt+
−π

π 
 iεxeixθ

1−εe
 dθ+
R

ε 
 tx−1e−ixπ

1+t
 dt.
Let R→∞ and ε→0 so that the first and third integrals tend to zero and the second and fourth combine to give (25) for 0 < x < 1. The full result follows by analytic continuation.

1.4.5  Double-contour integral

We have seen that it is possible to replace the integral for Γ(z) along a half line by a contour integral which converges for all values of z. A similar process can be carried out for the beta integral.

Let P be any point between 0 and 1. We have the following Pochhammer's extension of the beta integral:

(1+,0+,1−,0−)

P 
tx−1(1−t)y−1dt =  −4π2eπi(x+y)

Γ(1−x)Γ(1−y)Γ(x+y)
 .
Here the contour starts at P, encircles the point 1 in the positive (counterclockwise) direction, returns to P, then encircles the origin in the positive direction, and returns to P. The 1−,0− indicates that now the path of integration is in the clockwise direction, first around 1 and then 0. The formula is proved by the same method as Hankel's formula. Notice that it is true for any complex x and y: both sides are entire functions of x and y.

2  Hypergeometric functions

2.1  Introduction

In this Lecture we give the definition and main properties of the Gauss (F=2F1) hypergeometric function and shortly mention its generalizations, the pFq generalized and pφq basic (or q-) hypergeometric functions.

Almost all of the elementary functions of mathematics and some not very elementary, like the error function erf(x) and dilogarithm function Li2(x), are special cases of the hypergeometric functions, or they can be expressed as ratios of hypergeometric functions.

We will first derive Euler's fractional integral representation for the Gauss hypergeometric function F, from which many identities and transformations will follow. Then we talk about hypergeometric differential equation, as a general linear second order differential equation having three regular singular points. We derive contiguous relations satisfied by the function F. Finally, we explain the Barnes approach to the hypergeometric functions and Barnes-Mellin contour integral representation for the function F.

2.2  Definition

Directly from the definition of a hypergeometric series ∑cn, on factorizing the polynomials in n, we obtain
 cn+1

cn
=  (n+a1)(n+a2)…(n+ap)x

(n+b1)(n+b2)…(n+bq)(n+1)
 .
Hence, we can get a more explicit definition.

Definition 1 The (generalized) hypergeometric series is defined by the following series representation
pFq



a1,…,ap
b1,…, bq
;x



=

n=0 
 (a1)n…(ap)n

(b1)n…(bq)n
   xn

n!
 .

Sometimes, we will use other notation: pFq(a1,…, ap;b1,…,bq;x).

If we apply the ratio test to determine the convergence of the series,

 cn+1

cn

 |x|np−q−1(1+|a1|/n)…(1+|ap|/n)

|(1+1/n)(1+b1/n)…(1+bq/n)|
 ,
then we get the following theorem.

Theorem 2 The series pFq(a1,…, ap; b1,…,bq;x) converges absolutely for all x if p ≤ q and for |x| < 1 if p=q+1, and it diverges for all x ≠ 0 if p > q+1 and the series does not terminate.

PROOF. It is clear that |cn+1/cn|→0 as n→∞ if p ≤ q. For p=q+1, limn→∞ |cn+1/cn|=|x|, and for p > q+1, |cn+1/cn|→∞ as n→ ∞. This proves the theorem.

[¯]

The case |x|=1 when p=q+1 is of interest. Here we have the following conditions for convergence.

Theorem 3 The series q+1Fq(a1,…, aq+1; b1,…,bq;x) with |x|=1 converges absolutely if ℜ (∑bi−∑ai) > 0. The series converges conditionally if x=e ≠ 1 and 0 ≥ ℜ (∑bi−∑ai) > −1 and the series diverges if ℜ (∑bi−∑ai) ≤ −1.

PROOF. Notice that the shifted factorial can by expressed as a ratio of two gamma functions:
(x)n=  Γ(x+n)

Γ(x)
 .
By the definition of the gamma function

lim
n→∞ 
   Γ(n+x)

Γ(n+y)
 ny−x=  Γ(x)

Γ(y)
 
lim
n→∞ 
   (x)n

(y)n
 ny−x =  Γ(x)

Γ(y)
 ·   Γ(y)

Γ(x)
=1.
The coefficient of nth term in q+1Fq
 (a1)n…(aq+1)n

(b1)n…(bq)nn!
 ∼ 

Γ(bi)


Γ(ai)
 n∑a − ∑b −1
as n→∞. The statements about absolute convergence and divergence follow immediately. The part of the theorem concerning conditional convergence can by proved by summation by parts.

[¯]

The 2F1 series was studied extensively by Euler, Pfaff, Gauss, Kummer and Riemann. Examples:
log(1+x)=x 2F1



1,1
2
;−x



;

tan−1x=x 2F1



1/2,1
3/2
;−x2



;

sin−1x=x 2F1



1/2,1/2
3/2
;x2



;

(1−x)−a=1F0



a
;x



;

sinx=x 0F1



3/2
;−x2/4



;

cosx=0F1



1/2
;−x2/4



;

ex0F0



;x



.
The next set of examples uses limits:
ex=
lim
b→∞ 
2F1



1,b
1
;x/b



;

coshx=
lim
a,b→∞ 
2F1



a,b
1/2
;x2/(4ab)



;

1F1



a
c
;x



=
lim
b→∞ 
2F1



a,b
c
;x/b



;

0F1



c
;x



=
lim
a,b→∞ 
2F1



a,b
c
;x/(ab)



.
The example of log(1−x)=−x 2F1(1,1;2;x) shows that though the series converges for |x| < 1, it has a continuation as a single-valued function in the complex plane from which a line joining 1 to ∞ is deleted. This describes the general situation; a 2F1 function has a continuation to the complex plane with branch points at 1 and ∞.

Definition 4 The (Gauss) hypergeometric function 2F1(a,b;c;x) is defined by the series


n=0 
   (a)n(b)n

(c)nn!
 xn
for |x| < 1, and by continuation elsewhere.

2.3  Euler's integral representation

Theorem 5 If ℜc > ℜb > 0, then
2F1



a,b
c
;x



=  Γ(c)

Γ(b)Γ(c−b)
 
1

0 
tb−1(1−t)c−b−1(1−xt)−adt
(26)
in the x plane cut along the real axis from 1 to ∞. Here it is understood that argt=arg(1−t)=0 and (1−xt)−a has its principal value.

PROOF. Suppose |x| < 1. Expand (1−xt)−a by the binomial theorem
 Γ(c)

Γ(b)Γ(c−b)


n=1 
 (a)n

n!
 xn
1

0 
tn+b−1(1−t)c−b−1dt.
Since for ℜb > 1, ℜ(c−b) > 1 and |x| < 1 the series


n=0 
Un(t),       Un(t)=xn   (a)n

n!
 tb+n−1(1−t)c−b−1
converges uniformly with respect to t ∈ [0,1], we are able to interchange the order of integration and summation for these values of b, c and x.

Now, use the beta integral to prove the result for |x| < 1. Since the integral is analytic in the cut plane, the theorem holds for x in this region as well; also we apply the analytic continuation with respect to b and c in order to arrive at the conditions announced in the formulation of the theorem.

[¯]

Hence we have obtained the analytic continuation of F, as a function of x, outside the unit disc, but only when ℜc > ℜb > 0. It is important to note that we view 2F1(a,b;c;x) as a function of four complex variables a, b, c, and x instead of just x. It is easy to see that [ 1/Γ(c)] 2F1(a,b;c;x) is an entire function of a, b, c if x is fixed and |x| < 1, for in this case the series converges uniformly in every compact domain of the a, b, c space.

Gauss found evaluation of the series in the point 1.

Theorem 6 For ℜ(c−a−b) > 0


n=0 
 (a)n(b)n

(c)nn!
=2F1



a,b
c
;1



=  Γ(c)Γ(c−a−b)

Γ(c−a)Γ(c−b)
 .

PROOF. Let x→ 1 in Euler's integral for 2F1. Then when ℜc > ℜb > 0 and ℜ(c−a−b) > 0 we get
 Γ(c)

Γ(b)Γ(c−b)
 
1

0 
tb−1 (1−t)c−a−b−1dt =  Γ(c)Γ(c−a−b)

Γ(c−a)Γ(c−b)
 .
The condition ℜc > ℜb > 0 may be removed by continuation.

[¯]

Corollary 7 (Chu-Vandermonde)
2F1



−n,b
c
;1



=  (c−b)n

(c)n
 ,    n=0,1,2,….

2.4  Two functional relations

The hypergeometric function satisfies a great number of relations. The most simple and obvious is the symmetry a↔ b. Let us prove two more relations
2F1



a,b
c
;x



=(1−x)−a2F1



a,c−b
c
;  x

x−1




       (Pfaff),
(27)

2F1



a,b
c
;x



=(1−x)c−a−b2F1



c−a,c−b
c
;x



       (Euler).
(28)
First relation is proved through the change of variable t=1−s in Euler's integral formula. The second relation follows by using the first relation twice.

The right-hand series in Pfaff's transformation converges for |x/(x−1)| < 1. This condition is implied by ℜx < 1/2; so we have a continuation of the series 2F1(a,b;c;x) to this region by Pfaff's formula.

Now, let us rewrite Euler's transformation as
(1−x)a+b−c2F1



a,b
c
;x



=2F1



c−a,c−b
c
;x



.
Equate the coefficient of xn on both sides to get
n

j=0 
 (a)j(b)j(c−a−b)n−j

j!(c)j(n−j)!
=  (c−a)n(c−b)n

n!(c)n
 .
Rewrite this as:

Theorem 8 (Pfaff-Saalschütz)
3F2



−n,a,b
c,1+a+b−c−n
;1



=  (c−a)n(c−b)n

(c)n(c−a−b)n
 .

The Pfaff-Saalschütz identity can be written as
(c)n(c+a+b)n  3F2



−n,−a,−b
c,1−a−b−c−n
;1



= (c+a)n(c+b)n.
This is a polynomial identity in a, b, c. Dougall (1907) took the view that both sides of this equation are polynomials of degree n in a. Therefore, the identity is true if both sides are equal for n+1 distinct values of a. By the same method he proved a more general identity:
7F6







a,1+  1

2
 a,−b,−c,−d,−e,−n
 1

2
 a, 1+a+b, 1+a+c, 1+a+d, 1+d+e, 1+a+n
;1








=  (1+a)n(1+a+b+c)n(1+a+b+d)n(1+a+c+d)n

(1+a+b)n(1+a+c)n(1+a+d)n(1+a+b+c+d)n
 ,
where 1+2a+b+c+d+e+n=0 and n is a positive integer. Taking the limit n→∞ we get
5F4







a,1+  1

2
 a,−b,−c,−d
 1

2
 a, 1+a+b, 1+a+c, 1+a+d
;1








=  Γ(1+a+b)Γ(1+a+c)Γ(1+a+d)Γ(1+a+b+c+d)

Γ(1+a)Γ(1+a+b+c)Γ(1+a+b+d)Γ(1+a+c+d)
when ℜ(a+b+c+d+1) > 0. Now, take d=−a/2 to get Dixon's summation formula
3F2



a,−b,−c
1+a+b, 1+a+c
;1



=
Γ(1+  1

2
a)Γ(1+a+b)Γ(1+a+c)Γ(1+  1

2
a+b+c)

Γ(1+a)Γ(1+  1

2
+b)Γ(1+  1

2
a+c)Γ(1+  1

2
a+b+c)
 .

2.5  Contour integral representations

A more general integral representation for the 2F1 hypergeometric function is the loop integral defined by
2F1



a,b
c
;x



=  Γ(c)Γ(1+b−c)

2πiΓ(b)
 
(1+)

0 
tb−1(t−1)c−b−1(1−xt)−adt,    ℜb > 0.

The contour starts and terminates at t=0 and encircles the point t=1 in the positive direction. The point 1/x should be outside the contour. The many-valued functions of the integrand assume their principal branches: arg(1−xt) tends to zero when x→0, and argt, arg(t−1) are zero at the point where the contour cuts the real positive axis (at the right of 1). Observe that no condition on c is needed, whereas in (26) we need ℜ(c−b) > 0. The proof of the above representation runs as for (26), with the help of the corresponding loop integral for the beta function.

Alternative representation involves a contour encircling the point 0:
2F1



a,b
c
;x



=  Γ(c)Γ(1−b)

2πiΓ(c−b)
 
(0+)

1 
(−t)b−1(1−t)c−b−1(1−xt)−adt,    ℜc > ℜb.

Using the double-loop (or Pochhammer's) contour integral one can derive the following representation
 1

Γ(c)
 2F1



a,b
c
;x



=−  e−i πc

4Γ(b)Γ(c−b)sinπb sinπ(c−b)
 

×
(1+,0+,1−,0−)

P 
tb−1(1−t)c−b−1(1−xt)−adt.
Here we have following conditions: |arg(1−x)| < π, argt=arg(1−t)=0 at the starting point P of the contour, and (1−xt)−a=1 when x=0. Note that there are no conditions on a, b, or c.

2.6  The hypergeometric differential equation

Let us introduce the differential operator ϑ = x   d/dx. We have
ϑ(ϑ+c−1) xn = n(n+c−1) xn.
Hence
ϑ(ϑ+c−1)   2F1(a,b;c;x) = x(ϑ+a)(ϑ+b)   2F1(a,b;c;x) .
In explicit form it reads
x(1−x)F′′+[c−(a+b+1)x)F′−abF=0,
(29)

F=F(a,b;c;x)=2F1(a,b;c;x).
This is the hypergeometric differential equation, which was given by Gauss.

It is easy to show that, in addition to F(a,b;c;x), a second solution of (29) is given by
x1−cF(a−c+1,b−c+1; 2−c;x).
When c=1 it does not give a new solution, but, in general, the second solution of (29) appears to be of the form
PF(a,b;c;x)+Qx1−cF(a−c+1,b−c+1; 2−c;x),
(30)
where P and Q are independent of x.

Next we observe that with the help of (29) and (30) we can express a hypergeometric function with argument 1−x or 1/x in terms of functions with argument x. For example, when in (29) we introduce a new variable x′=1−x we obtain a hypergeometric differential equation, but now with parameters a, b and a+b−c+1. Hence, besides the solutions in (30) we have F(a,b;a+b−c+1;1−x) as a solution as well. Any three solutions have to be linearly dependent. Therefore we get
F(a,b;a+b−c+1;1−x) = PF(a,b;c;x)+Qx1−cF(a−c+1,b−c+1; 2−c;x).
To find P and Q we substitute z=0 and z=1. If we also use Pfaff's and Euler's transformations we can get the following list of relations:
F(a,b;c;x)
=
AF(a,b;a+b−c+1;1−x)
+B(1−x)c−a−bF(c−a,c−b;c−a−b+1;1−x)
(31)
=
C(−x)−aF(a,1−c+a;1−b+a;1/x)
+D(−x)−bF(b,1−c+b;1−a+b;1/x)
(32)
=
C(1−x)−aF(a,c−b;a−b+1;1/(1−x))
+D(1−x)−bF(b,c−a;b−a+1;1/(1−x))
(33)
=
Ax−aF(a,a−c+1;a+b−c+1;1−1/x)
+Bxa−c(1−x)c−a−bF(c−a,1−a;c−a−b+1;1−1/x).
(34)
Here
A=  Γ(c)Γ(c−a−b)

Γ(c−a)Γ(c−b)
 ,      B=  Γ(c)Γ(a+b−c)

Γ(a)Γ(b)
 ,

C=  Γ(c)Γ(b−a)

Γ(b)Γ(c−a)
 ,      D=  Γ(c)Γ(a−b)

Γ(a)Γ(c−b)
 .
Since the Pfaff's formula (27) gives a continuatiuon of 2F1 from |x| < 1 to ℜx < [ 1/2], then (31) gives the continuation to ℜx > [ 1/2] cut along the real axis from x=1 to x=∞. The cut comes from the branch points of the factor (1−x)c−a−b. Analogously, (32) holds when |arg(−x)| < π; (33) holds when |arg(1−x)| < π; (34) holds when |arg(1−x)| < π and |argx| < π.

2.7  The Riemann-Papperitz equation

The hypergeometric differential equation (29) for the function 2F1 has three regular singular points, at 0, 1 and ∞ with exponents 0, 1−c; 0, c−a−b; and a, b respectively. Its Riemann symbol has the following form:
F=P





0
1
0
0
a
1−c
c−a−b
b
    x





.
In fact, this equation is a generic equation that has only three regular singularities.

Theorem 9 Any homogeneous linear differential equation of the second order with at most three singularities, which are regular singular points, can be transformed into the hypergeometric differential equation (29).

PROOF. Let us only sketch the proof. First we consider the equation
 d2f

dz2
+p(z)   df

dz
+q(z)f=0
and assume that it has only three finite regular singular points ξ, η and ζ with the exponents (α12), (β12) and (γ12). Then we find that such equation can be always brought into the form
f′′+
 1−α1−α2

z−ξ
+  1−β1−β2

z−η
+  1−γ1−γ2

z−ζ

f′
(35)


 α1α2

(z−ξ)(η−ζ)
+  β1β2

(z−η)(ζ−ξ)
+  γ1γ2

(z−ζ)(ξ−η)

    (ξ−η)(η−ζ)(ζ−ξ)

(z−ξ)(z−η)(z−ζ)
 f=0.
Next we introduce the following fractional linear transformation:
x=  (ζ−η)(z−ξ)

(ζ−ξ)(z−η)
 ,
and also a `gauge-transformation' of the function f:
F=x−α1(1−x)−γ1 f.
This transformation changes singularities to 0, 1 and ∞. The exponents in these points are
(0,α2−α1),       (0,γ2−γ1),      (α111, α121).
It is easy to check that we arrive at the hypergeometric differential equation (29) for the function F(x) with the following parameters:
a=α111,      b=α121,      c=1+α1−α2.

[¯]

Equation (35) is called the Riemann-Papperitz equation.

2.8  Barnes' contour integral for F(a,b;c;x)

The pair of Mellin transformations (direct and inverse) is defined by
F(s)=


0 
xs−1f(x)dx,      f(x)=  1

2πi

c+i∞

c−i∞ 
x−sF(s)ds.
It is true for some class of functions. For example,
Γ(s)=


0 
xs−1e−xdx,   e−x=  1

2πi

c+i∞

c−i∞ 
x−sΓ(s)ds,      c > 0.
This can be proved by Cauchy's residue theorem. Take a rectangular contour L with vertices c±iR, c−(N+[ 1/2])±iR, where N is a positive integer. The poles of Γ(s) inside this contour are at 0,1,…,N and the residues are (−1)j/j!. Now let R and N tend to ∞.

The Mellin transform of the hypergeometric function is



0 
xs−12F1



a,b
c
;−x



dx =  Γ(c)

Γ(a)Γ(b)
 Γ(s)Γ(a−s)Γ(b−s)

Γ(c−s)
 .

Theorem 10
 Γ(a)Γ(b)

Γ(c)
  2F1



a,b
c
;x



=  1

2πi

i∞

−i∞ 
 Γ(a+s)Γ(b+s)Γ(−s)

Γ(c+s)
 (−x)sds,
|arg(−x)| < π. The path of integration is curved, if necessary, to separate the poles s=−a−n, s=−b−n, from the poles s=n, where n is an integer ≥ 0. (Such a contour can always be drawn if a and b are not negative integers.)

3  Orthogonal polynomials

3.1  Introduction

In this lecture we talk about general properties of orthogonal polynomials and about classical orthogonal polynomials, which appear to be hypergeometric orthogonal polynomials. One way to link the hypergeometric function to orthogonal polynomials is through a formula of Jacobi. Multiply the hypergeometric equation by xc−1(1−x)a+b−c and write it in the following form
 d

dx
 [x(1−x)xc−1(1−x)a+b−cy′] = abxc−1(1−x)a+b−cy.
From
 d

dx
 2F1



a,b
c
;x



=  ab

c
 2F1



a+1,b+1
c+1
;x



,
by induction,
 d

dx
 [xk(1−x)kMy(k)] = (a+k−1)(b+k−1)xk−1(1−x)k−1My(k−1),
where M=xc−1(1−x)a+b−c. Then
 dk

dxk
 [xk(1−x)kMy(k)] = (a)k(b)k My.
Substitute
y(k)=  (a)k(b)k

(c)k
 2F1



a+k,b+k
c+k
;x



,
to get
 dk

dxk
 



xk(1−x)kM2F1



a+k,b+k
c+k
;x







=(c)k2F1



a,b
c
;x



.
Put b=−n, k=n, then
2F1



−n,a
c
;x



=  x1−c(1−x)c+n−a

(c)n
   dn

dxn
 [xc+n−1(1−x)a−c].
This is Jacobi's formula.

Set x=(1−y)/2, c=α+1, and a=n+α+β+1 to get
2F1



−n,n+α+β+1
α+1
;  1−y

2




=  (1−y)−α(1+y)−β

(α+1)n2n
   dn

dxn
 [(1−y)n+α(1+y)n+β].
(36)

Definition 1 The Jacobi polynomial of degree n is defined by
P(α,β)n(x):=  (α+1)n

n!
2F1



−n,n+α+β+1
α+1
;  1−x

2




.
(37)

Its orthogonality relation is as follows:

+1

−1 
Pn(α,β)(x)Pm(α,β)(x)(1−x)α(1+x)βdx

=  2α+β+1Γ(n+α+1)Γ(n+β+1)

(2n+α+β+1)Γ(n+α+β+1)n!
 δmn.
Formula (36) is called Rodrigues formula for Jacobi polynomials.

3.2  General orthogonal polynomials

Consider the linear space P of polynomials of the real variable x with real coefficients.

A set of orthogonal polynomials is defined by the interval (a,b) and by the measure dμ(x)=w(x)dx of orthogonality. The positive function w(x), with the property that

b

a 
w(x)xkdx < ∞,       ∀k=0,1,2,…,
is called the weight function.

Definition 2 We say that a sequence of polynomials {pn(x)}0, where pn(x) has exact degree n, is orthogonal with respect to the weight function w(x) if

b

a 
pn(x)pm(x)w(x)dx=hnδmn.

Theorem 3 A sequence of orthogonal polynomials {pn(x)} satisfies the three-term recurrence relation
pn+1(x)=(Anx+Bn)pn(x)−Cnpn−1(x)        for    n ≥ 0,
where we set p−1(x)=0. Here An, Bn, and Cn are real constants, n=0,1,2,…, and An−1AnCn > 0, n=1,2,…. If the highest coefficient of pn(x) is kn, then
An=  kn+1

kn
 ,        Cn+1=  An+1

An
   hn+1

hn
 .

An important consequence of the recurrence relation is the Christoffel-Darboux formula.

Theorem 4 Suppose that pn(x) are normalized so that hn=1. The
n

m=0 
pm(y)pm(x) =  kn

kn+1
   pn+1(x)pn(y)−pn+1(y)pn(x)

x−y
 .

Corollary 5 pn+1′(x)pn(x)−pn+1(x)pn′(x) > 0        ∀x.

3.3  Zeros of orthogonal polynomials

Theorem 6 Suppose that {pn(x)} is a sequence of orthogonal polynomials with respect to the weight function w(x) on the interval [a,b]. Then pn(x) has n simple zeros in [a,b].

Theorem 7 The zeros of pn(x) and pn+1(x) separate each other.

3.4   Gauss quadrature

Theorem 8 There are positive numbers λ1,…,λn such that for every polynomial f(x) of degree at most 2n−1

b

a 
f(x)w(x)dx= n

j=1 
λjf(xj),
where xj, j=1,…,n, are zeros of the polynomial pn(x) from the set of polynomials orthogonal with respect to the weight function w(x), and λj have the form
λj=
b

a 
 pn(x)w(x)dx

p′n(xj)(x−xj)
 .

3.5  Classical orthogonal polynomials

The orthogonal polynomials associated with the names of Jacobi, Gegenbauer, Chebyshev, Legendre, Laguerre and Hermite are called the classical orthogonal polynomials. The following properties are characteristic of the classical orthogonal polynomials:

(i) the family {pn′} is also an orthogonal system;

(ii) pn satisfies a second order linear differential equation
A(x)y′′+B(x)y′+λny=0,
where A and B do not depend on n and λn does not depend on x;

(iii) there is a Rodrigues formula of the form
pn(x)=  1

Knw(x)
 
 d

dx

n

 
[w(x)Xn(x)],
where X is a polynomial in x with coefficients not depending on n, and Kn does not depend on x.

As was said earlier all classical orthogonal polynomials are, in fact, hypergeometric polynomials, in the sense that they can be expressed in terms of the hypergeometric function. With the Jacobi polynomials expressed by (37), all the other appear to be either partial cases or limits from hypergeometric function to confluent hypergeometric function.

Gegenbauer polynomials:


Cnγ(x)=  (2γ)n

(γ+  1

2
)n
 Pn(γ−[ 1/2],γ−[ 1/2])(x),
(38)

Chebyshev polynomials:
Tn(x)=  n!

(  1

2
)n
 Pn(−[ 1/2],−[ 1/2])(x),
(39)

Legendre polynomials:
Pn(x)=Pn(0,0)(x),
(40)

Laguerre polynomials:
Lαn(x)=



n+α
n




  1F1(−n;α+1;x);
(41)

Hermite polynomials:
H2n(x)
=
 (−1)n(2n)!

n!
 1F1
−n;  1

2
;x2
,
(42)
H2n+1(x)
=
 (−1)n(2n+1)!

n!
 2x 1F1
−n;  3

2
;x2
.
(43)

3.6  Hermite polynomials

Hermite polynomials are orthogonal on (−∞,+∞) with the e−x2 as the weight function. This weight function is its own Fourier transform:
e−x2=  1




π
 


−∞ 
e−t2e2ixtdt.
(44)
Hermite polynomials can be defined by the Rodrigues formula:
Hn(x)=(−1)nex2   dne−x2

dxn
 .
It is easy to check that Hn(x) is a polynomial of degree n.

If we repeatedly differentiate (44) we get
 dne−x2

dxn
=  (2i)n




π
 


−∞ 
e−t2tne2ixtdt.
Hence
Hn(x)=  (−2i)nex2




π
 


−∞ 
e−t2tne2ixtdt.
It is easy now to prove the orthogonality property,



−∞ 
e−x2 Hn(x)Hm(x)dx = 2nn!

 

π
 
δmn,
using the Rodrigues formula and integrating by parts.

The Hermite polynomials have a simple generating function


n=0 
 Hn(x)

n!
 rn=e2xr−r2.
(45)
Recurrence relation has the form
H2n+1(x)−2xHn(x)+2nHn−1(x)=0.
From the integral representation we can derive the Poisson Kernel for the Hermite polynomials


n=0 
 Hn(x)Hn(y)

2nn!
 rn = (1−r2)−1/2e[2xyr−(x2+y2)r2]/(1−r2).
The following integral equation for |r| < 1 can be derived from the Poisson kernel by using orthogonality:
 1




π
 


−∞ 
 e[2xyr−(x2+y2)r2]/(1−r2)




1−r2
 Hn(y) dy=Hn(x)rn.
Let r→ i and we have, at least formally,
 1





 


−∞ 
eixye−y2/2Hn(y) dy=ine−x2Hn(x).
Hence, e−x2Hn(x) is an eigenfunction of the Fourier transform with eigenvalue in. This can be proved by using the Rodrigues formula for Hn(x).

4  Separation of variables and special functions

4.1  Introduction

This lecture is about some applications of special functions. It will also give some answer to the question: Where special functions come from? Special functions usually appear when solving linear partial differential equations (PDEs), like heat equation, or when solving spectral problems arising in quantum mechanics, like finding eigenfuntions of a Schödinger operator. Many equations of this kind, including many PDEs of mathematical physics can be solved by the method of Separation of Variables (SoV). We will give an introduction to this very powerful method and also will see how it fits into the theory of special functions.

Definition 1 Separation of Variables M is a transformation which brings a function ψ(x1,…,xn) of many variables into a factorized form
M: ψ→ ϕ1(y1)·…·ϕn(yn).

Functions ϕj(yj) are usually some known special functions of one variable. The transformation M could be a change of variables from {x} to {y}, but could also be an integral transform. Usually the function ψ satisfies a simple linear PDE.

4.2  SoV for the heat equation

Let the complex valued function q(x,t) satisfy the heat equation
iqt+qxx=0,       x,t ∈ [0,∞).
(46)

q(x,0)=q1(x),       q(0,t)=q2(t),
where q1(x) and q2(t) are given functions decaying sufficiently fast for large x and large t, respectively.

Divide (46) by q(x,t) and rewrite it as
iqt=k2q,       qxx=−k2q,
which is a separation of variables, since there is a factorized solution of last two equations:
qk(x,t)=e−ikx−ik2t.
Notice that this gives solution to (46) ∀k ∈ C. Because our equation is linear, the following integral is also a solution to the equation (46)
q(x,t)=


L 
e−ikx−ik2tρ(k)dk,
(47)
where L is some contour in the complex k-plane, and the function ρ(k) (`spectral data') can be expressed in terms of certain integral transforms of q1(x) and q2(t), in order to satisfy the initial data.

This is just a simple demonstration of the method of Separation of Variables, also called Ehrenpreis principle when applied to such kind of problems. It is interesting to note that all solutions of (46) can be given by (47) with the appropriate choice of the contour L and the function ρ(x).

4.3  SoV for a quantum problem

Now consider another simple problem that comes from quantum mechanics, namely: the linear spectral problem for the stationary Schödinger operator describing bound states of the 2-dimensional harmonic oscillator. That is, consider ψ(x1,x2) ∈ L2(R2) which is an eigenfunction of the following differential operator:
H=−
 ∂2

∂x12
+  ∂2

∂x22

+x12+x22,

Hψ(x1,x2)=hψ(x1,x2).
This problem can be solved by the straightforward application of the method of SoV without any intermediate transformations. We get
H1ψ =
 ∂2

∂x12
+x12
ψ(x1,x2)=h1ψ(x1,x2),

H2ψ =
 ∂2

∂x22
+x22
ψ(x1,x2)=h2ψ(x1,x2),

h1+h2=h.
Then
ψ(x1,x2)=ψ(x1)ψ(x2),
where

 ∂2

∂xi2
+xi2
ψ(xi)=hiψ(xi).
Square-integrable solution is expressed in terms of the Hermite polynomials
ψ(xi) ∈ L2(R)⇔ ψ(xi) = e−xi2/2Hn(xi),    h=2ni+1,    ni=0,1,2,….
Notice that
[H1,H2]=0.
Hence, we get the basis in L2(R2) of the form
ψn1n2(x1,x2)=e−(x12+x22)/2Hn1(x1)Hn2(x2),

h=2(n1+n2)+2.
The functions {ψn1n2} constitute an orthogonal set of functions in R2



−∞ 
ψn1n2(x1,x2m1m2(x1,x2) = an1n2δn1m1δn2m2.
Every function f ∈ L2(R2) can be decomposed into a series with respect to these basis functions
f(x1,x2)=

m,n=0 
fmnψmn(x1,x2).

4.4  SoV and integrability

As we have seen, SoV can provide a very constructive way to finding general solutions of some PDE's. Of course, a PDE in question should possess some unique property to allow application of the above technique. This extra quality is called integrability. In very rough terms it means existence of several commuting operators, like H1 and H2 above. This new notion will become clearer in the next example.

Now, we can give slightly more precise definition of SoV, when applied to the spectral problems like the one in the previous subsection.

Definition 2 SoV is a transformation of the multi-dimensional spectral problem into a set of 1-dimensional spectral problems.

4.5  Another SoV for the quantum problem

It might be surprising at first, but SoV is not unique. To demonstrate this, let us construct another solution of the same oscillator problem.

Consider functions Θ(u):
Θ(u) =  x12

u−a
+  x22

u−b
− 1 = −  (u−u1)(u−u2)

(u−a)(u−b)
,
(48)
where u,a,b ∈ R are some parameters, and u1, u2 are zeros of Θ(u). Taking the residues at u=a and u=b in both sides of (48), we have
x12 =  (u1−a)(u2−a)

b−a
,    x22 =  (u1−b)(u2−b)

a−b
.
(49)
The variables u1 and u2 are called elliptic coordinates in R2, because by definition they satisfy the equation
 x12

u−a
+  x22

u−b
=1
(50)
with the roots u1,u2 given by the equations
u1 + u2 = a + b + x12 + x22,    u1u2 = ab + bx12 + ax22.

Let all variables satisfy the inequalities
a < u1 < b < u2 < ∞.
(51)

Now introduce the functions
ψλ(x1,x2):=x1k1x2k2e−(x12+x22)/2 n

i=1 
Θ(λi),
where ki=0,1 and λ1,…,λnR are indeterminates.

Theorem 3 Functions ψλ(x1,x2) are eigenfunctions of the operator H iff {λi} satisfy the following algebraic equations


j ≠ i 
 1

λi−λj
 1

2
+
k1/2+  1

4

λi−a
+
k2/2+  1

4

λi−b
=0,    i=1,…,n.
(52)
Parameters λi have the following properties (generalized Stieltjes theorem): i) they are simple, ii) they are placed along the real axis inside the intervals (51), iii) they are the critical points of the function |G|,
G(λ1,…,λn) = exp(−  1

2
1+…+λn))
× n

p=1 
p−a)k1/2+1/4p−b)k2/2+1/4

r > p 
r−λp).
(53)

This Theorem gives another basis for the oscillator problem. In this case we have two commuting operators G1 and G2:
G1=−  ∂2

∂x12
+x12+  1

a−b
 
x1  ∂

∂x2
−x2  ∂

∂x1

2

 
,

G2=−  ∂2

∂x22
+x22  1

a−b
 
x1  ∂

∂x2
−x2  ∂

∂x1

2

 
,

[G1,G2]=0,
which are diagonal on this basis. Notice that
H=G1+G2.

?

5  Integrable systems and special functions

5.1  Introduction

Separation of variables (SoV) for linear partial differential operators of two varaibles Dx1,x2 can be defined by the following procedure. Assume that by some transformation we could transform the operator Dx1,x2 into the form:
Dx1,x2→ Dy1,y2=  1

ϕ1(y1)−ϕ2(y2)
 (D(1)y1 −D(2)y2),
(54)
where ϕi(yi) are some functions of one variable and D(i)yi are some ordinary differential operators. It could be done by changing variables (coordinate transform) {x}→ {y}, but it could also involve an integral transform. Then, we can introduce another operator Gy1,y2 such that
D(1)y1 − ϕ1(y1) Dy1,y2=Gy1,y2 = D(2)y2 − ϕ2(y2) Dy1,y2.
Notice, that D and G commute
[D,G]=0.
The operator G that appeared in the procedure of SoV is called operator of constant of separation.

The above definition is easily expanded to the more variable case. Essential step is to keep bringing the operator into the `separable form' (54) that allows to introduce more and more operators of constants of separation. If one could break the operator down to a set of one-variable operators, then separation of variables is done. This obviously requires that the number of operators G will be equal to the number of variables minus 1, and also that they commute between themselves and with the operator D. The latter condition defines an integrable system. So, we can say that necessary condition for an operator to be separable is that it can be supplemented by a full set of mutually commuting operators (G), or in other words the operator has to belong to an integrable family of operators.

As we have seen in the previous Lecture, special functions of one variable often appear when one separates variables in linear PDEs in attempt to find a general solution in terms of a large set of factorized partial (or separated) solutions. Usually, the completeness of the set of separated solutions can also be proved, so that we can indeed expand any solution of our equation, which is a multi-variable `special function', into a basis of separated (one-variable) special functions.

There are two aspects of this procedure. First, the separated functions of one variable will satisfy ODEs, so that we can, in principle, `classify' the initial multi-variable special function by the procedure of separation of variables and by the type of the obtained ODEs. It is clear that when some regularity conditions are set on the class of allowed transformations which used in a SoV, one should expect a good correspondence between complexity of both functions, the multivariable and any of the corresponding one-variable special functions.

In the example of isotropic harmonic oscillator, we had a trivial separation of variables first (in Cartesian coordinates), which gave us a basis as a product of Hermite polynomials. Hence, we might conclude that the operator H is, in a sense, a two-dimensional analogue of the hypergeometric differential operator, because one of its separated bases is given in terms of hypergeometric functions (Hermite polynomials).

Curiously, the second separation, in elliptic coordinates, led to the functions of the Heun type, which is beyond the hypergeometric class. The explanation of this seeming contradiction is that the operator H is `degenerate' in the sense that it separates in many coordinate systems. To avoid this degeneracy, one could disturb this operator by adding some additional terms that will break one separation, but will still allow the other.

Therefore generically, if an operator can be separated, it usually separates by a unique transformation, leading to a unique set of separated special functions of one variable.

The second aspect of the problem is understanding what are the sufficient conditions of separability? Or, which integrable operators can be separated and which can not? A very close question is: what class of transformations should be allowed when trying to separate an operator?

In order to demonstrate the point about a class of transformations, take a square of the Laplace operator:

 ∂2

∂x12
+  ∂2

∂x22

2

 
.
Of course, this operator is integrable, although there is a Theorem saying that one can not separate this operator by any coordinate transform {x}→ {y}:
y1=y1(x1,x2),       y2=y2(x1,x2).
It means that although one can find some partial solutions in the factorized form, they will never form a basis. There is no statement like that if one allows integral transforms, which means that this operator might be still separable in a more general sense.

It is interesting to note that the operators that are separable through a change of coordinates, although being very important in applications, constitute a very small subclass of all separable operators. Correspondingly, the class of special functions of many variables, which is related to integrable systems, is much larger then its subclass that is reducible to one-variable special functions by a coordinate change of variables.

Below we give an example of the integrable system that can not be separated by a coordinate change of variables, but is neatly separable by a special integral transform.

5.2  Calogero-Sutherland system

In this subsection, an integral operator M is constructed performing a separation of variables for the 3-particle quantum Calogero-Sutherland (CS) model. Under the action of M the CS eigenfunctions (Jack polynomials for the root system A2) are transformed to the factorized form φ(y1)φ(y2), where φ(y) is a trigonometric polynomial of one variable expressed in terms of the 3F2 hypergeometric series. The inversion of M produces a new integral representation for the A2 Jack polynomials.

The set of commuting differential operators defining the integrable system called (3-particle) Calogero-Sutherland model is generated by the following partial differential operators
H1
=
−i(∂1+∂2+∂3),
H2
=
−(∂12+∂13+∂23) −g(g−1)(sin−2q12+sin−2q13+sin−2q23),
H3
=
i∂123 +ig(g−1)(sin−2q23 ∂1 +sin−2q13 ∂2+sin−2q12 ∂3),

qij=qi−qj,       ∂i=  ∂

∂qi
,
or, by the equivalent set, acting on Laurent polynomials in variables tj=e2iqj, j=1,2,3:
~
H
 

1 
=
−i(∂1+∂2+∂3)
~
H
 

2 
=
−(∂12+∂13+∂23)
g[cotq12(∂1−∂2)+cotq13(∂1−∂3) +cotq23(∂2−∂3)]
−4g2
~
H
 

3 
=
i∂123
−ig[cotq12(∂1−∂2)∂3 +cotq13(∂1−∂3)∂2 +cotq23(∂2−∂3)∂1 ]
+2ig2[ (1+cotq12cotq13)∂1 +(1−cotq12cotq23)∂2 +(1+cotq13cotq23)∂3 ]
the vacuum function being
Ω(

q
 
)=|sinq12sinq13sinq23|g.
(55)
Their eigenvectors Ψn, resp. Jn,
Ψn(

q
 
)=Ω(

q
 
)Jn(

q
 
),
(56)
are parametrized by the triplets of integers {n1 ≤ n2 ≤ n3} ∈ Z3, the corresponding eigenvalues being
h1=2(m1+m2+m3),     h2=4(m1m2+m1m3+m2m3),     h3=8m1m2m3,
(57)
where,
m1=n1−g,        m2=n2,        m3=n3+g.
(58)

5.3  Integral transform

We will denote the separating operator acting on Ψn as K, and the one acting on the Jack polynomials Jn as M.

To describe both operators, let us introduce the following notation.
x1=q1−q3,        x2=q2−q3,        Q=q3,

x±=x1±x2,        y±=y1±y2.

We shall study the action of K locally, assuming that q1 > q2 > q3 and hence x+ > x.

The operator K:Ψ(q1,q2,q3)→~Ψ(y1,y2;Q) is defined as an integral operator
~
Ψ
 
(y1,y2;Q)=
y+

y 
dξ K(y1,y2;ξ) Ψ
 y+

2
+Q,  y+−ξ

2
+Q,Q
(59)
with the kernel
K


sin
  ξ+y

2

sin
  ξ−y

2

sin
  y+

2

sin
  y+−ξ

2


siny1siny2sinξ



g−1


 
(60)
where κ is a normalization coefficient to be fixed later. It is assumed in (59) and (60) that y < x=ξ < y+=x+. The integral converges when g > 0 which will always be assumed henceforth.

The motivation for such a choice of K takes its origin from considering the problem in the classical limit (g→∞) where there exists effective prescription for constructing a separation of variables for an integrable system from the poles of the so-called Baker-Akhiezer function.

Theorem 1 Let HkΨn1n2n3=hkΨn1n2n3. Then the function ~Ψn=KΨn satisfies the differential equations
Q
~
Ψ
 

n 
=0,       Yj
~
Ψ
 

n 
=0,     j=1,2
(61)
where
Q=−i∂Q−h1,
(62)

Yj=i∂yj3+h1yj2−i
h2+3  g(g−1)

sin2yj

yj

h3+  g(g−1)

sin2yj
h1 +2ig(g−1)(g−2)  cosyj

sin3 yj

.
(63)

The proof is based on the following proposition.

Proposition 2 The kernel K satisfies the differential equations
[−i∂Q−H1*]K=0,


i∂3yj+H1*2yj −i
H2*+  3g(g−1)

sin2 yj

yj

H3*+H1*  g(g−1)

sin2 yj
+2ig(g−1)(g−2)  cosyj

sin3 yj


K=0,
where Hn* is the Lagrange adjoint of Hn

ϕ(q)(Hψ)(q) dq=
(H*φ)(q)ψ(q) dq

H*1
=
i(∂q1+∂q2+∂q3),
H*2
=
−∂q1q2−∂q1q3−∂q2q3 −g(g−1)[sin−2q12+sin−2q13+sin−2q23],
H*3
=
−i∂q1q2q3−ig(g−1) [sin−2q23 ∂q1+sin−2q13 ∂q2 +sin−2q12q3].

The proof is given by a direct, though tedious calculation.

To complete the proof of the theorem 5.1, consider the expressions Qn and Yjn using the formulas (59) and (60) for K. The idea is to use the fact that Ψn is an eigenfunction of Hk and replace hkΨn by HkΨn. After integration by parts in the variable ξ the operators Hk are replaced by their adjoints Hk* and the result is zero by virtue of proposition 5.2.

The following theorem gives the separation of variables.

Theorem 3 The function ~Ψn1n2n3 is factorized
~
Ψ
 

n1n2n3 
(y1,y2;Q) = eih1Qψn1n2n3(y1n1n2n3(y2).
(64)
The factor ψn(y) allows further factorization
ψn(y)=(siny)2gφn(y)
(65)
where φn(y) is a Laurent polynomial in t=e2iy
φn(y)= n3

k=n1 
tk ck(

n
 
;g).
(66)
The coefficients ck(n;g) are rational functions of k, nj and g. Moreover, φn(y) can be expressed explicitely in terms of the hypergeometric function 3 F2 as
φn(y)=tn1(1−t)1−3g3 F2



a1,a2,a3
b1,b2
;t



(67)
where
aj=n1−n4−j+1−(4−j)g,        bj=aj+g.
(68)

Note that, by virtue of the theorem 5.1, the function ~Ψn(y1,y2;Q) satisfies an ordinary differential equation in each variable. Since Qf=0 is a first order differential equation having a unique, up to a constant factor, solution f(Q)=eih1Q, the dependence on Q is factorized. However, the differential equations Yjψ(yj)=0 are of third order and have three linearly independent solutions. To prove the theorem 5.3 one needs thus to study the ordinary differential equation

i∂y3+h1y2−i
h2+3  g(g−1)

sin2y

y

h3+  g(g−1)

sin2y
h1 +2ig(g−1)(g−2)  cosy

sin3 y


ψ = 0.
(69)
and to select its special solution corresponding to ~Ψ.

The proof will take several steps. First, let us eliminate from Ψ and ~Ψ the vacuum factors Ω, see (56), and, respectively
~
Ψ
 
(y1,y2;Q)=ω(y1)ω(y2)
~
J
 
(y1,y2;Q),        ω(y)=sin2gy.
(70)

Conjugating the operator K with the vacuum factors
M=ω1−1ω2−1KΩ:J→
~
J
 
(71)
we obtain the integral operator
~
J
 
(y1,y2;Q)=
y+

y 
dξ M(y1,y2;ξ) J
 y+

2
+Q,  y+−ξ

2
+Q,Q
(72)
with the kernel
M(y1,y2;ξ)=K(y1,y2;ξ)

 y+

2
+Q,  y+−ξ

2
+Q,Q

ω(y1)ω(y2)
= κsinξ

sin
 ξ+y

2

sin
 ξ−y

2


g−1

 

sin
 y+

2

sin
 y+−ξ

2


2g−1

 

[ siny1siny2]3g−1
.
(73)

Proposition 4 Let S be a trigonometric polynomial in qj, i.e. Laurent polynomial in tj=e2iqj, which is symmetric w.r.t. the transpositon q1↔ q2. Then ~S=MS is a trigonometric polynomial symmetric w.r.t.  y1↔ y2.

5.4  Separated equation

To complete the proof of the theorem 5.3 we need to learn more about the separated equation (69).

Eliminating from ψ the vacuum factor ω(y)=sin2gy via the substitution ψ(y)=φ(y)ω(y) one obtains
[i∂y3+(h1+6igcoty)∂y2
+(−i(h2+12g2)+4gh1coty+3ig(3g−1)sin−2y)∂y
+(−(h3+4g2h1)−2ig(h2+4g2)coty +g(3g−1)h1sin−2y)]φ = 0.
(74)

The change of variable t=e2iy brings the last equation to the Fuchsian form:
[∂t3+w1t2+w2t+w3]φ = 0
(75)
where
w1
=
3(g−1)+  1

2
h1

t
+  6g

t−1
,
w2
=
(3g2−3g+1)+  1

2
(2g−1)h1+  1

4
h2

t2
+  3g(3g−1)

(t−1)2
 g(9(g−1)+2h1)

t(t−1)
,
w3
=
g3+  1

2
g2h1+  1

4
gh2+  1

8
h3

t3
+
 1

2
g((h2+4g2)(t−1)−(3g−1)h1)

t2(t−1)2
.

The points t=0,1,∞ are regular singularities with the exponents
t ∼ 1   
ϕ ∼ (t−1)μ   
μ ∈ {−3g+2,−3g+1,0}
t ∼ 0   
ϕ ∼ tρ   
ρ ∈ {n1,n2+g,n3+2g}
t ∼ ∞   
ϕ ∼ t−σ   
−σ ∈ {n1−2g,n2−g,n3}

The equation (75) is reduced by the substitution φ(t)=tn1(1−t)1−3gf(t) to the standard 3F2 hypergeometric form
[t∂t(t∂t+b1−1)(t∂t+b2−1) −t(t∂t+a1)(t∂t+a2)(t∂t+a3)]f=0,
(76)
the parameters a1, a2, a3, b1, b2 being given by the formulas (68) which read
a1=n1−n3+1−3g,        a2=n1−n2+1−2g,        a3=1−g,

b1=n1−n3+1−2g,        b2=n1−n2+1−g.

Proposition 5 Let the parameters hk be given by (57), (58) for a triplet of integers {n1 ≤ n2 ≤ n3} and g ≠ 1,0,−1,−2,…. Then the equation (75) has a unique, up to a constant factor, Laurent-polynomial solution
φ(t)= n3

k=n1 
tk ck(

n
 
;g),
(77)
the coefficients ck(n;g) being rational functions of k, nj and g.

Proof. Consider first the hypergeometric series for 3F2 which converges for |t| < 1. Using for aj and bj the expressions (68) one notes that aj+1=bj+n4−j−n3−j and therefore
 (aj+1)k

(bj)k
=  (bj+k)n4−j−n3−j

(bj)n4−j−n3−j
.

The expression
 (a2)k(a3)k

(b1)k(b2)k
=  (b1+k)n3−n2(b2+k)n2−n1

(b1)n3−n2(b2)n2−n1
= Pn3−n1(k)
is thus a polynomial in k of degree n3−n1. So we have
3F2(a1,a2,a3;b1,b2;t) =

k=0 
 (a1)k

k!
Pn3−n1(k)tk
from which it follows that
3F2(a1,a2,a3;b1,b2;t) =
~
P
 

n3−n1 
(t)(1−t)3g−1
where ~Pn3−n1(t) is a polynomial of degree n3−n1 in t.

To prove now the proposition 5.5 it is sufficient to notice that the hypergeometric series 3F2(a1,a2,a3;b1,b2;t) satisfies the same equation (76) as f(t) and therefore the Laurent polynomial Fn(t) constructed above satisfies the equation (75). The uniqueness follows from the fact that all the linearly independent solutions to (75) are nonpolynomial which is seen from the characteristic exponents.

Now everything is ready to finish the proof of the theorem 5.3. Since the function ~Jn1 n2 n3 (y1, y2; Q) satisfies (74) in variables y1,2 and is a Laurent polynomial it inevitably has the factorized form
~
J
 

n1n2n3 
(y1,y2;Q) = eih1Qφn1n2n3(y1n1n2n3(y2)
(78)
by virtue of the proposition 5.5.

5.5  Integral representation for Jack polynomials

The formula (78) presents an interesting opportunity to construct a new integral representation of the Jack polynomial Jn in terms of the 3F2 hypergeometric polynomials φn(y) constructed above. To achieve this goal, it is necessary to invert explicitely the operator M:J→~J.

It is possible to show that this problem (after changing variables) is equivalent to the problem of finding an inverse of the following integral transform:
~
s
 
)=
η+

η 
   (ξ−η)g−1

Γ(g)
s(ξ)
(79)
which is known as Riemann-Liouville integral of fractional order g. Its inversion is formally given by changing sign of g
s(ξ)=
ξ+

ξ 
   (η−ξ)−g−1

Γ(−g)
~
s
 
)
(80)
and is called fractional differentiation operator. We will not give details of this calculation, just giving the final result. The formula for M−1:~J→ J is
J(x+,x;Q)=
x+

x 
dy  \checkM(x+,x;y)
~
J
 
(x+,y;Q)
(81)

\checkM=\checkκ
siny
sin
 x++y

2

sin
 x+−y

2


3g−1

 


sin
 y+x

2

sin
 y−x

2


g+1

 
[sinx1sinx2]2g−1
(82)
where
\checkκ =  Γ(2g)

2Γ(−g)Γ(3g)
.
(83)
The operators M (and M−1) are normalized by M: 1→ 1.

For the kernel of K−1 we have respectively
\checkK=\checkκ
sing xsiny
sin
 x++y

2

sin
 x+−y

2


g−1

 


sin
 y+x

2

sin
 y−x

2


g+1

 
[sinx1sinx2]g−1
.
(84)

The formulas (78), (81), (82) provide a new integral representation for Jack polynomial Jn of three variables in terms of the 3F2 hypergeometric polynomials φn(y). It is remarkable that for positive integer g the operators K−1, M−1 become differential operators of order g. In particular, for g=1 we have K−1=∂/∂y.

References

[1]
George E. Andrews, Richard Askey, and Ranjan Roy. Special functions. Cambridge University Press, Cambridge, 1999.

[2]
Willard Miller, Jr. Lie Theory and Special Functions. Academic Press, New York, 1968. Mathematics in Science and Engineering, Vol. 43.

[3]
N. Ja. Vilenkin. Special Functions and the Theory of Group Representations. American Mathematical Society, Providence, R. I., 1968. Translated from the Russian by V. N. Singh. Translations of Mathematical Monographs, Vol. 22.

Index (showing section)

(Gauss) hypergeometric function, 2.2
(generalized) hypergeometric series, 2.2

Chebyshev polynomials, 3.5
Christoffel-Darboux formula, 3.2
classical orthogonal polynomials, 3.5

fractional linear transformation, 2.7

gamma function, 1.2
Gauss quadrature, 0.0, 3.4
Gegenbauer polynomials, 3.5

Hermite polynomials, 3.5
hypergeometric differential equation, 2.6
hypergeometric function, 1.1
hypergeometric series, 1.1

integrable system, 5.1
integral transform, 4.1, 5.1

Jacobi polynomial, 3.1
Laguerre polynomials, 3.5
Legendre polynomials, 3.5

Rodrigues formula, 3.5

Separation of Variables, 4.1
SoV, 4.4



File translated from TEX by TTH, version 3.13.
On 22 May 2003, 12:56.