Bessel, Airy and Coulomb wave functions

The Bessel, Airy and Coulomb wave functions are confluent hypergeometric functions: each of them is a combination of ${}_0F_1$ or ${}_1F_1$ series near the origin and of $U$ functions with exponential prefactors near infinity. FLINT evaluates them with the general methods of Hypergeometric functions, plus a few specializations: the cancellation estimate that decides the working precision, propagated error bounds based on derivative bounds, dedicated error bounds for the Airy functions, and interval Newton iteration for the zeros of the Airy functions.

Bessel functions#

For all complex $\nu$ and $z$,

with the regularized ${}_0F_1$ for negative integer $\nu$, and the asymptotic expansions come from

The prefactors $A_\pm$ simplify to $(\pm i)^{n}(2\pi z)^{-1/2}$ when $\nu + \tfrac12 = n$ is an integer and to $\exp(\mp i [(2\nu+1)/4]\pi)(2\pi z)^{-1/2}$ when $\operatorname{Re}(z) > 0$; the general formula is used otherwise. The functions of the second kind follow from

for non-integer $\nu$. At integer $\nu = n$, $K_n$ is computed as a limit with power series in $\nu$ (acb_hypgeom_bessel_k_0f1_series), and $Y_n(z) = -\frac{2}{\pi}(i^n K_n(iz) + [\log(iz) - \log z] J_n(z))$.

Choice of method. acb_hypgeom_bessel_j uses the convergent series for $|z| < 8$, and the asymptotic expansion for $|z| \ge 2^{64}$ or $2|z| \ge p$; the asymptotic expansion of $U^*$ reaches about $(1 + \log 2)/\log 2 \cdot |z| \approx 2.44|z|$ bits for small $\nu$, and the factor 2 is a safety margin. Between these regions the convergent series is used with a working precision raised by the estimated cancellation. The terms of the ${}_0F_1$ series for $J_\nu(x + yi)$ grow to about $e^{|z|}$ while the function has size about $e^{|y|}$, so about $(|z| - |y|)/\log 2$ bits are lost; no precision is added when the terms do not cancel (purely imaginary $z$ with real $\nu > 0$ or integer $\nu$, or $\nu > |z|^2/4$). The modified functions $I_\nu$ and $K_\nu$ use the series for $|z| < 16$ and for $2|z| < p$, and the asymptotic expansion otherwise.

Propagated error. For inexact $z$ and integer $\nu$, $J_\nu$ is evaluated at the midpoint and the error is bounded with acb_hypgeom_bessel_j_deriv_bound, which uses for example $|J_n^{(d)}(z)| \le e^{|\operatorname{Im} z|}$; this avoids the large overestimation that ball arithmetic produces when a wide ball is pushed through a series with cancellation. When the resulting enclosure is still poor, it is intersected with the direct bound for $|J_\nu(z)|$.

Large orders and integration. For real $\nu$ and $z$, if the result has less than half the target accuracy (typically because $\nu$ is large compared to $z$ and the precision), acb_hypgeom_bessel_i and acb_hypgeom_bessel_k retry with rigorous numerical integration of an integral representation (arb_hypgeom_bessel_i_integration, arb_hypgeom_bessel_k_integration).

Power series. The Taylor expansion of $J_\nu$ at a point $z_0 \ne 0$ follows from $J_\nu(z_0)$, $J'_\nu(z_0) = \frac12(J_{\nu-1}(z_0) - J_{\nu+1}(z_0))$ and a five-term recurrence for the higher coefficients derived from the Bessel differential equation $z^2 y'' + z y' + (z^2 - \nu^2) y = 0$. This is implemented generically in gr_bessel_j_jet and _gr_poly_bessel_j_series, which arb_hypgeom_bessel_j_series and acb_hypgeom_bessel_j_series call; the recurrence divides by $z_0^2$, so $z_0$ must be bounded away from 0 unless it is exactly 0 and $\nu$ is an integer.

Airy functions#

The Airy functions are entire solutions of $y'' = zy$. Near the origin FLINT sums

with $c_1 = \operatorname{Ai}(0) = 3^{-2/3}/\Gamma(2/3)$ and $c_2 = -\operatorname{Ai}'(0) = 3^{-1/3}/\Gamma(1/3)$ (cached constants), by rectangular splitting in powers of $z^3$ with word-size coefficients (acb_hypgeom_airy_direct). With $w = z^3/9$ the terms are bounded by $3|w|^n/((n-1)!)^2$, which gives the tail bound. Expressing the Airy functions through Bessel functions of order $\pm 1/3$ is avoided, since such formulas hold only in sectors (the Stokes phenomenon).

Asymptotic expansions#

The error bounds are based on Olver [NIST2012, §9.7]. For $|\arg z| < \pi$ and $\zeta = \frac23 z^{3/2}$,

For $n \ge 1$ the remainders are bounded by

with $\chi(n) = \sqrt\pi\,\Gamma(\frac n2 + 1)/\Gamma(\frac n2 + \frac12)$.

For $\operatorname{Bi}$ when $z$ is roughly in the right half-plane, the connection formulas

are combined with the expansions, which after collecting the roots of unity gives

where the upper signs are valid for $-\pi/3 < \arg(z) < \pi$ and the lower signs for $-\pi < \arg(z) < \pi/3$; the derivative is analogous. In the left half-plane the formulas

with $z_1 = -z e^{\pi i/3}$ and $z_2 = -z e^{-\pi i/3}$, reduce to arguments with $|\arg(z_1)|, |\arg(z_2)| < \pi$ when $|\arg(-z)| < 2\pi/3$, and collecting roots of unity gives, with $\zeta = \frac23 (-z)^{3/2}$,

The differentiated formulas are analogous. acb_hypgeom_airy_asymp implements these formulas, and acb_hypgeom_airy_bound evaluates them with $n = 1$ in mag_t arithmetic to get cheap upper bounds for the four functions on a ball.

Choice of method#

acb_hypgeom_airy uses the asymptotic expansion when $|z| \ge 4$ and the terms decrease below $2^{-p}$ (the number of terms comes from an estimate of $\log_2 |u(k)/\zeta^k|$), and the series otherwise. For the series it estimates in double precision the magnitude of the largest term, about $\exp(2|z^3/9|^{1/2})$, and the magnitude of the functions from the leading asymptotic term, and adds the difference to the working precision (this ignores the extra cancellation near the zeros). The series is evaluated at the midpoint of $z$; the error propagated from the radius $r$ is bounded by $r$ times the bounds of acb_hypgeom_airy_bound for the derivatives, using $\operatorname{Ai}''(z) = z \operatorname{Ai}(z)$.

Zeros#

arb_hypgeom_airy_zero computes the $n$-th real zeros $a_n, a'_n, b_n, b'_n$ of $\operatorname{Ai}, \operatorname{Ai}', \operatorname{Bi}, \operatorname{Bi}'$ from the asymptotic expansions [NIST2012, §9.9] $a_n \sim -T(\tfrac38\pi(4n-1))$, $b_n \sim -T(\tfrac38\pi(4n-3))$, $a'_n \sim -U(\tfrac38\pi(4n-3))$, $b'_n \sim -U(\tfrac38\pi(4n-1))$, where $T(t) = t^{2/3}(1 + \frac{5}{48}t^{-2} - \cdots)$ and $U(t) = t^{2/3}(1 - \frac{7}{48}t^{-2} + \cdots)$. Pittaluga and Sacripante [PS1991] proved that two consecutive truncations of the expansions for $a_n$ and $b_n$ bracket the zero; for $a'_n$ and $b'_n$ the truncated expansion gives a starting interval that is verified by an interval Newton step. For $n \le 10$ tabulated 48-bit enclosures are used, since the expansions are poor for the first zero. The enclosure is then refined by interval Newton steps $x \mapsto m - f(m)/f'(X)$ with doubling precision, using $f'' = xf$ for the derivatives.

Coulomb wave functions#

The Coulomb wave functions $F_\ell(\eta, z)$, $G_\ell(\eta, z)$ and $H^\pm_\ell(\eta,z) = G_\ell \pm i F_\ell$ solve $y'' + (1 - 2\eta/z - \ell(\ell+1)/z^2) y = 0$. With the normalization constant

(computed through $\log\Gamma(1+\ell \pm i\eta)$ so that it is analytic in $\ell$ and $\eta$), the regular solution is

and all four functions are expressed through $U(1+\ell \pm i\eta, 2\ell+2, \mp 2iz)$ with exponential and power prefactors. The formulas, with the branch cuts placed on the negative real $z$ axis, follow the connection formulas collected in [Fungrim] (see also [DYF1999; Gas2018; Mic2007] and [NIST2012, chapter 33]); different formulas are used in the upper and lower half-planes and on the real line so that each function is computed without unnecessary cancellation, and $H^\pm$ are computed directly rather than as $G \pm iF$, which can cancel. The $U$ functions are evaluated with the asymptotic expansion when $|z|$ is large enough and with acb_hypgeom_u otherwise. The jets (acb_hypgeom_coulomb_jet) use $F$ and $G$ at $\ell$ and $\ell + 1$ for the first derivative and the differential equation for the higher coefficients.

Use in FLINT#

• Bessel functions: acb_hypgeom_bessel_j, acb_hypgeom_bessel_y, acb_hypgeom_bessel_jy, acb_hypgeom_bessel_i, acb_hypgeom_bessel_k and their scaled versions, with the methods acb_hypgeom_bessel_j_0f1, acb_hypgeom_bessel_j_asymp, acb_hypgeom_bessel_i_0f1, acb_hypgeom_bessel_i_asymp, acb_hypgeom_bessel_k_0f1, acb_hypgeom_bessel_k_asymp; the real wrappers in arb_hypgeom; the generic gr_bessel_j_jet and _gr_poly_bessel_j_series.
• ${}_0F_1$: acb_hypgeom_0f1 uses the Bessel connection formulas for its asymptotic expansion.
• Airy functions: acb_hypgeom_airy, acb_hypgeom_airy_direct, acb_hypgeom_airy_asymp, acb_hypgeom_airy_bound, acb_hypgeom_airy_jet, the series versions and arb_hypgeom_airy_zero.
• Coulomb wave functions: acb_hypgeom_coulomb, acb_hypgeom_coulomb_jet, acb_hypgeom_coulomb_series and the real versions.

@misc{flint-doc-algorithms-special-functions-bessel-airy,
  author       = {{The FLINT team}},
  title        = {{Bessel, Airy and Coulomb wave functions}},
  howpublished = {FLINT: Fast Library for Number Theory, documentation, version 3.7.0-dev},
  year         = {2026},
  url          = {https://flintlib.org/doc/algorithms/special-functions/bessel-airy.html},
  note         = {Generated 2026-11-02}
}