Higher Derivatives of Airy Functions and of their Products

The problem of evaluation of higher derivatives of Airy functions in a closed form is investigated. General expressions for the polynomials which have arisen in explicit formulae for these derivatives are given in terms of particular values of Gegenbauer polynomials. Similar problem for products of Airy functions is solved in terms of terminating hypergeometric series.


Introduction
Explicit formulae for higher derivatives is a usual part of investigation of special functions in mathematical physics [10,11,12,22,24,25]. A collection of these results can be found in [6,18].
For any solution y(x) of a homogeneous differential equation of second order y + p(x)y + q(x)y = 0 one can obtain that y = −qy − py , y = (pq − q )y + p 2 − p − q y , y IV = −p 2 q + pq + 2p q + q 2 − q y + −p 3 + 3pp + 2pq − p − 2q y , etc. However, it would be quite difficult to get the explicit formula in general: y (n) = P n (x; p, q)y + Q n (x; p, q)y , because successful finding of coefficients P n (x; p, q) and Q n (x; p, q) in a closed form depends on many circumstances. The main purpose of this paper is to obtain general formulae for n-th derivatives of Airy functions, Ai(x) and Bi(x). Both functions satisfy the same equation y = xy.
(1.1) Therefore, Ai (n) (x) = P n (x)Ai(x) + Q n (x)Ai (x), Table 1. First 16 polynomials P n (x) and Q n (x). n P n (x) Q n (x) 0 1 0 x 3 + 4 6x 7 9x 2 x 3 + 10 8 x 4 + 28x 12x 2 9 16x 3 + 28 x 4 + 52x 10 x 5 + 100x 2 20x 3 + 80 11 25x 4 + 280x x 5 + 160x 2 12 x 6 + 260x 3 + 280 30x 4 + 600x 13 36x 5 + 1380x 2 x 6 + 380x 3 + 880 14 x 7 + 560x 4 + 3640x 42x 5 + 2520x 2 15 49x 6 + 4760x 3 + 3640 x 7 + 770x 4 + 8680x where P n (x) and Q n (x) are some polynomials and the index n corresponds to the derivative order but not the polynomials degree. By differentiating the first equation of (1.2), Ai (n+1) (x) = P n (x) + xQ n (x) Ai(x) + P n (x) + Q n (x) Ai (x) = P n+1 (x)Ai(x) + Q n+1 (x)Ai (x), we have two differential difference relations P n+1 (x) = P n (x) + xQ n (x), Q n+1 (x) = P n (x) + Q n (x) (1.3) with initial conditions P 0 (x) = 1 and Q 0 (x) = 0 which help to determine P n (x) and Q n (x) for any natural n (see the Table 1). It is quite possible that the problem of evaluation of polynomials P n (x) and Q n (x) was formulated in XIX century when G.B. Airy introduced the function which is now denoted Ai(x) (see [1,2]). However, as far as we know, the first solution of the problem has been published by Maurone and Phares in 1979 in terms of double finite sums containing factorials and gamma function ratio [17]. We rewrite it using binomial coefficients and Pochhammer symbols , where x is the integral part of x, δ = 0, 1, 2, and k 1 = 0, 1 = 0, m 1 = 0, k 2 = 1, 2 = 1, m 2 = 0, if δ = 0, k 1 = 1, 1 = 2, m 1 = 1, k 2 = 0, 2 = 0, m 2 = 0, if δ = 1, Recently, the problem of evaluation of P n (x) and Q n (x) has been investigated by Laurenzi in [15] where the solution has been written in terms of Bell polynomials. More general problem of evaluation of higher derivatives of Bessel and Macdonald functions of arbitrary order has been solved by Brychkov in [7]. Corresponding polynomials have been written as finite sums containing products of terminated hypergeometric series 3 F 2 and 2 F 3 but they look quite cumbersome for the case of Ai(x), namely, And third, the Wronskian of f (x) and g(x) is equal to 1, that helps to simplify polynomials P n (x) and Q n (x) as solutions of the system (2.1) Substitution of the series expansions of Airy atoms and their derivatives into (2.2) yields where both series on k are naturally terminating, i.e., the index of summation runs over all values that produce non-zero summands.
Let us define γ m by the formula Then (2.3) can be rewritten as follows Reducing the product of Pochhammer symbols to Euler's beta function and applying one of the function integral definitions Here, ω = e 2πi/3 andω = e −2πi/3 are the cubic roots of unity, and an overline means complex conjugation.
Changing the variable in the last integral, t → 1/t, we get another expression for γ m : Then the differences in square brackets in (2.4) can be written as Expansions (2.4) show that exponents of the factor (1 + ωt) in both integrals are nonnegative.
The integral for γ m (n − 1, 0) − γ m (n − 1, n) looks a little simpler than for γ m (n, n) − γ m (n, 1), so we begin with it. Integrating the analytic function (1+z) 3m+1−n /(1+z 3 ) m+1 over the contour Then the integral of our interest can be written as It is worth noting that the first term on the right side makes a zero contribution to the difference Using the generating function, we prove that the right part of (2.7) vanishes for n = 0, 1, . . . , 2m. Let s ∈ R. Then Changing the variable ζ = (z + a)/(1 + bz), where parameters a and b will be chosen later, one gets With the aim to simplify the integrand, we set b = −i and a = e πi/6 + s. Then Therefore, Re F (s) = 0, and the expansion of Q n (x) in (2.4) can be simplified by throwing out the terms with n < 2m + 1. As result, the function is actually a polynomial.
The double inequality n−1 is quite restrictive for the summation index. In particular, it leads immediately to zero polynomials when n is equal to 0 or 2.
Now we need to evaluate the derivative in (2.8). As before, we use the generating function approach. Let m, n 0 and s ∈ R. Then assuming that is a small positive number and |s| < min |z|= |z(z+i)| = (1− ) for the geometric series convergence. Roots of the equation z 2 + iz − s = 0 are z ± = −i(1 ± √ 1 − 4s )/2, and z − is the only pole of the integrand within the contour |z| = . Therefore, Here, we use the notation proposed in [14]. Namely, if A(z) is any power series k a k z k , then [[z k ]]A(z) denotes the coefficient of z k in A(z). In our view, this notation is more convenient to manipulate power series than usual analytic description, Returning to (2.8) and renaming the coefficients for brevity, we havẽ n + 2m + 1 n · 2 F 1 n + 2m + 3 2 , n + 2m + 2 2 ; 2m + 3 2 − 1 3 .

Using Euler's transformation
and one of hypergeometric representations of Gegenbauer's polynomials [21, equation (7.3.1.202)], we can simplify the coefficients tõ (2.11) and obtain required formulae for both polynomial families, first for Q n (x), then for P n (x) = Q n+1 (x) − Q n (x). We write them in the same style, shifting the index in Q n (x) The last expression in (2.11) follows from the familiar generating function for Gegenbauer poly- Of course, all the coefficientsg m,n−2m andg m,n−2m −g m,n−2m−1 are positive but this is easier to deduce from (1.3) than from (2.11). It is worth noting, that the termg m,n−2m−1 in (2.13) vanishes for the case n = 2m since (2.10) holds. Left and right inequalities on the summation index in (2.12) and (2.13) reveal the different structures of analytic expressions of the polynomials for small and large values of x. Namely, for x ∼ 0 the expressions depend on n (mod 3): 14) while for x → ∞ they depend on n (mod 2): In addition, it may be worth noting that the methods of finding these expansions are essentially different: (2.15) follows directly from (2.12) and (2.13), while for proving (2.14) the easiest way is to apply (2.3).

Hypergeometric series and difference equations
The initial statement of the problem (1.2) and our final expressions for polynomials P n (x) and Q n (x) lead to some corollaries whose relation to the Airy functions does not look evident a priori. Below we consider two of them. The first helps to find the Gauss hypergeometric function special values. The second is connected with difference equations of third and fourth orders.

Some special values of the hypergeometric function of Gauss
Power series expansions of polynomials P n (x) and Q n (x) help to find values of the Gauss hypergeometric function in some special cases. The simplest of them is obtained equating the terms Q 3n+1 (0) in (2.12) and (2.14): This formula is not new (see Exercise 15 in [25,Chapter 14]) and can be extended to real or complex a's by replacing n by −2a [21, equation (7.3.9.15)]: As is typical in the theory of hypergeometric functions, the extension from n to a is valid due to Carlson's theorem [23, Section 5.8.1]: If f (z) is regular and of the form O(e c|z| ), where c < π, for Re z 0, and f (z) = 0 for z = 0, 1, 2, . . ., then f (z) = 0 identically. (See [4,5,8] for many examples of application of Carlson's theorem to hypergeometric series.) Next two hypergeometric function special values are related to coefficients of Q 3n+3 (x) and Q 3n+2 (x): By replacing n by −2a in the first relation and by 1 − 2a in the second, we get In order to use the coefficients of P n (x), we need to do some preliminary work. Substitution of (2.9) in both terms of the differenceg m,n −g m,n−1 and the binomial identity Then, by considering the first nonzero coefficients of P 3n (x) and P 3n+2 (x), P 3n (0) = 3 n 1 3 n = n!{g n,n −g n,n−1 }, we get the following special values of the Gauss hypergeometric function: It is quite clear that this approach provides a way for evaluation of in a closed form for any n ∈ Z, while, of course, the expressions will become more and more cumbersome as |n| increases.

Difference equations of third and fourth orders
Relations (2.2) help to find generating functions for polynomials P n (x) and Q n (x): Since f (x) and g(x) satisfy the equation (1.1), then P(x, t) and Q(x, t) being the functions on t are solutions of the differential equation If a solution of this equation is expanded in a power series then the series substitution into the equation leads to a recurrence relation for its coefficients, functions Y n (x): As result, polynomials P n (x) and Q n (x) are two of three solutions of (3.2), corresponding to different initial conditions: It is not evident how to find the third independent solution of (3.2), say and how this relates with Airy functions' derivatives. Investigation of the polynomials Z n (x) shows that where coefficients λ m,n satisfy the relation mλ m,n = (3m − n)λ m−1,n−2 + nλ m−1,n−3 , m 1.
We mention also two difference equations for the coefficients of polynomials P n (x) and Q n (x). These equations can be found considering the Laplace transform of the shifted Airy function, We will use a formal power series approach for simplicity, while there is, of course, a more rigorous but tedious way based on an asymptotic behaviour of Ai(t) for large t's. Assuming that p → +∞, the asymptotic expansion of L (p, x) can be written in two ways. The first follows from (1.2) On the other hand, the function L (p, x) satisfies a differential equation of first order d dp With initial condition where the last notation is due to Aspnes [3], the solution is Integration by parts yields Then, in general, L (p, x) has the form with coefficients µ k , ν k ,μ k ,ν k depending on x. Substituting (3.6) in (3.5), we obtain the same recurrence relations for the pairs (µ k , ν k ) and (μ k ,ν k ) but different initial conditions (µ 0 , ν 0 ) = (1, 0), (µ 1 , ν 1 ) = (0, 1), (μ 0 ,ν 0 ) = (0, 1), (μ 1 ,ν 1 ) = (0, 0). Now we need to reduce (3.6) to the form of (3.4). Transforming the first series of (3.6), we have For the second series of (3.6), transformation is the same. As result, we obtain expansions of P n (x) and Q n (x) depending on parity of n where we add an explicit dependence of the coefficients on x. All the coefficients satisfy difference equations of fourth order based on (3.7): µ k ,μ k : y k+4 = (2k + 3)(2k + 6)y k+1 + (2k + 2)(2k + 6)xy k , (3.8) ν k ,ν k : y k+4 = (2k + 4)(2k + 7)y k+1 + (2k + 2)(2k + 6)xy k . Unfortunately, two other pairs of independent solutions of (3.8) and (3.9) are yet remain unknown, while a relation of one of the pairs (say,μ k andν k ) to polynomials (3.3) seems quite predictable.

Higher derivatives of Ai 2 (x), Ai(x)Bi(x) and Bi 2 (x)
In a similar fashion to the evaluation of higher derivatives of Airy functions considered in Section 2, let us investigate the same problem for the products of these functions. It is quite evident that the problem is reduced to a determination of polynomials R n (x), S n (x) and T n (x) satisfying the following three equations where R 0 (x) = 1, S 0 (x) = 0 and T 0 (x) = 0. By taking the derivative of any equation above, we get the system of differential difference relations which help to determine R n (x), S n (x) and T n (x) for any n (see the Table 2). As before, it is better to use Airy atoms for finding the polynomials than Airy functions. By rewriting the system (4.1) as follows and applying the Wronskian W [f, g] to the matrix inversion, it is easy to obtain the solution Since the hypergeometric description of f 2 (x), f (x)g(x) and g 2 (x) is known and simple [6, equation (8.3.2.38)] the initial problem is reduced to manipulations of power series. Combining all three formulae of (4.4) together, (here and below, we use the symbol ' ' for separation of subparts) we write n-th derivatives of f 2 (x), f (x)g(x) and g 2 (x) in the same manner: and substitute all the series into (4.3).
We begin with the case of T n (x) and employ the same approach as in Section 2: To simplify the expression of T n (x), we shift n by 2: Applying Euler's beta integral we transform the inner series into an integral representation γ m (n, δ) = 2 π m! sin π 6 (1 + 2δ) As result, , and (4.6) is reduced to the following expression: In particular, T 0 (x) = T 1 (x) = T 3 (x) ≡ 0 because the set of possible m's values is empty. As before, we find h m,n using the generating function Then In spite of a terminating form of the hypergeometric series, the expression (4.8) does not look a perfect because the argument of the series is located in the exterior of the unit disk. Having in mind the results discussed in Section 3, we come back to the finite sum over k and reverse the order of summation that gives Its presence in (4.11) means that the second term in braces should be removed if n = 0. Now we briefly discuss two corollaries of the above formulae -a difference equation of third order and special values of the function 3 F 2 .
The case R 3n+2 (0) = 12 n+1 1 6 n+1 gives (4.18) Finally, it is worth noting that R n (x), S n (x) and T n (x) can be expressed in terms of P n (x) and Q n (x). Since we have that However, it seems difficult to obtain the expansions (4.9)-(4.11) by using formulae (2.12) and (2.13) directly.

Concluding remarks
In this paper, we have found higher derivatives of Airy functions and their products in a closed form. Explicit formulae for the polynomials which are contained in these derivatives help to obtain special values of hypergeometric functions 2 F 1 and 3 F 2 . We did not try to give a complete solution of this problem and considered only a few simple cases. In our view, a construction of the general solution will not require significant efforts. Moreover, combining this approach with the ideas presented in [13], one can find more general values of 3 F 2 -function depending on two parameters. For example (cf. with (4.16) and (4.17)), Nevertheless, there are some problems which are yet remain open and look quite difficult. For example, how to find the solution Z n (x) of the equation (3.2) and to reveal its relation to Airy functions.
The next problem is connected with zeros of all above polynomials. Let us define reduced polynomials removing a zero at the origin and an excessive cubicity: where δ = {0 1 2}. For example, Numerical investigations show that zeros of all reduced polynomials are real (therefore, negative) and simple but we could not prove this and find an asymptotical behaviour of zeros and extreme points of the polynomials.
In conclusion, it is interesting to note that negative counterparts of the coefficients (2.11) appear in some integral transforms of Airy functions. For example, where a > 0, b < 2 2/3 a, Re( 2 2/3 a − ωb) > 0, and is a periodic function with period 2 (see Fig. 1). A closer scrutiny of the curve reveals that with a hope thatτ (a) = τ (a). However, the numerical plotting procedure reveals that for all a. As result, one can obtain the formula F(a) = sin π a − 5 6 sin π 2a − 5 6 sin π a − 1 3 sin π a − 2 which coincides with equation (54) (see Fig. 2, where F(a) is shown for a ∈ [−1.5, 2.5]). Now, "when you have satisfied yourself that the theorem is true, you start proving it" [20,Chapter 5]. At this stage, it is convenient to replace the parameter a by a complex z: sin πz + π 6 .
We will prove that F(z) = F 0 (z) for all z in the region of joint analyticity by fitting F(z) − F 0 (z) to conditions of Carlson's theorem.
First of all, we cannot use the left half-plane, Re z 0, since F(z) is meromorphic there. Let us consider the right half-plane, Re z 0, where, as we will see later, F(z) is holomorphic. We need to define a sequence of positive numbers {z n } with z n → ∞ such that F(z n ) can be found easily in order to examine the validity of the relation F(z n ) = F 0 (z n ). The simplest way is to choose z n 's for which the hypergeometric series terminates: z = 1 6 and 3 2 − 3z = 0, −1, −2, . . .. Thus, one can obtain the following sequences: z n = n + 1 6 ,z n = n + 1 2 ,z n = n + 5 6 , n 0.
In conclusion, it is interesting to note that for proving identities (4.15)-(4.18) the numerical part of the method above leads to the same constant, −2, as in (A.1), while other components of final expressions are rather different.