On the complex zeros of Airy and Bessel functions and those of their derivatives

We study the distribution of zeros of general solutions of the Airy and Bessel equations in the complex plane. Our results characterize the patterns followed by the zeros for any solution, in such a way that if one zero is known it is possible to determine the location of the rest of zeros.


Introduction
Airy and Bessel functions are important examples of special functions satisfying second order lineal ODEs. Their zeros appear in a great number of applications in physics and engineering. Particularly, the complex zeros are quantities appearing in some problems of quantum physics [1], electromagnetism [2] and wave propagation and scattering [12,4].
Information is available regarding the distribution of zeros of some particular solutions of these differential equations. For example, the zeros of Airy functions Ai(z), Bi(z) and imAi(z) − Bi(z) (m real) are studied in [9] (see also [3]) together with the zeros of the Bessel functions J n (z) and Y n (z) of integer order (see also [11]) and the zeros of semi-integer order for combinations cos αJ ν (z)−sin αY ν (z) are discussed in [2]. The zeros of Hankel functions H (1) ν (z) and H (2) ν (z) for real order ν, which are also solutions of the Bessel equation, are analyzed in [1]. However, up to date there is no available description of the distribution of zeros for general solutions of the Airy and Bessel equations.
We perform this analysis by considering two complementary approaches. In the first place, a qualitative picture of the possible patterns of the complex zeros is obtained by assuming that the Liouville-Green approximation holds. This is the starting point of the numerical method of [14]. In the second place, the use of asymptotics (both of Poincaré type and uniform asymptotics for large orders in the case of Bessel functions) will provide more detailed and quantitative information. These results will characterize the possible patterns followed by the zeros for any solution, leading to the development of methods which are able to compute safely and accurately all the zeros in a given region and having as only input data the location of just one zero (and also the order for the case of the Bessel equation).

Liouville-Green approximation. Stokes and anti-Stokes lines.
An apparently naive method for studying the distribution of the complex zeros of second order ODEs y ′′ (z) + A(z)y(z) = 0 (1) consist in considering that A(z) is "locally constant" in the sense we next describe. We start by considering the trivial case of A(z) constant. Then the general solution of (1) reads y(z) = C sin A(z) (z − ψ) , and the zeros are over the line z = ψ + e −i ϕ 2 λ, λ ∈ R, ϕ = arg A(z).
In other words, writing z = u + iv we have that the zeros are over an integral line of dv du = − tan(ϕ/2). (2) and those curves are also given by (2). These are the so-called anti-Stokes lines (ASLs) in the LG approximation. Similarly, the curve defined by ℜ z z (0) gives the Stokes line (SL) passing through z 0 . It is a well-known fact that when SLs and ASLs intersect they do so perpendicularly. However, two ASLs (or SLs) can not intersect except at the zeros and poles of A(z). This fact singularizes the role of the ASLs (or SLs) emerging from the zeros (also called turning points) and the poles. The ASLs (or SLs) are called principal when they emerge from the zeros of A(z). In particular, if A(z) has a zero of multiplicity m at z = z 0 such that A(z) = a(z − z 0 ) m (1 + O(z − z 0 )) as z → z 0 , m ∈ N, then m + 2 principal ASLs emerge from z 0 in the directions arg(z − z 0 ) = 1 m + 2 (− arg(a) + 2kπ), k = 0, . . . , m + 1.
Principal ASLs (SLs) divide the complex plane in different disjoint domains such that any ASL (SL) is either inside one of these domains or is a principal ASL (SL) itself. In the literature these domains are called anti-Stokes and Stokes domains.
As discussed in [14], the qualitative picture provided by the ASLs gives generally a quite accurate picture of how the complex zeros may be distributed in the complex plane. For the zeros of y ′ (z), as we will see, the approximation that they are on ASL curves given by (2) (and equivalently 4) will be also in general sufficiently accurate. Furthermore, the numerical method for complex zeros of y(z) given in [14] is also successful for the zeros of y ′ (z) with only a simple change in the iteration function of the fixed point method.

Zeros of Airy functions
Next we analyze the distribution of the zeros of Airy functions. We begin with a description of the ASLs and SLs, followed by an analysis of the distribution of the zeros for the general solution cos αAi(z) + sin αBi(z); this analysis will be important in the description of the zeros of general Bessel functions cos αJ ν (z) − sin αY ν (z).

Anti-Stokes lines for the Airy equation in the LG approximation
For the Airy equation y ′′ (z) − A(z)y(z) = 0, with A(z) = −z, the ASL passing through a point z i is given in the LG approximation by This is the curve ℜ(z 3/2 − z 3/2 i ) = 0 (z = re iθ , z i = r i e iθi ). For z i = 0 we obtain the three principal ASLs emerging from the origin in the directions arg z = −π, ±π/3. The remaining ASLs (not having z i = 0 in the path) have the expression with θ running in three possible different intervals, and therefore with ASLs given by r(θ)e iθ in three different sectors In each sector, the curve (7) approaches asymptotically the boundaries of the sector as θ approaches the boundaries of the intervals, and it passes through the point z = r 0 e 2πij/3 . The curves (7) can be related to lines parallel to the positive real axis by a simple transformations. This relation will be interesting when we later obtain asymptotic approximations. We define where the factor 3/2 inside the parenthesis is introduced for later convenience. Then it is straightforward to verify that the values e iπ/3 z(ξ) with and δ a fixed value different from zero, represent the curve (7) with and with θ ∈ [0, π/3) if δ < 0 and θ ∈ (π/3, 2π/3] if δ > 0; similarly, e −iπ/3 z(ξ) gives a curve (7) with θ ∈ (−π/3, 0] if δ > 0 and θ ∈ [−2π/3, π/3] if δ < 0. Finally, e ±iπ z(ξ) gives the curves (7) for θ ∈ (π, 2π/3) and θ ∈ (−2π/3, −π).
Using this combination and letting α be in general complex we can study the distribution of zeros for any solution of the Airy equation. Indeed, given a solution y(z) = βAi(z) + γBi(z) the solution A(arctan(γ/β), z) is proportional to y(z) except when γ/β = ±i; in this last case the zeros of y(z) can be obtained from those of A(α, z) by letting ℑα → ∓∞.
According to the description of the ASLs, we expect that zeros can appear which asymptotically tend to one or several of the rays arg z = π, ±π/3. Later, when we analyze the zeros of Bessel functions, the zeros in the directions arg z = ±π/3 will be particularly important. We analyze the distribution of zeros in the three directions, starting with the negative real axis.

Zeros approaching asymptotically the negative real axis
For any real α, the function A(α, z) has an infinite number of negative real zeros, and for any finite non real value of α, there are zeros approaching the negative real axis, but not on the axis. To see this, we write Ai(−z) ∼ π −1/2 z −1/4 (cos ϕ R(ξ) + sin ϕ S(ξ)) , Bi(−z) ∼ π −1/2 z −1/4 (− sin ϕ R(ξ) + cos ϕ S(ξ)) , with where R(ξ) and S(ξ) have the Poincaré-type expansions [10, 9.7.9, 9.7.11] in a sector | arg z| ≤ 2π/3−δ, with u m numbers defined in [10, 9.7(i)]. The equality A(α, −z) = 0 gives Using the leading terms in (14) we have cot ϕ ∼ tan α and inverting this relation we have the estimation for the zeros as asymptotic approximation for large positive k (ℜξ < 0 does not make sense because (14) is valid for | arg ξ| ≤ π − ǫ). This corresponds to values of −z approaching the negative real axis from above if ℑα > 0 and from below of ℑα < 0 as k → +∞; the approximation (16) gives zeros which are on an ASL given by the LG approximation (see section 2.1).
Eq. (16) gives the dominant term in the asymptotic expansion for the negative real zeros of A(α, z), more terms can be obtained by inverting (15) using additional terms of the expansions (12). For later use, it will be interesting to invert a more general relation than (15), and we consider the asymptotic inversion of If f (α) = ±i and using (12) we obtain that The case f (α) = ±i would correspond to the zeros of Ai(−z) ± iBi(−z), but considering [10, 9.2.11] we have Ai(−z) ± iBi(−z) = 2e ∓πi/3 Ai(ze ±iπ/3 ), which does not have zeros as ℜz → +∞ (all the zeros of Ai(z) are real and negative); we have in this case: Inverting (18) and using the expansions (14) we have and following the steps [8, section 7.6.1] we conclude that the zeros of (18) as ℜz → +∞ are given by where T (t) has an asymptotic expansion given in [10, 9.2.11] and additional terms can be found in [3]. The first terms of this expansion for large t are Considering the particular case f (α) = − tan α, we have that the zeros of A(α, z) approaching the negative real axis as ℜz → −∞ are given by Particular cases are the real zeros of Ai(z) (α = 0) and Bi(z) (α = π/2), which correspond to Eqs. 9.9.6 and 9.9.10 of [10].
Remark 1 Observe that the smallest possible value of k in (22) may vary depending on the value of α. For instance, comparing with Eqs. 9.9.6 and 9.9.10 of [10] it is clear that k = 1, 2, ... for α = 0, π/2; but for α = −π/2 (which algo gives the zeros of Bi(z)) we should take k = 0, 1, .... In order to be sure that the smallest possible positive value of t k is meaningful we should make a detailed counting of zeros, using the argument principle as in [9]; we will not consider this analysis here.
The analysis of the zeros of the first derivative is very similar and leads to analogous results. One would start with Eqs. 9.7.10 and 9.7.12 of [10], which we can write as Ai ′ (−z) ∼ z 1/4 √ π (cos(ϕ − π/2)R(ξ) + sin(ϕ − π/2)S(ξ)), where R(ξ) and S(ξ) have the Poincaré-type expansions in a sector | arg z| ≤ 2π/3 − δ, with v m numbers defined in [10, 9.7(i)]. Then, the equation Therefore, the first approximation for the zeros of the derivatives can be obtained from those of the function by replacing χ by χ + π/2. In other words, in first approximation the location of the nodal curve with respect to the axis is the same for the zeros and the zeros of the derivative, as it only depends on ℑχ. This is consistent with the previous analysis based on the assumption that A(z) is locally constant: the zeros of a given solution and those of the derivative lie on the same ASL under this approximation.
For this reason, we will concentrate mainly on the zeros of the solutions of the ODEs, and not so much on the zeros of the derivative. The results we will enunciate, except for specific asymptotic estimates, can be applied both for the function and its derivative. For instance, Theorem 1 is also true for the derivative and we have: Theorem 2 Up to arbitrary multiplicative factors, Ai(ze ±iπ/3 ) are the only solutions of the Airy equation without zeros of its first derivative as ℜz → −∞ The same type of annotations can be made with respect to the results appearing in subsequent sections.

Zeros tending to the rays arg z = ±π/3 and general results
If α = 0, there is also an infinite number of zeros tending to the ray arg z = π/3 (and similarly for arg z = −π/3). For α = 0 we are in the case of the Airy function Ai(z), which only has negative real zeros.
Let us see, depending on α, when there is a string of zeros above the ray arg z = π/3 and when it lies below this ray.

Also, as a consequence
Theorem 4 Up to arbitrary multiplicative factors, A(±π/6, z) are the only solutions of the Airy equation with zeros on the three principal anti-Stokes lines (the rays arg z = 0, ±π/3).
On the other hand, similarly as in the case of Theorem 1, it is easy to see that when in (29) we have 1 − 2e 2iα = ±1 (α = 0 or ℑα → +∞), there are no zeros of A(α, ze iπ/3 ) as ℜz → +∞, and this corresponds to the fact that Ai(z) and Ai(ze −iπ/3 ) are, up to multiplicative factors, the only solutions without zeros approaching the ray arg z = π/3. This is summarized in the following theorem: Theorem 5 All the solutions of the Airy equation have zeros tending asymptotically to the three rays arg z = −π, ±π/3 with the only exception (up to multiplicative factors) of Ai(z), Ai(ze iπ/3 ) and Ai(ze −iπ/3 ) which only have zero on one ray (arg z = π, arg z = π/3 and arg z = −π/3 respectively).
The same arguments are true for the zeros of A ′ (α, ze ±iπ/3 ) = 0 and the results of this section can be easily adapted to these zeros. The starting point would be to replace (29) and (30).
The only thing left to describe is the position of the zeros with respect to the principal ASLs when the zeros are not on the rays themselves and to describe some asymptotic approximations. We already did this for the ray arg z = π, and now we complete the analysis for the other two principal ASLs.
From (29) we conclude that the zeros of A(α, ze iπ/3 ) as ℜz → +∞ are given by (20) with Setting the right-hand side to zero in (20) gives the first approximation Therefore, there are zeros of A(α, ze iπ/3 ) above (below) the positive real axis if |1 − e −2iα | is greater (smaller) than 1. This means that the zeros of A(α, z) approach the ASL arg z = π/3 from above if |1 − e −2iα | is greater than 1 and from below if it is smaller. The first approximation for the zeros of zeros of A(α, ze −iπ/3 ) as ℜz → +∞ can be obtained in a similar way. Now which gives the first approximation and the zeros approach the ASL arg z = −π/3 from below (above) if |1 − e −2iα | is greater (smaller) than 1.
For the zeros of the derivative similar approximations hold, with χ replaced by χ + π/2, and the location of the zeros of the derivative with respect to the axis is the same.
From the discussion of section (2.1), we see that the values ofz = e iπ/3 z = e iπ/3 3 2 ξ 2/3 corresponding to (32) (that is, the approximations for the zeros of A(α, z) approaching arg z = π/3) lie exactly over the ASL (in the LG approximation) given by (7), with where the angle θ in (7) is such that θ Similarly, the first approximations for the zeros of A(α, z) approaching arg z = −π/3 lie exactly over the LG-ASL with and The same is true for the zeros of the first derivative.
Similarly, A(α, z) has zeros approaching the ray arg z = −π/3 which are given by For the zeros of the derivative A ′ (α, z) we have zeros at where k ∈ Z with ℜt k > 0 (see Remark 1). U (t) has the asymptotic expansion [10, 9.9.19]; the first terms of the expansion are

Zeros of Bessel functions
For the zeros of Bessel functions we follow a similar scheme as for Airy functions with one important difference. As before, we start by analyzing the structure of the Stokes and anti-Stokes lines and then a more detailed analysis is performed using asymptotics. But differently from the Airy case, asymptotics for large z will be not enough to obtain a complete picture of the distribution of zeros, and we will also need to consider uniform asymptotic approximations for large order [9] in order to describe some of the zeros.
As happened for Airy functions, the picture of the zeros of the first derivative is very similar to that of the zeros of the function, and their relative position with respect to the anti-Stokes lines the same. We are not carrying a detailed description of the zeros of the first derivative, and we accept that the correspondence can be made very easily.

Stokes and anti-Stokes lines for the Bessel equation in the LG approximation
We consider the differential equation satisfied by the Riccati-Bessel functions Let us focus on the case of Bessel functions of real orders |ν| > 1/2. The differential equation has two simple turning points at z ± = ± ν 2 − 1/4 and because sign(A ′ (z ± )) = sign(z ± ) the principal ASLs emerge from z ± at angles 0, ±2π/3 for z + and π, ±π/3 for z − .
Writing the Riccati-Bessel equation w ′′ (z) + (1 − λ 2 /z 2 )w(z) = 0, λ = ν 2 − 1/4 in the variable η = z/λ, the transformed equation has coefficient A(η) = λ 2 (1 − 1/η 2 ). For real λ the ASLs emerging from the η = ±1 which are given by ℑ z 1 A(u)du = 0 can be written as in [14,Eq. (36)] or equivalently as the curve in the complex plane F = 1, where The principal lines are shown in Fig. 1. The eye-shaped region cuts the imaginary axis at η = ±ic where c is the real root of The rest of ASLs are also given by F = C with C a constant; for 0 < C < 1 we have closed curves inside the eye-shaped region and for C > 1 a curve for ℑz > 0 with horizontal asymptotes as ℜz → ±∞ and its complex conjugated curve for ℑz < 0. F < 1 gives the interior of the eye-shaped region and F > 1 the exterior. According to the scheme of ASLs depicted in Fig. 1 we expect that several types of zeros may appear (in the LG approximation): 1. A string of zeros as ℜz → +∞ and with asymptote ℑz = d for some real value d. Only one string (or none) of this type can be expected, either with positive or negative d. 3. Zeros either inside or above/below the eye shaped region. If there are zeros inside the eye-shaped region, they are located on a closed curve given by F = C < 1 for some positive C. If the zeros are outside (F > 1), they may appear to lie on a same curve as the zeros as ℜz → ±∞. Also, zeros on the principal ASL forming the eye-shaped region are possible. This is coherent with the asymptotic behavior for the zeros of Bessel functions of n integer (see [11,10.21(ix)]), but also for other cases described in the literature like Hankel functions of real order [1] or the Bessel functions J ν (z) and Y ν (z) of semi-integer order [2]. For instance, the Hankel function H (1) ν (z) has zeros close to the eye-shaped region for ℑz < 0 (for an illustration for integer order see [11, Fig. 10.21.4]) and, depending on the value of ν (and as we will later see) an additional string of zeros below the cut [1]. On the other hand the Bessel function Y n (z) of integer order has all the three types of zeros, with the zeros for ℑz < 0 complex conjugated of those for ℑz > 0; in fact, this is true for Y ν (z) with real orders and more generally also for cos αJ ν (z) − sin αY ν (z) for some values of α, as we will later see.
In the following, we will provide information on the location of the zeros for general solutions of the Bessel function of real orders, by considering the solutions We can restrict the analysis to ν ≥ 0, because using the connections formulas of [11, 10.4], one readily sees that We will consider arg z ∈ (−π, π], that is, we will stay in the principal Riemann sheet; the zeros in other Riemann sheets can be expressed in terms of zeros in the principal Riemann sheet using continuation formulas [11, 10.11]. We discuss in some detail the case of real α, but letting α be complex results for general combinations with complex coefficients will become available. It is in particular convenient to analyze the limiting cases ℑα → ±∞, which correspond to the zeros of Hankel functions H (1) (z) and H (2) (z).

Zeros of Hankel functions
Hankel functions H (1) (z) and H (2)  That these functions have no zeros as ℜ(z) → +∞ is obvious from its asymptotic behaviour [11, 10.17.5-6]. That the rest of solutions have zeros as ℜz → +∞ follows from this same asymptotic estimations, by considering the combination (43) with general complex α. We have that, as ℜz → +∞ the equation C ν (α, z) = 0 implies that and inverting we have the estimation for the zeros Therefore, for any α ∈ C we have zeros that run parallel to the positive real axis with asymptote ℑz = −ℑα; only if ℑα → ±∞ these zeros do not exist (case of the Hankel functions). Eq. (46) is the starting point for MacMahon asymptotic expansions [11, 10.21.19]; in [8, Example 7.9] additional terms in the expansion where given. The same expansions are valid in the general case with β in [8,7.33] replaced by β = −α + ν 2 − 1 4 + k π (this parameter is denoted as a in [11]).
Regarding the zeros close to the eye-shaped region, H ν (z) has zeros on the lower part (ℑz < 0) but not on the upper part, and the contrary happens with H (2) ν (z). These, together with J ν (z), which only has real zeros, are (up to multiplicative factors) the only solutions of the Bessel equation that do not have zeros both with positive and negative imaginary part in |ℜz| < ν for sufficiently large ν. We postpone this analysis for section 3.3.3.
Finally, H ν (z) may have zeros below the negative real axis, depending on the value of ν, while it does not have zeros over the negative real axis.
In the first place, is easy to see that there are no zeros over the negative real axis by using the continuation formula [11, 10.11.5]: and, after inverting, we see that the zeros of H ν (ze −πi ) = 0 can be estimated as where p = 0 if if cos νπ < 0 and p = 1 if cos νπ > 0. We observe that ℑz s > 0 if |2 cos νπ| > 1, which happens when {ν} ∈ [0, 1/3) ∪ (2/3, 1), where {ν} = ν − ⌊ν⌋ is the fractional part of ν. Only for these values does H ν (z) no longer has zeros below the negative imaginary axis in the principal Riemann sheet because, as {ν} becomes larger than 1/3 (or smaller than 2/3) these zeros are on the next Riemann sheet (arg(e −iπ z k ) < −π). For an analysis of the trajectories of this type of zeros depending on the value of ν we refer to [1].
In [1] it is also shown that the zeros below the negative real axis all cross the negative real axis precisely at {ν} = 1/3, 2/3. A way to see this is by using (48) and writing the Hankel functions in terms of J ν (z) and Y ν (z); with this we see that H Clearly A = 0, B = 0 and because J ν (z) and Y ν (z) are real for real positive z and do not have common zeros, the previous equality for z real and positive implies that A/B ∈ R, which only occurs when | cos νπ| = 1/2 and then {ν} = 1/3, 2/3. Therefore we conclude that only if {ν} / ∈ [1/3, 2/3] there are zeros of H (1) (z) below the imaginary axis and with ℑz s → − 1 2 log |2 cos νπ| as s → +∞ (ℜz s → −∞). The same is true for the zeros of H (1) ′ (z).

Zeros of
Here we provide a description of the distribution of zeros for general cylinder functions C ν (α, z) = cos αJ ν (z) − sin αY ν (z), where α may be complex. We will analyze in more detail some particular cases, like the case of real α, but the general results provide a complete picture on the distribution of zeros.

Zeros parallel to the positive real axis as ℜz → +∞ and MacMahon-type expansions
As explained before, for any α there is a string of zeros running parallel to the real axis as ℜz → +∞ which are approximately given by (46). Only as ℑα → ±∞ these zeros do not exist (case of the Hankel functions). We briefly summarize how asymptotic expansions can be obtained for these zeros; the same procedure can be applied to the zeros running parallel to the negative real axis (if any), that we will consider in the next section. Inverting C ν (α, z) = 0 is equivalent to inverting Y ν (z)/J ν (z) = cot α; we write χ = π/2 − α and then we consider the inversion of as ℜz → +∞. We can let χ be complex in general and therefore obtain expansions for general combinations of Bessel functions (with Hankel functions as the limits ℑχ → ±∞). Now, we write [11, 10.4] J ν (z) = 2 πz where ω = z − νπ 2 − π 4 . P(ν, z) and Q(ν, z) have asymptotic expansions in inverse powers of z for | arg z| < π − δ; P(ν, With this, the equation (52) becomes and inverting this relation gives Neglecting the last term (Q/P = O(z −1 )) we have our first approximation, and by re-substitution we can generate as many terms as needed of the asymptotic expansion for large s (see [8,Example 7.9]). This gives the MacMahon expansions where p i (µ) are certain polynomials of µ of degree i. The parameter β is given by In particular, if we are considering the zeros of C ν (α, z) = cos αJ ν (z) − sin αY ν (z) as ℜz → ∞, then χ = π/2 − α and And for α = 0 and α = π/2 we are in the cases of the positive zeros of J ν (z) and Y ν (z), and the expansions correspond to [11, 10.21.19] (where β is denoted as a). For any other real α, we have positive real zeros, while for complex α the zeros as ℜz → +∞ have asymptote ℑz = −ℑα.
The procedure for inverting (52) can be used to provide MacMahon-type expansions also for the case of the zeros parallel to the negative real axis (if any) just by replacing χ by its corresponding value.
For the zeros of the derivatives the analysis is very similar. We would start with [11, 10.17.9-10], which we can write as (cos(ω + π/2)P (ν, z) − sin(ω + π/2)Q(ν, z)), where, as before ω = z − νπ 2 − π 4 .P (ν, z) andQ(ν, z) have asymptotic expansions in inverse powers of z for | arg z| < π − δ;P (ν, With this, the equation (52) becomes The situation is analogous to (54) but with χ replaced by χ − π/2. It is then simple to write down the corresponding expansion by using the MacMahon expansions [11, 10.21.20] Because it is the imaginary part of the zeros which controls the relative position with respect to the real axis, the results for the zeros of C(α, z) and its first derivative are the same in this sense (because χ is shifted by a real amount); we will not consider again explicitly the case of the first derivative.

Zeros parallel to the branch cut as ℜz → −∞
Similarly as happened for the Hankel functions, for estimating the zeros running parallel to the branch cut as ℜz → −∞, we need to use continuation formulas so that the expansions (59) are applicable.
Using the first approximation (64) we have and there is an infinite number of zeros of C ν (α, ze imπ ) = 0 for ℜz > 0 with asymptote ℑz = a. Now, if we take m = 1 and it turns out that a < 0 then there is an infinite number of zeros of C ν (α, ze imπ ) with z below the positive real axis; therefore, C ν (z) will have an infinite number of zeros above the negative real axis. Contrary, if a > 0 the zeros of C α,ν (z) would be on the next Riemann sheet after turning with angle π (because then arg(e iπ z k ) > π).
Let us consider the particular case of α real, α ∈ [0, π). Observe that for real parameters sign(a) = sign(B). Therefore, we have zeros over the branch cut (and also the conjugated zeros below the branch cut) only if B < 0. Then the set of values for which there exist zeros above (and below) the branch cut is: We note that as, α varies, the zeros of C ν (z) can be at the branch cut only if B = 0, because C(α, e imπ z) = 0 implies that is real for real and positive z. Therefore, if the zeros parallel to the branch cut as ℜz → −∞ cross the negative real axis as α varies, they will all cross this axis for a same value of α. The zeros lie over the negative real axis in the following cases:
For the particular case of integer order, there exist zeros over and below the branch cut for any α = 0, with asymptotes as ℜz → −∞ given by And for α = 0, the result 10.21.46 of [11] is reproduced.
For the more general case of complex parameters, there is a string of zeros of C ν (α, z) running above the negative real with asymptote if a < 0. Similarly (taking m = −1 in the previous analysis), there is a string of zeros running below the negative real with asymptote The same is true for the zeros of the first derivative.

Zeros close to the eye-shaped region or inside the eye-shaped region
The zeros close to the eye-shaped region are not so simple to analyze. Poincaré asymptotics is of no use here, but uniform asymptotics for large ν gives a clear picture. In 1954, Olver developed powerful uniform asymptotic expansions for Bessel functions of large order [9], uniformly valid forz = z/ν ∈ (0, +∞) as the order goes to infinity. These are expansions with Airy functions as main approximants, and they read [11, 10.20.4-5]: with ζ the solution of the differential equation that is infinitely differentiable on the interval 0 <z < ∞, includingz = 1, and continued analytically to the complex plane cut along the negative real axis. From these approximations, we observe that for large enough ν the zeros of Bessel functions C ν (α, z) = cos αJ ν (z) − sin αY ν (z), z = νz, are given, in first approximation, by the zeros of A(α, ν 2/3 ζ) = cos αAi(ν 2/3 ζ) + sin αBi(ν 2/3 ζ). Therefore, for ν sufficiently large, we can infer the distribution of the zeros of Bessel functions from that of Airy functions that we studied before.
Similarly, we can deduce the distribution of zeros of derivative C ′ ν (α, z) from that of A ′ (α, ν 2/3 ζ) (see [10.20.7-8] [11]). And because the zeros of Airy functions and of their first derivative lie on the same curves (in first approximation), the discussion for the zeros of C ′ ν (α, z) will be similar to that of C ν (α, z).
In this section we will not be concerned with the analysis of the asymptotic expansions of the zeros for large ν. Our main goal is to determine when the zeros related to the eye-shaped region exist and to describe the curve (ASL) where they lie. In particular, we will determine the intersection of this ASLs with the imaginary axis, and also with the real axis for curves inside the eye-shaped region. As we discuss later, this is an important piece of information for the numerical computation of these zeros, particularly the intersection with the imaginary axis.
It is important to bear in mind the correspondence between the values ofz and those of ζ ( Figure 2). Observe that the curves BP 1 E 1 and BP 2 E 2 in thez-plane (the limits of the eye shaped region) correspond to the line segments respectively. Therefore, for analyzing the zeros close to the eye-shaped region we need to consider the zeros of Airy functions as arg z → ±π/3, |z| → ∞ (section 2.2.2). We consider the zeros approaching arg z = π/3, which give the zeros corresponding to the lower part of the eye-shaped region; the description for the upper part is analogous but considering the direction arg z = −π/3. We note the similarity between the domains determined by the ASLs in section 3.1 and the different regions in thez-domain of Fig. 2. The only difference is that the eye-shaped curve defined from (41) and setting F (η) = 1 is given in terms of η = z/ ν 2 − 1/2, while the equation for the curve E 1 P 1 BP 2 E 2 of Fig.2 is given in terms ofz = z/ν (F (z) = 1). In practice, those curves are very similar if ν is not small, and we expect that the information we will obtain from asymptotics is consistent with the LG approximation. The eye-shaped domain inside the curve E 1 P 1 BP 2 E 2 , which will be denoted as K, is with F given by (41).
In the ζ plane, we will have zeros below (above) the ray arg ζ = π/3 if |1 − e −2iα | is smaller (greater) than 1 (section 2.2.2). If the zeros are below the ray arg ζ = π/3, this gives a finite number of zeros because the boundary given by the curve D 2 E 2 in the ζ-domain would be reached by the string of Airy zeros as |z| becomes large; this corresponds to zeros inside the eye-shaped region and is in agreement with the picture given by the LG approximation. If the zeros are above the ray arg z = π/3 then, in the LG approximation for Airy functions, we have an ASL with an infinite number of zeros both as arg ζ → π/3 and as arg ζ → π; as arg ζ → π, in the z variable this would give a string of zeros asymptotically parallel to the positive real axis while, as arg ζ → π/3, we have a string of zeros below the negative real axis. This however, is an approximate picture and, as we have already shown, the string of zeros as arg ζ → π/3 may be absent in the principal Riemann sheet; it is important to bear in mind that in the ζ-domain there is a branch cut given by the ray arg z = π/3 from the point E 2 to infinity (there is also the cut arg z = −π/3 from the point E 1 to infinity).
we can put in correspondence θ andz as follows, where we have used that cos 3θ 2 > 0. When θ = 0, ζ becomes real and in (0, 1) and the corresponding value ofz gives the cut of the curve containing the zeros inside the eye-shaped region with the real axis. We have that the curve of zeros of C ν (α, z) cuts the positive real x-axis at νx 0 withx 0 the positive real root of and it also cuts the negative real axis at −νx 0 . As θ increases from θ = 0, the curve in thez plane moves towards the negative imaginary real axis until a value of θ in (0, π/3) is reached such that becomes purely imaginary. Setting z = −iỹ 0 in (75) the real part gives: In a similar way, it is easy to check that in the case |1 − e −2iα | > 1, the cut with the imaginary axisz = −iỹ 0 is also given by the solution of (77). For the case |1 − e −2iα | = 1 the solution of this equation (g(ỹ 0 ) = 0) gives the cut of the boundary F (z) = 1 with the imaginary axis (−iy 0 ). For the Airy-type zeros with positive imaginary part (corresponding to arg ζ → −π/3) the same equations but replacing α by −α give the cut with the axis.
Summarizing, the cuts of the curve containing the Airy-type zeros of C ν (α, z) with the imaginary axis are jiνỹ j , j = ±1, withỹ j the (positive) solution of We observe that the function g(y) is monotonic and g(0 + ) = +∞, g(+∞) = −∞, therefore the value of the cut with the positive (or the negative) imaginary axis, as expected, can be any positive real value because 1 2ν log |1 − e ±2iα | can take any possible real values for complex α. However, as ν → +∞ for a fixed α all the zeros tend to cluster on the boundary of the domain K.
The approximation of the point of intersection of the curve containing the Airy-type zeros with the imaginary axis given by (78) turns out to be accurate, not only for large ν but also for small ν. It typically gives the imaginary part of the zero closest to the imaginary axis with 2 − 3 digits for ν > 4 and improving as the order becomes larger.

20
In [14], an algorithm for computing complex zeros of special functions satisfying second order ODEs y ′′ (z) + A(z)y(z) = 0 was constructed which was based on the assumption that the coefficient A(z) is "locally constant". The idea of the method is to follow the ASLs in the Liouville-Green approximation by taking steps in the direction of the ASLs, and to apply a fourth order fixed point iteration x n+1 = g(x n ) with iteration function in order to compute the zeros. Once a first zero z 1 is computed by iterating (79) with an starting value conveniently chosen, the method takes a stepẑ = ±π/ A(z 1 ), where the sign is preferably chosen in such a way that the step is taken in the direction of decreasing |A(z)|. If A(z) was constant,ẑ would be another zero. In generalẑ will not be another zero, but it will be a value which gives convergence to the next zero in the ASL by iterating this value with (79). Once a second zero z 2 is computed by iteratingẑ with (79) we would take another stepẑ = ±π/ A(z 2 ), with the sign chosen as before, and so on. The main difficulty consists in obtaining good initial values for computing the first zero in each ASL. Although the good global convergence properties of (79) are such that it is possible to build a method not using such estimations, it is in any case convenient to have sharper first estimates in order to improve the speed of computation of the first zero. On the other hand, it is useful to know in advance which are the approximate ASLs where the zeros lie in order to avoid testing regions where there are no zeros. This information is contained in the present paper. Furthermore, if a zero of the solution is known, with this information it is possible to determine the rest of zeros (or a finite amount of them) reliably.

Algorithm for Airy functions
Consider first the case of the Airy functions and let suppose that we are interested in computing the zeros of a solution with a zero at z 0 . Such solution is proportional to y(z) = Bi(z 0 )Ai(z) − Ai(z 0 )Bi(z 0 ); then, taking we have that the zeros of y(z) are the zeros of A(α, z).
For computing the zeros of A(α, z) we can consider initial values given by the asymptotic approximations, which will be more accurate as |z| is larger, and then take the steps (80) in the direction of decreasing |z|. In fact, the asymptotic approximations are not really necessary for our method and taking initial values close to the principal ASLs is enough (though using the expansions for estimating a first zero is convenient).
For the zeros of the derivative the same procedure can be considered, but replacing (79) by a fixed point giving convergence to the zeros of A ′ (α, z). A possibility is to consider g(z) = z + 1 which gives convergence to the zeros of y ′ (z). Contrary to (79), which has order of convergence 4, (82) has order 2, like Newton's method. This fixed point method was discussed in [6] for the case of real zeros and it was proved to be globally convergent when A(z) is monotonic [13,Theorem 3.2]. A better possibility consists in applying the idea of [7] for computing zeros of the derivatives of solutions of ODEs: take the derivative of the differential equation and eliminate y(z) by using the same differential equation; this gives a second order ODE for w(z) = y ′ (z); transform to normal form and compute the corresponding fourth order fixed point method. A straightforward computation gives the following fourth order fixed point method for computing the zeros of A ′ (α, z):

Algorithm for Bessel functions
A simple but effective algorithm for computing complex zeros of Bessel functions C ν (α, z) was given in [14]. The method can be largely improved by using the information given in the present paper. In addition, as in the case of Airy functions, we can design a method that, given a zero z 0 of a solution of the ODE is able to compute the rest of zeros (a finite amount of them) just by considering that α = arctan(J ν (z 0 )/Y ν (z 0 )). As discussed before, there are always zeros as ℜz → +∞ except for Hankel functions. We compute these zeros starting from large ℜz and use (80) with the minus sign. We stop when ℜz < ν is reached (when the region of Airy-type zeros is reached). The MacMahon expansions can be used for the estimation of a first large zero. A simpler and also very effective starting value is z = L − iℑα (see Eqs. (56) and (58)), with L a positive value that can be chosen at will (depending on how many zeros are wanted).
With respect to the zeros running parallel to the negative axis as ℜz → −∞, the information given in section 3.3.2 can be used to determine when there are zeros above and/or below the negative real axis (see also section 3.2 for the case of Hankel functions). When there are zeros, it is possible to compute them starting from large |z| and in the direction of decreasing |z| (using (80 with the plus sign) until ℜz > −ν; the starting value can be obtained from McMahon-type expansions, with a first approximation given by (64) with m = ±1. It is simple and equally effective to consider a value −L − ia for zeros above the branch cut and −L − ib below the branch cut, with a and b given in (67) and (68).
Finally, for computing the Airy-type zeros, we can give good starting values using the information of section 3.3.3. We compute the intersection of the ASLs containing these zeros with the imaginary axis solving the equations (78). After computing one of such intersections, we iterate with (79) to compute a first zero; then we move to the right to compute subsequent zeros using the steps (80) with plus sign; and starting again from this first zero we move to the left with the minus sign. In both cases, we stop when |ℜz| > ν is reached or when the imaginary part of z becomes very small or changes sign (this last criterion is needed for the zeros inside the eye-shaped region).
For the zeros of the first derivative, the same scheme works, but the fixed point method (79) has to be replaced. We can use the fixed point iteration developed for the real zeros of γC ν (α, z) + zC ′ ν (α, z) in [7], which is of fourth order (also for complex zeros). Taking γ = 0 the fixed point iteration reads: Maple and Fortran codes implementing these methods are under construction [5].