Computation of a numerically satisfactory pair of solutions of the differential equation for conical functions of non-negative integer orders

We consider the problem of computing satisfactory pairs of solutions of the differential equation for Legendre functions of non-negative integer order $\mu$ and degree $-\frac12+i\tau$, where $\tau$ is a non-negative real parameter. Solutions of this equation are the conical functions ${\rm{P}}^{\mu}_{-\frac12+i\tau}(x)$ and ${Q}^{\mu}_{-\frac12+i\tau}(x)$, $x>-1$. An algorithm for computing a numerically satisfactory pair of solutions is already available when $-1<x<1$ (see \cite{gil:2009:con}, \cite{gil:2012:cpc}).In this paper, we present a stable computational scheme for a real valued numerically satisfactory companion of the function ${\rm{P}}^{\mu}_{-\frac12+i\tau}(x)$ for $x>1$, the function $\Re\left\{e^{-i\pi \mu} {{Q}}^{\mu}_{-\frac{1}{2}+i\tau}(x) \right\}$. The proposed algorithm allows the computation of the function on a large parameter domain without requiring the use of extended precision arithmetic.

− 1 2 +iτ (x), x > −1. An algorithm for computing a numerically satisfactory pair of solutions is already available when −1 < x < 1 (see [6], [5]). In this paper, we present a stable computational scheme for a real valued numerically satisfactory companion of the function P µ − 1 2 +iτ (x) for x > 1, the function ℜ e −iπµ Q µ − 1 2 +iτ (x) . The proposed algorithm allows the computation of the function on a large parameter domain without requiring the use of extended precision arithmetic.

Introduction
Conical or Mehler functions appear in a large number of applications in applied physics, particle physics (related to the amplitude for Yukawa potential scattering) or cosmology, among others. These functions also appear when solving the Laplace's problem in spherical coordinates for two intersecting cones [9].
The conical functions are solutions of the second order differential equation In applications, one encounters values of the conical functions where the order µ is integer positive (or zero), so we will consider µ = m ∈ Z + .

Numerically satisfactory pairs of solutions
In the interval −1 < x < 1, a numerically satisfactory pair of solutions (real valued) of the equation (1.1) is P −µ − 1 2 +iτ (x) and P −µ − 1 2 +iτ (−x), µ ∈ ℜ, τ ∈ ℜ + [3]. The Wronskian for this pair of solutions is: . (2.1) When µ = m ∈ Z, it should be noted that the functions P m The definition of the function P m − 1 2 +iτ (x) in terms of the Gauss hypergeometric function 2 F 1 is given in Eq. (3.3). In [6] and [5] we have described an algorithm for computing the conical functions P m − 1 2 +iτ (x) for x > −1; therefore the problem of computing a numerically satisfactory pair of solutions of the differential equation for conical functions in the interval −1 < x < 1 can be considered as already solved. The algorithm is based on the use of x 10 8 different methods of computation, depending on the range of the parameters: quadrature methods, recurrence relations and uniform asymptotic expansions in terms of elementary functions or in terms of modified Bessel function K ia (x) and its derivative K ′ ia (x). When x > 1, a real-valued satisfactory pair of solutions is P m . As an example, Fig. 1 shows the graphs of P 10 − 1 2 +i5 (x) and Q 10 − 1 2 +i5 (x). As in the previous case, an algorithm for computing the function P m − 1 2 +iτ (x) for x > 1 is already given in [6] and [5]. However, no computational schemes can be found in the literature for computing the numerically satisfactory companion solution of P m − 1 2 +iτ (x), the function Q m − 1 2 +iτ (x). Our purpose is to give an algorithm for computing this function.

Recurrence relations
Three-term recurrence relations y n+1 + b n y n + a n y n−1 = 0, (3.1) are useful methods of computation when two starting values are available for starting the recursive process. Usually, the direction of application of the recursion can not be chosen arbitrarily, and the conditioning of the computation of a given solution fixes the direction. The conical function P m −1/2+iτ (x) satisfies a three-term recurrence relation, for x > 1, where we have adopted the following definition of the function P m This is the definition used in [5] for x > −1. Note the difference of the sign in the second term of the recurrence relation (3.2) in comparison with Eq. 14.10.6 in [3], which reflects the fact that the definition of our function P m − 1 2 +iτ (x) and the function P m −1/2+ıτ (x) ( [2], [3]) differ by a factor (−1) m : The conical functions P m − 1 2 +iτ (x) are monotonic in the interval (1, x c ) and oscillating in (x c , +∞), where x c = 1 + β 2 /β and β = τ /m. As an example, for the parameters of the functions plotted in Fig. 1, the transition point is found at x c ≈ 2.24. The stability analysis based on Perron's theorem discussed in [6], revealed that the function P m − 1 2 +iτ (x) was the minimal solution of the recurrence relation (3.2). Hence, for this function backward recursion will be generally stable for x > 1. However, and similar for other special functions, recurrence relations in the oscillatory regime of the conical functions (x > x c > 1) are not ill conditioned in either backward and forward directions; thus, both recursions are possible.
It is very simple to show that the function Q m − 1 2 +iτ (x) also satisfies the three-term recurrence relation given in (3.2), and this function is a dominant solution of (3.2) (the minimal solution, if it exists, is unique). In this case, the stable direction of application of the recurrence relation is with increasing m, although the same comment as for the function P m We next consider the problem of computing two starting values ( Q 0

Power series expansions by using the hypergeometric functions
From the many representations of the conical functions in terms of the Gauss hypergeometric functions we choose one representation that can be used for small values of x − 1 and one for large values of x (and that one can also be used for large values of τ ).

Expansions valid near the point x = 1 and moderate τ
We use the representation (see [1,Eqn. (32), page 131]) We use power series in z, and because we want these expansion for µ = 0 and µ = 1, we need to find the limit for the representation in (4.1). First we write The limit for µ → 0 gives where ψ(α) = Γ ′ (α)/Γ(α). We obtain the expansion where w and z are given in (4.2).
which is real when τ is real, and the computation easily follows from recursion. The quantity ψ(ν + 1) is the only complex term to consider in more detail. For ψ(k + 1) we can use the recursion ψ(k + 1) = ψ(k) For the computation of ψ(α) (with α complex) we use the asymptotic expansion valid for α → ∞ in |ph α| < π. We use this expansion if |α| ≥ 12 with 8 terms of the series (or less), and we use the backward recurrence relation ψ(α) = ψ(α + 1) − 1/α for smaller values of |α|. We only need an algorithm for ℜα ≥ 1 2 . For Q 1 ν (x) we can also use a limiting procedure, but it is more convenient to use the derivative of Q 0 ν (x), because where This gives the expansion (4.12) When computing the expansion for τ large, it is convenient to use the following asymptotic expansion for the ratio of two gamma functions [8]: (4.14) In the present case, we have and the first coefficients are The representation in (4.12) becomes and We can obtain these quantities from the recurrence relations with u 0 (τ ) = 1, v 0 (τ ) = 0, w 0 (τ ) = 1. We can write (4.12) with a single complex representation, by writing which gives where (4.23)

Expansions in terms of the Kummer U −functions
First we explain how the method works for the Gauss hypergeometric function. We take the integral representation and using the transformation u = 1 − e −t , we write it in the form We assume that ω is large and that z may be large as well. In that case α will be small. The easiest approach is to expand The functions Φ n can be expressed in terms of the confluent hypergeometric function U (a, c, z). We have Φ n = (b) n α n+b−a U (n + b, n + b + 1 − a, αω) = (b) n ω a−n−b U (a, a + 1 − n − b, αω).

(5.7)
The functions Φ n can be obtained by using recurrence relations for the U −function 1 , or from integrating by parts in (5.6). We have For numerical aspects of such recursions, we refer to [7] and [4, §4.5.1]. The expansion in (5.5) has an asymptotic character for large values of ω, uniformly with respect to α. The second line of (5.7) gives the integral representation uniformly with respect to α ≥ 0, when n is such that convergence at infinity of the integral is guaranteed if α = 0, that is, if n > ℜ(a − b).

Expansion of
For the representation of the conical Q−function in (4.10) we have a = 1 2 +µ, b = 1 2 − µ, c = 1, ω = iτ . This gives where α = ln z+1 z , z is defined in (4.11), and (5.12) 14) The first f k coefficients 2 are given by The Φ k can be written in terms of the Hankel functions. The first one is The representation in terms of of the Hankel functions is convenient, because it is easy to separate real and imaginary parts by using

Numerical tests of the expansions
In order to test the range of validity of the expansions, we have first compared the results with the real part of the value obtained in the computation of Q 0 − 1 2 +iτ (x) using the command LegendreQ in Maple with 50 digits. Fig. 2 shows, as a function of x and for three different values of τ , the accuracy obtained in the computation of Q 0 − 1 2 +iτ (x) using the power series expansion given in (4.6). We have used the expansion with 0 ≤ k ≤ 8. As can be seen, the expansion provides full accuracy for values of x very close to 1 but the accuracy worsens as τ increases. As an example, the number of terms of the expansion has to be increased to more than 100 in order to obtain an accuracy better 10 −14 when computing Q 0 − 1 2 +i200 (1.05). The same test has been performed using the power series of (4.12) for the computation of Q 0 − 1 2 +iτ (x). The results are shown in Fig. 3. In this case, as expected the expansion behaves better as τ increases but the accuracy worsen as x approaches to 1.
A test for the expansion in terms of the Kummer U -function 3 is shown in Fig. 4: the relative errors obtained in the computation of Q 0  with starting values those given in (5.16) and (5.17). As can be seen, the expansion provides an accuracy better than 10 −8 (single precision) even for values of the parameter τ lower than 5. Two values of the argument x are chosen (x = 1.1, 100) in order to illustrate the validity of the expansion for large and small (larger than 1) values of x.
A test independent of the Maple procedure LegendreQ has also been considered. It consists in testing that the values computed using the expansion are consistent with the three-term recurrence relation (3.2) for Q m − 1 2 +iτ (x). We first compute Q 0 The deviations from 1 of the left-hand side of (6.1) for m = 1 are shown in Fig. 5. The results are consistent with those obtained when testing the computed value of Q 0 − 1 2 +iτ (x). Another test involving the application of the three-term recurrence relation can be seen in Fig. 6. In this figure we show as a function of m, the relative error obtained in the computation of Q m − 1 2 +iτ (x) by applying the Figure 5: Test of the recurrence relation given in (6.1) for m = 1 and x = 1.1, 100 using the asymptotic expansion (5.13) with n = 8.

Computational scheme
From the results obtained in the previous section, a stable computational scheme for evaluating the function Q m 2. For x close to 1 and moderate/large values of τ , compute Q m − 1 2 +iτ (x) using the expansion for large τ given in Sect. 5.1. Figure 6: Test of the recurrence relation given in (6.1) for τ fixed (τ = 50) and x = 1.00001, 10, 100, 500 using the asymptotic expansion (5.13) with n = 8. For higher accuracies, the regions of application of the methods may need some adjustments.