The Airy functions, denoted Ai(z) and Bi(z), are a pair of linearly independent solutions to the second-order linear differential equation w''(z) - z w(z) = 0, known as Airy's equation.[1] These entire functions, named after the Britishastronomer and mathematicianGeorge Biddell Airy, were first introduced in his 1838 paper analyzing the intensity of light near optical caustics within the framework of the wave theory of light.[2] Ai(z) is the principal solution, characterized by its oscillatory behavior for negative real arguments and exponential decay as z \to +\infty along the positive real axis, while Bi(z) exhibits exponential growth in that direction, making Ai(z) suitable for bounded physical problems. Their values at the origin are Ai(0) ≈ 0.355028 and Bi(0) ≈ 0.614927, with corresponding derivatives Ai'(0) ≈ -0.258819 and Bi'(0) ≈ 0.448289.[1]Asymptotic expansions provide key insights into their large-argument behavior: for |\arg z| \leq \pi - \delta, Ai(z) ∼ \frac{e^{-\zeta}}{2 \sqrt{\pi} z^{1/4}} \sum_{k=0}^\infty (-1)^k \frac{u_k}{\zeta^k}, where \zeta = \frac{2}{3} z^{3/2} and the u_k are coefficients defined recursively; similarly, for |\arg z| \leq \frac{1}{3}\pi - \delta, Bi(z) ∼ \frac{e^{\zeta}}{\sqrt{\pi} z^{1/4}} \sum_{k=0}^\infty \frac{u_k}{\zeta^k}.[3] These expansions, valid uniformly in specified sectors, underpin uniform approximations near turning points in various boundary-value problems.[3]In mathematics, Airy functions serve as prototypes for studying special functions, with integral representations such as Ai(z) = \frac{1}{\pi} \int_0^\infty \cos\left(\frac{t^3}{3} + z t\right) \, dt for real z > 0,[4] and connections to other functions like Bessel functions via limits or transformations. They also appear in the theory of linear differential equations with turning points and in asymptotic analysis of integrals.[5]Physically, Airy functions model phenomena across classical and quantum domains, including light intensity near caustics and rainbows in optics, electromagnetic diffraction and radiowave propagation, stability in viscous fluid flows via the Orr-Sommerfeld equation, nonlinear wave propagation such as solitons in shallow water, and quantum scattering in linear potentials as solutions to the Schrödinger equation.[2] Their role in semiclassical approximations extends to rainbow scattering and tunneling effects in quantum mechanics.[2]
Mathematical Foundations
Differential Equation
The Airy differential equation is a second-order linear ordinary differential equation of the formw''(z) - z w(z) = 0,where primes denote differentiation with respect to the argument z. Equivalently, in terms of a variable x, it can be written as y''(x) - x y(x) = 0. This equation belongs to the class of Sturm-Liouville problems, where the eigenvalue equation takes the form -u''(x) + x u(x) = \lambda u(x), and it features a turning point at x = 0 where the nature of the solutions transitions from oscillatory to exponential behavior.[1][6]The two linearly independent solutions to this equation are the Airy function of the first kind, \operatorname{Ai}(z), and the Airy function of the second kind, \operatorname{Bi}(z). The function \operatorname{Ai}(z) decays exponentially for large positive real z, while \operatorname{Bi}(z) grows exponentially in that regime, making \operatorname{Ai}(z) particularly useful for bounded solutions in certain physical contexts. These solutions are entire functions of z, analytic everywhere in the complex plane.[1][7]The functions are normalized such that at z = 0,\operatorname{Ai}(0) = \frac{1}{3^{2/3} \Gamma\left( \frac{2}{3} \right)}, \quad \operatorname{Ai}'(0) = -\frac{1}{3^{1/3} \Gamma\left( \frac{1}{3} \right)},and\operatorname{Bi}(0) = \frac{1}{3^{1/6} \Gamma\left( \frac{2}{3} \right)}, \quad \operatorname{Bi}'(0) = \frac{3^{1/6}}{\Gamma\left( \frac{1}{3} \right)},where \Gamma denotes the gamma function. This normalization ensures consistency across applications and computational implementations.[1]The Airy equation arises naturally in turning-point problems for ordinary differential equations, particularly those modeling systems with a linear potential, such as a particle under constant force (e.g., gravitational acceleration). In such cases, the Schrödinger equation or similar wave equations reduce to this form after a change of variables, capturing the transition between classically allowed and forbidden regions.[1][8][6]
Integral Representations
The Airy functions can be expressed through various integral representations, which offer explicit formulas for computation and facilitate the derivation of their properties, distinct from their definition via the differential equation. These representations are particularly useful for complex arguments and enable analysis in different regions of the complex plane.A primary contour integral representation for the Airy function of the first kind is\mathrm{Ai}(z) = \frac{1}{2\pi i} \int_{\infty e^{-\pi i /3}}^{\infty e^{\pi i /3}} \exp\left(\frac{t^3}{3} - z t\right) \, dt,where the Hankel contour starts at infinity along the ray \arg t = -\pi/3, proceeds to the origin while avoiding it from above, encircles the origin counterclockwise, and returns to infinity along \arg t = \pi/3. This form converges for all finite complex z and defines \mathrm{Ai}(z) as an entire function.[4]For the Airy function of the second kind,\mathrm{Bi}(z) = \frac{1}{\pi} \left[ \int_{-\infty}^{(0+)} \exp\left( -\frac{t^3}{3} + z t \right) \, dt + \int_{0+}^{\infty} \sin\left( \frac{t^3}{3} + z t \right) \, dt \right],though a complex contour variant is\mathrm{Bi}(z) = \frac{1}{\pi} \left[ \int_{-\infty e^{\pi i /3}}^{(0+)} + \int_{-\infty e^{-\pi i /3}}^{(0+)} \right] \exp\left( \frac{t^3}{3} - z t \right) \, dt,with paths along the indicated rays to a small circle around the origin. These contours ensure convergence across the complex plane, with the choice of rays determined by the phase of z to dampen the exponential growth in the integrand. For instance, in sectors where |\arg z| < \pi/3, the integrals along rays at \pm 2\pi/3 may be used for better numerical stability.[4]For real arguments x, simpler oscillatory and exponential integrals apply. The Airy function of the first kind is\mathrm{Ai}(x) = \frac{1}{\pi} \int_0^\infty \cos\left( \frac{t^3}{3} + x t \right) \, dt,which converges conditionally due to the oscillatory nature at large t, and can be accelerated via integration by parts or other methods for computation. Similarly,\mathrm{Bi}(x) = \frac{1}{\pi} \int_0^\infty \left[ \exp\left( -\frac{t^3}{3} + x t \right) + \sin\left( \frac{t^3}{3} + x t \right) \right] \, dt.These real forms are valid for all real x, with the exponential term in \mathrm{Bi}(x) diverging as x \to +\infty, reflecting its growth, while the sine term contributes to oscillatory behavior for x < 0. Convergence for the real integrals holds absolutely for x > 0 in the exponential parts and conditionally otherwise.[4]
Properties for Real Arguments
Oscillatory and Exponential Behavior
The Airy functions Ai(x) and Bi(x) exhibit distinct behaviors for positive and negative real arguments. For x > 0, Ai(x) is positive and monotonically decreases to zero as x increases, decaying exponentially with a leading asymptotic form of \operatorname{Ai}(x) \sim \frac{\exp\left(-\frac{2}{3} x^{3/2}\right)}{2\sqrt{\pi} x^{1/4}}. In contrast, Bi(x) is positive and monotonically increases to infinity as x increases, growing exponentially with a leading asymptotic form of \operatorname{Bi}(x) \sim \frac{\exp\left(\frac{2}{3} x^{3/2}\right)}{\sqrt{\pi} x^{1/4}}.For x < 0, both Ai(x) and Bi(x) display oscillatory behavior with amplitudes that increase as |x| grows. The oscillations have a frequency proportional to |x|^{1/2} and a phase involving \frac{2}{3} |x|^{3/2}, leading to progressively more rapid and larger-amplitude waves as x decreases. Specifically, Ai(x) can be expressed as M(x) \sin \theta(x) and Bi(x) as M(x) \cos \theta(x) for x ≤ 0, where M(x) is the modulus that increases with |x|, and the phase \theta(x) satisfies \theta(x) \sim \frac{\pi}{4} + \frac{2}{3} (-x)^{3/2} as x \to -\infty. This oscillatory pattern begins near x = 0 and intensifies toward negative infinity.A key relation governing these functions is their Wronskian, defined as W{\operatorname{Ai}(x), \operatorname{Bi}(x)} = \operatorname{Ai}(x) \operatorname{Bi}'(x) - \operatorname{Ai}'(x) \operatorname{Bi}(x) = \frac{1}{\pi}, which remains constant for all real x and underscores their linear independence.Numerical evaluations highlight these behaviors: at x = 0, \operatorname{Ai}(0) \approx 0.355028 and \operatorname{Bi}(0) \approx 0.614927, both positive. The first zero of Ai(x) occurs at approximately x \approx -2.338, marking the onset of its oscillations, after which Ai(x) crosses zero repeatedly with growing amplitude.Qualitatively, plots of Ai(x) and Bi(x) along the real line show Ai(x) starting at about 0.355 at x=0, decaying smoothly to near zero for large positive x, while crossing the axis multiple times for x < 0 with waves of increasing height and frequency. Bi(x) starts at about 0.615 at x=0, rises steeply for x > 0, and for x < 0 oscillates similarly to Ai(x) but phase-shifted by \pi/2, maintaining the same envelope of growing amplitude. These patterns can be computed using integral representations for practical evaluation.
Zeros and Orthogonality
The Airy function \mathrm{Ai}(x) has infinitely many zeros, all of which are real and located on the negative real axis, denoted a_k for k = 1, 2, \dots. The first few are a_1 \approx -2.338, a_2 \approx -4.088, and a_3 \approx -5.521. For large n, these zeros follow the asymptotic approximation a_n \approx -\left( \frac{3\pi (4n-1)}{8} \right)^{2/3}.[9]The Airy function \mathrm{Bi}(x) has infinitely many real zeros, all on the negative real axis, with the first few given by b_1 \approx -1.174, b_2 \approx -3.271, and b_3 \approx -4.831. Additionally, \mathrm{Bi}(x) possesses infinitely many complex zeros occurring in conjugate pairs within the sectors \frac{\pi}{3} < |\arg z| < \frac{\pi}{2}.[9]Airy functions satisfy orthogonality relations that underpin their use in expansions and spectral analysis. For discrete orthogonality, consider boundary value problems governed by the Airy equation in Sturm-Liouville form, such as -y'' + x y = \lambda y on [0, \infty) with the decaying condition as x \to \infty and Dirichlet condition y(0) = 0. The eigenvalues are \lambda_k = -a_k, and the eigenfunctions \mathrm{Ai}(x + a_k) are orthogonal in L^2(0, \infty), with squared norms \int_0^{\infty} [\mathrm{Ai}(x + a_k)]^2 \, dx = \frac{1}{[\mathrm{Ai}'(a_k)]^2}.[6][10]Any sufficiently regular function f(x) on [0, \infty) can be expanded as f(x) = \sum_{k=1}^{\infty} c_k \mathrm{Ai}(x + a_k), where the coefficients are c_k = [\mathrm{Ai}'(a_k)]^2 \int_0^{\infty} f(x) \mathrm{Ai}(x + a_k) \, dx.[6]These properties position the Airy functions centrally within Sturm-Liouville theory for the Airy operator -\frac{d^2}{dx^2} + x, facilitating eigenvalue expansions, Green's functions, and solutions to related differential equations.[6][10]
Asymptotic Approximations
Large Positive Arguments
For large positive arguments x \to +\infty, the Airy functions \operatorname{Ai}(x) and \operatorname{Bi}(x) exhibit distinct exponential behaviors, with \operatorname{Ai}(x) decaying rapidly and \operatorname{Bi}(x) growing exponentially, reflecting their roles as recessive and dominant solutions to the Airy differential equation, respectively.[3] These asymptotics are derived from the method of steepest descent applied to the integral representations of the functions and are essential for approximating solutions in quantum mechanics and wave propagation problems where x represents a scaled energy or potential parameter far from the classical turning point.[3]The leading asymptotic expansion for \operatorname{Ai}(x) is\operatorname{Ai}(x) \sim \frac{1}{2} \pi^{-1/2} x^{-1/4} \exp\left(-\zeta\right),where \zeta = \frac{2}{3} x^{3/2}, capturing the dominant exponential decay modulated by a slowly varying prefactor.[3] Similarly, for \operatorname{Bi}(x),\operatorname{Bi}(x) \sim \pi^{-1/2} x^{-1/4} \exp\left(\zeta\right),highlighting the exponential growth.[3] The exponential terms \exp(\pm \zeta) dominate the behavior, with the x^{-1/4} factor providing the subdominant algebraic correction arising from the WKB approximation near the turning point scaled to large x.[3]Higher-order terms extend these expansions as asymptotic series:\operatorname{Ai}(x) \sim \frac{\exp(-\zeta)}{2 \sqrt{\pi} x^{1/4}} \sum_{k=0}^{\infty} (-1)^k \frac{u_k}{\zeta^k}, \quad \operatorname{Bi}(x) \sim \frac{\exp(\zeta)}{\sqrt{\pi} x^{1/4}} \sum_{k=0}^{\infty} \frac{u_k}{\zeta^k},where the coefficients u_k satisfy the recurrence relation u_0 = 1 and u_{k} = \frac{(6k-5)(6k-3)(6k-1)}{216k(2k-1)} u_{k-1} for k \geq 1.[3] These series are divergent but provide successively better approximations when truncated optimally, emphasizing the exponential dominance over oscillatory phases, which are absent for positive arguments.[3]The Airy functions themselves offer a uniform approximation near the turning point x=0 via their inherent scaling, where the asymptotic expansions transition smoothly but lose accuracy as x \to 0^+ since \zeta \to 0 and the series terms grow.[3] These asymptotic expansions hold as x \to +\infty, uniformly in the sector |\arg x| \leq \pi - \delta for small \delta > 0, with relative error bounds of O(1/\zeta) for the leading term in \operatorname{Ai}(x) and a more refined bound incorporating \chi-function estimates for \operatorname{Bi}(x), ensuring reliability in physical applications for moderately large x \gtrsim 1.[3]
Large Negative Arguments
For large negative real arguments x \to -\infty, the Airy functions \mathrm{Ai}(x) and \mathrm{Bi}(x) exhibit oscillatory behavior, contrasting with their exponential decay and growth for positive arguments. The leading-order asymptotic expansions are given by\mathrm{Ai}(x) \sim \pi^{-1/2} |x|^{-1/4} \sin\left( \xi + \frac{\pi}{4} \right),\mathrm{Bi}(x) \sim \pi^{-1/2} |x|^{-1/4} \cos\left( \xi + \frac{\pi}{4} \right),where \xi = \frac{2}{3} |x|^{3/2}. These expansions hold uniformly in sectors of the complex plane excluding small neighborhoods around the anti-Stokes lines, with higher-order terms involving recursive coefficients u_k that refine the approximation.[3]The phase shift of \pi/4 in these expressions originates from the turning point analysis at x=0 in the WKB (Wentzel–Kramers–Brillouin) approximation to the Airy differential equation w''(x) = x w(x), where the oscillatory solution in the classically allowed region (x < 0) connects to the evanescent solution across the turning point with this characteristic shift to ensure matching.[11]The amplitude factor |x|^{-1/4} varies slowly compared to the rapidly increasing phase \xi, resulting in oscillations whose frequency grows as |x|^{1/2} while the envelope decreases gradually, producing increasingly rapid but diminishing waves as |x| enlarges. This behavior aligns with the WKB plane wave approximation \exp(\pm i \int^x \sqrt{-t} \, dt) in the allowed region away from the turning point, where the Airy functions provide the exact uniform transition without singularity at x=0.[11]In the complex plane, these oscillatory forms apply in sectors where the real part of the phase is negative, but crossing Stokes lines—rays from the origin at \arg z = 0, \pm 2\pi/3—triggers a switch to exponential dominance, with the subdominant oscillatory contribution multiplying by a Stokes multiplier (often i or -i) to maintain analytic continuation of the function. This Stokes phenomenon ensures the asymptotic representations remain valid across sectors by adjusting the relative contributions of exponentially growing and decaying terms.[12]
Extension to Complex Arguments
Analytic Continuation
The Airy functions \operatorname{Ai}(z) and \operatorname{Bi}(z), originally defined for real arguments through integral representations or asymptotic behaviors, admit a unique analytic continuation to the entire complex plane, where they become entire functions satisfying the differential equation w''(z) = z w(z) at every finite point.\] As entire functions, they are holomorphic everywhere in $\mathbb{C}$, with no poles or other singularities in the finite plane.\[ This continuation preserves their single-valued nature, distinguishing them from multi-valued functions that require branch cuts for definition across the complex domain.$$]The explicit form of the analytic continuation is provided by the Maclaurin series expansions around z = 0, which converge for all complex z and thus define the functions globally. For \operatorname{Ai}(z),[
\operatorname{Ai}(z) = \operatorname{Ai}(0) \sum_{n=0}^{\infty} \frac{1 \cdot 4 \cdots (3n-2)}{(3n)!} z^{3n} + \operatorname{Ai}'(0) \sum_{n=0}^{\infty} \frac{2 \cdot 5 \cdots (3n-1)}{(3n+1)!} z^{3n+1},
where the leading coefficients are $\operatorname{Ai}(0) = 3^{-2/3} / \Gamma(2/3)$ and $\operatorname{Ai}'(0) = -3^{-1/3} / \Gamma(1/3)$, both involving Gamma functions for precise evaluation.$$\] A similar power series holds for $\operatorname{Bi}(z)$,
\[
\operatorname{Bi}(z) = \operatorname{Bi}(0) \sum_{n=0}^{\infty} \frac{1 \cdot 4 \cdots (3n-2)}{(3n)!} z^{3n} + \operatorname{Bi}'(0) \sum_{n=0}^{\infty} \frac{2 \cdot 5 \cdots (3n-1)}{(3n+1)!} z^{3n+1},with \operatorname{Bi}(0) = 3^{-1/6} / \Gamma(2/3) and \operatorname{Bi}'(0) = 3^{1/6} / \Gamma(1/3).\] These series, derived from the differential equation and initial conditions, ensure the continuation is rigorous and uniform in compact subsets of $\mathbb{C}$.\[An alternative path for continuing \operatorname{Bi}(z) leverages the analyticity of \operatorname{Ai}(z) via the connection formula\operatorname{Bi}(z) = e^{-\pi i / 6} \operatorname{Ai}\left( z e^{-2 \pi i / 3} \right) + e^{\pi i / 6} \operatorname{Ai}\left( z e^{2 \pi i / 3} \right),which expresses \operatorname{Bi} in terms of rotated arguments of \operatorname{Ai}, valid for all complex z on the principal branch with |\arg z| \leq \pi.\] No Riemann surface is required for either function, as their entire character implies they are single-valued without branches in the finite plane.\[ At infinity, however, both exhibit an essential singularity, reflecting their order $3/2 growth behavior in certain sectors.$$]
Branch Cuts and Plots
The Stokes lines of the Airy functions in the complex plane are the rays emanating from the origin where \arg z = \pm \pi/3 and \arg z = \pi. These lines delineate sectors in which the relative dominance of exponential terms in the asymptotic expansions changes, marking transitions between regions of distinct qualitative behavior for large |z|.[12][13]In the principal sector |\arg z| < \pi/3, \mathrm{Ai}(z) decays exponentially as |z| increases, while \mathrm{Bi}(z) grows exponentially in the same region. Beyond these Stokes lines, in the sectors \pi/3 < |\arg z| < \pi, both functions transition to oscillatory behavior, characterized by complex exponentials with imaginary exponents leading to wave-like patterns. Anti-Stokes lines, perpendicular to the Stokes lines and defined where the phase \zeta = (2/3) z^{3/2} has constant imaginary part, further outline the loci of equal oscillatory phase.[12][13]Contour plots of |\mathrm{Ai}(z)| vividly demonstrate these regional differences, displaying compact regions of rapid exponential decay confined to |\arg z| < \pi/3 and extending outward into oscillatory zones with alternating ridges and valleys of magnitude in the adjacent sectors. Similarly, contour plots of \mathrm{Re}(\mathrm{Ai}(z)) reveal smooth, undulating patterns that intensify across the Stokes lines, underscoring the shift from monotonic decay to periodic variation without any discontinuities, as \mathrm{Ai}(z) is an entire function. Plots for \mathrm{Bi}(z) mirror this but inverted in the principal sector, with growth amplifying the contrast between regions.Although \mathrm{Ai}(z) and \mathrm{Bi}(z) are entire functions with no intrinsic branch cuts in the complex plane, certain integral representations or asymptotic forms for \mathrm{Bi}(z) introduce apparent branch cuts along the negative real axis (\arg z = \pi) or the rays \arg z = \pm 2\pi/3 to ensure single-valuedness of the phase factor in the contour integrals. These effective cuts align with secondary Stokes lines and guide the analytic continuation process.[14][15]Numerical evaluation of Airy functions in the complex domain is challenging due to exponential growth in \mathrm{Bi}(z) within |\arg z| < \pi/3 and severe cancellations from rapid oscillations elsewhere, which can lead to loss of precision. Computations typically start in the stable sector |\arg z| < \pi/3 using power series or integrals, then employ connection formulas to propagate values across Stokes lines into other regions while avoiding paths that traverse areas of numerical instability.For instance, along the imaginary ray z = i y with real y > 0 (\arg z = \pi/2), \mathrm{Ai}(i y) resides in an oscillatory sector and exhibits damped sinusoidal variations, with |\mathrm{Ai}(i y)| decreasing initially before modulating with amplitude growing like y^{-1/4} for large y. This behavior contrasts sharply with the monotonic decay along positive real arguments.[16]
Relations to Other Special Functions
Bessel and Modified Bessel Functions
The Airy functions Ai(z) and Bi(z) for positive real arguments z > 0 can be expressed in terms of modified Bessel functions of the first kind, I_{\nu}(\zeta), where \zeta = \frac{2}{3} z^{3/2}. Specifically,[
\operatorname{Ai}(z) = \frac{1}{3} \sqrt{z} \left[ I_{-1/3}(\zeta) - I_{1/3}(\zeta) \right],\operatorname{Bi}(z) = \sqrt{\frac{z}{3}} \left[ I_{-1/3}(\zeta) + I_{1/3}(\zeta) \right].
These representations arise from a change of variables that transforms the Airy differential equation $w''(z) - z w(z) = 0$ into the modified Bessel equation of order $\pm 1/3$.[](https://dlmf.nist.gov/9.6)
For negative real arguments, letting $x > 0$, the Airy functions exhibit oscillatory behavior and relate to ordinary [Bessel functions](/page/Bessel_function) of the first kind, $J_{\nu}(\eta)$, where $\eta = \frac{2}{3} x^{3/2}$. In particular,
\operatorname{Ai}(-x) = \frac{\sqrt{x}}{3} \left[ J_{1/3}(\eta) + J_{-1/3}(\eta) \right],
with a corresponding expression for $\operatorname{Bi}(-x)$ obtained by replacing the sum with an appropriate difference involving the same [Bessel functions](/page/Bessel_function). This connection follows from the same [transformation](/page/Transformation) of the Airy [equation](/page/Equation), now yielding the standard Bessel [equation](/page/Equation) $t^2 u''(t) + t u'(t) + (t^2 - \nu^2) u(t) = 0$ for $\nu = \pm 1/3$.[](https://dlmf.nist.gov/9.6)
The Airy functions emerge as special cases of solutions to the Bessel equation through this substitution, where the independent variable is scaled as $t \propto z^{3/2}$ and the dependent variable as $w(z) = \sqrt{z} \, u(t)$. In the limit as the Bessel order $\nu$ approaches $1/3$ (or $-1/3$), the solutions align precisely with Ai(z) and Bi(z), providing a unified framework for deriving properties of the Airy functions from the more general Bessel theory.[](https://dlmf.nist.gov/9.6)
These relations offer practical advantages, as the well-established asymptotic expansions of Bessel and modified Bessel functions directly imply those for the Airy functions, facilitating analysis in regions of large |z|.[](https://dlmf.nist.gov/9.6)
### Hypergeometric and Other Connections
The Airy function of the first kind, $\operatorname{Ai}(z)$, admits a power series expansion around $z = 0$ that can be compactly expressed using the confluent hypergeometric function ${}_0F_1(; b; w)$:
\operatorname{Ai}(z) = \frac{1}{3^{2/3} \Gamma\left(\frac{2}{3}\right)} {}_0F_1\left(; \frac{2}{3}; \frac{z^3}{9}\right) - \frac{z}{3^{1/3} \Gamma\left(\frac{1}{3}\right)} {}_0F_1\left(; \frac{4}{3}; \frac{z^3}{9}\right).
Similarly, the Airy function of the second kind, $\operatorname{Bi}(z)$, has the representation
\operatorname{Bi}(z) = \frac{1}{3^{1/6} \Gamma\left(\frac{2}{3}\right)} {}_0F_1\left(; \frac{2}{3}; \frac{z^3}{9}\right) + \frac{3^{1/6} z}{\Gamma\left(\frac{1}{3}\right)} {}_0F_1\left(; \frac{4}{3}; \frac{z^3}{9}\right).
These expressions follow directly from the Maclaurin series expansions of the Airy functions, where the terms correspond to the defining series of the ${}_0F_1$ function.[](https://dlmf.nist.gov/9.4)
The Airy functions also connect to the Kummer confluent hypergeometric function of the first kind, $M(a, b, w) = {}_1F_1(a; b; w)$, and Tricomi's confluent hypergeometric function of the second kind, $U(a, b, w)$. For positive real $z$, define the auxiliary variable $\zeta = \frac{2}{3} z^{3/2}$. Then,
\operatorname{Ai}(z) = 3^{-1/6} \pi^{-1/2} \zeta^{2/3} e^{-\zeta} U\left(\frac{5}{6}, \frac{5}{3}, 2\zeta\right),
which employs Tricomi's $U$ function to capture the decaying behavior for large positive arguments. For $\operatorname{Bi}(z)$,
\operatorname{Bi}(z) = \frac{1}{3^{1/6} \Gamma\left(\frac{2}{3}\right)} e^{-\zeta} M\left(\frac{1}{6}, \frac{1}{3}, 2\zeta\right) + \frac{3^{5/6}}{2^{2/3} \Gamma\left(\frac{1}{3}\right)} \zeta^{2/3} e^{-\zeta} M\left(\frac{5}{6}, \frac{5}{3}, 2\zeta\right).
These forms bridge the Airy equation to the confluent hypergeometric equation through limiting processes and are useful for asymptotic analysis and numerical evaluation in regions where the power series may converge slowly.[](https://dlmf.nist.gov/9.6)[](https://dlmf.nist.gov/13.2)
In the limit of the order parameter, the parabolic cylinder function $D_\nu(w)$ reduces to the Airy function as $\nu \to -\frac{1}{3}$. Specifically, a scaled version of $D_{-\frac{1}{3}}(w)$ approaches $\operatorname{Ai}\left(\frac{w}{2^{1/3}}\right)$ up to normalization constants, reflecting the confluence of the Weber parabolic cylinder differential equation to the Airy equation when the quadratic potential term dominates. This connection facilitates uniform asymptotic expansions for parabolic cylinder functions near turning points, akin to Airy-type approximations in WKB theory.[](https://arxiv.org/abs/1207.6861)
The Scorer functions $\operatorname{Gi}(z)$ and $\operatorname{Hi}(z)$, solutions to the inhomogeneous Airy equation $w''(z) - z w(z) = -\frac{1}{\pi}$, relate linearly to the homogeneous Airy functions via
\operatorname{Gi}(z) + \operatorname{Hi}(z) = \operatorname{Bi}(z).
They can also be expressed as
\operatorname{Gi}(z) = \frac{1}{3} \operatorname{Bi}(z) - \frac{z^2}{2\pi} {}_1F_2\left(1; \frac{4}{3}, \frac{5}{3}; \frac{z^3}{9}\right),
where the ${}_1F_2$ is a [generalized hypergeometric function](/page/Generalized_hypergeometric_function). These relations stem from integral representations of the Scorer functions, such as $\operatorname{Gi}(z) = \frac{1}{\pi} \int_0^\infty \sin\left(zt + \frac{t^3}{3}\right) dt$, which overlap conceptually with Airy integrals but include an inhomogeneous source term.[](https://dlmf.nist.gov/9.12)
## Integral Transforms
### Fourier Transform
The Fourier transform of the Airy function $\mathrm{Ai}(x)$ takes a remarkably simple form, $\hat{\mathrm{Ai}}(k) = e^{i k^3 / 3}$, under the [angular frequency](/page/Angular_frequency) convention defined by\hat{f}(k) = \int_{-\infty}^{\infty} f(x) , e^{-i k x} , dx.This expression holds in the sense of tempered distributions, reflecting the oscillatory decay of $\mathrm{Ai}(x)$ for $x < 0$ and exponential decay for $x > 0$.[](https://people.tamu.edu/~fulling/m412/f07/airywkb.pdf)
To derive this, apply the [Fourier transform](/page/Fourier_transform) to the Airy differential equation $w''(x) - x w(x) = 0$. The transform of the second derivative yields $-k^2 \hat{w}(k)$, while multiplication by $x$ in the spatial domain corresponds to $i \frac{d}{dk}$ in the [frequency domain](/page/Frequency_domain), giving the transformed equation-k^2 \hat{w}(k) = i \frac{d \hat{w}}{dk}.Rearranging produces the [first-order](/page/First-order) [ordinary differential equation](/page/Ordinary_differential_equation)\frac{d \hat{w}}{dk} = i k^2 \hat{w}(k),with solution $\hat{w}(k) = C \exp\left(i \int_0^k s^2 \, ds \right) = C \exp\left( i k^3 / 3 \right)$. The constant $C = 1$ is selected to match the normalization of $\mathrm{Ai}(x)$, ensuring [exponential decay](/page/Exponential_decay) for large positive arguments. Alternatively, [contour integration](/page/Contour_integration) techniques can evaluate the inverse transform directly, confirming the result via the method of stationary phase in the oscillatory regime.[](https://people.tamu.edu/~fulling/m412/f07/airywkb.pdf)
The inverse relation expresses $\mathrm{Ai}(x)$ as\mathrm{Ai}(x) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \exp\left( i k x + i \frac{k^3}{3} \right) , dk,where the factor $1/(2\pi)$ arises from the asymmetric [convention](/page/Convention) (forward transform without $1/(2\pi)$, inverse with it). This [Fourier](/page/Fourier) integral representation aligns with the real oscillatory form\mathrm{Ai}(x) = \frac{1}{\pi} \int_0^{\infty} \cos\left( \frac{t^3}{3} + x t \right) , dt,obtained by taking the real part and exploiting the even nature of the integrand.[](https://people.tamu.edu/~fulling/m412/f07/airywkb.pdf)
Properties of the transform leverage standard [Fourier](/page/Fourier) theorems adapted to this cubic phase. The [convolution theorem](/page/Convolution_theorem) states that if $h(x) = (\mathrm{Ai} * g)(x) = \int_{-\infty}^{\infty} \mathrm{Ai}(x - u) g(u) \, du$, then $\hat{h}(k) = e^{i k^3 / 3} \hat{g}(k)$, enabling efficient computation of convolutions involving Airy kernels in linear systems. [Parseval's theorem](/page/Parseval's_theorem) applies in a distributional [sense](/page/Sense), relating integrals of $\mathrm{Ai}(x)$ products to those over the unit-magnitude [phase factor](/page/Phase_factor), though the non-$L^2$ integrability of $\mathrm{Ai}(x)$ requires careful [interpretation](/page/Interpretation) via limits or principal values. Different [normalization](/page/Normalization) conventions, such as the unitary transform with $1/\sqrt{2\pi}$ factors in both directions, rescale the result to $\hat{\mathrm{Ai}}(k) = \sqrt{2\pi} \, e^{i k^3 / 3}$. In higher dimensions, radial or separable extensions appear in contexts like 2D [diffraction](/page/Diffraction) patterns, where the 1D transform serves as a building block for cylindrical or spherical symmetries.[](https://people.tamu.edu/~fulling/m412/f07/airywkb.pdf)
### Laplace Transform Applications
The Laplace transform of the Airy function of the first kind, $\mathcal{L}\{\Ai(t)\}(p) = \int_0^\infty \Ai(t) e^{-p t} \, dt$, admits a closed-form expression in terms of confluent hypergeometric functions:
e^{-p^3/3} \left[ \frac{1}{3} - \frac{p \cdot {}_1F_1\left( \frac{1}{3}; \frac{4}{3}; \frac{p^3}{3} \right)}{3^{4/3} \Gamma\left( \frac{4}{3} \right)} + \frac{p^2 \cdot {}_1F_1\left( \frac{2}{3}; \frac{5}{3}; \frac{p^3}{3} \right)}{3^{5/3} \Gamma\left( \frac{5}{3} \right)} \right],
valid for $p \in \mathbb{C}$ in a suitable region ensuring convergence. This representation facilitates analytical manipulation and inversion in applications. For the Airy function of the second kind, $\Bi(t)$, the standard Laplace transform from 0 to $\infty$ diverges [due](/page/A_due) to the asymptotic growth $\Bi(t) \sim \frac{1}{\sqrt{\pi t^{1/4}}} \exp\left( \frac{2}{3} t^{3/2} \right)$ for large positive $t$, which outpaces the exponential decay $e^{-p t}$ for any fixed $\Re p > 0$. However, the transform of $\Bi(-t)$ converges and is given by
\frac{1}{\sqrt{3}} e^{p^3/3} \left[ \frac{\Gamma\left( \frac{2}{3}, \frac{p^3}{3} \right)}{\Gamma\left( \frac{2}{3} \right)} - \frac{\Gamma\left( \frac{1}{3}, \frac{p^3}{3} \right)}{\Gamma\left( \frac{1}{3} \right)} \right],
for $\Re p > 0$, where $\Gamma(s, x)$ denotes the upper incomplete gamma function; this form accounts for the oscillatory behavior of $\Bi(-t)$.
These transforms prove invaluable for solving initial value problems of the form $y''(x) = x y(x) + f(x)$, $y(0) = a$, $y'(0) = b$, which model systems with linear potentials, such as in quantum mechanics or wave propagation. Applying the Laplace transform yields the first-order equation $Y'(s) + s^2 Y(s) = \mathcal{L}\{f\}(s) + s a + b$, where $Y(s) = \mathcal{L}\{y\}(s)$. The integrating factor is $e^{s^3/3}$, leading to
Y(s) = e^{-s^3/3} \int_s^\infty e^{u^3/3} \left( \mathcal{L}{f}(u) + u a + b \right) , du,
assuming the boundary condition $Y(\infty) = 0$ for the physically relevant decaying solution. For $f \equiv 0$, the inverse transform recovers linear combinations of $\Ai(x)$ and $\Bi(x)$ matching the initial conditions, with $\Ai(0) = 3^{-2/3} / \Gamma(2/3)$ and $\Ai'(0) = -3^{-1/3} / \Gamma(1/3)$. For general $f$, the solution involves a convolution integral whose inversion may require numerical methods or series expansions, but the Laplace domain simplifies handling of initial conditions and forcing terms.
Numerical computation of these transforms leverages the hypergeometric or incomplete gamma representations for moderate $|p|$, with asymptotic expansions (e.g., via [integration by parts](/page/Integration_by_parts) on the defining [integral](/page/Integral)) for large $|p|$ to achieve high precision. For instance, the leading asymptotic for $\mathcal{L}\{\Ai(t)\}(p)$ as $|p| \to \infty$ in $|\arg p| < \pi/3$ is $\frac{e^{-p^3/3}}{3 p^{1/3}}$, derived from the integral form of $\Ai(t)$. Software implementations often employ adaptive quadrature on the original [integral](/page/INTEGRAL) $\int_0^\infty \Ai(t) e^{-p t} \, dt$, truncated at a point where $\Ai(t)$ decays sufficiently, ensuring accuracy for engineering and scientific simulations.[](https://dlmf.nist.gov/9.10)
## Physical and Scientific Applications
### Quantum Tunneling and Mechanics
In quantum mechanics, the Airy function arises as the exact solution to the time-independent Schrödinger equation for a particle subject to a linear potential, such as that induced by a constant force $F$. The equation takes the form
-\frac{\hbar^2}{2m} \frac{d^2 \psi}{dx^2} + F x \psi = E \psi,
where $\hbar$ is the reduced Planck's constant, $m$ is the particle mass, and $E$ is the energy. By introducing the scaled variable $z = \left( \frac{2m F}{\hbar^2} \right)^{1/3} \left( x - \frac{E}{F} \right)$, the equation reduces to the standard Airy equation $\frac{d^2 \psi}{dz^2} - z \psi = 0$, with the wavefunction given by $\psi(x) \propto \mathrm{Ai}(z)$.[](https://www.classe.cornell.edu/~dlr/teaching/p443/LinearPotential.pdf)[](https://www.sciencedirect.com/science/article/pii/S2211379721009384)
The Airy function $\mathrm{Ai}(z)$ describes both bound and continuum states depending on boundary conditions. For bound states in a potential well terminated by an infinite wall (e.g., a triangular well), the energies are determined by the zeros of $\mathrm{Ai}(z)$, providing quantized levels. In the classically forbidden region where $x > E/F$ (assuming $F > 0$), the wavefunction exhibits [exponential decay](/page/Exponential_decay), reflecting quantum tunneling. The asymptotic behavior of $\mathrm{Ai}(z)$ for large positive $z$ is $\mathrm{Ai}(z) \sim \frac{1}{2 \sqrt{\pi} z^{1/4}} \exp\left( -\frac{2}{3} z^{3/2} \right)$, which yields the tunneling probability through a linear barrier via the [WKB approximation](/page/WKB_approximation) or exact matching.[](https://bingweb.binghamton.edu/~suzuki/QM_Graduate/WKB_approximation_I.pdf)[](https://pubs.aip.org/aip/jap/article/106/11/113701/900334/Tunneling-escape-time-from-a-semiconductor-quantum)
This framework applies to the linear Stark effect in atomic systems, where an external electric field $F = e \mathcal{E}$ (with $e$ the [electron](/page/Electron) charge and $\mathcal{E}$ the field strength) perturbs hydrogen-like atoms, leading to wavefunctions involving Airy functions and [energy](/page/Energy) shifts proportional to $\mathcal{E}$. In gravitational quantum states, ultracold [neutron](/page/Neutron)s in Earth's field $F = m [g](/page/M&G)$ ( $[g](/page/M&G)$ the [acceleration due to gravity](/page/Acceleration_due_to_gravity)) form bound states described by Airy functions, as verified in [interferometry](/page/Interferometry) experiments measuring quantized [energy](/page/Energy) levels and [wave packet](/page/Wave_packet) dynamics. Recent realizations include the generation of propagating neutron Airy beams, observing the function's zeros critical for these states.[](https://link.springer.com/article/10.1007/s40094-014-0140-x)[](https://arxiv.org/pdf/1207.2953)[](https://link.aps.org/doi/10.1103/PhysRevLett.134.153401)
Airy wave packets, solutions propagating without [diffraction](/page/Diffraction) in linear potentials, have been experimentally realized in the 2020s with cold atomic gases and [semiconductor](/page/Semiconductor) structures. In atom lasers from Bose-Einstein condensates, gravitational caustics produce Airy-like [matter wave](/page/Matter_wave) profiles, enabling studies of tunneling and [interference](/page/Interference) in controlled linear fields. In [semiconductor](/page/Semiconductor) quantum wells under tilted [electric fields](/page/Electric_Fields), Airy functions model [electron](/page/Electron) tunneling escape times, with recent chip-scale experiments enhancing precision in quantum devices.[](https://www.nature.com/articles/s41467-021-27555-3)[](https://pubs.aip.org/aip/jap/article/106/11/113701/900334/Tunneling-escape-time-from-a-semiconductor-quantum)
### Optics and Diffraction
In optics, the Airy functions play a key role in describing diffraction phenomena, particularly through uniform approximations that capture wave behavior near focal points or turning points in ray optics. In Fraunhofer diffraction through a circular aperture of radius $a$, the far-field intensity pattern is given by $ I(\theta) \propto \left[ \frac{2 J_1(k a \sin \theta)}{k a \sin \theta} \right]^2 $, where $ J_1 $ is the first-order Bessel function of the first kind, $ k = 2\pi / \lambda $ is the wavenumber, and $ \theta $ is the angular deviation from the optical axis.[](https://opg.optica.org/josaa/abstract.cfm?uri=josaa-31-5-914) This pattern forms the Airy disk, the central bright spot surrounded by concentric rings, which represents the fundamental limit of resolution for imaging systems like telescopes. The radius of the Airy disk to its first minimum is approximately $ 1.22 \lambda / D $, where $ D = 2a $ is the aperture diameter, establishing the scale for diffraction-limited performance.[](https://www.edmundoptics.com/knowledge-center/application-notes/imaging/limitations-on-resolution-and-contrast-the-airy-disk/) Near the focus, uniform approximations using Airy functions provide a more accurate description of the wave field by accounting for the turning point where rays converge, extending the validity of semiclassical methods beyond simple geometric optics.[](https://dlmf.nist.gov/9.16)
Airy functions also arise directly in solutions to the paraxial wave equation for beam propagation in media with a linear refractive index profile, such as $ n(x) = n_0 + \alpha x $, where $ n_0 $ is the base index and $ \alpha $ characterizes the gradient. In this case, the paraxial Helmholtz equation reduces to the Airy differential equation, yielding exact solutions of the form involving $ \mathrm{Ai}(s) $, where $ s $ is a scaled transverse coordinate incorporating the propagation distance and gradient strength. These solutions describe guided modes that maintain beam integrity over extended distances, contrasting with free-space diffraction. Such configurations are realized in engineered graded-index materials, enabling controlled light routing in integrated optics.[](https://opg.optica.org/josaa/abstract.cfm?uri=josaa-31-5-914)
Finite-energy Airy beams, truncated versions of ideal Airy solutions to the paraxial equation, exhibit non-diffracting [propagation](/page/Propagation) with parabolic trajectories and self-healing properties, where the [main lobe](/page/Main_lobe) reconstructs after obstructions. The transverse profile is typically $ \psi(x,0) = \mathrm{Ai}\left( s - \frac{(x/w)^2}{4} \right) \exp\left( i \frac{(x/w)^2}{2} \right) $, with $ s $ a scaling parameter, $ w $ a width, and the [exponential](/page/Exponential) providing finite energy; during [propagation](/page/Propagation), the [intensity](/page/Intensity) follows $ I(x,z) \propto \left| \mathrm{Ai}\left( s - \frac{(x - \beta z^2 / 4)^2}{4} \right) \right|^2 \exp(-\text{attenuation terms}) $, where $ \beta $ governs acceleration. These beams were first theoretically proposed and experimentally observed in [2007](/page/2007), revolutionizing beam shaping.[](https://opg.optica.org/ol/abstract.cfm?uri=ol-32-8-979)[](https://link.aps.org/doi/10.1103/PhysRevLett.99.213901)
Recent applications leverage Airy beams for advanced optical manipulation and fabrication. In [optical tweezers](/page/Optical_tweezers), focused Airy beams create multiple stable traps along curved paths, enabling precise 3D control of microparticles with reduced diffraction blurring, as demonstrated in high-numerical-aperture setups post-2020.[](https://opg.optica.org/ao/abstract.cfm?uri=ao-50-1-43)[](https://pubs.acs.org/doi/10.1021/acsphotonics.2c01986) For [lithography](/page/Lithography) and micro-machining, their self-healing and non-diffracting nature improves resolution in beam scanning, allowing sub-wavelength patterning in photoresists with enhanced uniformity, particularly in metasurface-generated configurations since 2021.[](https://www.mdpi.com/2304-6732/12/8/794)
### Caustics and Wave Phenomena
In classical wave theory, Airy functions provide a uniform asymptotic description of wave fields near fold caustics, where ray paths coalesce and geometric optics predicts infinite intensity. The simplest such singularity, known as the A_2 catastrophe, is modeled by the diffraction integral whose exact solution is the Airy function Ai(z). For a wave field near a fold caustic, the uniform approximation takes the form
u \sim \mathrm{Ai}\left( -\left(\frac{3}{2} \phi \right)^{2/3} \zeta \right) \pi^{1/2} \left(\frac{3}{2} \phi \right)^{1/6},
where $\phi$ represents the phase difference across the coalescing stationary points, and $\zeta$ is the scaled coordinate perpendicular to the caustic. This expression captures the transition from oscillatory behavior on the illuminated side ($\zeta < 0$) to exponential decay in the shadow ($\zeta > 0$), with the Airy function's oscillations accounting for [interference](/page/Interference) fringes adjacent to the caustic.
A prominent application arises in the formation of [rainbow](/page/Rainbow) caustics from [light](/page/Light) scattering by spherical droplets, such as raindrops. Airy theory approximates the scattered [intensity](/page/Intensity) near the rainbow angle by treating the phase stationary points as a fold caustic, where the Airy function's oscillatory tail explains the supernumerary [interference](/page/Interference) fringes inside the primary [rainbow](/page/Rainbow). This model, refined from geometric [optics](/page/Optics), predicts a finite maximum [intensity](/page/Intensity) at the caustic edge, with the fringe spacing inversely proportional to the droplet radius, enabling [quantitative analysis](/page/Quantitative_analysis) of [rainbow](/page/Rainbow) patterns observed in nature.
In acoustics and [seismology](/page/Seismology), Airy functions model wave amplification and [phase](/page/Phase) shifts near [caustics](/page/Caustic) formed by Earth's heterogeneous structure, such as those from the core-mantle boundary. Seismic waves crossing such folds exhibit enhanced amplitudes described by the [Airy integral](/page/Integral), with the uniform [approximation](/page/Approximation) yielding a [phase](/page/Phase) advance of $\pi/4$ at the caustic and oscillatory precursors that match observed triplicated waveforms in long-period seismograms. This approach extends beyond simple geometric predictions, incorporating [diffraction](/page/Diffraction) to resolve shadow zones and illuminated regions accurately.[](https://pubs.geoscienceworld.org/ssa/bssa/article-abstract/71/5/1445/118219/Amplitude-and-phase-shift-due-to-caustics)
Catastrophe optics formalizes these phenomena, identifying the Airy integral as the canonical diffraction pattern for the A_2 (fold) singularity in generic ray families. Developed as part of Thom's classification, this framework unifies the morphologies of caustics—such as cusps and swallows-tails—with their wave corrections, where the Airy function provides the leading-order uniform approximation across the singularity. Seminal analyses demonstrate how [structural stability](/page/Structural_stability) ensures these patterns appear universally in optical, acoustic, and elastic wave propagation near simple folds.[](https://doi.org/10.1016/S0079-6638(08)70009-4)
Modern applications leverage these properties in high-resolution imaging and simulations. In ultrasound, Airy-based models describe focusing at fold caustics in nonlinear acoustic fields, enabling precise control of shock wave amplification for therapeutic beamforming without singularities.[](https://pubs.aip.org/asa/jasa/article/127/4/2129/897993/A-detailed-analysis-about-penumbra-caustics) Similarly, 2020s wave-optics simulations of gravitational lensing incorporate Airy functions to resolve diffraction near fold caustics in binary lens systems, predicting interference patterns in lensed gravitational wave signals from black hole mergers that deviate from ray-tracing approximations.[](https://arxiv.org/abs/2504.10128)
### Probabilistic and Statistical Uses
In random matrix theory, the Airy kernel arises as the limiting kernel for the determinantal point process describing eigenvalue correlations at the soft edge of the spectrum in Gaussian Orthogonal Ensemble (GOE) and Gaussian Unitary Ensemble (GUE) models. Specifically, for the edge scaling limit, the kernel takes the form
K(x,y) = \frac{\mathrm{Ai}(x) \mathrm{Ai}'(y) - \mathrm{Ai}'(x) \mathrm{Ai}(y)}{x - y},
where $\mathrm{Ai}$ denotes the Airy function of the first kind and $\mathrm{Ai}'$ its derivative; this structure governs the local statistics of eigenvalues near the largest one after appropriate rescaling by $n^{2/3}$, with $n$ the matrix dimension. The resulting fluctuations of the largest eigenvalue, centered and scaled, converge to the Tracy-Widom distribution, a cornerstone for understanding extreme value statistics in these ensembles. For the GUE case ($\beta=2$), the cumulative distribution function is given by
F_2(s) = \exp\left( -\int_s^\infty q(r)^2 , dr \right),
where $q(r)$ solves the Painlevé II equation $q''(r) = r q(r) + 2 q(r)^3$ with the Hastings-McLeod boundary condition $q(r) \sim \mathrm{Ai}(r)$ as $r \to +\infty$.
Beyond spectral theory, Airy functions appear in the asymptotics of combinatorial objects, particularly in the distribution of the longest increasing subsequence in a random permutation of $\{1, \dots, n\}$. The length $L_n$, after centering at $2\sqrt{n}$ and scaling by $n^{1/6}$, converges in distribution to the Tracy-Widom law $F_2$, linking discrete enumeration problems to continuous random matrix limits via integrable structures like the Robinson-Schensted-Knuth correspondence. Similar Airy-driven edge asymptotics govern fluctuations in random tilings, such as domino tilings of the Aztec diamond, where boundary height variations follow the Airy process, reflecting universal scaling in two-dimensional growth models.
In [statistical mechanics](/page/Statistical_mechanics), the Airy function underpins the fixed point of the Kardar-Parisi-Zhang (KPZ) [universality class](/page/Universality_class), which describes interface growth and directed polymers in random media. The one-point height distribution at the KPZ fixed point coincides with the Tracy-Widom $F_2$, while multi-point correlations are captured by the Airy sheet—a random surface whose finite-dimensional marginals involve Airy processes, ensuring scale-invariant fluctuations across models like the totally asymmetric simple exclusion process or polynuclear growth in the long-time limit. This framework unifies one- and higher-dimensional stochastic growth phenomena through exact solvability via the replica method or Macdonald processes.[](https://projecteuclid.org/journals/acta-mathematica/volume-227/issue-1/The-KPZ-fixed-point/10.4310/ACTA.2021.v227.n1.a3.pdf)
Recent developments (2020–2025) extend these probabilistic roles to non-Hermitian random matrices, where deformed Airy kernels interpolate between Hermitian edge statistics and bulk behaviors in elliptic Ginibre ensembles, enabling universal descriptions of [complex](/page/Complex) eigenvalue clustering near spectral boundaries; analogous Tracy-Widom-like laws emerge for extreme real eigenvalues in real non-Hermitian settings, with applications to stability analysis in disordered systems. In quantum chaos contexts, Airy-mediated edge fluctuations align chaotic [billiard](/page/Billiard) spectra with random matrix predictions, informing recent studies of non-equilibrium dynamics in open [quantum systems](/page/Quantum-Systems).[](https://doi.org/10.1007/s00220-023-04751-3)
## Historical Context
### Discovery and Early Development
The Airy functions originated in the work of British astronomer and mathematician [George Biddell Airy](/page/George_Biddell_Airy), who derived them in 1838 while studying optical [caustics](/page/Caustic), particularly those associated with rainbows. In his seminal paper "On the Intensity of Light in the Neighbourhood of a [Caustic](/page/Caustic)," published in the Transactions of the [Cambridge](/page/Cambridge) Philosophical Society, Airy addressed the limitations of [geometrical optics](/page/geometrical_optics) in explaining the intensity distribution near caustics, such as the bright edge and supernumerary fringes observed in rainbows. He modeled the problem by considering light propagation in a medium where the [refractive index](/page/refractive_index) varies linearly with position near the turning point of rays, leading to a second-order [ordinary differential equation](/page/ordinary_differential_equation) whose solutions are expressed as integrals now known as Airy integrals.[](https://www.ams.org/publicoutreach/feature-column/fc-2015-02)[](https://dlmf.nist.gov/9.16)
Airy's mathematical formalization involved approximating the ray paths in a stratified medium with a linear density gradient for the [refractive index](/page/Refractive_index), transforming the wave [equation](/page/Equation) into a form amenable to integral representation. This approach provided a more accurate description of [light intensity](/page/Light_intensity) oscillations beyond the [caustic](/page/Caustic), capturing effects like the interference patterns in rainbow supernumeraries that earlier ray-tracing methods by Descartes and [Newton](/page/Newton) could not. Airy himself performed initial numerical computations and tabulations of these integrals specifically tailored to rainbow geometry, enabling quantitative predictions of [fringe](/page/Fringe) spacing and brightness.[](https://www.philiplaven.com/p8a.html)
Shortly after, in the late [1840s](/page/1840s), G.G. Stokes extended Airy's work by providing asymptotic expansions for the Airy integral using [divergent series](/page/Divergent_series) and expressing it in terms of [Bessel functions](/page/Bessel_function) of fractional orders ±1/3, which facilitated computations for large arguments and highlighted connections to cylindrical wave solutions developed by Friedrich Wilhelm Bessel in the [1820s](/page/1820s). Though Airy was the first to isolate this particular integral form for linear caustics, Stokes' analysis was pivotal for its application in diffraction theory. Initially termed the "Airy integral" in reference to its oscillatory representation, the functions received their modern designation as "Airy functions" only after the [1920s](/page/1920s), as they gained prominence in diverse mathematical contexts.[](https://www.ams.org/publicoutreach/feature-column/fc-2015-02)[](https://doras.dcu.ie/19203/)
### Key Mathematical Contributions
In the 1920s, G. N. Watson provided a systematic treatment of the asymptotic expansions of [Bessel functions](/page/Bessel_function) and their connections to the Airy functions, particularly through integral representations and limiting cases involving fractional orders such as 1/3 and -1/3. In his seminal treatise, Watson derived rigorous asymptotic forms for [Bessel functions](/page/Bessel_function) of large argument or order, showing how the Airy functions emerge as uniform approximations near turning points in the [differential equation](/page/Differential_equation), with explicit error bounds using methods like steepest descent and Hankel's contour integrals. These developments established the Airy functions as a key bridge between [Bessel functions](/page/Bessel_function) and more general asymptotic theory for linear second-order [differential equation](/page/Differential_equation)s.[](https://archive.org/details/treatiseontheory00watsuoft)
During the [1940s](/page/1940s), amid wartime demands for accurate numerical computations in physics and [engineering](/page/Engineering), researchers computed extensive tables of Airy functions and their zeros to support applications in wave propagation and [boundary layer](/page/Boundary_layer) theory. V. A. Fock calculated tables of Airy functions for complex arguments, providing values essential for solving turning-point problems in [quantum mechanics](/page/Quantum_mechanics) and [aerodynamics](/page/Aerodynamics), with computations extending to high precision for oscillatory and [exponential](/page/Exponential) behaviors. These efforts, often conducted under resource constraints, laid the [foundation](/page/Foundation) for reliable numerical evaluation of Airy zeros and integrals, influencing subsequent handbook compilations.[](https://ntrs.nasa.gov/api/citations/20050028467/downloads/20050028467.pdf)
By the 1950s, the Airy functions were formally classified within the canon of [special functions](/page/Special_functions) in authoritative handbooks published by the [American Mathematical Society](/page/American_Mathematical_Society) and the National Bureau of Standards, reflecting their growing role in [asymptotic analysis](/page/Asymptotic_analysis) and differential equations. The [AMS](/page/A!MS) series included dedicated sections on Airy integrals and zeros, with tables for real and complex arguments up to 10 decimal places for the first several roots. Concurrently, L. J. Slater established explicit links between Airy functions and confluent hypergeometric functions through series representations, such as Ai(z) expressed as a [limit](/page/Limit) of the Kummer function M(a, b, w) with appropriate scaling, enabling new summation formulas and analytic continuations. These contributions integrated the Airy functions into the broader hypergeometric framework, facilitating computations and theoretical extensions.
Painlevé analysis reveals a deep connection between the Airy equation and the first [Painlevé transcendent](/page/Painlevé_transcendent), where the Airy functions arise as a linear scaling limit of the nonlinear [Painlevé I equation](/page/Painlevé_I_equation) y'' = 6y² + x. Specifically, by rescaling y = ε u and x = ε^{-2/3} ζ with ε → 0, the nonlinear term vanishes, recovering the Airy equation u'' = ζ u, which governs the asymptotic behavior of [Painlevé I](/page/Painlevé_I) solutions near the origin or in tronquée forms. This relation underscores the integrability properties of both, with Airy functions providing exact solutions in the perturbative regime.[](https://www.math.ucdavis.edu/~tracy/selectedPapers/2010s/CV92.pdf)
In recent decades up to 2025, enhancements in numerical libraries have significantly improved the computation of Airy functions, incorporating high-precision algorithms for zeros, integrals, and complex arguments. The GNU Scientific Library (GSL) version 2.8 (released May 2024) and later includes optimized routines for Airy functions using continued fractions and power series, with updates in the 2020s focusing on vectorized implementations for parallel computing. Similarly, SciPy's special module, updated through releases like 1.16.3 (October 2025), employs Cephes-based methods with rational approximations for real arguments in [-10, 10] and asymptotic expansions elsewhere, achieving machine epsilon accuracy. These advancements support applications in scientific computing, while ongoing research explores Airy functions in integrable systems, such as exact solutions to nonlinear evolution equations like the modified KdV via Airy ansätze. Additionally, 2020s studies have applied Airy functions to model dispersive waves in integrable hierarchies, revealing new soliton structures through scaling limits.[](https://www.gnu.org/software/gsl/doc/html/specfunc.html)[](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.airy.html)[](https://www.researchgate.net/publication/361714966_Airy_function_solution_of_two_classes_of_nonlinear_evolution_equations)