PG course on
SPECIAL FUNCTIONS AND THEIR SYMMETRIES

Vadim KUZNETSOV

Course Outline

1  Gamma and Beta functions

1.1  Introduction

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

Definition 1 A hypergeometric series is a series ∑_n=0^∞ a_n with a_n+1/a_n 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 XX^th 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.

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!.

Several other useful forms of the beta integral can be derived by a change of variables. For example, take t=sin²θ 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 l₁(t) and l₂(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 I_x,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 ∫_−∞⁽⁰⁺⁾. 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 z^x−1 has its principal value. Thus

−2πi= ⌠ ⌡ π −π   iRxeixθ 1−Reiθ  dθ+ ⌠ ⌡ ε R   tx−1eixπ 1+t  dt+ ⌠ ⌡ −π π   iεxeixθ 1−εeiθ  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

Almost all of the elementary functions of mathematics and some not very elementary, like the error function erf(x) and dilogarithm function Li₂(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 ∑c_n, 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: _pF_q(a₁,…, a_p;b₁,…,b_q;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 _pF_q(a₁,…, a_p; b₁,…,b_q;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 |c_n+1/c_n|→0 as n→∞ if p ≤ q. For p=q+1, lim_n→∞ |c_n+1/c_n|=|x|, and for p > q+1, |c_n+1/c_n|→∞ 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+1F_q(a₁,…, a_q+1; b₁,…,b_q;x) with |x|=1 converges absolutely if ℜ (∑b_i−∑a_i) > 0. The series converges conditionally if x=e^iθ ≠ 1 and 0 ≥ ℜ (∑b_i−∑a_i) > −1 and the series diverges if ℜ (∑b_i−∑a_i) ≤ −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+1F_q

(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 ₂F₁ 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      ;

ex= 0F0      − − ;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 ₂F₁(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 ₂F₁ function has a continuation to the complex plane with branch points at 1 and ∞.

Definition 4 The (Gauss) hypergeometric function ₂F₁(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 ₂F₁(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)] ₂F₁(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 ₂F₁. 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 ₂F₁(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 xⁿ 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 ₂F₁ 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 ₂F₁ 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 ₂F₁ 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 (α₁,α₂), (β₁,β₂) and (γ₁,γ₂). 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),      (α1+β1+γ1, α1+β2+γ1).

It is easy to check that we arrive at the hypergeometric differential equation (29) for the function F(x) with the following parameters:

a=α1+β1+γ1,      b=α1+β2+γ1,      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

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 {p_n(x)}₀^∞, where p_n(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 {p_n(x)} satisfies the three-term recurrence relation

pn+1(x)=(Anx+Bn)pn(x)−Cnpn−1(x)        for    n ≥ 0,

where we set p₋₁(x)=0. Here A_n, B_n, and C_n are real constants, n=0,1,2,…, and A_n−1A_nC_n > 0, n=1,2,…. If the highest coefficient of p_n(x) is k_n, 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 p_n(x) are normalized so that h_n=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 p_n+1′(x)p_n(x)−p_n+1(x)p_n′(x) > 0        ∀x.

3.3  Zeros of orthogonal polynomials

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

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

3.4   Gauss quadrature

Theorem 8 There are positive numbers λ₁,…,λ_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 x_j, j=1,…,n, are zeros of the polynomial p_n(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 {p_n′} is also an orthogonal system;

(ii) p_n 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 K_n 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

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 √ 2π   ⌠ ⌡ ∞ −∞  eixye−y2/2Hn(y) dy=ine−x2Hn(x).

Hence, e^−x2H_n(x) is an eigenfunction of the Fourier transform with eigenvalue iⁿ. This can be proved by using the Rodrigues formula for H_n(x).

4  Separation of variables and special functions

4.1  Introduction

Definition 1 Separation of Variables M is a transformation which brings a function ψ(x₁,…,x_n) of many variables into a factorized form

M: ψ→ ϕ1(y1)·…·ϕn(yn).

Functions ϕ_j(y_j) 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 q₁(x) and q₂(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 q₁(x) and q₂(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 ψ(x₁,x₂) ∈ L²(R²) 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 L²(R²) 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 R²

⌠ ⌡ ∞ −∞  ψn1n2(x1,x2)ψm1m2(x1,x2) = an1n2δn1m1δn2m2.

Every function f ∈ L²(R²) 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 H₁ and H₂ 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 u₁, u₂ 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 u₁ and u₂ are called elliptic coordinates in R², because by definition they satisfy the equation

x12 u−a +  x22 u−b =1

(50)

with the roots u₁,u₂ 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 k_i=0,1 and λ₁,…,λ_n ∈ R are indeterminates.

Theorem 3 Functions ψ_→λ(x₁,x₂) 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/4(λp−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 G₁ and G₂:

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 D_x1,x2 can be defined by the following procedure. Assume that by some transformation we could transform the operator D_x1,x2 into the form:

Dx1,x2→ Dy1,y2=  1 ϕ1(y1)−ϕ2(y2)  (D(1)y1 −D(2)y2),

(54)

where ϕ_i(y_i) are some functions of one variable and D⁽ⁱ⁾_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 G_y1,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 A₂) are transformed to the factorized form φ(y₁)φ(y₂), where φ(y) is a trigonometric polynomial of one variable expressed in terms of the ₃F₂ hypergeometric series. The inversion of M produces a new integral representation for the A₂ 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 = −(∂1∂2+∂1∂3+∂2∂3) −g(g−1)(sin−2q12+sin−2q13+sin−2q23), H3 = i∂1∂2∂3 +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 t_j=e^2iq_j, j=1,2,3:

~ H   1  = −i(∂1+∂2+∂3) ~ H   2  = −(∂1∂2+∂1∂3+∂2∂3) g[cotq12(∂1−∂2)+cotq13(∂1−∂3) +cotq23(∂2−∂3)] −4g2 ~ H   3  = i∂1∂2∂3 −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. J_→n,

Ψ→n( → q   )=Ω( → q   )J→n( → q   ),

(56)

are parametrized by the triplets of integers {n₁ ≤ n₂ ≤ n₃} ∈ Z³, 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 J_→n 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 q₁ > q₂ > q₃ and hence x₊ > x₋.

The operator K:Ψ(q₁,q₂,q₃)→~Ψ(y₁,y₂;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 H_kΨ_n1n2n3=h_kΨ_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+h1∂yj2−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 H_n^* is the Lagrange adjoint of H_n

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

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