Stokes phenomenon for the prolate spheroidal wave equation

Applied Numerical Mathematics 60 (2010) 1309–1319

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Stokes phenomenon for the prolate spheroidal wave equation F. Fauvet a , J.-P. Ramis b , F. Richard-Jung c,∗ , J. Thomann a a b c

Mathématiques – IRMA, Université de Strasbourg, 67084 Strasbourg cedex, France Laboratoire Emile Picard, Université Paul Sabatier, 31062 Toulouse, France Laboratoire Jean Kuntzmann, Université Joseph Fourier, BP 53, 38041 Grenoble cedex, France

a r t i c l e

i n f o

Article history: Received 21 December 2009 Received in revised form 10 May 2010 Accepted 20 May 2010 Available online 25 May 2010 Keywords: Stokes matrices Ordinary differential equations Asymptotics Summability Computer algebra Galois differential theory

a b s t r a c t The prolate spheroidal wave functions can be characterized as the eigenfunctions of a differential operator of order 2:



t2 − τ 2



ϕ  + 2t ϕ  + σ 2 t 2 ϕ = μϕ .

In this article we study the formal solutions of this equation in the neighborhood of the singularities (the regular ones ±τ , and the irregular one, at infinity) and perform some numerical experiments on the computation of Stokes matrices and monodromy, using formal/numerical algorithms we developed recently in the Maple package Desir. This leads to the following conjecture: the series appearing in the formal solutions at infinity, depending on the parameter μ, are in general divergent; they become convergent for some particular values of the parameter, corresponding exactly to the eigenvalues of the prolate operator. We give the proof of this result and its interpretation in terms of differential Galois groups. © 2010 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction The main aim of this paper is to present a new result on the Stokes phenomena at infinity of prolate spheroidal wave equations in relation with the spectral theory for these equations. More precisely, we will prove that the power series expansions appearing in the asymptotic expansions at infinity are in general divergent but become convergent for some particular values of the spectral parameter μ corresponding exactly to the spectral values of the corresponding spheroidal operator. In [21], David Slepian wrote the following about spheroidal wave functions: Most of us feel that there is something deeper here than we currently understand – that there is a way of viewing these problems more abstractly that will explain their elegant solution in a more natural and profound way, so that these nice results will not appear so much as a lucky accident. We think that our main result is a first step in the direction suggested by D. Slepian. Today there is an abstract, natural and profound way to understand linear differential equations with polynomial coefficients, the Riemann–Hilbert theory and its generalization to the irregular case (cf. [25], Malgrange [13]). It is based on the properties of solutions in the complex domain, in particular on their monodromy and Stokes phenomena. (In the Fuchsian

*

Corresponding author. E-mail addresses: [email protected] (F. Fauvet), [email protected] (J.-P. Ramis), [email protected] (F. Richard-Jung), [email protected] (J. Thomann). 0168-9274/$30.00 © 2010 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2010.05.010

1310

F. Fauvet et al. / Applied Numerical Mathematics 60 (2010) 1309–1319

case, monodromy is sufficient but in the irregular case it is necessary to “add” the Stokes phenomena, anyway in all cases it is necessary to “leave” the real domain.) Our feeling is that it is possible to go in the direction suggested by D. Slepian using this generalized version of Riemann–Hilbert theory and its transformation under Fourier transform. (This opinion could be confirmed by the recent success of Riemann–Hilbert methods in special functions theory.) Our main result was first guessed by the second author, then this intuition was confirmed by experiments using a mixed formal-numerical package (DESIR) developed for the study of general algebraic linear differential equations in the complex domain and finally we found a proof of the result. In the present paper we will begin by the statement of the main theorem and its proof, afterwards we will “reinforce” the proof using the DESIR package. This step is important because it is an example of what could be done in more general cases: using the DESIR package it is possible to study the monodromy and Stokes phenomena not only for the solutions of the spheroidal wave equations but more generally for the complex solutions of an arbitrary algebraic Sturm–Liouville equation defined on the real line (theoretically there is a relation between the eigenvalues for the solutions of such an equation and the Stokes phenomena at infinity but it is more complicated than in the spheroidal case). DESIR package allows one to work in the complex domain and with general classes of equations, for the particular case of spheroidal functions on the real line there exists numerical algorithms more efficient than DESIR algorithms. Since the 80s a lot of algorithmic work has been done (in Grenoble) by the team of J. Della Dora in the domain of Computer Algebra and Differential equations, among numerous projects one was named DESIR, an acronym for Differential Equations Singularities Irregular and Regular [8,16]. From the beginning this was developed in collaboration with the team of the second author (in Strasbourg) on the theoretical counterpart: k-summability, Gevrey asymptotics and Stokes phenomena. The name DESIR is now given to a package in MAPLE [20], which can be divided into three parts. – The first part contains the formal algorithms for computing a basis of formal solutions for a linear homogeneous ODE with polynomial coefficients, in particular in the neighborhood of its singularities in the complex plane. From the theoretical point of view, these formal solutions are completely understood thanks in particular to the theory of Gevrey asymptotics k and multi-summability [3,17,14]. From the numerical point of view the theory leads to some algorithms of summation of divergent series generalizing the well-known Borel–Laplace transform [22,23], allowing to compute actual solutions on sectors whose vertex is a singular point. – The second part of the DESIR package collects some functions for the numerical computation of k-summable series and the graphical visualization of the solutions in the complex plane [18,19]. – The third part of the package was developed recently by F. Fauvet, F. Richard-Jung and J. Thomann. It deals with Stokes phenomenon at an irregular singular point and it contains a mixed formal/numerical method for the computation of the Stokes matrices [10–12]. In Section 1 of the paper we will recall definition and basic properties of the spheroidal wave functions. In Section 2, using DESIR code, we will describe solutions at the singular points, including infinity. The main theorem will be stated and proved in Section 3. In Section 4, we will explain how to use DESIR code to “rediscover” the main results. Finally in Section 5 we will “translate” the main theorem in terms of differential Galois groups. 2. The prolate spheroidal wave functions: definition and basic properties The literature on these objects is large, from the mathematical point of view as much as from the applied point of view. For this brief presentation, we refer to [26] or [9]. The prolate spheroidal wave functions, (PSWFs) {ϕn,σ ,τ }, constitute an orthonormal basis of the space of σ -bandlimited functions on the real line, i.e. functions whose Fourier transforms have support on the interval [−σ , σ ]. The PSWFs are maximally concentrated on the interval [−τ , τ ] (as precised below) and depend on positive parameters σ and τ . They can be characterized by one of the following properties:

• as the eigenfunctions of a differential operator arising from a Helmholtz equation on a prolate spheroid:





τ 2 − t 2 ϕn,σ ,τ − 2t ϕn ,σ ,τ − σ 2 t 2 ϕn,σ ,τ = μn,σ ,τ ϕn,σ ,τ ;

• as the maximum energy concentration of a σ -bandlimited function on the interval [−τ , τ ]; that is, ϕ0,σ ,τ is the function of total energy 1 (= ϕ0,σ ,τ 2 ) such that

τ

   f (t )2 dt

−τ

is maximized, ϕ0,σ ,τ , etc.;

ϕ1,σ ,τ is the function with the maximum energy concentration among those functions orthogonal to

F. Fauvet et al. / Applied Numerical Mathematics 60 (2010) 1309–1319

1311

• as the eigenfunctions of an integral operator with kernel arising from the sinc function S (t ) = sin(π t )/π t:

σ π

τ −τ





σ ϕn,σ ,τ (x) S (t − x) dx = λn,σ ,τ ϕn,σ ,τ (t ). π

For our study, we concentrate on the first characterization in terms of eigenfunctions of a differential operator of order 2. 3. Solutions in the complex plane 3.1. Classification of the singularities in the complex plane, formal solutions Let us consider the differential equation:

L σ ,τ (ϕ ) = μϕ ,



where L σ ,τ (ϕ ) = t 2 − τ 2



ϕ  + 2t ϕ  + σ 2t 2 ϕ .

Here the parameters σ , τ and μ can be any complex numbers. We will give the first results in the general case, in order to see how the parameters influence the singularities, and the shape of the solutions near the singularities. This equation has only two singularities at finite distance ±τ . In the neighborhood of each of these points, we have a basis of solutions constituted of a regular (holomorphic at the singularity) function f and a solution of the form f (t ) log(t ± τ ) + g (t ), where g is also regular [15]. For example, in the neighborhood of the point τ , we can use the Maple package Desir to obtain:



t(x) − τ = x, y(x) = 1 −



(−μ + σ 2 τ 2 )x 2τ

  + O x2 ,

    (−1 − 2μ + 2σ 2 τ 2 )x (−μ + σ 2 τ 2 )x + t(x) − τ = x, y(x) = log(x) 1 − + O x2 log(x) . 2τ

2τ

In general, in the neighborhood of a point t 0 , the solutions are computed using the rational Newton algorithm [4], this means that we obtain “generalized formal solutions”: they are parametrized by a new variable x and have the form (t (x) − t 0 = Λxr , y (x) = exp( Q (1/x))xλ Φ(x)). In this expression, Λ is a complex number (useful in order to reduce the algebraic extension needed), r is a positive integer, called the ramification, Q is a polynomial without constant term, λ is a complex number, called the exponent and, xλ Φ(x) is the regular part of the formal solution. Φ(x) is a polynomial in log(x) with power series coefficients. In our example, as τ is a regular singularity, the parametrization is only a translation t (x) − τ = x, and Q = 0. For all values of the parameters, except σ = 0, the equation has an irregular singularity at infinity, with a basis of formal solutions constituted of:



t(x) =

1 x

, y(x) = e −

_

RootOf(σ 2 + Z 2 ) x



x 1−

1 RootOf(σ 2 + _Z 2 )(μ − σ 2 τ 2 )x 2

σ2

  + O x2



.

Remark that we have above a condensed shape representing two solutions, corresponding to the two possible values of the RootOf: ±iσ . The parametrization is now given by the change of variable t = 1/x, the ramification is trivial, and the two series which appear in these solutions are a priori divergent, but 1-summable in each direction but ±iσ R+ . Of course, we give here only the first terms of the series, but the following ones can be generated using a recurrence formula, so (theoretically) we can obtain as much terms as wanted in all the series appearing in the solutions. Our goal is now to have a more precise idea on the way these series diverge, depending on the parameter μ, and for this to compute the Stokes multipliers of the equation in this basis, expecting them to vanish for the known values of the eigenvalues. In the following sections, we assume that the parameters σ and τ are real positive numbers. Moreover, we specialize to particular values of the two parameters: σ = τ = 1. Indeed, the change of variable t = τ u transforms the operator L σ ,τ into the operator L σ τ ,1 , and the general case for σ can be treated in the same way as σ = 1. We consider thus the following differential operator:









D μ : y −→ t 2 − 1 y  + 2t y  + t 2 − μ y . D μ has two finite regular singularities, at 1 and −1. In the neighborhood of the point 1, we find a basis of solutions:

f μ (t ), where

f μ (t ) log(t − 1) + ϕμ (t ),

1312

F. Fauvet et al. / Applied Numerical Mathematics 60 (2010) 1309–1319



μ

f μ (t ) = 1 +

+

−

2



1

11 288



 (t − 1) +

2

7

μ+

288

288

1

2

16

11

−

1

μ − μ− 4

μ2 +

1 288



1 16

 (t − 1)2 



μ3 (t − 1)3 + O (t − 1)4 ,

and



ϕμ (t ) =

   1 1 3 − μ (t − 1) + μ+ − μ2 (t − 1)2

1 2

2

16

16

    17 31 47 2 11 3 + − − μ+ μ − μ (t − 1)3 + O (t − 1)4 . 216

864

432

864

In the neighborhood of the point −1, we find a basis:

g μ (t ),

g μ (t ) log(t + 1) + ψμ (t ),

where



1

g μ (t ) = 1 +

2



+ −

−

11 288

μ



 (t + 1) +

2

μ−

7

16

11

+

288

1

288

1

1

4

16

μ2 − μ − 1

2

μ −

 3

288

μ

 (t + 1)2

  (t + 1)3 + O (t + 1)4 ,

and



   1 1 3 2 ψμ (t ) = μ − μ+ − μ (t + 1)2 (t + 1) + 1

 +

2

2

17

31

+

216

864

μ−

16

47 432

μ2 +

16

11 864







μ3 (t + 1)3 + O (t + 1)4 .

All the series f μ , ϕμ , g μ , ψμ are convergent, with a radius of convergence at least 2. D μ has an irregular singularity at ∞, with a basis of formal solutions constituted of:

yˆ 1 (t ) = yˆ 2 (t ) =

e −it

 1−

t

e

it 

t

1+

1 i( 2

μ − 1) t

1 i( 2

μ − 1) t

+

+

− 18 μ2 + 12 μ +

− 18

t2 1 2

μ + μ+ 2

t2

1 8

1 8

+

−

1 i(11 48

1 i(11 48

μ + 7 − 11μ2 + μ3 ) t3

μ + 7 − 11μ + μ ) 2

t3

3

  1 +O 4 , 

+O

t

1 t4

 .

Thus the formal monodromy at infinity is trivial. This fact could be verified by calling our function monodromy implemented in the Desir package, provided that a numerical value is assigned to the parameter μ. 3.2. Stokes phenomenon Now we want to compute the Stokes matrices in the previous basis. Let ˆf 1 be the series e i t yˆ 1 and ˆf 2 the series e −i t yˆ 2 .1 Both are divergent. The first one is 1-summable in all directions but −iR+ (its Borel transform has just one singularity at −2i), the second one is 1-summable in all directions but iR+ (its Borel , defined transform has just one singularity at 2i). By summation in a direction θ in ]− π2 , π2 [, we obtain actual solutions y − i defined in the sector in the sector ]−π , π [, and by summation in a direction θ in ] π2 , 32π [, we obtain actual solutions y + i + ]0, 2π [. The functions y − and y are defined together on the sector ] 0 , π [ : the Stokes matrix associated to the Stokes ray i i π is the transition matrix which expresses the change of basis 2









π

− + + 2 y− 1 , y2 = y1 , y2 S .

+ − + + As y − 1 = y 1 , this matrix is an upper triangular matrix. Moreover, as y 2 − y 2 is a multiple of y 1 , it is of the form:

S

1

π 2



=

1 0

α 1



.

Of course, all the formal series yˆ 1 , yˆ 2 , ˆf 1 , ˆf 2 depend on

μ, even if it is not explicit in their name.

F. Fauvet et al. / Applied Numerical Mathematics 60 (2010) 1309–1319

1313

For similar reasons, we obtain a Stokes matrix associated to the Stokes ray − π2 , of the form:

S



− π2

=

1

β

0 1



.

We now state our main result: Theorem 1. Let μ be in C. The following properties are equivalent: 1.

2

μ is an eigenvalue of the differential operator L 1,1 = (t 2 − 1) dtd 2 + 2t dtd + t 2 ; more precisely, there exists a solution to the d 2 (t − 1) dy + t 2 y = μ y on [−1, 1], with the boundary conditions: y (−1), y (1) finite2 ; corresponding Sturm–Liouville problem dt dt

2. the series f μ and g μ are entire functions;

3. the series ˆf 1 and ˆf 2 are convergent; 4. the Stokes phenomenon is trivial (the Stokes matrices in the basis [ yˆ 1 , yˆ 2 ] are identity); 5. the monodromy around [−1, 1] is trivial. Some parts of this result, concerning essentially the properties of the series f μ and g μ in relation with the existence of the PSWFs, are already known: we recall them to have a complete description of the behavior of the solutions near the singularities. The original part of the result concerns the divergent series yˆ 1 and yˆ 2 and the Stokes phenomenon at infinity. Proof. We begin with the easy points.

• Properties (3) and (4) above are equivalent. It is obvious (if the Stokes phenomenon is trivial, the Borel transforms of ˆf i have no singularities, that means that ˆf i are convergent). • Properties (4) and (5) are equivalent. Indeed, as the formal monodromy at infinity is trivial, the monodromy around π π the segment [−1, 1] is conjugated to the product of the two matrices S 2 and S − 2 , namely S

π 2

S

− π2



=

1 + αβ

α

β

1



.

And this matrix is the identity matrix if and only if α = β = 0. • Property (5) implies property (2). Indeed, if the monodromy around the segment [−1, 1] is trivial, then the holomorphic function at 1, f μ , can be continued at −1, and has no other singularity elsewhere, thus is entire. Of course, for similar reasons, g μ is also entire.

• Property (2) implies property (1). Indeed, if f μ is an entire function, we consider its restriction y μ to the segment

1 [−1, 1] and L 1,1 as an operator on C ∞ ([−1, 1], C), equipped with the classical inner product  f , g = −1 f (t ) g (t ) dt. Then L 1,1 ( y μ ) = μ y μ . Moreover, using integration by parts, we verify that









L 1,1 ( y μ ), y μ = y μ , L 1,1 ( y μ ) ,

¯ ) y μ , y μ = 0, and so that (μ − μ to the value μ.

μ is real. This means that ( y μ ) is an eigenfunction of the operator L 1,1 associated

The part of the proof, which is more difficult, is the converse, i.e. that property (1) implies property (5). So assume that μ is an eigenvalue of the operator L 1,1 . The first lemma gives the relation between a corresponding eigenfunction and the solutions of D μ in the complex variable. Lemma 1. An eigenfunction, associated to the eigenvalue μ, is the restriction on R of an entire function f , solution of D μ . Proof. Let y be an eigenfunction. The behavior of any solution on ]−1, 1[ is of type φ − log + ψ − in a neighborhood of −1 and φ + log + ψ + in a neighborhood of 1. Using the boundary conditions we get, for the solution y, φ + = φ − = 0. The functions ψ − and ψ + are the restrictions of analytic functions respectively on complex neighborhoods of −1 and 1, therefore y is the restriction of an entire function f , and f is a solution of D μ . 2 The continuation of the proof consists in exhibiting another solution, linearly independent of f , that has such a shape that makes it possible to explicitly compute the action of the monodromy around the segment [−1, 1], and to check its triviality. Lemma 2. The function f is even or odd.

2

This property implies that

μ is real.

1314

F. Fauvet et al. / Applied Numerical Mathematics 60 (2010) 1309–1319

Proof. It is known that the eigenspaces of L 1,1 , associated to the eigenvalues μ, are of dimension 1. As the operator D μ is invariant under the transformation, which changes t into −t, the function g (x) = f (−x) is also an eigenfunction, thus there exists a scalar ν such that g = ν f . But then the characterization of the prolate functions in terms of maximization of the L 2 norm on [−τ , τ ] imposes that ν = ±1 and we are done. 2 We give the end of the proof, assuming that f is even (the proof will be mostly the same if f is odd). Considering the behavior of f in the neighborhood of −1, we see that necessarily f is a multiple of g μ . For similar reasons, f is also a multiple of f μ . So in the neighborhood of −1, besides the solution f , we know another solution h1 (t ) = f (t ) log(t + 1) + u (t ), with u holomorphic at −1. We consider the behavior of h1 in the neighborhood of 1:





h1 (t ) = af μ (t ) + b f μ (t ) log(t − 1) + ϕμ (t ) , thus

u (t ) = bf μ (t ) log(t − 1) + af μ (t ) + bϕμ (t ) − f (t ) log(t + 1) = b f (t ) log(t − 1) + v (t ), where v is holomorphic at 1, by construction. But we have also v (t ) = u (t ) − b f (t ) log(t − 1), which proves that v is holomorphic at −1. And so, v is an entire function. Now, h1 (−t ) is also a solution of D μ (remember that D μ is invariant under the change of variable, which changes t into −t). But

h1 (−t ) = f (−t ) log(−t + 1) + u (−t ), this means that h2 (t ) = f (t ) log(t − 1) + u (−t ) is also a solution (the difference is ±iπ f , which is a solution). Finally:

h1 (t ) − h2 (t ) = f (t ) log

t+1 t−1

+ u (t ) − u (−t )

  t+1 = 1 − b f (t ) log + v (t ) − v (−t ) ± iπ b f (t ). t−1 1 Thus, as announced, we obtain another (linearly independent) solution h = f (t ) log tt + −1 + g (t ), where g is an entire function. The action of a single loop around the segment [−1, 1] leaves both f and g unchanged, as they are entire, and so the monodromy acts trivially on h, because the contributions of 2iπ in the numerator and denominator neutralize one another. The monodromy operator is the identity. 2

4. Experimental results on Stokes phenomenon and monodromy In this section, we make a report on the numerical experiments we have done to confirm the initial idea on the vanishing of the Stokes multipliers for the eigenvalues of the prolate operator, before the proof of Theorem 1. 4.1. First Stokes matrices In order to compute numerically the Stokes multiplier α , we apply on this example the formal/numerical algorithms presented in [12] and implemented in the Desir package. We only recall that these algorithms are based on the analysis of the singularities of the Borel transform of the divergent series, and this for a large class of differential equations of single rank k (even if k is not 1): the method works under the hypotheses that the Borel transforms don’t have irregular singularities and don’t have many singularities aligned on a half line issued from the origin, but it allows in the Borel plane any polar, ramified or logarithmic singularities. The main function is StokesMatrices. It takes two arguments: a list (basis) of formal solutions, and a list of two integers, fixing the order of the Padé approximants used to sum the Borel transform of the divergent series. The output is a description of the Stokes phenomenon as a list of objects [arg, M_arg]: arg is the argument of a Stokes ray and M_arg is the corresponding Stokes matrix. The first experiment consists in computing the Stokes matrices for the first eigenvalue: we can find its approximate value for example in [1] (our equation corresponds to m = 0, c = 1), μ0 = 0.319. The result is the following (resu is a list containing yˆ 1 and yˆ 2 ):

> StokesMatrices(subs(mu = 319/1000, resu), [8, 8]);   π 1 0 , − , 0.1095087321 10−23 − 0.4342357172 10−6 I 1 2   π 1 −0.1095087321 10−23 − 0.4342357172 10−6 I , . 2

0

1

Thus we obtain matrices which are very close to identity.

F. Fauvet et al. / Applied Numerical Mathematics 60 (2010) 1309–1319

1315

Fig. 1.

4.2. A curve α (μ) The second experiment consists in computing the coefficient α for many values of μ. We observe that the real part of α is numerically very small (less than 10−9 ) and we draw the curve (α ) in function of μ, for μ ∈ ]0, 7]. We obtain the following curve (see Fig. 1). Again, we find in [1] the values of the second and third eigenvalues: the coefficient α annihilates for the first three eigenvalues.

μ1 = 2.593084 and μ2 = 6.533471. It seems that

4.3. Monodromy around [−1, 1] The third experiment consists in computing the monodromy around [−1, 1]. Let us explain how to do this with our formal/numerical tools, available in the package Desir. Let us denote by y 1 and y 2 two series (linearly independent) solutions in the neighborhood of the origin. They are convergent, with a radius of convergence at least 1. Using the connect function, we perform the numerical connection at the point −1/2 with the basis [ g μ , g μ log + ψμ ], that is, we obtain a numerical approximation of the coefficients a−1 , b−1 such that:





y 1 (−1/2) = a−1 g μ (−1/2) + b−1 g μ (−1/2) log(1/2) + ψμ (−1/2) , and a numerical approximation of the coefficients c −1 , d−1 such that:





y 2 (−1/2) = c −1 g μ (−1/2) + d−1 g μ (−1/2) log(1/2) + ψμ (−1/2) . Let M be the matrix



a −1 c −1

b −1 d −1



,

ˆ the matrix of formal monodromy around −1 in the basis [ g μ , g μ log + ψμ ]: and M

ˆ = M



1 2iπ 0 1



.

Let y˜ i the function obtained by analytic continuation of y i along a path turning anticlockwise around −1. Then ( y 1 , y 2 ) = ˆ M , and we put M −1 = M −1 M ˆ M. ( y˜ 1 , y˜ 2 ) M −1 M In the same manner, we obtain a matrix N of connection between the basis [ y 1 , y 2 ] and the basis [ f μ , f μ log + ψμ ],

ˆ N (the formal monodromy is the same). and we express the matrix of monodromy around 1 by: M 1 = N −1 M Then we obtain the monodromy around [−1, 1] by multiplying: M 1 M −1 . On the following graph, we measure the distance between M 1 M −1 and Id by drawing  M 1 M −1 − Id∞ (the sup norm is the maximum of the coefficients of the matrix in absolute value) (see Fig. 2). This plot, albeit less striking than the previous one, shows that the monodromy is trivial for the first three eigenvalues.

1316

F. Fauvet et al. / Applied Numerical Mathematics 60 (2010) 1309–1319

Fig. 2.

Fig. 3.

4.4. Connection between 1 and −1 The last experiment consists in computing the connection between the basis [ f μ , f μ log + ϕμ ] and the basis [ g μ , g μ log + ψμ ]. Indeed, for the particular values of the parameter μ which correspond to the eigenvalues, the differential operator D μ admits an entire function solution (see Lemma 1 below). This solution is necessarily a multiple of f μ in the neighborhood of 1 and a multiple of g μ in the neighborhood of −1. In other words, this means that both functions f μ and g μ are in fact entire functions, or still that the function g μ (or at least a multiple of g μ ) is the analytic continuation of the function f μ . This last fact can be verified numerically in the following way. We consider the solution f μ and express it in the basis [ g μ , g μ log +ψμ ], using the function connect with the connection point at 0. So we obtain two (real) coefficients aμ and bμ such that f μ (t ) = aμ g μ (t ) + bμ ( g μ (t ) log(t + 1) + ψμ (t )). Now we draw the coefficients aμ (in grey) and bμ (in black) in function of μ for μ ∈ ]0, 14] (see Fig. 3).

F. Fauvet et al. / Applied Numerical Mathematics 60 (2010) 1309–1319

1317

Again, we see that the coefficient bμ vanishes for the first four eigenvalues (the numerical values of these eigenvalues, found in [1], are μ0 = 0.319, μ1 = 2.593084, μ2 = 6.533471, μ3 = 12.514462). Moreover, we observe that the corresponding coefficient aμ is ±1 (alternatively). 5. The prolate spheroidal equation and differential Galois theory We will now translate our main result in terms of differential Galois theory. For the necessary definitions and results we refer the reader to the articles of M. Singer and J. Martinet–J.P. Ramis in [24] or, for more details, to [25]. Let be μ ∈ C. The local differential Galois group G ∞ (μ) of the differential equation D μ y = 0 at infinity (that is relatively to the field C({1/t })) is a complex linear algebraic group topologically (in Zariski sense) generated by the formal monodromy (at infinity), the exponential torus and the Stokes multipliers. Here the formal monodromy is trivial, therefore G ∞ (μ) is Zariski connected. The exponential torus is isomorphic to the multiplicative group C∗ . Let G (μ) be the global differential Galois group (that is relatively to the field C(t )) of the differential equation D μ y = 0, the Wronskian differential equation is (t 2 − 1) w  + 2t w = 0, it admits a rational solution t 2 1−1 , therefore G (μ) is isomorphic to a subgroup of SL2 (C). The group G (μ) is connected: the local Galois groups at −1, 1, ∞ are connected and the global group is topologically generated by the local groups. Theorem 2. Let be μ ∈ C, the five equivalent conditions of Theorem 1 are equivalent to the following: 6. The local differential Galois group at infinity G ∞ (μ) of the differential equation D μ y = 0 is solvable. If these conditions are satisfied, then G ∞ (μ) is isomorphic to the multiplicative group C∗ , therefore it is abelian. If not, G ∞ (μ) is isomorphic to SL2 (C). Proof. The change of variables t → −t changes yˆ 1 into − yˆ 2 , yˆ 2 into − yˆ 1 , it exchanges the directions π2 and − π2 , and the summations y + and y − (i = 1, 2), therefore we have, on the sector ]π , 2π [: i i









π

+ − − 2 y+ 2 , y1 = y2 , y1 S ,

as we have also, on the same sector:









π

+ − − −2 y+ , 1 , y2 = y1 , y2 S

we get α = β . If the five conditions of Theorem 1 are satisfied, then G ∞ (μ) is isomorphic to the exponential torus, that is to C∗ , therefore it is abelian and a fortiori solvable. Conversely, if G ∞ (μ) is solvable, then it is possible to put it in triangular matrix form. As it contains the Stokes multipliers, one of them is the identity, but then the other is also the identity, therefore condition (4) is satisfied. The group G ∞ (μ) is connected and isomorphic to a subgroup of SL2 (C), if G ∞ (μ) is not solvable, then it is necessarily isomorphic to SL2 (C) [25]. 2 This result is an answer to a question asked by Alain Connes to J.P. Ramis about a possible interpretation of the eigenvalues of prolate spheroidal equations in differential Galois terms. It generalizes what happens for the harmonic oscillator and the hydrogen atom. In these cases eigenvalues correspond to the existence of orthogonal polynomials (Hermite and Laguerre polynomials), therefore to the vanishing of one of the Stokes multipliers, and this is equivalent to the solvability of the connected component of the differential groups. (For a study of some relations between spectra of Schrödinger equations and differential Galois theory, cf. [2].) However, if in the case of Laguerre and Hermite it is possible to find the eigenvalues using the global Galois group, it is not the case for the prolate spheroidal equation: for all μ ∈ C, the global differential Galois group is isomorphic to SL2 (C). We already know this result when μ is not an eigenvalue (G ∞ (μ) ⊂ G (μ)), it remains to prove it in the eigenvalues cases. The group G (μ) is connected and it is isomorphic to a subgroup of SL2 (C), therefore G (μ) = SL2 (C) if and only if G (μ) is solvable, or equivalently if and only if there exists a solution y of D μ y = 0 such that

y y

∈ C(t ). f

In an eigenvalue case, such an y must be collinear to f μ , we will prove that the meromorphic function f μ is not a μ rational function. t t In the neighborhood of infinity in C, we have f μ (t ) = cos (a + O( t12 )) + sin (b + O( t12 )) and f μ (t ) = − sint t (a + O( t12 )) + t t2 cos t (b t2

− a + O( t12 )) or f μ (t ) = to the parity (with a = 0).

sin t (a t

+ O( t12 )) +

 (t ) fμ sin t +O(1/t ) = − cos t +O(1/t ) f μ (t )  (t ) fμ sign of f (t ) is the sign μ

Then, for t ∈ R \ Z π2 , we have and |t | sufficiently large, the

cos t (b t2

zeroes or poles and it is not a rational function.

+ O( t12 )) and f μ (t ) =

or

 (t ) fμ f μ (t )

=

cos t +O(1/t ) , sin t +O(1/t )

cos t (a t

+ O( t12 )) −

sin t (a t2

+ b + O( t12 )), according

according to the parity. Therefore, for t ∈ π4 + Zπ

of − tan t = ±1 (or cotan t = ±1). Hence

 fμ fμ

has an infinite number of

1318

F. Fauvet et al. / Applied Numerical Mathematics 60 (2010) 1309–1319

6. Conclusion There are some natural questions in relation with the present article, we plan to go back to the subject in a future work. – Fourier transform d By a Fourier transform ( dt → iξ , t → i ddξ , up to a choice of normalization), the prolate differential operator L σ ,τ ∈ C[t , dtd ]

(“in t”) is changed into the differential operator L τ ,σ ∈ C[ξ, ddξ ] (“in ξ ”), cf. [26]. This is related to the bispectrality of prolate spheroidal differential equation. Traditionally this is studied in the real domain, but it is interesting to understand what is happening in the complex domain, using the results of the present article and Malgrange dictionary (cf. Chapter XII of [13]). The prolate differential operators are of exponential type in Malgrange sense (their Newton polygons admit only a nonzero slope and this slope is 1). For such operators there is a good dictionary exchanging “solutions at finite distance + monodromy” and “solutions at infinity + Stokes phenomena”. The operator L τ ,τ , is “invariant by Fourier transform” (cf. [6], 3.3, page 356), then it is possible to give a Fourier transform interpretation of the results of the present article. One important point is that, in the case of L 1,1 , f μ (resp. g μ ), interpreted as a power series in the variable (t − 1) (resp. (t + 1)) is the Borel-transform of −ie it yˆ 1 (resp. −ie −it yˆ 2 ) interpreted as a power series in the variable − ti . – Eigenvalues as zeroes of entire functions

We consider the Stokes matrix

1 0



α (μ) 1

associated to L 1,1 y = μ y. The function

its zeroes are the eigenvalues. It would be interesting to know the growth of (one can compare with the case of the quartic oscillator, cf. [7]).

μ ∈ C → α (μ) is an entire function and α at infinity and if it is a resurgent function

– Numerical study of the convergent series at ∞ (eigenvalue case), in comparison with Hermite case In the case of the Hermite differential equation: y  (x) + ay (x) − xy  (x) = 0, one series solution turns into a polynomial for values of a in N∗ . This situation has been studied by J.L. Callot in [5] with techniques of non-standard analysis, and it is shown that so called “duck solutions” appear. It will be interesting to do an analog of this study for the prolate wave equation. – Degenerate cases The “degenerate cases” τ = 1, σ = 0 and τ = 0, σ = 1 of prolate differential operators L σ ,τ correspond respectively to Legendre differential equations (t 2 − 1) y  + 2t y  = n(n + 1) y and Bessel differential equations t 2 y  + 2t y  + t 2 y = n2 y. It would be interesting to understand what happens to the structures studied above in the corresponding limits ((1, σ ) → (1, 0) and (τ , 1) → (0, 1)), in particular in relation with the study of the Helmoltz P.D.E. when the spheroid tends to a sphere. One can also study the Hermite–Weber approximation of spheroidal functions (L τ ,τ , τ → +∞). Acknowledgements J.P. Ramis thanks A. Connes who introduced him to prolate spheroidal functions and asked the question which motivated this work. This work was done with the support of the ANR Galois. References [1] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, Dover Publications, New York, 1964. [2] P.B. Acosta Humánez, Galoisian approach to supersymmetric quantum mechanics, PhD Dissertation, Departament de Matematica Aplicada I, Universitat Politécnica de Catalunya, April 20, 2009. [3] W. Balser, Formal Power Series and Linear Systems of Meromorphic Ordinary Differential Equations, Springer-Verlag, 1999, 299 pp. [4] M. Barkatou, Rational Newton algorithm for computing formal solutions of linear differential equations, in: Proceedings of ISSAC’88, Rome, Italy, ACM Press, 1988, pp. 183–195. [5] J.L. Callot, Bifurcations du portrait de phase pour des équations différentielles linéaires du second ordre ayant pour type l’équation d’Hermite, Thèse, IRMA Strasbourg, 24 juin 1981. [6] A. Connes, M. Marcolli, Noncommutative Geometry, Quantum Fields and Motives, A.M.S. Colloquium Publications, vol. 55, A.M.S., 2007. [7] E. Delabaere, F. Pham, Unfolding the quartic oscillator, Ann. Phys. 261 (1997) 180–218. [8] J. Della Dora, C. Di Crescenzo, E. Tournier, An algorithm to obtain formal solutions of a linear homogeneous differential equation at an irregular singular point, in: J. Calmet (Ed.), EUROSAM 82, in: Lecture Notes in Computer Science, vol. 144, Springer-Verlag, Berlin/Heidelberg, 1982, p. 273. [9] P.E. Falloon, Theory and computation of spheroidal harmonics with general arguments, Thesis, University of Western Australia, 2001. [10] F. Fauvet, J. Thomann, Formal and numerical calculations with resurgent functions, Numer. Algorithms 40 (4) (December 2005), online at http://dx.doi.org/10.1007/s11075-005-5326-5. [11] F. Fauvet, F. Richard-Jung, J. Thomann, Algorithms for the splitting of formal series; applications to alien differential calculus, in: The Proceedings of Transgressive Computing 2006, a Conference in Honor of Jean Della Dora, Granada, Spain, April 2006.

F. Fauvet et al. / Applied Numerical Mathematics 60 (2010) 1309–1319

1319

[12] F. Fauvet, F. Richard-Jung, J. Thomann, Automatic computation of Stokes matrices, Numer. Algorithms 50 (2) (February 2009), online at http://dx.doi.org/10.1007/s11075-008-9223-6. [13] B. Malgrange, Équations différentielles à coefficients polynomiaux, Progress in Mathematics, Birkhäuser, 1991. [14] J. Martinet, J.P. Ramis, Elementary acceleration and multisummability. I, Ann. Inst. H. Poincaré Phys. Théor. 54 (4) (1991) 331–401. [15] J. Meixner, F.W. Schafke, Mathieusche Funktionen und Sphäroidfunktionen, Springer-Verlag, 1954. [16] E. Pflügel, On the latest version of DESIR-II, Theor. Comput. Sci. 187 (1–2) (1997) 81–86. [17] J.P. Ramis, Séries divergentes et théories asymptotiques (Divergent series and asymptotic theories), Bull. Soc. Math. France 121 (1993), Panoramas et Synthèses, suppl., 74 pp. (in French). [18] F. Richard-Jung, Représentation graphique de solutions d’équations différentielles dans le champ complexe, Thèse de Doctorat, Université L. Pasteur, septembre 1988. [19] F. Richard-Jung, An implicit differential equation with Maple and IRONDEL, in: The Proceedings of Computer Algebra in Scientific Computing 2004, Saint Petersburg, Russia, July 2004, pp. 383–397. [20] F. Richard-Jung, The DESIR package, software demonstration, in: ISSAC’09, Seoul, July 2009, see also http://www-ljk.imag.fr/CASYS/LOGICIELS/ desir2009.html. [21] D. Slepian, Some comments on Fourier analysis, uncertainty and modeling, SIAM Rev. 25 (1983) 379–393. [22] J. Thomann, Resommation des séries formelles solutions d’équations différentielles linéaires ordinaires du second ordre dans le champ complexe au voisinage de singularités irrégulières, Numer. Math. 58 (1990) 503–535. [23] J. Thomann, Procédés formels et numériques de sommation de séries solutions d’équations différentielles, Expo. Math. 13 (1995) 223–246. [24] E. Tournier (Ed.), Computer Algebra and Differential Equations, Lecture Note Series, vol. 193, Cambridge University Press, 1994. [25] M. van der Put, M. Singer, Galois Theory of Linear Differential Equations, Springer-Verlag, Berlin, 2003. [26] G. Walter, T. Soleski, A friendly method of computing prolate spheroidal wave functions and wavelets, Comp. Harmonic Anal. 19 (2005) 432–443.